Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-04-09 02:02:20

0001 #include "SimG4Core/Application/interface/StackingAction.h"
0002 #include "SimG4Core/Application/interface/TrackingAction.h"
0003 #include "SimG4Core/Notification/interface/MCTruthUtil.h"
0004 #include "SimG4Core/Notification/interface/TrackInformation.h"
0005 #include "SimG4Core/Notification/interface/CMSSteppingVerbose.h"
0006 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0007 #include "SimG4Core/Physics/interface/CMSG4TrackInterface.h"
0008 
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 
0011 #include "G4VProcess.hh"
0012 #include "G4EmProcessSubType.hh"
0013 #include "G4LogicalVolumeStore.hh"
0014 #include "G4RegionStore.hh"
0015 #include "Randomize.hh"
0016 #include <CLHEP/Units/SystemOfUnits.h>
0017 #include "G4VSolid.hh"
0018 #include "G4TransportationManager.hh"
0019 #include "G4GammaGeneralProcess.hh"
0020 #include "G4LossTableManager.hh"
0021 
0022 StackingAction::StackingAction(const edm::ParameterSet& p, const CMSSteppingVerbose* sv) : steppingVerbose(sv) {
0023   trackNeutrino = p.getParameter<bool>("TrackNeutrino");
0024   killHeavy = p.getParameter<bool>("KillHeavy");
0025   killGamma = p.getParameter<bool>("KillGamma");
0026   kmaxGamma = p.getParameter<double>("GammaThreshold") * CLHEP::MeV;
0027   kmaxIon = p.getParameter<double>("IonThreshold") * CLHEP::MeV;
0028   kmaxProton = p.getParameter<double>("ProtonThreshold") * CLHEP::MeV;
0029   kmaxNeutron = p.getParameter<double>("NeutronThreshold") * CLHEP::MeV;
0030   killDeltaRay = p.getParameter<bool>("KillDeltaRay");
0031   limitEnergyForVacuum = p.getParameter<double>("CriticalEnergyForVacuum") * CLHEP::MeV;
0032   maxTrackTime = p.getParameter<double>("MaxTrackTime") * CLHEP::ns;
0033   maxTrackTimeForward = p.getParameter<double>("MaxTrackTimeForward") * CLHEP::ns;
0034   maxZCentralCMS = p.getParameter<double>("MaxZCentralCMS") * CLHEP::m;
0035   maxTrackTimes = p.getParameter<std::vector<double> >("MaxTrackTimes");
0036   maxTimeNames = p.getParameter<std::vector<std::string> >("MaxTimeNames");
0037   deadRegionNames = p.getParameter<std::vector<std::string> >("DeadRegions");
0038   savePDandCinAll = p.getUntrackedParameter<bool>("SaveAllPrimaryDecayProductsAndConversions", true);
0039   savePDandCinTracker = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInTracker", false);
0040   savePDandCinCalo = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInCalo", false);
0041   savePDandCinMuon = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInMuon", false);
0042   saveFirstSecondary = p.getUntrackedParameter<bool>("SaveFirstLevelSecondary", false);
0043 
0044   gRusRoEnerLim = p.getParameter<double>("RusRoGammaEnergyLimit") * CLHEP::MeV;
0045   nRusRoEnerLim = p.getParameter<double>("RusRoNeutronEnergyLimit") * CLHEP::MeV;
0046 
0047   gRusRoEcal = p.getParameter<double>("RusRoEcalGamma");
0048   gRusRoHcal = p.getParameter<double>("RusRoHcalGamma");
0049   gRusRoMuonIron = p.getParameter<double>("RusRoMuonIronGamma");
0050   gRusRoPreShower = p.getParameter<double>("RusRoPreShowerGamma");
0051   gRusRoCastor = p.getParameter<double>("RusRoCastorGamma");
0052   gRusRoWorld = p.getParameter<double>("RusRoWorldGamma");
0053 
0054   nRusRoEcal = p.getParameter<double>("RusRoEcalNeutron");
0055   nRusRoHcal = p.getParameter<double>("RusRoHcalNeutron");
0056   nRusRoMuonIron = p.getParameter<double>("RusRoMuonIronNeutron");
0057   nRusRoPreShower = p.getParameter<double>("RusRoPreShowerNeutron");
0058   nRusRoCastor = p.getParameter<double>("RusRoCastorNeutron");
0059   nRusRoWorld = p.getParameter<double>("RusRoWorldNeutron");
0060 
0061   gRusRoZDC = p.getParameter<double>("RusRoZDCGamma");
0062   gRusRoHGcal = p.getParameter<double>("RusRoHGcalGamma");
0063   nRusRoZDC = p.getParameter<double>("RusRoZDCNeutron");
0064   nRusRoHGcal = p.getParameter<double>("RusRoHGcalNeutron");
0065 
0066   if (gRusRoEnerLim > 0.0 && (gRusRoEcal < 1.0 || gRusRoHcal < 1.0 || gRusRoMuonIron < 1.0 || gRusRoPreShower < 1.0 ||
0067                               gRusRoCastor < 1.0 || gRusRoWorld < 1.0)) {
0068     gRRactive = true;
0069   }
0070   if (nRusRoEnerLim > 0.0 && (nRusRoEcal < 1.0 || nRusRoHcal < 1.0 || nRusRoMuonIron < 1.0 || nRusRoPreShower < 1.0 ||
0071                               nRusRoCastor < 1.0 || nRusRoWorld < 1.0)) {
0072     nRRactive = true;
0073   }
0074 
0075   if (p.exists("TestKillingOptions")) {
0076     killInCalo = (p.getParameter<edm::ParameterSet>("TestKillingOptions")).getParameter<bool>("KillInCalo");
0077     killInCaloEfH = (p.getParameter<edm::ParameterSet>("TestKillingOptions")).getParameter<bool>("KillInCaloEfH");
0078     edm::LogWarning("SimG4CoreApplication")
0079         << " *** Activating special test killing options in StackingAction \n"
0080         << " *** Kill secondaries in Calorimetetrs volume = " << killInCalo << "\n"
0081         << " *** Kill electromagnetic secondaries from hadrons in Calorimeters volume= " << killInCaloEfH;
0082   }
0083   m_trackInterface = CMSG4TrackInterface::instance();
0084 
0085   initPointer();
0086 
0087   edm::LogVerbatim("SimG4CoreApplication")
0088       << "StackingAction initiated with"
0089       << " flag for saving decay products in "
0090       << " Tracker: " << savePDandCinTracker << " in Calo: " << savePDandCinCalo << " in Muon: " << savePDandCinMuon
0091       << " everywhere: " << savePDandCinAll << "\n  saveFirstSecondary"
0092       << ": " << saveFirstSecondary << " Tracking neutrino flag: " << trackNeutrino
0093       << " Kill Delta Ray flag: " << killDeltaRay << " Kill hadrons/ions flag: " << killHeavy
0094       << " MaxZCentralCMS = " << maxZCentralCMS / CLHEP::m << " m"
0095       << " MaxTrackTimeForward = " << maxTrackTimeForward / CLHEP::ns << " ns";
0096 
0097   if (killHeavy) {
0098     edm::LogVerbatim("SimG4CoreApplication") << "StackingAction kill protons below " << kmaxProton / CLHEP::MeV
0099                                              << " MeV, neutrons below " << kmaxNeutron / CLHEP::MeV << " MeV and ions"
0100                                              << " below " << kmaxIon / CLHEP::MeV << " MeV";
0101   }
0102   killExtra = killDeltaRay || killHeavy || killInCalo || killInCaloEfH;
0103 
0104   edm::LogVerbatim("SimG4CoreApplication") << "StackingAction kill tracks with "
0105                                            << "time larger than " << maxTrackTime / CLHEP::ns << " ns ";
0106   numberTimes = maxTimeNames.size();
0107   if (0 < numberTimes) {
0108     for (unsigned int i = 0; i < numberTimes; ++i) {
0109       edm::LogVerbatim("SimG4CoreApplication")
0110           << "          MaxTrackTime for " << maxTimeNames[i] << " is " << maxTrackTimes[i] << " ns ";
0111       maxTrackTimes[i] *= CLHEP::ns;
0112     }
0113   }
0114   if (limitEnergyForVacuum > 0.0) {
0115     edm::LogVerbatim("SimG4CoreApplication")
0116         << "StackingAction LowDensity regions - kill if E < " << limitEnergyForVacuum / CLHEP::MeV << " MeV";
0117     printRegions(lowdensRegions, "LowDensity");
0118   }
0119   if (deadRegions.size() > 0.0) {
0120     edm::LogVerbatim("SimG4CoreApplication") << "StackingAction Dead regions - kill all secondaries ";
0121     printRegions(deadRegions, "Dead");
0122   }
0123   if (gRRactive) {
0124     edm::LogVerbatim("SimG4CoreApplication")
0125         << "StackingAction: "
0126         << "Russian Roulette for gamma Elimit(MeV)= " << gRusRoEnerLim / CLHEP::MeV << "\n"
0127         << "                 ECAL Prob= " << gRusRoEcal << "\n"
0128         << "                 HCAL Prob= " << gRusRoHcal << "\n"
0129         << "             MuonIron Prob= " << gRusRoMuonIron << "\n"
0130         << "            PreShower Prob= " << gRusRoPreShower << "\n"
0131         << "                HGCAL Prob= " << gRusRoHGcal << "\n"
0132         << "                  ZDC Prob= " << gRusRoZDC << "\n"
0133         << "               CASTOR Prob= " << gRusRoCastor << "\n"
0134         << "                World Prob= " << gRusRoWorld;
0135   }
0136   if (nRRactive) {
0137     edm::LogVerbatim("SimG4CoreApplication")
0138         << "StackingAction: "
0139         << "Russian Roulette for neutron Elimit(MeV)= " << nRusRoEnerLim / CLHEP::MeV << "\n"
0140         << "                 ECAL Prob= " << nRusRoEcal << "\n"
0141         << "                 HCAL Prob= " << nRusRoHcal << "\n"
0142         << "             MuonIron Prob= " << nRusRoMuonIron << "\n"
0143         << "            PreShower Prob= " << nRusRoPreShower << "\n"
0144         << "                HGCAL Prob= " << nRusRoHGcal << "\n"
0145         << "                  ZDC Prob= " << nRusRoZDC << "\n"
0146         << "               CASTOR Prob= " << nRusRoCastor << "\n"
0147         << "                World Prob= " << nRusRoWorld;
0148   }
0149 
0150   if (savePDandCinTracker) {
0151     edm::LogVerbatim("SimG4CoreApplication") << "StackingAction Tracker regions: ";
0152     printRegions(trackerRegions, "Tracker");
0153   }
0154   if (savePDandCinCalo) {
0155     edm::LogVerbatim("SimG4CoreApplication") << "StackingAction Calo regions: ";
0156     printRegions(caloRegions, "Calo");
0157   }
0158   if (savePDandCinMuon) {
0159     edm::LogVerbatim("SimG4CoreApplication") << "StackingAction Muon regions: ";
0160     printRegions(muonRegions, "Muon");
0161   }
0162   worldSolid = G4TransportationManager::GetTransportationManager()
0163                    ->GetNavigatorForTracking()
0164                    ->GetWorldVolume()
0165                    ->GetLogicalVolume()
0166                    ->GetSolid();
0167 }
0168 
0169 G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack(const G4Track* aTrack) {
0170   // G4 interface part
0171   G4ClassificationOfNewTrack classification = fUrgent;
0172   const int pdg = aTrack->GetDefinition()->GetPDGEncoding();
0173   const int abspdg = std::abs(pdg);
0174   auto track = const_cast<G4Track*>(aTrack);
0175   const G4VProcess* creatorProc = aTrack->GetCreatorProcess();
0176 
0177   if (creatorProc == nullptr && aTrack->GetParentID() != 0) {
0178     edm::LogWarning("StackingAction::ClassifyNewTrack")
0179         << " TrackID=" << aTrack->GetTrackID() << " ParentID=" << aTrack->GetParentID() << " "
0180         << aTrack->GetDefinition()->GetParticleName() << " Ekin(MeV)=" << aTrack->GetKineticEnergy();
0181   }
0182   if (aTrack->GetKineticEnergy() < 0.0) {
0183     edm::LogWarning("StackingAction::ClassifyNewTrack")
0184         << " TrackID=" << aTrack->GetTrackID() << " ParentID=" << aTrack->GetParentID() << " "
0185         << aTrack->GetDefinition()->GetParticleName() << " Ekin(MeV)=" << aTrack->GetKineticEnergy() << " creator "
0186         << creatorProc->GetProcessName();
0187   }
0188   // primary
0189   if (creatorProc == nullptr || aTrack->GetParentID() == 0) {
0190     if (!trackNeutrino && (abspdg == 12 || abspdg == 14 || abspdg == 16 || abspdg == 18)) {
0191       classification = fKill;
0192     } else if (worldSolid->Inside(aTrack->GetPosition()) == kOutside) {
0193       classification = fKill;
0194     } else {
0195       MCTruthUtil::primary(track);
0196     }
0197   } else {
0198     // secondary
0199     const G4Region* reg = aTrack->GetVolume()->GetLogicalVolume()->GetRegion();
0200     const double time = aTrack->GetGlobalTime();
0201 
0202     // definetly killed tracks
0203     if (aTrack->GetTrackStatus() == fStopAndKill) {
0204       classification = fKill;
0205     } else if (!trackNeutrino && (abspdg == 12 || abspdg == 14 || abspdg == 16 || abspdg == 18)) {
0206       classification = fKill;
0207 
0208     } else if (std::abs(aTrack->GetPosition().z()) >= maxZCentralCMS) {
0209       // very forward secondary
0210       if (time > maxTrackTimeForward) {
0211         classification = fKill;
0212       } else {
0213         const G4Track* mother = m_trackInterface->getCurrentTrack();
0214         MCTruthUtil::secondary(track, *mother, 0);
0215       }
0216 
0217     } else if (isItOutOfTimeWindow(reg, time)) {
0218       // time window check
0219       classification = fKill;
0220 
0221     } else {
0222       // potentially good for tracking
0223       const double ke = aTrack->GetKineticEnergy();
0224       G4int subType = (nullptr != creatorProc) ? creatorProc->GetProcessSubType() : 0;
0225 
0226       LogDebug("SimG4CoreApplication") << "##StackingAction:Classify Track " << aTrack->GetTrackID() << " Parent "
0227                                        << aTrack->GetParentID() << " " << aTrack->GetDefinition()->GetParticleName()
0228                                        << " Ekin(MeV)=" << ke / CLHEP::MeV << " subType=" << subType << " ";
0229 
0230       // kill tracks in specific regions
0231       if (isThisRegion(reg, deadRegions)) {
0232         classification = fKill;
0233       }
0234       if (classification != fKill && ke <= limitEnergyForVacuum && isThisRegion(reg, lowdensRegions)) {
0235         classification = fKill;
0236 
0237       } else if (classification != fKill) {
0238         // very low-energy gamma
0239         if (pdg == 22 && killGamma && ke < kmaxGamma) {
0240           classification = fKill;
0241         }
0242 
0243         // specific track killing - not for production
0244         if (killExtra && classification != fKill) {
0245           if (killHeavy && classification != fKill) {
0246             if (((pdg / 1000000000 == 1) && (((pdg / 10000) % 100) > 0) && (((pdg / 10) % 100) > 0) &&
0247                  (ke < kmaxIon)) ||
0248                 ((pdg == 2212) && (ke < kmaxProton)) || ((pdg == 2112) && (ke < kmaxNeutron))) {
0249               classification = fKill;
0250             }
0251           }
0252 
0253           if (killDeltaRay && classification != fKill && aTrack->GetParentID() > 0 &&
0254               G4TrackToParticleID::isGammaElectronPositron(aTrack)) {
0255             classification = fKill;
0256           }
0257 
0258           if (killInCalo && classification != fKill && isThisRegion(reg, caloRegions)) {
0259             classification = fKill;
0260           }
0261           if (killInCaloEfH && classification != fKill) {
0262             int pdgMother = std::abs(m_trackInterface->getCurrentTrack()->GetDefinition()->GetPDGEncoding());
0263             if ((pdg == 22 || abspdg == 11) && pdgMother != 11 && pdgMother != 22 && isThisRegion(reg, caloRegions)) {
0264               classification = fKill;
0265             }
0266           }
0267         }
0268 
0269         // Russian roulette && MC truth
0270         if (classification != fKill) {
0271           const G4Track* mother = m_trackInterface->getCurrentTrack();
0272           int flag = 0;
0273           if (savePDandCinAll) {
0274             flag = isItPrimaryDecayProductOrConversion(subType, *mother);
0275           } else {
0276             if ((savePDandCinTracker && isThisRegion(reg, trackerRegions)) ||
0277                 (savePDandCinCalo && isThisRegion(reg, caloRegions)) ||
0278                 (savePDandCinMuon && isThisRegion(reg, muonRegions))) {
0279               flag = isItPrimaryDecayProductOrConversion(subType, *mother);
0280             }
0281           }
0282           if (saveFirstSecondary && 0 == flag) {
0283             flag = isItFromPrimary(*mother, flag);
0284           }
0285 
0286           // Russian roulette
0287           if (2112 == pdg || 22 == pdg) {
0288             double currentWeight = aTrack->GetWeight();
0289 
0290             if (1.0 >= currentWeight) {
0291               double prob = 1.001;
0292               double elim = 0.0;
0293 
0294               // neutron
0295               if (nRRactive && pdg == 2112) {
0296                 elim = nRusRoEnerLim;
0297                 if (reg == regionEcal) {
0298                   prob = nRusRoEcal;
0299                 } else if (reg == regionHcal) {
0300                   prob = nRusRoHcal;
0301                 } else if (reg == regionMuonIron) {
0302                   prob = nRusRoMuonIron;
0303                 } else if (reg == regionPreShower) {
0304                   prob = nRusRoPreShower;
0305                 } else if (reg == regionHGcal) {
0306                   prob = nRusRoHGcal;
0307                 } else if (reg == regionZDC) {
0308                   prob = nRusRoZDC;
0309                 } else if (reg == regionCastor) {
0310                   prob = nRusRoCastor;
0311                 } else if (reg == regionWorld) {
0312                   prob = nRusRoWorld;
0313                 }
0314 
0315                 // gamma
0316               } else if (gRRactive && pdg == 22) {
0317                 elim = gRusRoEnerLim;
0318                 if (reg == regionEcal || reg == regionPreShower) {
0319                   if (rrApplicable(aTrack, *mother)) {
0320                     if (reg == regionEcal) {
0321                       prob = gRusRoEcal;
0322                     } else {
0323                       prob = gRusRoPreShower;
0324                     }
0325                   }
0326                 } else {
0327                   if (reg == regionHcal) {
0328                     prob = gRusRoHcal;
0329                   } else if (reg == regionMuonIron) {
0330                     prob = gRusRoMuonIron;
0331                   } else if (reg == regionHGcal) {
0332                     prob = gRusRoHGcal;
0333                   } else if (reg == regionZDC) {
0334                     prob = gRusRoZDC;
0335                   } else if (reg == regionCastor) {
0336                     prob = gRusRoCastor;
0337                   } else if (reg == regionWorld) {
0338                     prob = gRusRoWorld;
0339                   }
0340                 }
0341               }
0342               if (prob < 1.0 && aTrack->GetKineticEnergy() < elim) {
0343                 if (G4UniformRand() < prob) {
0344                   track->SetWeight(currentWeight / prob);
0345                 } else {
0346                   classification = fKill;
0347                 }
0348               }
0349             }
0350           }
0351           if (classification != fKill) {
0352             MCTruthUtil::secondary(track, *mother, flag);
0353           }
0354           LogDebug("SimG4CoreApplication")
0355               << "StackingAction:Classify Track " << aTrack->GetTrackID() << " Parent " << aTrack->GetParentID()
0356               << " Type " << aTrack->GetDefinition()->GetParticleName() << " Ekin=" << ke / CLHEP::MeV
0357               << " MeV from process subType=" << subType << " as " << classification << " Flag: " << flag;
0358         }
0359       }
0360     }
0361   }
0362   if (nullptr != steppingVerbose) {
0363     steppingVerbose->stackFilled(aTrack, (classification == fKill));
0364   }
0365   return classification;
0366 }
0367 
0368 void StackingAction::NewStage() {}
0369 
0370 void StackingAction::PrepareNewEvent() {}
0371 
0372 void StackingAction::initPointer() {
0373   // prepare region vector
0374   const unsigned int num = maxTimeNames.size();
0375   maxTimeRegions.resize(num, nullptr);
0376 
0377   // Russian roulette
0378   const std::vector<G4Region*>* rs = G4RegionStore::GetInstance();
0379 
0380   for (auto& reg : *rs) {
0381     const G4String& rname = reg->GetName();
0382     if ((gRusRoEcal < 1.0 || nRusRoEcal < 1.0) && rname == "EcalRegion") {
0383       regionEcal = reg;
0384     }
0385     if ((gRusRoHcal < 1.0 || nRusRoHcal < 1.0) && rname == "HcalRegion") {
0386       regionHcal = reg;
0387     }
0388     if ((gRusRoMuonIron < 1.0 || nRusRoMuonIron < 1.0) && rname == "MuonIron") {
0389       regionMuonIron = reg;
0390     }
0391     if ((gRusRoPreShower < 1.0 || nRusRoPreShower < 1.0) && rname == "PreshowerRegion") {
0392       regionPreShower = reg;
0393     }
0394     if ((gRusRoCastor < 1.0 || nRusRoCastor < 1.0) && rname == "CastorRegion") {
0395       regionCastor = reg;
0396     }
0397     if ((gRusRoZDC < 1.0 || nRusRoZDC < 1.0) && rname == "ZDCRegion") {
0398       regionZDC = reg;
0399     }
0400     if ((gRusRoHGcal < 1.0 || nRusRoHGcal < 1.0) && rname == "HGCalRegion") {
0401       regionHGcal = reg;
0402     }
0403     if ((gRusRoWorld < 1.0 || nRusRoWorld < 1.0) && rname == "DefaultRegionForTheWorld") {
0404       regionWorld = reg;
0405     }
0406 
0407     // time limits
0408     for (unsigned int i = 0; i < num; ++i) {
0409       if (rname == (G4String)(maxTimeNames[i])) {
0410         maxTimeRegions[i] = reg;
0411         break;
0412       }
0413     }
0414     //
0415     if (savePDandCinTracker &&
0416         (rname == "BeamPipe" || rname == "BeamPipeVacuum" || rname == "TrackerPixelSensRegion" ||
0417          rname == "TrackerPixelDeadRegion" || rname == "TrackerDeadRegion" || rname == "TrackerSensRegion" ||
0418          rname == "FastTimerRegionBTL" || rname == "FastTimerRegionETL" || rname == "FastTimerRegionSensBTL" ||
0419          rname == "FastTimerRegionSensETL")) {
0420       trackerRegions.push_back(reg);
0421     }
0422     if (savePDandCinCalo && (rname == "HcalRegion" || rname == "EcalRegion" || rname == "PreshowerSensRegion" ||
0423                              rname == "PreshowerRegion" || rname == "APDRegion" || rname == "HGCalRegion")) {
0424       caloRegions.push_back(reg);
0425     }
0426     if (savePDandCinMuon && (rname == "MuonChamber" || rname == "MuonSensitive_RPC" || rname == "MuonIron" ||
0427                              rname == "Muon" || rname == "MuonSensitive_DT-CSC")) {
0428       muonRegions.push_back(reg);
0429     }
0430     if (rname == "BeamPipeOutside" || rname == "BeamPipeVacuum") {
0431       lowdensRegions.push_back(reg);
0432     }
0433     for (auto& dead : deadRegionNames) {
0434       if (rname == (G4String)(dead)) {
0435         deadRegions.push_back(reg);
0436       }
0437     }
0438   }
0439 }
0440 
0441 bool StackingAction::isThisRegion(const G4Region* reg, std::vector<const G4Region*>& regions) const {
0442   bool flag = false;
0443   for (auto& region : regions) {
0444     if (reg == region) {
0445       flag = true;
0446       break;
0447     }
0448   }
0449   return flag;
0450 }
0451 
0452 int StackingAction::isItPrimaryDecayProductOrConversion(const int stype, const G4Track& mother) const {
0453   int flag = 0;
0454   auto motherInfo = static_cast<const TrackInformation*>(mother.GetUserInformation());
0455   // Check whether mother is a primary
0456   if (motherInfo->isPrimary()) {
0457     if (stype == fDecay) {
0458       flag = 1;
0459     } else if (stype == fGammaConversion) {
0460       flag = 2;
0461     }
0462   }
0463   return flag;
0464 }
0465 
0466 bool StackingAction::rrApplicable(const G4Track* aTrack, const G4Track& mother) const {
0467   auto motherInfo = static_cast<const TrackInformation*>(mother.GetUserInformation());
0468 
0469   // Check whether mother is gamma, e+, e-
0470   const int genID = motherInfo->genParticlePID();
0471   return (22 != genID && 11 != std::abs(genID));
0472 }
0473 
0474 int StackingAction::isItFromPrimary(const G4Track& mother, int flagIn) const {
0475   int flag = flagIn;
0476   if (flag != 1) {
0477     auto ptr = static_cast<const TrackInformation*>(mother.GetUserInformation());
0478     if (ptr->isPrimary()) {
0479       flag = 3;
0480     }
0481   }
0482   return flag;
0483 }
0484 
0485 bool StackingAction::isItOutOfTimeWindow(const G4Region* reg, const double& t) const {
0486   double tofM = maxTrackTime;
0487   for (unsigned int i = 0; i < numberTimes; ++i) {
0488     if (reg == maxTimeRegions[i]) {
0489       tofM = maxTrackTimes[i];
0490       break;
0491     }
0492   }
0493   return (t > tofM);
0494 }
0495 
0496 void StackingAction::printRegions(const std::vector<const G4Region*>& reg, const std::string& word) const {
0497   for (unsigned int i = 0; i < reg.size(); ++i) {
0498     edm::LogVerbatim("SimG4CoreApplication")
0499         << " StackingAction: " << word << "Region " << i << ". " << reg[i]->GetName();
0500   }
0501 }