CMS 3D CMS Logo
CMSG4CheckOverlap.cc
Go to the documentation of this file.
1 #include "SimG4Core/Geometry/interface/CMSG4CheckOverlap.h"
2 #include "SimG4Core/Geometry/interface/CMSG4RegionReporter.h"
3 #include "SimG4Core/Geometry/interface/CustomUIsession.h"
4 
5 #include "FWCore/MessageLogger/interface/MessageLogger.h"
6 
7 #include "G4GDMLParser.hh"
8 #include "G4GeomTestVolume.hh"
9 #include "G4LogicalVolume.hh"
10 #include "G4LogicalVolumeStore.hh"
11 #include "G4PhysicalVolumeStore.hh"
12 #include "G4GeometryManager.hh"
13 #include "G4Region.hh"
14 #include "G4RegionStore.hh"
15 #include "G4Element.hh"
16 #include "G4ElementTable.hh"
17 #include "G4Material.hh"
18 #include "G4MaterialTable.hh"
19 #include "G4ProductionCutsTable.hh"
20 #include "G4MaterialCutsCouple.hh"
21 #include "G4SystemOfUnits.hh"
22 #include "G4VPhysicalVolume.hh"
23 #include "G4UnitsTable.hh"
24 #include "G4ios.hh"
25 #include "globals.hh"
26 
27 #include <string>
28 #include <vector>
29 #include <iostream>
30 #include <iomanip>
31 
32 CMSG4CheckOverlap::CMSG4CheckOverlap(const edm::ParameterSet& p,
33  std::string& regFile,
34  CustomUIsession* session,
35  G4VPhysicalVolume* world) {
36  bool mat = p.getParameter<bool>("MaterialFlag");
37  const std::string& ss = p.getParameter<std::string>("OutputBaseName");
38  if (ss.empty()) {
39  edm::LogWarning("SimG4CoreGeometry")
40  << "CMSG4CheckOverlap: OutputFileBaseName is not provided - no check is performed";
41  return;
42  }
43  if (mat) {
44  const std::string sss = "Materials_" + ss + ".txt";
45  std::ofstream fout(sss.c_str(), std::ios::out);
46  if (fout.fail()) {
47  edm::LogWarning("SimG4CoreGeometry")
48  << "CMSG4CheckOverlap: file <" << sss << "> is not opened - no report provided";
49  } else {
50  edm::LogVerbatim("SimG4CoreGeometry") << "CMSG4CheckOverlap: output file <" << sss << "> is opened";
51  makeReportForMaterials(fout);
52  fout.close();
53  }
54  }
55 
56  bool reg = p.getParameter<bool>("RegionFlag");
57  if (reg) {
58  const std::string qqq = (regFile.empty()) ? ss : regFile;
59  const std::string sss = "Regions_" + qqq + ".txt";
60  CMSG4RegionReporter rrep;
61  rrep.ReportRegions(sss);
62  }
63 
64  bool geom = p.getParameter<bool>("GeomFlag");
65  if (geom) {
66  const std::string sss = "Geometry_" + ss + ".txt";
67  std::ofstream fout(sss.c_str(), std::ios::out);
68  if (fout.fail()) {
69  edm::LogWarning("SimG4CoreGeometry")
70  << "CMSG4CheckOverlap: file <" << sss << "> is not opened - no report provided";
71  } else {
72  edm::LogVerbatim("SimG4CoreGeometry") << "CMSG4CheckOverlap: output file <" << sss << "> is opened";
73  session->sendToFile(&fout);
74  makeReportForGeometry(fout, world);
75  session->stopSendToFile();
76  fout.close();
77  }
78  }
79  bool oFlag = p.getParameter<bool>("OverlapFlag");
80  if (oFlag) {
81  const std::string sss = "Overlaps_" + ss + ".txt";
82  std::ofstream fout(sss.c_str(), std::ios::out);
83  if (fout.fail()) {
84  edm::LogWarning("SimG4CoreGeometry")
85  << "CMSG4CheckOverlap: file <" << sss << "> is not opened - no report provided";
86  } else {
87  edm::LogVerbatim("SimG4CoreGeometry") << "CMSG4CheckOverlap: output file <" << ss << "> is opened";
88  session->sendToFile(&fout);
89  makeReportForOverlaps(fout, p);
90  session->stopSendToFile();
91  fout.close();
92  }
93  }
94 }
95 
96 void CMSG4CheckOverlap::makeReportForMaterials(std::ofstream& fout) {
97  int nelm = G4Element::GetNumberOfElements();
98  int nmat = G4Material::GetNumberOfMaterials();
99  G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
100  int ncouples = theCoupleTable->GetTableSize();
101  fout << "====================================================================="
102  << "\n";
103  fout << "NumberOfElements " << nelm << "\n";
104  fout << "NumberOfMaterials " << nmat << "\n";
105  fout << "NumberOfCouples " << ncouples << "\n";
106  fout << "====================================================================="
107  << "\n";
108  fout << "ElementsDump:"
109  << "\n";
110  G4ElementTable* elmtab = G4Element::GetElementTable();
111  fout << *elmtab;
112  fout << "====================================================================="
113  << "\n";
114  G4MaterialTable* mattab = G4Material::GetMaterialTable();
115  fout << "MaterialsDump:"
116  << "\n";
117  //fout << *mattab << "\n";
118  for (int i = 0; i < nmat; ++i) {
119  fout << "### Material " << i << " " << ((*mattab)[i])->GetNumberOfElements() << " elements\n";
120  fout << (*mattab)[i] << "\n";
121  }
122  fout << "====================================================================="
123  << "\n";
124  fout << "MaterialsCutsCoupleDump:"
125  << "\n";
126  const std::vector<G4double>* gcut = theCoupleTable->GetEnergyCutsVector(0);
127  const std::vector<G4double>* ecut = theCoupleTable->GetEnergyCutsVector(1);
128  const std::vector<G4double>* pcut = theCoupleTable->GetEnergyCutsVector(2);
129  const std::vector<G4double>* icut = theCoupleTable->GetEnergyCutsVector(3);
130  for (int i = 0; i < ncouples; ++i) {
131  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
132  const G4ProductionCuts* aCut = couple->GetProductionCuts();
133  fout << "Index : " << i << " used in the geometry : ";
134  if (couple->IsUsed()) {
135  fout << "Yes\n";
136  } else {
137  fout << "No \n";
138  }
139  fout << " Material : " << couple->GetMaterial()->GetName() << "\n";
140  fout << " Range cuts : "
141  << " gamma " << G4BestUnit(aCut->GetProductionCut(0), "Length") << " e- "
142  << G4BestUnit(aCut->GetProductionCut(1), "Length") << " e+ "
143  << G4BestUnit(aCut->GetProductionCut(2), "Length") << " proton "
144  << G4BestUnit(aCut->GetProductionCut(3), "Length") << "\n";
145  fout << " Energy thresholds : ";
146  fout << " gamma " << G4BestUnit((*gcut)[i], "Energy") << " e- " << G4BestUnit((*ecut)[i], "Energy")
147  << " e+ " << G4BestUnit((*pcut)[i], "Energy") << " proton " << G4BestUnit((*icut)[i], "Energy") << "\n";
148  }
149  fout << "======================================================================"
150  << "\n";
151 }
152 
153 void CMSG4CheckOverlap::makeReportForGeometry(std::ofstream& fout, G4VPhysicalVolume* world) {
154  const G4RegionStore* regs = G4RegionStore::GetInstance();
155  const G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
156  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
157  int numPV = pvs->size();
158  int numLV = lvs->size();
159  int nreg = regs->size();
160  fout << "====================================================================="
161  << "\n";
162  fout << "NumberOfRegions " << nreg << "\n";
163  fout << "NumberOfLogicalVolumes " << numLV << "\n";
164  fout << "NumberOfPhysicalVolumes " << numPV << "\n";
165  fout << "====================================================================="
166  << "\n";
167  G4GeometryManager::GetInstance()->CloseGeometry(true, true, world);
168  fout << "====================================================================="
169  << "\n";
170 }
171 
172 void CMSG4CheckOverlap::makeReportForOverlaps(std::ofstream& fout, const edm::ParameterSet& p) {
173  std::vector<std::string> nodeNames = p.getParameter<std::vector<std::string>>("NodeNames");
174  std::string PVname = p.getParameter<std::string>("PVname");
175  std::string LVname = p.getParameter<std::string>("LVname");
176  double tolerance = p.getParameter<double>("Tolerance") * CLHEP::mm;
177  int nPoints = p.getParameter<int>("Resolution");
178  bool verbose = p.getParameter<bool>("Verbose");
179  bool regionFlag = p.getParameter<bool>("RegionFlag");
180  bool gdmlFlag = p.getParameter<bool>("gdmlFlag");
181  int nPrints = p.getParameter<int>("ErrorThreshold");
182  int level = p.getParameter<int>("Level");
183  int depth = p.getParameter<int>("Depth");
184 
185  const G4RegionStore* regStore = G4RegionStore::GetInstance();
186 
187  G4LogicalVolume* lv;
188  const G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
189  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
190  unsigned int numPV = pvs->size();
191  unsigned int numLV = lvs->size();
192  unsigned int nn = nodeNames.size();
193 
194  std::vector<G4String> savedgdml;
195 
196  fout << "====================================================================="
197  << "\n";
198  fout << "CMSG4OverlapCheck is initialised with " << nodeNames.size() << " nodes; "
199  << " nPoints= " << nPoints << "; tolerance= " << tolerance / mm << " mm; verbose: " << verbose << "\n"
200  << " RegionFlag: " << regionFlag << " PVname: " << PVname << " LVname: " << LVname << "\n"
201  << " Nlv= " << numLV << " Npv= " << numPV << "\n";
202  fout << "====================================================================="
203  << "\n";
204 
205  if (0 < nn) {
206  for (unsigned int ii = 0; ii < nn; ++ii) {
207  if (nodeNames[ii].empty() || "world" == nodeNames[ii] || "World" == nodeNames[ii]) {
208  nodeNames[ii] = "DDDWorld";
209  fout << "### Check overlaps for DDDWorld "
210  << "\n";
211  G4VPhysicalVolume* pv = pvs->GetVolume("DDDWorld");
212  G4GeomTestVolume test(pv, tolerance, nPoints, verbose);
213  test.SetErrorsThreshold(nPrints);
214  test.TestRecursiveOverlap(level, depth);
215  } else if (regionFlag) {
216  fout << "---------------------------------------------------------------"
217  << "\n";
218  fout << "### Check overlaps for G4Region Node[" << ii << "] : " << nodeNames[ii] << "\n";
219  G4Region* reg = regStore->GetRegion((G4String)nodeNames[ii]);
220  if (!reg) {
221  fout << "### NO G4Region found - EXIT"
222  << "\n";
223  return;
224  }
225  std::vector<G4LogicalVolume*>::iterator rootLVItr = reg->GetRootLogicalVolumeIterator();
226  unsigned int numRootLV = reg->GetNumberOfRootVolumes();
227  fout << " " << numRootLV << " Root Logical Volumes in this region"
228  << "\n";
229 
230  for (unsigned int iLV = 0; iLV < numRootLV; ++iLV, ++rootLVItr) {
231  // Cover each root logical volume in this region
232  lv = *rootLVItr;
233  fout << "### Check overlaps for G4LogicalVolume " << lv->GetName() << "\n";
234  for (unsigned int i = 0; i < numPV; ++i) {
235  if (((*pvs)[i])->GetLogicalVolume() == lv) {
236  G4String pvname = ((*pvs)[i])->GetName();
237  G4bool isNew = true;
238  for (unsigned int k = 0; k < savedgdml.size(); ++k) {
239  if (pvname == savedgdml[k]) {
240  isNew = false;
241  break;
242  }
243  }
244  if (!isNew) {
245  fout << "### Check overlaps for PhysVolume " << pvname << " is skipted because was already done"
246  << "\n";
247  continue;
248  }
249  savedgdml.push_back(pvname);
250  fout << "### Check overlaps for PhysVolume " << pvname << "\n";
251  // gdml dump only for 1 volume
252  if (gdmlFlag) {
253  G4GDMLParser gdml;
254  gdml.Write(pvname + ".gdml", (*pvs)[i], true);
255  }
256  G4GeomTestVolume test(((*pvs)[i]), tolerance, nPoints, verbose);
257  test.SetErrorsThreshold(nPrints);
258  test.TestRecursiveOverlap(level, depth);
259  }
260  }
261  }
262  } else {
263  fout << "### Check overlaps for PhysVolume Node[" << ii << "] : " << nodeNames[ii] << "\n";
264  G4VPhysicalVolume* pv = pvs->GetVolume((G4String)nodeNames[ii]);
265  G4GeomTestVolume test(pv, tolerance, nPoints, verbose);
266  test.SetErrorsThreshold(nPrints);
267  test.TestRecursiveOverlap(level, depth);
268  }
269  }
270  }
271  if (!PVname.empty()) {
272  fout << "----------- List of PhysVolumes by name -----------------"
273  << "\n";
274  for (unsigned int i = 0; i < numPV; ++i) {
275  if (PVname == ((*pvs)[i])->GetName()) {
276  fout << " ##### PhysVolume " << PVname << " [" << ((*pvs)[i])->GetCopyNo()
277  << "] LV: " << ((*pvs)[i])->GetLogicalVolume()->GetName()
278  << " Mother LV: " << ((*pvs)[i])->GetMotherLogical()->GetName()
279  << " Region: " << ((*pvs)[i])->GetLogicalVolume()->GetRegion()->GetName() << "\n";
280  fout << " Translation: " << ((*pvs)[i])->GetObjectTranslation() << "\n";
281  fout << " Rotation: " << ((*pvs)[i])->GetObjectRotationValue() << "\n";
282  if (gdmlFlag) {
283  G4GDMLParser gdml;
284  gdml.Write(PVname + ".gdml", (*pvs)[i], true);
285  }
286  }
287  }
288  }
289  if (!LVname.empty()) {
290  fout << "---------- List of Logical Volumes by name ------------------"
291  << "\n";
292  for (unsigned int i = 0; i < numLV; ++i) {
293  if (LVname == ((*lvs)[i])->GetName()) {
294  G4int np = ((*lvs)[i])->GetNoDaughters();
295  fout << " ##### LogVolume " << LVname << " " << np << " daughters"
296  << " Region: " << ((*lvs)[i])->GetRegion()->GetName() << "\n";
297  fout << *(((*lvs)[i])->GetSolid()) << "\n";
298  for (G4int j = 0; j < np; ++j) {
299  G4VPhysicalVolume* pv = ((*lvs)[i])->GetDaughter(j);
300  if (pv) {
301  fout << " PV: " << pv->GetName() << " [" << pv->GetCopyNo() << "]"
302  << " type: " << pv->VolumeType() << " multiplicity: " << pv->GetMultiplicity()
303  << " LV: " << pv->GetLogicalVolume()->GetName() << "\n";
304  fout << " Translation: " << pv->GetObjectTranslation() << "\n";
305  fout << " Rotation: " << pv->GetObjectRotationValue() << "\n";
306  fout << *(pv->GetLogicalVolume()->GetSolid()) << "\n";
307  }
308  }
309  }
310  }
311  }
312  fout << "---------------- End of overlap checks ---------------------"
313  << "\n";
314 }
315 
316 CMSG4CheckOverlap::~CMSG4CheckOverlap() {}
personalPlayback.level
level
Definition: personalPlayback.py:22
mps_fire.i
i
Definition: mps_fire.py:355
MessageLogger.h
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
np
int np
Definition: AMPTWrapper.h:43
CMSG4CheckOverlap.h
CMSG4RegionReporter.h
CMSG4CheckOverlap::makeReportForMaterials
void makeReportForMaterials(std::ofstream &fout)
Definition: CMSG4CheckOverlap.cc:96
contentValuesCheck.ss
ss
Definition: contentValuesCheck.py:33
test
Definition: SmallWORMDict.h:13
ctpps_dqm_sourceclient-live_cfg.test
test
Definition: ctpps_dqm_sourceclient-live_cfg.py:7
CMSG4RegionReporter
Definition: CMSG4RegionReporter.h:6
relativeConstraints.geom
geom
Definition: relativeConstraints.py:72
CustomUIsession.h
verbose
static constexpr int verbose
Definition: HLTExoticaSubAnalysis.cc:25
dqmdumpme.k
k
Definition: dqmdumpme.py:60
LEDCalibrationChannels.depth
depth
Definition: LEDCalibrationChannels.py:65
CustomUIsession::stopSendToFile
void stopSendToFile()
Definition: CustomUIsession.h:22
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
g4SimHits_cfi.LVname
LVname
Definition: g4SimHits_cfi.py:69
edm::LogWarning
Definition: MessageLogger.h:141
CMSG4CheckOverlap::makeReportForOverlaps
void makeReportForOverlaps(std::ofstream &fout, const edm::ParameterSet &p)
Definition: CMSG4CheckOverlap.cc:172
edm::ParameterSet
Definition: ParameterSet.h:36
groupFilesInBlocks.fout
fout
Definition: groupFilesInBlocks.py:162
MetAnalyzer.pv
def pv(vc)
Definition: MetAnalyzer.py:7
tolerance
const double tolerance
Definition: HGCalGeomParameters.cc:27
edm::LogVerbatim
Definition: MessageLogger.h:297
CMSG4CheckOverlap::~CMSG4CheckOverlap
~CMSG4CheckOverlap()
Definition: CMSG4CheckOverlap.cc:316
FSQDQM_cfi.pvs
pvs
Definition: FSQDQM_cfi.py:12
g4SimHits_cfi.PVname
PVname
Definition: g4SimHits_cfi.py:68
groupFilesInBlocks.nn
nn
Definition: groupFilesInBlocks.py:150
g4SimHits_cfi.gdmlFlag
gdmlFlag
Definition: g4SimHits_cfi.py:61
CustomUIsession
Definition: CustomUIsession.h:13
relativeConstraints.empty
bool empty
Definition: relativeConstraints.py:46
CustomUIsession::sendToFile
void sendToFile(std::ofstream *)
Definition: CustomUIsession.cc:35
CMSG4CheckOverlap::makeReportForGeometry
void makeReportForGeometry(std::ofstream &fout, G4VPhysicalVolume *world)
Definition: CMSG4CheckOverlap.cc:153
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
cuy.ii
ii
Definition: cuy.py:590
CMSG4CheckOverlap::CMSG4CheckOverlap
CMSG4CheckOverlap(edm::ParameterSet const &p, std::string &regFile, CustomUIsession *, G4VPhysicalVolume *world)
Definition: CMSG4CheckOverlap.cc:32
CMSG4RegionReporter::ReportRegions
void ReportRegions(const std::string &ss)
Definition: CMSG4RegionReporter.cc:18
Generated for CMSSW Reference Manual by doxygen 1.8.16