Extended checks on geometry

SimG4Core/Application/interface/RunManagerMT.h

@@ -54,9 +54,12 @@ namespace HepPDT {
  * (acting as the Geant4 master thread), and there should be exactly
  * one instance of it. 
  */
+class RunManagerMTWorker;
 
 class RunManagerMT 
 {
+  friend class RunManagerMTWorker;
+
 public:
   RunManagerMT(edm::ParameterSet const & p);
   ~RunManagerMT();

SimG4Core/Application/python/g4SimHits_cfi.py

@@ -53,6 +53,8 @@
     G4CheckOverlap = cms.PSet(
         Tolerance = cms.untracked.double(0.0),
         Resolution = cms.untracked.int32(10000),
+        RegionFlag = cms.untracked.bool(True),  # if true - selection by G4Region name
+        gdmlFlag = cms.untracked.bool(True),  # if true - dump gdml file
         NodeNames = cms.vstring('World')
     ),
     G4Commands = cms.vstring(),

SimG4Core/Application/src/ParametrisedEMPhysics.cc

@@ -34,6 +34,7 @@
 #include "G4MuonNuclearProcess.hh"
 #include "G4MuonVDNuclearModel.hh"
 
+#include "G4EmParameters.hh"
 #include "G4EmProcessOptions.hh"
 #include "G4PhysicsListHelper.hh"
 #include "G4SystemOfUnits.hh"
@@ -44,10 +45,28 @@ ParametrisedEMPhysics::ParametrisedEMPhysics(std::string name,
 					     const edm::ParameterSet & p) 
   : G4VPhysicsConstructor(name), theParSet(p) 
 {
-  theEcalEMShowerModel = 0;
-  theEcalHadShowerModel = 0;
-  theHcalEMShowerModel = 0;
-  theHcalHadShowerModel = 0;
+  theEcalEMShowerModel = nullptr;
+  theEcalHadShowerModel = nullptr;
+  theHcalEMShowerModel = nullptr;
+  theHcalHadShowerModel = nullptr;
+
+  // bremsstrahlung threshold and EM verbosity
+  G4EmParameters* param = G4EmParameters::Instance();
+  G4int verb = theParSet.getUntrackedParameter<int>("Verbosity",0);
+  param->SetVerbose(verb);
+
+  G4double bremth = theParSet.getParameter<double>("G4BremsstrahlungThreshold")*GeV; 
+  param->SetBremsstrahlungTh(bremth);
+
+  bool fluo = theParSet.getParameter<bool>("FlagFluo");
+  param->SetFluo(fluo);
+
+  edm::LogInfo("SimG4CoreApplication") 
+    << "ParametrisedEMPhysics::ConstructProcess: bremsstrahlung threshold Eth= "
+    << bremth/GeV << " GeV" 
+    << "\n                                         verbosity= " << verb
+    << "  fluoFlag: " << fluo; 
+
 }
 
 ParametrisedEMPhysics::~ParametrisedEMPhysics() {
@@ -157,18 +176,9 @@ void ParametrisedEMPhysics::ConstructProcess() {
       }
     }
   }
-  // bremsstrahlung threshold and EM verbosity
-  G4EmProcessOptions opt;
-  G4int verb = theParSet.getUntrackedParameter<int>("Verbosity",0);
-  opt.SetVerbose(verb - 1);
-
-  G4double bremth = theParSet.getParameter<double>("G4BremsstrahlungThreshold")*GeV; 
-  edm::LogInfo("SimG4CoreApplication") 
-    << "ParametrisedEMPhysics::ConstructProcess: bremsstrahlung threshold Eth= "
-    << bremth/GeV << " GeV"; 
-  opt.SetBremsstrahlungTh(bremth);
 
   // Russian roulette and tracking cut for e+-
+  G4EmProcessOptions opt;
   const G4int NREG = 6; 
   const G4String rname[NREG] = {"EcalRegion", "HcalRegion", "MuonIron",
 				"PreshowerRegion","CastorRegion",
@@ -219,10 +229,10 @@ void ParametrisedEMPhysics::ConstructProcess() {
   if(fluo) {
     G4VAtomDeexcitation* de = new G4UAtomicDeexcitation();
     G4LossTableManager::Instance()->SetAtomDeexcitation(de);
-    de->SetFluo(true);
   }
   // enable muon nuclear (valid option for Geant4 10.0pX only)
-  bool munuc = theParSet.getParameter<bool>("FlagMuNucl");
+  // bool munuc = theParSet.getParameter<bool>("FlagMuNucl");
+  bool munuc = false;
   if(munuc) {
     G4MuonNuclearProcess* muNucProcess = new G4MuonNuclearProcess();
     muNucProcess->RegisterMe(new G4MuonVDNuclearModel());

SimG4Core/Application/src/RunManager.cc

@@ -300,7 +300,7 @@ void RunManager::initG4(const edm::EventSetup & es)
   }
 
   if("" != m_WriteFile) {
-    G4GDMLParser gdml(new G4GDMLReadStructure(), new CMSGDMLWriteStructure());
+    G4GDMLParser gdml;
     gdml.Write(m_WriteFile, world->GetWorldVolume(), true);
   }
 

SimG4Core/Application/src/RunManagerMT.cc

@@ -71,7 +71,7 @@ RunManagerMT::RunManagerMT(edm::ParameterSet const & p):
       m_G4Commands(p.getParameter<std::vector<std::string> >("G4Commands")),
       m_fieldBuilder(nullptr)
 {    
-  m_currentRun = 0;
+  m_currentRun = nullptr;
   G4RunManagerKernel *kernel = G4MTRunManagerKernel::GetRunManagerKernel();
   if(!kernel) m_kernel = new G4MTRunManagerKernel();
   else {
@@ -96,13 +96,19 @@ void RunManagerMT::initG4(const DDCompactView *pDD, const MagneticField *pMF,
 			  const HepPDT::ParticleDataTable *fPDGTable)
 {
   if (m_managerInitialized) return;
+
+  edm::LogInfo("SimG4CoreApplication") 
+    << "RunManagerMT: start initialisation of geometry";
 
   // DDDWorld: get the DDCV from the ES and use it to build the World
   G4LogicalVolumeToDDLogicalPartMap map_;
   m_world.reset(new DDDWorld(pDD, map_, m_catalog, false));
   m_registry.dddWorldSignal_(m_world.get());
 
   // setup the magnetic field
+  edm::LogInfo("SimG4CoreApplication") 
+    << "RunManagerMT: start initialisation of magnetic field";
+
   if (m_pUseMagneticField)
     {
       const GlobalPoint g(0.,0.,0.);
@@ -119,6 +125,9 @@ void RunManagerMT::initG4(const DDCompactView *pDD, const MagneticField *pMF,
     }
 
   // Create physics list
+  edm::LogInfo("SimG4CoreApplication") 
+    << "RunManagerMT: create PhysicsList";
+
   std::unique_ptr<PhysicsListMakerBase>
     physicsMaker(PhysicsListFactory::get()->create(
       m_pPhysics.getParameter<std::string> ("type")));
@@ -161,15 +170,14 @@ void RunManagerMT::initG4(const DDCompactView *pDD, const MagneticField *pMF,
     throw SimG4Exception("G4RunManagerKernel initialization failed!"); 
   }
 
-  if (m_StorePhysicsTables)
-    {
-      std::ostringstream dir;
-      dir << m_PhysicsTablesDir << '\0';
-      std::string cmd = std::string("/control/shell mkdir -p ")+m_PhysicsTablesDir;
-      if (!std::ifstream(dir.str().c_str(), std::ios::in))
-        G4UImanager::GetUIpointer()->ApplyCommand(cmd);
-      m_physicsList->StorePhysicsTable(m_PhysicsTablesDir);
-    }
+  if (m_StorePhysicsTables) {
+    std::ostringstream dir;
+    dir << m_PhysicsTablesDir << '\0';
+    std::string cmd = std::string("/control/shell mkdir -p ")+m_PhysicsTablesDir;
+    if (!std::ifstream(dir.str().c_str(), std::ios::in))
+      G4UImanager::GetUIpointer()->ApplyCommand(cmd);
+    m_physicsList->StorePhysicsTable(m_PhysicsTablesDir);
+  }
 
   initializeUserActions();
 
@@ -183,7 +191,7 @@ void RunManagerMT::initG4(const DDCompactView *pDD, const MagneticField *pMF,
 
   // geometry dump
   if("" != m_WriteFile) {
-    G4GDMLParser gdml(new G4GDMLReadStructure(), new CMSGDMLWriteStructure());
+    G4GDMLParser gdml;
     gdml.Write(m_WriteFile, m_world->GetWorldVolume(), true);
   }
 
@@ -227,8 +235,6 @@ void RunManagerMT::terminateRun() {
   m_userRunAction->EndOfRunAction(m_currentRun);
   delete m_userRunAction;
   m_userRunAction = 0;
-  // delete m_currentRun;
-  //m_currentRun = 0;
   if(m_kernel && !m_runTerminated) {
     m_kernel->RunTermination();
     m_runTerminated = true;

SimG4Core/Geometry/BuildFile.xml

@@ -3,6 +3,7 @@
 <use   name="Geometry/Records"/>
 <use   name="FWCore/MessageLogger"/>
 <use   name="FWCore/Utilities"/>
+<use   name="xerces-c"/>
 <use   name="expat"/>
 <export>
   <lib   name="1"/>

SimG4Core/Geometry/src/G4CheckOverlap.cc

@@ -2,11 +2,15 @@
 
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 
+#include "G4RegionStore.hh"
+#include "G4Region.hh"
 #include "G4PhysicalVolumeStore.hh"
 #include "G4VPhysicalVolume.hh"
+#include "G4LogicalVolume.hh"
 #include "G4GeomTestVolume.hh"
 #include "globals.hh"
 #include "G4SystemOfUnits.hh"
+#include "G4GDMLParser.hh"
 
 #include <string>
 
@@ -18,27 +22,76 @@ G4CheckOverlap::G4CheckOverlap(const edm::ParameterSet &p) {
     = p.getUntrackedParameter<double>("Tolerance", 0.0)*CLHEP::mm;
   int nPoints = p.getUntrackedParameter<int>("Resolution", 10000);
   bool verbose = p.getUntrackedParameter<bool>("Verbose", true);
+  bool regionFlag = p.getUntrackedParameter<bool>("RegionFlag", true);
+  bool gdmlFlag = p.getUntrackedParameter<bool>("gdmlFlag", false);
   int nPrints = p.getUntrackedParameter<int>("ErrorThreshold", 1);
   int level = p.getUntrackedParameter<int>("Level", 0);
   int depth = p.getUntrackedParameter<int>("Depth", -1);
 
+  const G4RegionStore* regStore =  G4RegionStore::GetInstance();
+
+  G4LogicalVolume* lv;
+  const G4PhysicalVolumeStore * pvs = G4PhysicalVolumeStore::GetInstance();
+  unsigned int numPV = pvs->size();
+
   unsigned int nn = nodeNames.size();
-  G4cout << "Initialised with " 
+  G4cout << "G4OverlapCheck is initialised with " 
 	 << nodeNames.size() << " nodes; " << " nPoints= " << nPoints
 	 << "; tolerance= " << tolerance/mm << " mm; verbose: " 
-	 << verbose << G4endl;
-  const G4PhysicalVolumeStore * pvs = G4PhysicalVolumeStore::GetInstance();
+	 << verbose << " RegionFlag: " << regionFlag << G4endl;
+
   if(0 < nn) { 
     for (unsigned int ii=0; ii<nn; ++ii) {
-      if("" == nodeNames[ii] || "world" == nodeNames[ii] 
-	 || "World" == nodeNames[ii] ) { 
+      if("" == nodeNames[ii] || "world" == nodeNames[ii] || "World" == nodeNames[ii] ) { 
 	nodeNames[ii] = "DDDWorld"; 
+	G4cout << "### Check overlaps for DDDWorld " << G4endl; 
+	G4VPhysicalVolume* pv = pvs->GetVolume("DDDWorld");
+	G4GeomTestVolume test(pv, tolerance, nPoints, verbose);
+	test.SetErrorsThreshold(nPrints);
+	test.TestRecursiveOverlap(level, depth);
+      } else if(regionFlag) {
+        G4cout << "---------------------------------------------------------------" << G4endl;
+	G4cout << "### Check overlaps for G4Region Node[" << ii << "] : " << nodeNames[ii] 
+	       << G4endl; 
+        G4Region* reg = regStore->GetRegion((G4String)nodeNames[ii]);
+        if(!reg) {
+	  G4cout << "### NO G4Region found - EXIT" << G4endl;
+          return; 
+	}
+	std::vector<G4LogicalVolume*>::iterator rootLVItr
+	  = reg->GetRootLogicalVolumeIterator();
+	unsigned int numRootLV = reg->GetNumberOfRootVolumes();
+
+	for(unsigned int iLV=0; iLV < numRootLV; ++iLV, ++rootLVItr ) {
+	  // Cover each root logical volume in this region
+	  lv = *rootLVItr;
+	  G4cout << "### Check overlaps for G4LogicalVolume " << lv->GetName() << G4endl; 
+          for(unsigned int i=0; i<numPV; ++i) {
+            if(((*pvs)[i])->GetLogicalVolume() == lv) {
+	      G4cout << "### Check overlaps for PhysVolume  " <<  ((*pvs)[i])->GetName()
+		     << G4endl;
+	      // gdml dump only for 1 volume
+              if(gdmlFlag) {
+		G4GDMLParser gdml;
+		gdml.Write(((*pvs)[i])->GetName()+".gdml", (*pvs)[i], true);
+                gdmlFlag = false;
+	      }
+
+	      G4GeomTestVolume test(((*pvs)[i]), tolerance, nPoints, verbose);
+	      test.SetErrorsThreshold(nPrints);
+	      test.TestRecursiveOverlap(level, depth);
+	    }
+	  }
+	}
+
+      } else {
+	G4cout << "### Check overlaps for PhysVolume Node[" << ii << "] : " << nodeNames[ii] 
+	       << G4endl; 
+	G4VPhysicalVolume* pv = pvs->GetVolume((G4String)nodeNames[ii]);
+	G4GeomTestVolume test(pv, tolerance, nPoints, verbose);
+	test.SetErrorsThreshold(nPrints);
+	test.TestRecursiveOverlap(level, depth);
       }
-      G4cout << "Check overlaps for Node[" << ii << "] : " << nodeNames[ii] << G4endl; 
-      G4VPhysicalVolume* pv = pvs->GetVolume((G4String)nodeNames[ii]);
-      G4GeomTestVolume test(pv, tolerance, nPoints, verbose);
-      test.SetErrorsThreshold(nPrints);
-      test.TestRecursiveOverlap(level, depth);
     }
   }
 }