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