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
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
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
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 }
