working version of DRC monolothic in k4geo

example/SteeringFile_IDEA_o1_v03.py

@@ -20,7 +20,7 @@
 ## Outputfile from the simulation: .slcio, edm4hep.root and .root output files are supported
 SIM.outputFile = "testIDEA_o1_v03.root"
 ## Physics list to use in simulation
-SIM.physicsList = None
+SIM.physicsList = "FTFP_BERT"
 ## Verbosity use integers from 1(most) to 7(least) verbose
 ## or strings: VERBOSE, DEBUG, INFO, WARNING, ERROR, FATAL, ALWAYS
 SIM.printLevel = 3
@@ -177,7 +177,8 @@
 ##     default filter for calorimeter sensitive detectors;
 ##     this is applied if no other filter is used for a calorimeter
 ##
-SIM.filter.calo = "edep0"
+# note: do not turn on the calo filter, otherwise all optical photons will be killed!
+SIM.filter.calo = ""
 
 ##  list of filter objects: map between name and parameter dictionary
 SIM.filter.filters = {
@@ -251,7 +252,7 @@
 ################################################################################
 
 ##  direction of the particle gun, 3 vector
-SIM.gun.direction = (1.0, 1.0, 1.0)
+SIM.gun.direction = (1.0,0.1,0.1)
 
 ## choose the distribution of the random direction for theta
 ##
@@ -285,10 +286,10 @@
 SIM.gun.isotrop = False
 
 ## Maximal momentum when using distribution (default = 0.0)
-SIM.gun.momentumMax = 10000.0
+#SIM.gun.momentumMax = 10000.0
 
 ## Minimal momentum when using distribution (default = 0.0)
-SIM.gun.momentumMin = 0.0
+#SIM.gun.momentumMin = 0.0
 SIM.gun.multiplicity = 1
 SIM.gun.particle = "e-"
 
@@ -557,7 +558,7 @@ def Geant4Output2EDM4hep_DRC_plugin(dd4hepSimulation):
 ##     Set printlevel to DEBUG to see a printout of all range cuts,
 ##     but this only works if range cut is not "None"
 ##
-SIM.physics.rangecut = 0.7
+SIM.physics.rangecut = None
 
 ## Set of PDG IDs that will not be passed from the input record to Geant4.
 ##
@@ -621,7 +622,8 @@ def setupOpticalPhysics(kernel):
     opt = PhysicsList(kernel, 'Geant4OpticalPhotonPhysics/OpticalGammaPhys')
     opt.addParticleConstructor('G4OpticalPhoton')
     opt.VerboseLevel = 1
-    opt.BoundaryInvokeSD = True
+    # set BoundaryInvokeSD to true when using DRC wafer as the SD
+    # opt.BoundaryInvokeSD = True
     opt.enableUI()
     seq.adopt(opt)
 
@@ -653,7 +655,8 @@ def setupDRCFastSim(kernel):
     phys.adopt(ph)
     phys.dump()
 
-SIM.physics.setupUserPhysics(setupDRCFastSim)
+# turn-on fastsim if the skipScint option of the SD is set to false
+#SIM.physics.setupUserPhysics(setupDRCFastSim)
 
 ################################################################################
 ## Properties for the random number generator

plugins/DRCaloFastSimModel.h

@@ -82,42 +82,43 @@ class DRCFiberModel {
   DRCFiberModel()=default;
   ~DRCFiberModel()=default;
 
-  G4OpBoundaryProcess* pOpBoundaryProc{nullptr};
-  G4OpAbsorption* pOpAbsorption{nullptr};
-  G4OpWLS* pOpWLS{nullptr};
-  G4bool fProcAssigned{false};
-
-  FastFiberData mDataPrevious{FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0))};
-  FastFiberData mDataCurrent{FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0))};
-
-  G4int fSafety{2};
-  G4double mNtransport{0.};
-  G4double mTransportUnit{0.};
-  G4ThreeVector mFiberPos{G4ThreeVector(0)};
-  G4ThreeVector mFiberAxis{G4ThreeVector(0)};
-  G4bool fKill{false};
-  G4bool fTransported{false};
-  G4bool fSwitch{true};
-  G4int fVerbose{0};
+  G4OpBoundaryProcess* pOpBoundaryProc = nullptr;
+  G4OpAbsorption* pOpAbsorption = nullptr;
+  G4OpWLS* pOpWLS = nullptr;
+  G4bool fProcAssigned = false;
+
+  FastFiberData mDataPrevious = FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0));
+  FastFiberData mDataCurrent = FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0));
+
+  G4int fSafety = 1;
+  G4double mNtransport = 0.;
+  G4double mTransportUnit = 0.;
+  G4ThreeVector mFiberPos = G4ThreeVector(0);
+  G4ThreeVector mFiberAxis = G4ThreeVector(0);
+  G4bool fKill = false;
+  G4bool fTransported = false;
+  G4bool fSwitch = true;
+  G4int fVerbose = 0;
 
   G4bool checkTotalInternalReflection(const G4Track *track) {
     if (!fProcAssigned)
       setPostStepProc(track); // locate OpBoundaryProcess only once
 
+    G4int theStatus = pOpBoundaryProc->GetStatus();
+
+    if (fVerbose > 1) {
+      std::cout << "DRCFiberModel::checkTotalInternalReflection | TrackID = " << std::setw(4) << track->GetTrackID();
+      std::cout << " | G4OpBoundaryProcessStatus = " << std::setw(2) << theStatus;
+      std::cout << " | Track status = " << track->GetTrackStatus();
+      std::cout << " | StepLength = " << std::setw(9) << track->GetStepLength() << std::endl;
+    }
+
     if (track->GetTrackStatus() == fStopButAlive || track->GetTrackStatus() == fStopAndKill)
       return false;
 
     // accumulate step length
     mDataCurrent.AddStepLengthInterval(track->GetStepLength());
 
-    G4int theStatus = pOpBoundaryProc->GetStatus();
-
-    if (fVerbose > 1) {
-      G4cout << "DRCFiberModel::checkTotalInternalReflection | TrackID = " << std::setw(4) << track->GetTrackID();
-      G4cout << " | G4OpBoundaryProcessStatus = " << std::setw(2) << theStatus;
-      G4cout << " | StepLength = " << std::setw(9) << track->GetStepLength() << G4endl;
-    }
-
     // skip exceptional iteration with FresnelReflection
     if (theStatus == G4OpBoundaryProcessStatus::FresnelReflection)
       mDataCurrent.SetOpBoundaryStatus(theStatus);