CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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 
183  const G4RegionStore* regStore = G4RegionStore::GetInstance();
184 
185  G4LogicalVolume* lv;
186  const G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
187  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
188  unsigned int numPV = pvs->size();
189  unsigned int numLV = lvs->size();
190  unsigned int nn = nodeNames.size();
191 
192  std::vector<G4String> savedgdml;
193 
194  fout << "====================================================================="
195  << "\n";
196  fout << "CMSG4OverlapCheck is initialised with " << nodeNames.size() << " nodes; "
197  << " nPoints= " << nPoints << "; tolerance= " << tolerance / mm << " mm; verbose: " << verbose << "\n"
198  << " RegionFlag: " << regionFlag << " PVname: " << PVname << " LVname: " << LVname << "\n"
199  << " Nlv= " << numLV << " Npv= " << numPV << "\n";
200  fout << "====================================================================="
201  << "\n";
202 
203  if (0 < nn) {
204  for (unsigned int ii = 0; ii < nn; ++ii) {
205  if (nodeNames[ii].empty() || "world" == nodeNames[ii] || "World" == nodeNames[ii]) {
206  nodeNames[ii] = "DDDWorld";
207  fout << "### Check overlaps for DDDWorld "
208  << "\n";
209  G4VPhysicalVolume* pv = pvs->GetVolume("DDDWorld");
210  G4GeomTestVolume test(pv, tolerance, nPoints, verbose);
211  test.SetErrorsThreshold(nPrints);
212  test.TestOverlapInTree();
213  } else if (regionFlag) {
214  fout << "---------------------------------------------------------------"
215  << "\n";
216  fout << "### Check overlaps for G4Region Node[" << ii << "] : " << nodeNames[ii] << "\n";
217  G4Region* reg = regStore->GetRegion((G4String)nodeNames[ii]);
218  if (!reg) {
219  fout << "### NO G4Region found - EXIT"
220  << "\n";
221  return;
222  }
223  std::vector<G4LogicalVolume*>::iterator rootLVItr = reg->GetRootLogicalVolumeIterator();
224  unsigned int numRootLV = reg->GetNumberOfRootVolumes();
225  fout << " " << numRootLV << " Root Logical Volumes in this region"
226  << "\n";
227 
228  for (unsigned int iLV = 0; iLV < numRootLV; ++iLV, ++rootLVItr) {
229  // Cover each root logical volume in this region
230  lv = *rootLVItr;
231  fout << "### Check overlaps for G4LogicalVolume " << lv->GetName() << "\n";
232  for (unsigned int i = 0; i < numPV; ++i) {
233  if (((*pvs)[i])->GetLogicalVolume() == lv) {
234  G4String pvname = ((*pvs)[i])->GetName();
235  G4bool isNew = true;
236  for (unsigned int k = 0; k < savedgdml.size(); ++k) {
237  if (pvname == savedgdml[k]) {
238  isNew = false;
239  break;
240  }
241  }
242  if (!isNew) {
243  fout << "### Check overlaps for PhysVolume " << pvname << " is skipted because was already done"
244  << "\n";
245  continue;
246  }
247  savedgdml.push_back(pvname);
248  fout << "### Check overlaps for PhysVolume " << pvname << "\n";
249  // gdml dump only for 1 volume
250  if (gdmlFlag) {
251  G4GDMLParser gdml;
252  gdml.Write(pvname + ".gdml", (*pvs)[i], true);
253  }
254  G4GeomTestVolume test(((*pvs)[i]), tolerance, nPoints, verbose);
255  test.SetErrorsThreshold(nPrints);
256  test.TestOverlapInTree();
257  }
258  }
259  }
260  } else {
261  fout << "### Check overlaps for PhysVolume Node[" << ii << "] : " << nodeNames[ii] << "\n";
262  G4VPhysicalVolume* pv = pvs->GetVolume((G4String)nodeNames[ii]);
263  G4GeomTestVolume test(pv, tolerance, nPoints, verbose);
264  test.SetErrorsThreshold(nPrints);
265  test.TestOverlapInTree();
266  }
267  }
268  }
269  if (!PVname.empty()) {
270  fout << "----------- List of PhysVolumes by name -----------------"
271  << "\n";
272  for (unsigned int i = 0; i < numPV; ++i) {
273  if (PVname == ((*pvs)[i])->GetName()) {
274  fout << " ##### PhysVolume " << PVname << " [" << ((*pvs)[i])->GetCopyNo()
275  << "] LV: " << ((*pvs)[i])->GetLogicalVolume()->GetName()
276  << " Mother LV: " << ((*pvs)[i])->GetMotherLogical()->GetName()
277  << " Region: " << ((*pvs)[i])->GetLogicalVolume()->GetRegion()->GetName() << "\n";
278  fout << " Translation: " << ((*pvs)[i])->GetObjectTranslation() << "\n";
279  fout << " Rotation: " << ((*pvs)[i])->GetObjectRotationValue() << "\n";
280  if (gdmlFlag) {
281  G4GDMLParser gdml;
282  gdml.Write(PVname + ".gdml", (*pvs)[i], true);
283  }
284  }
285  }
286  }
287  if (!LVname.empty()) {
288  fout << "---------- List of Logical Volumes by name ------------------"
289  << "\n";
290  for (unsigned int i = 0; i < numLV; ++i) {
291  if (LVname == ((*lvs)[i])->GetName()) {
292  G4int np = ((*lvs)[i])->GetNoDaughters();
293  fout << " ##### LogVolume " << LVname << " " << np << " daughters"
294  << " Region: " << ((*lvs)[i])->GetRegion()->GetName() << "\n";
295  fout << *(((*lvs)[i])->GetSolid()) << "\n";
296  for (G4int j = 0; j < np; ++j) {
297  G4VPhysicalVolume* pv = ((*lvs)[i])->GetDaughter(j);
298  if (pv) {
299  fout << " PV: " << pv->GetName() << " [" << pv->GetCopyNo() << "]"
300  << " type: " << pv->VolumeType() << " multiplicity: " << pv->GetMultiplicity()
301  << " LV: " << pv->GetLogicalVolume()->GetName() << "\n";
302  fout << " Translation: " << pv->GetObjectTranslation() << "\n";
303  fout << " Rotation: " << pv->GetObjectRotationValue() << "\n";
304  fout << *(pv->GetLogicalVolume()->GetSolid()) << "\n";
305  }
306  }
307  }
308  }
309  }
310  fout << "---------------- End of overlap checks ---------------------"
311  << "\n";
312 }
CMSG4CheckOverlap::makeReportForGeometry
void makeReportForGeometry(std::ofstream &fout, G4VPhysicalVolume *world)
Definition: CMSG4CheckOverlap.cc:153
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
run_AlCaRecoTriggerBitsUpdateWorkflow.session
tuple session
Definition: run_AlCaRecoTriggerBitsUpdateWorkflow.py:57
mps_fire.i
i
Definition: mps_fire.py:428
tolerance
const double tolerance
Definition: HGCalGeomParameters.cc:29
cuy.ii
int ii
Definition: cuy.py:589
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
CustomUIsession::sendToFile
void sendToFile(std::ofstream *)
Definition: CustomUIsession.cc:27
CMSG4CheckOverlap::makeReportForMaterials
void makeReportForMaterials(std::ofstream &fout)
Definition: CMSG4CheckOverlap.cc:96
np
int np
Definition: AMPTWrapper.h:43
verbose
static constexpr int verbose
Definition: HLTExoticaSubAnalysis.cc:25
CMSG4CheckOverlap::makeReportForOverlaps
void makeReportForOverlaps(std::ofstream &fout, const edm::ParameterSet &p)
Definition: CMSG4CheckOverlap.cc:172
ctpps_dqm_sourceclient-live_cfg.test
test
Definition: ctpps_dqm_sourceclient-live_cfg.py:7
CMSG4RegionReporter::ReportRegions
void ReportRegions(const std::string &ss)
Definition: CMSG4RegionReporter.cc:18
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
gpuVertexFinder::nn
int32_t *__restrict__ nn
Definition: gpuClusterTracksDBSCAN.h:41
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
CMSG4CheckOverlap::CMSG4CheckOverlap
CMSG4CheckOverlap(edm::ParameterSet const &p, std::string &regFile, CustomUIsession *, G4VPhysicalVolume *world)
Definition: CMSG4CheckOverlap.cc:32
Generated for CMSSW Reference Manual by doxygen 1.8.5