Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-30 23:17:21

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