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