Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:24:54

0001 #include "SimG4Core/Geometry/interface/CMSG4CheckOverlap.h"
0002 #include "SimG4Core/Geometry/interface/CMSG4RegionReporter.h"
0003 #include "SimG4Core/Geometry/interface/CustomUIsession.h"
0004 
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 
0007 #include "G4GDMLParser.hh"
0008 #include "G4GeomTestVolume.hh"
0009 #include "G4LogicalVolume.hh"
0010 #include "G4LogicalVolumeStore.hh"
0011 #include "G4PhysicalVolumeStore.hh"
0012 #include "G4GeometryManager.hh"
0013 #include "G4Region.hh"
0014 #include "G4RegionStore.hh"
0015 #include "G4Element.hh"
0016 #include "G4ElementTable.hh"
0017 #include "G4Material.hh"
0018 #include "G4MaterialTable.hh"
0019 #include "G4ProductionCutsTable.hh"
0020 #include "G4MaterialCutsCouple.hh"
0021 #include "G4SystemOfUnits.hh"
0022 #include "G4VPhysicalVolume.hh"
0023 #include "G4UnitsTable.hh"
0024 #include "G4ios.hh"
0025 #include "globals.hh"
0026 
0027 #include <string>
0028 #include <vector>
0029 #include <iostream>
0030 #include <iomanip>
0031 
0032 CMSG4CheckOverlap::CMSG4CheckOverlap(const edm::ParameterSet& p,
0033                                      std::string& regFile,
0034                                      CustomUIsession* session,
0035                                      G4VPhysicalVolume* world) {
0036   bool mat = p.getParameter<bool>("MaterialFlag");
0037   const std::string& ss = p.getParameter<std::string>("OutputBaseName");
0038   if (ss.empty()) {
0039     edm::LogWarning("SimG4CoreGeometry")
0040         << "CMSG4CheckOverlap: OutputFileBaseName is not provided - no check is performed";
0041     return;
0042   }
0043   if (mat) {
0044     const std::string sss = "Materials_" + ss + ".txt";
0045     std::ofstream fout(sss.c_str(), std::ios::out);
0046     if (fout.fail()) {
0047       edm::LogWarning("SimG4CoreGeometry")
0048           << "CMSG4CheckOverlap: file <" << sss << "> is not opened - no report provided";
0049     } else {
0050       edm::LogVerbatim("SimG4CoreGeometry") << "CMSG4CheckOverlap: output file <" << sss << "> is opened";
0051       makeReportForMaterials(fout);
0052       fout.close();
0053     }
0054   }
0055 
0056   bool reg = p.getParameter<bool>("RegionFlag");
0057   if (reg) {
0058     const std::string qqq = (regFile.empty()) ? ss : regFile;
0059     const std::string sss = "Regions_" + qqq + ".txt";
0060     CMSG4RegionReporter rrep;
0061     rrep.ReportRegions(sss);
0062   }
0063 
0064   bool geom = p.getParameter<bool>("GeomFlag");
0065   if (geom) {
0066     const std::string sss = "Geometry_" + ss + ".txt";
0067     std::ofstream fout(sss.c_str(), std::ios::out);
0068     if (fout.fail()) {
0069       edm::LogWarning("SimG4CoreGeometry")
0070           << "CMSG4CheckOverlap: file <" << sss << "> is not opened - no report provided";
0071     } else {
0072       edm::LogVerbatim("SimG4CoreGeometry") << "CMSG4CheckOverlap: output file <" << sss << "> is opened";
0073       session->sendToFile(&fout);
0074       makeReportForGeometry(fout, world);
0075       session->stopSendToFile();
0076       fout.close();
0077     }
0078   }
0079   bool oFlag = p.getParameter<bool>("OverlapFlag");
0080   if (oFlag) {
0081     const std::string sss = "Overlaps_" + ss + ".txt";
0082     std::ofstream fout(sss.c_str(), std::ios::out);
0083     if (fout.fail()) {
0084       edm::LogWarning("SimG4CoreGeometry")
0085           << "CMSG4CheckOverlap: file <" << sss << "> is not opened - no report provided";
0086     } else {
0087       edm::LogVerbatim("SimG4CoreGeometry") << "CMSG4CheckOverlap: output file <" << ss << "> is opened";
0088       session->sendToFile(&fout);
0089       makeReportForOverlaps(fout, p);
0090       session->stopSendToFile();
0091       fout.close();
0092     }
0093   }
0094 }
0095 
0096 void CMSG4CheckOverlap::makeReportForMaterials(std::ofstream& fout) {
0097   int nelm = G4Element::GetNumberOfElements();
0098   int nmat = G4Material::GetNumberOfMaterials();
0099   G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
0100   int ncouples = theCoupleTable->GetTableSize();
0101   fout << "====================================================================="
0102        << "\n";
0103   fout << "NumberOfElements  " << nelm << "\n";
0104   fout << "NumberOfMaterials " << nmat << "\n";
0105   fout << "NumberOfCouples   " << ncouples << "\n";
0106   fout << "====================================================================="
0107        << "\n";
0108   fout << "ElementsDump:"
0109        << "\n";
0110   G4ElementTable* elmtab = G4Element::GetElementTable();
0111   fout << *elmtab;
0112   fout << "====================================================================="
0113        << "\n";
0114   G4MaterialTable* mattab = G4Material::GetMaterialTable();
0115   fout << "MaterialsDump:"
0116        << "\n";
0117   //fout << *mattab << "\n";
0118   for (int i = 0; i < nmat; ++i) {
0119     fout << "### Material " << i << "   " << ((*mattab)[i])->GetNumberOfElements() << " elements\n";
0120     fout << (*mattab)[i] << "\n";
0121   }
0122   fout << "====================================================================="
0123        << "\n";
0124   fout << "MaterialsCutsCoupleDump:"
0125        << "\n";
0126   const std::vector<G4double>* gcut = theCoupleTable->GetEnergyCutsVector(0);
0127   const std::vector<G4double>* ecut = theCoupleTable->GetEnergyCutsVector(1);
0128   const std::vector<G4double>* pcut = theCoupleTable->GetEnergyCutsVector(2);
0129   const std::vector<G4double>* icut = theCoupleTable->GetEnergyCutsVector(3);
0130   for (int i = 0; i < ncouples; ++i) {
0131     const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
0132     const G4ProductionCuts* aCut = couple->GetProductionCuts();
0133     fout << "Index : " << i << "  used in the geometry : ";
0134     if (couple->IsUsed()) {
0135       fout << "Yes\n";
0136     } else {
0137       fout << "No \n";
0138     }
0139     fout << " Material : " << couple->GetMaterial()->GetName() << "\n";
0140     fout << " Range cuts        : "
0141          << " gamma  " << G4BestUnit(aCut->GetProductionCut(0), "Length") << "    e-  "
0142          << G4BestUnit(aCut->GetProductionCut(1), "Length") << "    e+  "
0143          << G4BestUnit(aCut->GetProductionCut(2), "Length") << " proton "
0144          << G4BestUnit(aCut->GetProductionCut(3), "Length") << "\n";
0145     fout << " Energy thresholds : ";
0146     fout << " gamma  " << G4BestUnit((*gcut)[i], "Energy") << "    e-  " << G4BestUnit((*ecut)[i], "Energy")
0147          << "    e+  " << G4BestUnit((*pcut)[i], "Energy") << " proton " << G4BestUnit((*icut)[i], "Energy") << "\n";
0148   }
0149   fout << "======================================================================"
0150        << "\n";
0151 }
0152 
0153 void CMSG4CheckOverlap::makeReportForGeometry(std::ofstream& fout, G4VPhysicalVolume* world) {
0154   const G4RegionStore* regs = G4RegionStore::GetInstance();
0155   const G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
0156   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0157   int numPV = pvs->size();
0158   int numLV = lvs->size();
0159   int nreg = regs->size();
0160   fout << "====================================================================="
0161        << "\n";
0162   fout << "NumberOfRegions         " << nreg << "\n";
0163   fout << "NumberOfLogicalVolumes  " << numLV << "\n";
0164   fout << "NumberOfPhysicalVolumes " << numPV << "\n";
0165   fout << "====================================================================="
0166        << "\n";
0167   G4GeometryManager::GetInstance()->CloseGeometry(true, true, world);
0168   fout << "====================================================================="
0169        << "\n";
0170 }
0171 
0172 void CMSG4CheckOverlap::makeReportForOverlaps(std::ofstream& fout, const edm::ParameterSet& p) {
0173   std::vector<std::string> nodeNames = p.getParameter<std::vector<std::string>>("NodeNames");
0174   std::string PVname = p.getParameter<std::string>("PVname");
0175   std::string LVname = p.getParameter<std::string>("LVname");
0176   double tolerance = p.getParameter<double>("Tolerance") * CLHEP::mm;
0177   int nPoints = p.getParameter<int>("Resolution");
0178   bool verbose = p.getParameter<bool>("Verbose");
0179   bool regionFlag = p.getParameter<bool>("RegionFlag");
0180   bool gdmlFlag = p.getParameter<bool>("gdmlFlag");
0181   int nPrints = p.getParameter<int>("ErrorThreshold");
0182 
0183   const G4RegionStore* regStore = G4RegionStore::GetInstance();
0184 
0185   G4LogicalVolume* lv;
0186   const G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
0187   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0188   unsigned int numPV = pvs->size();
0189   unsigned int numLV = lvs->size();
0190   unsigned int nn = nodeNames.size();
0191 
0192   std::vector<G4String> savedgdml;
0193 
0194   fout << "====================================================================="
0195        << "\n";
0196   fout << "CMSG4OverlapCheck is initialised with " << nodeNames.size() << " nodes; "
0197        << " nPoints= " << nPoints << "; tolerance= " << tolerance / mm << " mm; verbose: " << verbose << "\n"
0198        << "               RegionFlag: " << regionFlag << "  PVname: " << PVname << "  LVname: " << LVname << "\n"
0199        << "               Nlv= " << numLV << "   Npv= " << numPV << "\n";
0200   fout << "====================================================================="
0201        << "\n";
0202 
0203   if (0 < nn) {
0204     for (unsigned int ii = 0; ii < nn; ++ii) {
0205       if (nodeNames[ii].empty() || "world" == nodeNames[ii] || "World" == nodeNames[ii]) {
0206         nodeNames[ii] = "DDDWorld";
0207         fout << "### Check overlaps for DDDWorld "
0208              << "\n";
0209         G4VPhysicalVolume* pv = pvs->GetVolume("DDDWorld");
0210         G4GeomTestVolume test(pv, tolerance, nPoints, verbose);
0211         test.SetErrorsThreshold(nPrints);
0212         test.TestOverlapInTree();
0213       } else if (regionFlag) {
0214         fout << "---------------------------------------------------------------"
0215              << "\n";
0216         fout << "### Check overlaps for G4Region Node[" << ii << "] : " << nodeNames[ii] << "\n";
0217         G4Region* reg = regStore->GetRegion((G4String)nodeNames[ii]);
0218         if (!reg) {
0219           fout << "### NO G4Region found - EXIT"
0220                << "\n";
0221           return;
0222         }
0223         std::vector<G4LogicalVolume*>::iterator rootLVItr = reg->GetRootLogicalVolumeIterator();
0224         unsigned int numRootLV = reg->GetNumberOfRootVolumes();
0225         fout << "      " << numRootLV << " Root Logical Volumes in this region"
0226              << "\n";
0227 
0228         for (unsigned int iLV = 0; iLV < numRootLV; ++iLV, ++rootLVItr) {
0229           // Cover each root logical volume in this region
0230           lv = *rootLVItr;
0231           fout << "### Check overlaps for G4LogicalVolume " << lv->GetName() << "\n";
0232           for (unsigned int i = 0; i < numPV; ++i) {
0233             if (((*pvs)[i])->GetLogicalVolume() == lv) {
0234               G4String pvname = ((*pvs)[i])->GetName();
0235               G4bool isNew = true;
0236               for (unsigned int k = 0; k < savedgdml.size(); ++k) {
0237                 if (pvname == savedgdml[k]) {
0238                   isNew = false;
0239                   break;
0240                 }
0241               }
0242               if (!isNew) {
0243                 fout << "### Check overlaps for PhysVolume  " << pvname << " is skipted because was already done"
0244                      << "\n";
0245                 continue;
0246               }
0247               savedgdml.push_back(pvname);
0248               fout << "### Check overlaps for PhysVolume  " << pvname << "\n";
0249               // gdml dump only for 1 volume
0250               if (gdmlFlag) {
0251                 G4GDMLParser gdml;
0252                 gdml.Write(pvname + ".gdml", (*pvs)[i], true);
0253               }
0254               G4GeomTestVolume test(((*pvs)[i]), tolerance, nPoints, verbose);
0255               test.SetErrorsThreshold(nPrints);
0256               test.TestOverlapInTree();
0257             }
0258           }
0259         }
0260       } else {
0261         fout << "### Check overlaps for PhysVolume Node[" << ii << "] : " << nodeNames[ii] << "\n";
0262         G4VPhysicalVolume* pv = pvs->GetVolume((G4String)nodeNames[ii]);
0263         G4GeomTestVolume test(pv, tolerance, nPoints, verbose);
0264         test.SetErrorsThreshold(nPrints);
0265         test.TestOverlapInTree();
0266       }
0267     }
0268   }
0269   if (!PVname.empty()) {
0270     fout << "----------- List of PhysVolumes by name -----------------"
0271          << "\n";
0272     for (unsigned int i = 0; i < numPV; ++i) {
0273       if (PVname == ((*pvs)[i])->GetName()) {
0274         fout << " ##### PhysVolume " << PVname << " [" << ((*pvs)[i])->GetCopyNo()
0275              << "]  LV: " << ((*pvs)[i])->GetLogicalVolume()->GetName()
0276              << " Mother LV: " << ((*pvs)[i])->GetMotherLogical()->GetName()
0277              << " Region: " << ((*pvs)[i])->GetLogicalVolume()->GetRegion()->GetName() << "\n";
0278         fout << "       Translation: " << ((*pvs)[i])->GetObjectTranslation() << "\n";
0279         fout << "       Rotation:    " << ((*pvs)[i])->GetObjectRotationValue() << "\n";
0280         if (gdmlFlag) {
0281           G4GDMLParser gdml;
0282           gdml.Write(PVname + ".gdml", (*pvs)[i], true);
0283         }
0284       }
0285     }
0286   }
0287   if (!LVname.empty()) {
0288     fout << "---------- List of Logical Volumes by name ------------------"
0289          << "\n";
0290     for (unsigned int i = 0; i < numLV; ++i) {
0291       if (LVname == ((*lvs)[i])->GetName()) {
0292         G4int np = ((*lvs)[i])->GetNoDaughters();
0293         fout << " ##### LogVolume " << LVname << "  " << np << " daughters"
0294              << " Region: " << ((*lvs)[i])->GetRegion()->GetName() << "\n";
0295         fout << *(((*lvs)[i])->GetSolid()) << "\n";
0296         for (G4int j = 0; j < np; ++j) {
0297           G4VPhysicalVolume* pv = ((*lvs)[i])->GetDaughter(j);
0298           if (pv) {
0299             fout << "   PV: " << pv->GetName() << " [" << pv->GetCopyNo() << "]"
0300                  << " type: " << pv->VolumeType() << "  multiplicity: " << pv->GetMultiplicity()
0301                  << " LV: " << pv->GetLogicalVolume()->GetName() << "\n";
0302             fout << "       Translation: " << pv->GetObjectTranslation() << "\n";
0303             fout << "       Rotation:    " << pv->GetObjectRotationValue() << "\n";
0304             fout << *(pv->GetLogicalVolume()->GetSolid()) << "\n";
0305           }
0306         }
0307       }
0308     }
0309   }
0310   fout << "---------------- End of overlap checks ---------------------"
0311        << "\n";
0312 }