copy cms-sw#43494

civanch · Dec 12, 2023 · 3ca07eb · 3ca07eb
1 parent ce34e99
commit 3ca07eb
Show file tree

Hide file tree

Showing 4 changed files with 28 additions and 9 deletions.
diff --git a/SimG4CMS/Forward/src/ZdcSD.cc b/SimG4CMS/Forward/src/ZdcSD.cc
@@ -122,18 +122,17 @@ double ZdcSD::getEnergyDeposit(const G4Step* aStep) {
 
   // preStepPoint information
   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
-  G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
-  std::string nameVolume = ForwardName::getName(currentPV->GetName());
 
   const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
   const G4ThreeVector& hit_mom = preStepPoint->GetMomentumDirection();
   G4double stepL = aStep->GetStepLength() / cm;
   G4double beta = preStepPoint->GetBeta();
   G4double charge = preStepPoint->GetCharge();
+  if (charge == 0.0) 
+    return 0.0;
 
   // theTrack information
   G4Track* theTrack = aStep->GetTrack();
-  G4String particleType = theTrack->GetDefinition()->GetParticleName();
   G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
 
 #ifdef EDM_ML_DEBUG
@@ -156,7 +155,10 @@ double ZdcSD::getEnergyDeposit(const G4Step* aStep) {
   // postStepPoint information
   G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
   G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume();
+  G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
+  std::string nameVolume = ForwardName::getName(currentPV->GetName());
   std::string postnameVolume = ForwardName::getName(postPV->GetName());
+  G4String particleType = theTrack->GetDefinition()->GetParticleName();
   edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: \n"
                                  << "  preStepPoint: " << nameVolume << "," << stepL << "," << stepE << "," << beta
                                  << "," << charge << "\n"
@@ -165,7 +167,7 @@ double ZdcSD::getEnergyDeposit(const G4Step* aStep) {
                                  << " Etot(GeV)= " << theTrack->GetTotalEnergy() / GeV;
 #endif
   const double bThreshold = 0.67;
-  if ((beta > bThreshold) && (charge != 0) && (nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber")) {
+  if (beta > bThreshold) {
 #ifdef EDM_ML_DEBUG
     edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit:  pass ";
 #endif
@@ -293,10 +295,12 @@ double ZdcSD::getEnergyDeposit(const G4Step* aStep) {
     // determine failure mode: beta, charge, and/or nameVolume
     if (beta <= bThreshold)
       edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: fail beta=" << beta;
-    if (charge == 0)
-      edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: fail charge=0";
+
+#ifdef EDM_ML_DEBUG
+    std::string nameVolume = ForwardName::getName(currentPV->GetName());
     if (!(nameVolume == "ZDC_EMFiber" || nameVolume == "ZDC_HadFiber"))
       edm::LogVerbatim("ForwardSim") << "ZdcSD::  getEnergyDeposit: fail nv=" << nameVolume;
+#endif
   }
 
   return NCherPhot;

diff --git a/SimG4Core/Application/interface/SteppingAction.h b/SimG4Core/Application/interface/SteppingAction.h
@@ -43,13 +43,15 @@ class SteppingAction : public G4UserSteppingAction {
 
   inline bool isInsideDeadRegion(const G4Region* reg) const;
   inline bool isOutOfTimeWindow(const G4Region* reg, const double& time) const;
+  inline bool isForZDC(const G4LogicalVolume* lv, int pdg) const;
 
   bool isLowEnergy(const G4LogicalVolume*, const G4Track*) const;
   void PrintKilledTrack(const G4Track*, const TrackStatus&) const;
 
   EventAction* eventAction_;
   const G4VPhysicalVolume *tracker, *calo;
   const CMSSteppingVerbose* steppingVerbose;
+  const G4LogicalVolume* m_CMStoZDC{nullptr};
   double theCriticalEnergyForVacuum;
   double theCriticalDensity;
   double maxTrackTime;
@@ -65,6 +67,7 @@ class SteppingAction : public G4UserSteppingAction {
 
   bool initialized;
   bool killBeamPipe;
+  bool m_CMStoZDCtransport;
   bool hasWatcher;
 
   std::vector<double> maxTrackTimes, ekinMins;
@@ -98,4 +101,8 @@ inline bool SteppingAction::isOutOfTimeWindow(const G4Region* reg, const double&
   return (time > tofM);
 }
 
+inline bool SteppingAction::isForZDC(const G4LogicalVolume* lv, int pdg) const {
+  return (m_CMStoZDCtransport && lv == m_CMStoZDC && (pdg == 22 || pdg == 2112));
+}
+
 #endif
diff --git a/SimG4Core/Application/python/g4SimHits_cfi.py b/SimG4Core/Application/python/g4SimHits_cfi.py
@@ -334,6 +334,7 @@
     ),
     SteppingAction = cms.PSet(
         common_maximum_time,
+        CMStoZDCtransport       = cms.bool(False),                 
         MaxNumberOfSteps        = cms.int32(20000),
         EkinNames               = cms.vstring(),
         EkinThresholds          = cms.vdouble(),
@@ -664,7 +665,9 @@
 ## Disable PPS from Run 3 PbPb runs
 ##
 from Configuration.Eras.Modifier_pp_on_PbPb_run3_cff import pp_on_PbPb_run3
-pp_on_PbPb_run3.toModify( g4SimHits, LHCTransport = False )
+pp_on_PbPb_run3.toModify(g4SimHits, LHCTransport = False, 
+                         SteppingAction = dict(
+                             CMStoZDCtransport = cms.bool(True) ) )
 
 ##
 ## Change ECAL time slices

diff --git a/SimG4Core/Application/src/SteppingAction.cc b/SimG4Core/Application/src/SteppingAction.cc
@@ -27,6 +27,7 @@ SteppingAction::SteppingAction(EventAction* e, const edm::ParameterSet& p, const
   if (0.0 < theCriticalEnergyForVacuum) {
     killBeamPipe = true;
   }
+  m_CMStoZDCtransport = (p.getParameter<bool>("CMStoZDCtransport"));
   theCriticalDensity = (p.getParameter<double>("CriticalDensity") * CLHEP::g / CLHEP::cm3);
   maxZCentralCMS = p.getParameter<double>("MaxZCentralCMS") * CLHEP::m;
   maxTrackTime = p.getParameter<double>("MaxTrackTime") * CLHEP::ns;
@@ -46,7 +47,7 @@ SteppingAction::SteppingAction(EventAction* e, const edm::ParameterSet& p, const
       << " MaxTrackTime = " << maxTrackTime / CLHEP::ns << " ns;"
       << " MaxZCentralCMS = " << maxZCentralCMS / CLHEP::m << " m"
       << " MaxTrackTimeForward = " << maxTrackTimeForward / CLHEP::ns << " ns"
-      << " MaxNumberOfSteps = " << maxNumberOfSteps;
+      << " MaxNumberOfSteps = " << maxNumberOfSteps << " ZDC: " << m_CMStoZDCtransport;
 
   numberTimes = maxTrackTimes.size();
   if (numberTimes > 0) {
@@ -157,8 +158,9 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) {
     const G4Region* theRegion = lv->GetRegion();
 
     // kill in dead regions
-    if (isInsideDeadRegion(theRegion))
+    if (isInsideDeadRegion(theRegion) && !isForZDC(lv, std::abs(theTrack->GetParticleDefinition()->GetPDGEncoding()))) {
       tstat = sDeadRegion;
+    }
 
     // kill out of time
     if (sAlive == tstat) {
@@ -245,6 +247,9 @@ bool SteppingAction::initPointer() {
     ekinVolumes.resize(numberEkins, nullptr);
     for (auto const& lvcite : *lvs) {
       const G4String& lvname = lvcite->GetName();
+      if (lvname == "CMStoZDC" || lvname == "cmsextent:CMStoZDC") {
+        m_CMStoZDC = lvcite;
+      }
       for (unsigned int i = 0; i < numberEkins; ++i) {
         if (lvname == (G4String)(ekinNames[i])) {
           ekinVolumes[i] = lvcite;