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