PrintGeomInfo/plugins/PrintMaterialBudgetInfo.cc

0001 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0002 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0003 #include "SimG4Core/Notification/interface/Observer.h"
0004 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0005
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0009 #include "DetectorDescription/Core/interface/DDCompactView.h"
0010 #include "DetectorDescription/Core/interface/DDFilter.h"
0011 #include "DetectorDescription/Core/interface/DDFilteredView.h"
0012 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
0013 #include "DetectorDescription/Core/interface/DDSplit.h"
0014 #include "DetectorDescription/Core/interface/DDValue.h"
0015
0016 #include "G4Run.hh"
0017 #include "G4PhysicalVolumeStore.hh"
0018 #include "G4LogicalVolumeStore.hh"
0019 #include "G4VPhysicalVolume.hh"
0020 #include "G4LogicalVolume.hh"
0021 #include "G4VSolid.hh"
0022 #include "G4Material.hh"
0023 #include "G4NavigationHistory.hh"
0024 #include "G4Track.hh"
0025 #include "G4VisAttributes.hh"
0026 #include "G4UserLimits.hh"
0027 #include "G4TransportationManager.hh"
0028 #include "G4UnitsTable.hh"
0029 #include "Randomize.hh"
0030
0031 #include <iostream>
0032 #include <fstream>
0033 #include <memory>
0034 #include <set>
0035 #include <vector>
0036 #include <string>
0037
0038 typedef std::map<G4VPhysicalVolume*, G4VPhysicalVolume*, std::less<G4VPhysicalVolume*> > mpvpv;
0039 typedef std::multimap<G4LogicalVolume*, G4VPhysicalVolume*, std::less<G4LogicalVolume*> > mmlvpv;
0040
0041 class PrintMaterialBudgetInfo : public SimWatcher,
0042                                 public Observer<const BeginOfJob*>,
0043                                 public Observer<const BeginOfRun*> {
0044 public:
0045   PrintMaterialBudgetInfo(edm::ParameterSet const& p);
0046   ~PrintMaterialBudgetInfo() override;
0047
0048 private:
0049   void update(const BeginOfJob* job) override {}
0050   void update(const BeginOfRun* run) override;
0051   void dumpHeader(std::ostream& out = G4cout);
0052   void dumpLaTeXHeader(std::ostream& out = G4cout);
0053   void dumpHierarchyLeaf(G4VPhysicalVolume* pv,
0054                          G4LogicalVolume* lv,
0055                          unsigned int leafDepth,
0056                          std::ostream& weightOut = G4cout,
0057                          std::ostream& texOut = G4cout);
0058   void printInfo(G4VPhysicalVolume* pv,
0059                  G4LogicalVolume* lv,
0060                  unsigned int leafDepth,
0061                  std::ostream& weightOut = G4cout,
0062                  std::ostream& texOut = G4cout);
0063   void dumpElementMassFraction(std::ostream& elementOut = G4cout);
0064   void dumpLaTeXFooter(std::ostream& out = G4cout);
0065
0066 private:
0067   std::string name;
0068   int nchar;
0069   mpvpv thePVTree;
0070   G4NavigationHistory fHistory;
0071   bool volumeFound;
0072   unsigned int levelFound;
0073   std::ofstream weightOutputFile;
0074   std::ofstream elementOutputFile;
0075   std::ofstream texOutputFile;
0076   std::vector<std::string> elementNames;
0077   std::vector<double> elementTotalWeight;
0078   std::vector<double> elementWeightFraction;
0079
0080   std::string stringLaTeXUnderscore(std::string stringname);
0081   std::string stringLaTeXSuperscript(std::string stringname);
0082 };
0083
0084 PrintMaterialBudgetInfo::PrintMaterialBudgetInfo(const edm::ParameterSet& p) {
0085   name = p.getUntrackedParameter<std::string>("Name", "*");
0086   nchar = name.find('*');
0087   name.assign(name, 0, nchar);
0088   G4cout << "PrintMaterialBudget selected volume " << name << G4endl;
0089   volumeFound = false;
0090   std::string weightFileName = name + ".weight";
0091   weightOutputFile.open(weightFileName.c_str());
0092   std::string elementFileName = name + ".element";
0093   elementOutputFile.open(elementFileName.c_str());
0094   std::string texFileName = name + "_table.tex";
0095   texOutputFile.open(texFileName.c_str());
0096   G4cout << "PrintMaterialBudget output file " << weightFileName << G4endl;
0097   G4cout << "PrintMaterialBudget output file " << elementFileName << G4endl;
0098   G4cout << "PrintMaterialBudget output file " << texFileName << G4endl;
0099   elementNames.clear();
0100   elementTotalWeight.clear();
0101   elementWeightFraction.clear();
0102 }
0103
0104 PrintMaterialBudgetInfo::~PrintMaterialBudgetInfo() {}
0105
0106 void PrintMaterialBudgetInfo::update(const BeginOfRun* run) {
0107   [[clang::suppress]] G4Random::setTheEngine(new CLHEP::RanecuEngine);
0108   // Physical Volume
0109   G4VPhysicalVolume* theTopPV =
0110       G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume();
0111   assert(theTopPV);
0112   // Logical Volume
0113   G4LogicalVolume* lv = theTopPV->GetLogicalVolume();
0114   unsigned int leafDepth = 0;
0115   // the first time fill the vectors of elements
0116   if (elementNames.empty() && elementTotalWeight.empty() && elementWeightFraction.empty()) {
0117     for (unsigned int iElement = 0; iElement < G4Element::GetNumberOfElements();
0118          iElement++) {  // first element in table is 0
0119       elementNames.push_back("rr");
0120       elementTotalWeight.push_back(0);
0121       elementWeightFraction.push_back(0);
0122     }
0123   }
0124   dumpHeader(weightOutputFile);
0125   dumpLaTeXHeader(texOutputFile);
0126   dumpHierarchyLeaf(theTopPV, lv, leafDepth, weightOutputFile, texOutputFile);
0127   dumpElementMassFraction(elementOutputFile);
0128   dumpLaTeXFooter(texOutputFile);
0129   //
0130 }
0131
0132 void PrintMaterialBudgetInfo::dumpHeader(std::ostream& out) {
0133   out << "Geom."
0134       << "\t"
0135       << "Volume"
0136       << "\t"
0137       << "\t"
0138       << "Copy"
0139       << "\t"
0140       << "Solid"
0141       << "\t"
0142       << "\t"
0143       << "Material"
0144       << "\t"
0145       << "Density"
0146       << "\t"
0147       << "\t"
0148       << "Mass"
0149       << "\t"
0150       << "\t" << G4endl;
0151   out << "Level"
0152       << "\t"
0153       << "Name"
0154       << "\t"
0155       << "\t"
0156       << "Number"
0157       << "\t"
0158       << "Name"
0159       << "\t"
0160       << "\t"
0161       << "Name"
0162       << "\t"
0163       << "\t"
0164       << "[g/cm3]"
0165       << "\t"
0166       << "\t"
0167       << "[g]    "
0168       << "\t"
0169       << "\t" << G4endl;
0170 }
0171
0172 void PrintMaterialBudgetInfo::dumpLaTeXHeader(std::ostream& out) {
0173   out << "\\begin{table}[h!]" << G4endl << "  \\caption{\\textsf {" << name << "} volume list.}" << G4endl
0174       << "  \\label{tab: " << name << "}" << G4endl << "  \\begin{center}" << G4endl << "    \\begin{tabular}{ccccccc}"
0175       << G4endl << "      \\hline" << G4endl;
0176   out << "      Geom."
0177       << "\t & "
0178       << "      Volume"
0179       << "\t & "
0180       << "      Copy"
0181       << "\t & "
0182       << "      Solid"
0183       << "\t & "
0184       << "      Material"
0185       << "\t & "
0186       << "      Density"
0187       << "\t & "
0188       << "      Mass"
0189       << "\t \\\\ " << G4endl;
0190   out << "      Level"
0191       << "\t & "
0192       << "      Name"
0193       << "\t & "
0194       << "      Number"
0195       << "\t & "
0196       << "      Name"
0197       << "\t & "
0198       << "      Name"
0199       << "\t & "
0200       << "                "
0201       << "\t & "
0202       << "                "
0203       << "\t \\\\ " << G4endl << "      \\hline\\hline" << G4endl;
0204 }
0205
0206 void PrintMaterialBudgetInfo::dumpLaTeXFooter(std::ostream& out) {
0207   out << "      \\hline" << G4endl << "    \\end{tabular}" << G4endl << "  \\end{center}" << G4endl << "\\end{table}"
0208       << G4endl;
0209 }
0210
0211 void PrintMaterialBudgetInfo::dumpHierarchyLeaf(
0212     G4VPhysicalVolume* pv, G4LogicalVolume* lv, unsigned int leafDepth, std::ostream& weightOut, std::ostream& texOut) {
0213   if (volumeFound && (leafDepth <= levelFound))
0214     return;
0215   if (volumeFound && (leafDepth > levelFound))
0216     printInfo(pv, lv, leafDepth, weightOut, texOut);
0217
0218   // choose mother volume
0219   std::string lvname = lv->GetName();
0220   lvname.assign(lvname, 0, nchar);
0221   if (lvname == name) {
0222     volumeFound = true;
0223     levelFound = leafDepth;
0224     printInfo(pv, lv, leafDepth, weightOut, texOut);
0225     texOut << "      \\hline" << G4endl;
0226   }
0227
0228   //----- Get LV daughters from list of PV daughters
0229   mmlvpv lvpvDaughters;
0230   std::set<G4LogicalVolume*> lvDaughters;
0231   int NoDaughters = lv->GetNoDaughters();
0232   while ((NoDaughters--) > 0) {
0233     G4VPhysicalVolume* pvD = lv->GetDaughter(NoDaughters);
0234     lvpvDaughters.insert(mmlvpv::value_type(pvD->GetLogicalVolume(), pvD));
0235     lvDaughters.insert(pvD->GetLogicalVolume());
0236   }
0237
0238   std::set<G4LogicalVolume*>::const_iterator scite;
0239   mmlvpv::const_iterator mmcite;
0240
0241   //----- Dump daughters PV and LV
0242   for (scite = lvDaughters.begin(); scite != lvDaughters.end(); scite++) {
0243     std::pair<mmlvpv::iterator, mmlvpv::iterator> mmER = lvpvDaughters.equal_range(*scite);
0244     //----- Dump daughters PV of this LV
0245     for (mmcite = mmER.first; mmcite != mmER.second; mmcite++)
0246       dumpHierarchyLeaf((*mmcite).second, *scite, leafDepth + 1, weightOut, texOut);
0247   }
0248 }
0249
0250 void PrintMaterialBudgetInfo::printInfo(
0251     G4VPhysicalVolume* pv, G4LogicalVolume* lv, unsigned int leafDepth, std::ostream& weightOut, std::ostream& texOut) {
0252   double density = lv->GetMaterial()->GetDensity();
0253   double weight = lv->GetMass(false, false);
0254
0255   std::string volumeName = lv->GetName();
0256   if (volumeName.size() < 8)
0257     volumeName.append("\t");
0258
0259   std::string solidName = lv->GetSolid()->GetName();
0260   if (solidName.size() < 8)
0261     solidName.append("\t");
0262
0263   std::string materialName = lv->GetMaterial()->GetName();
0264   if (materialName.size() < 8)
0265     materialName.append("\t");
0266
0267   //----- dump info
0268   weightOut << leafDepth << "\t" << volumeName << "\t" << pv->GetCopyNo() << "\t" << solidName << "\t" << materialName
0269             << "\t" << G4BestUnit(density, "Volumic Mass") << "\t" << G4BestUnit(weight, "Mass") << "\t" << G4endl;
0270   //
0271   texOut << "\t" << leafDepth << "\t & " << stringLaTeXUnderscore(volumeName) << "\t & " << pv->GetCopyNo() << "\t & "
0272          << stringLaTeXUnderscore(solidName) << "\t & " << stringLaTeXUnderscore(materialName) << "\t & "
0273          << stringLaTeXSuperscript(G4BestUnit(density, "Volumic Mass")) << "\t & "
0274          << stringLaTeXSuperscript(G4BestUnit(weight, "Mass")) << "\t \\\\ " << G4endl;
0275   //
0276   for (unsigned int iElement = 0; iElement < (unsigned int)lv->GetMaterial()->GetNumberOfElements(); iElement++) {
0277     // exclude Air in element weight fraction computation
0278     if (materialName.find("Air")) {
0279       std::string elementName = lv->GetMaterial()->GetElement(iElement)->GetName();
0280       double elementMassFraction = lv->GetMaterial()->GetFractionVector()[iElement];
0281       double elementWeight = weight * elementMassFraction;
0282       unsigned int elementIndex = (unsigned int)lv->GetMaterial()->GetElement(iElement)->GetIndex();
0283       elementNames[elementIndex] = elementName;
0284       elementTotalWeight[elementIndex] += elementWeight;
0285     }
0286   }
0287 }
0288
0289 void PrintMaterialBudgetInfo::dumpElementMassFraction(std::ostream& elementOut) {
0290   // calculate mass fraction
0291   double totalWeight = 0.0;
0292   double totalFraction = 0.0;
0293   for (unsigned int iElement = 0; iElement < (unsigned int)elementTotalWeight.size(); iElement++) {
0294     totalWeight += elementTotalWeight[iElement];
0295   }
0296   // calculate element mass fractions
0297   for (unsigned int iElement = 0; iElement < (unsigned int)elementTotalWeight.size(); iElement++) {
0298     elementWeightFraction[iElement] = elementTotalWeight[iElement] / totalWeight;
0299     totalFraction += elementWeightFraction[iElement];
0300   }
0301   // header
0302   elementOut << "Element"
0303              << "\t\t"
0304              << "Index"
0305              << "\t"
0306              << "Total Mass"
0307              << "\t"
0308              << "Mass Fraction "
0309              << "\t" << G4endl;
0310   // dump
0311   for (unsigned int iElement = 0; iElement < (unsigned int)elementTotalWeight.size(); iElement++) {
0312     if (elementNames[iElement] != "rr") {
0313       if (elementNames[iElement].size() < 8)
0314         elementNames[iElement].append("\t");
0315       elementOut << elementNames[iElement] << "\t" << iElement << "\t"
0316                  << G4BestUnit(elementTotalWeight[iElement], "Mass") << "\t" << elementWeightFraction[iElement]
0317                  << G4endl;
0318     }
0319   }
0320   elementOut << "\n\t\tTotal Weight without Air " << G4BestUnit(totalWeight, "Mass") << "\tTotal Fraction "
0321              << totalFraction << G4endl;
0322 }
0323
0324 std::string PrintMaterialBudgetInfo::stringLaTeXUnderscore(std::string stringname) {
0325   // To replace '\' with '\_' to compile LaTeX output
0326   std::string stringoutput;
0327
0328   for (unsigned int i = 0; i < stringname.length(); i++) {
0329     if (stringname.substr(i, 1) == "_") {
0330       stringoutput += "\\_";
0331     } else {
0332       stringoutput += stringname.substr(i, 1);
0333     }
0334   }
0335
0336   return stringoutput;
0337 }
0338
0339 std::string PrintMaterialBudgetInfo::stringLaTeXSuperscript(std::string stringname) {
0340   // To replace 'm3' with 'm$^3$' to compile LaTeX output
0341   std::string stringoutput = stringname.substr(0, 1);
0342
0343   for (unsigned int i = 1; i < stringname.length(); i++) {
0344     if (stringname.substr(i - 1, 1) == "m" && stringname.substr(i, 1) == "3") {
0345       stringoutput += "$^3$";
0346     } else {
0347       stringoutput += stringname.substr(i, 1);
0348     }
0349   }
0350
0351   return stringoutput;
0352 }
0353
0354 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0355 #include "FWCore/PluginManager/interface/ModuleDef.h"
0356
0357 DEFINE_SIMWATCHER(PrintMaterialBudgetInfo);