Application/src/Phase2SteppingAction.cc

0001 #include "SimG4Core/Application/interface/Phase2SteppingAction.h"
0002 #include "SimG4Core/Geometry/interface/DD4hep2DDDName.h"
0003 #include "SimG4Core/Notification/interface/TrackInformation.h"
0004 #include "SimG4Core/Notification/interface/CMSSteppingVerbose.h"
0005
0006 #include "G4LogicalVolumeStore.hh"
0007 #include "G4ParticleTable.hh"
0008 #include "G4PhysicalVolumeStore.hh"
0009 #include "G4RegionStore.hh"
0010 #include "G4UnitsTable.hh"
0011 #include <CLHEP/Units/SystemOfUnits.h>
0012
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/Utilities/interface/isFinite.h"
0015
0016 Phase2SteppingAction::Phase2SteppingAction(const CMSSteppingVerbose* sv,
0017                                            const edm::ParameterSet& p,
0018                                            bool hasW,
0019                                            bool dd4hep)
0020     : steppingVerbose(sv), hasWatcher(hasW), dd4hep_(dd4hep) {
0021   theCriticalEnergyForVacuum = (p.getParameter<double>("CriticalEnergyForVacuum") * CLHEP::MeV);
0022   if (0.0 < theCriticalEnergyForVacuum) {
0023     killBeamPipe = true;
0024   }
0025   theCriticalDensity = (p.getParameter<double>("CriticalDensity") * CLHEP::g / CLHEP::cm3);
0026   maxZCentralCMS = p.getParameter<double>("MaxZCentralCMS") * CLHEP::m;
0027   maxTrackTime = p.getParameter<double>("MaxTrackTime") * CLHEP::ns;
0028   maxTrackTimeForward = p.getParameter<double>("MaxTrackTimeForward") * CLHEP::ns;
0029   maxTrackTimes = p.getParameter<std::vector<double> >("MaxTrackTimes");
0030   maxTimeNames = p.getParameter<std::vector<std::string> >("MaxTimeNames");
0031   deadRegionNames = p.getParameter<std::vector<std::string> >("DeadRegions");
0032   maxNumberOfSteps = p.getParameter<int>("MaxNumberOfSteps");
0033   ekinMins = p.getParameter<std::vector<double> >("EkinThresholds");
0034   ekinNames = p.getParameter<std::vector<std::string> >("EkinNames");
0035   ekinParticles = p.getParameter<std::vector<std::string> >("EkinParticles");
0036   cmseName_ = (G4String)(p.getParameter<std::string>("CMSName"));
0037   trackerName_ = (G4String)(p.getParameter<std::string>("TrackerName"));
0038   caloName_ = (G4String)(p.getParameter<std::string>("CaloName"));
0039   btlName_ = (G4String)(p.getParameter<std::string>("BTLName"));
0040   cms2ZDCName_ = p.getParameter<std::string>("CMS2ZDCName");
0041
0042   edm::LogVerbatim("SimG4CoreApplication")
0043       << "Phase2SteppingAction:: KillBeamPipe = " << killBeamPipe
0044       << " CriticalDensity = " << theCriticalDensity * CLHEP::cm3 / CLHEP::g << " g/cm3\n"
0045       << "                 CriticalEnergyForVacuum = " << theCriticalEnergyForVacuum / CLHEP::MeV << " Mev;"
0046       << " MaxTrackTime = " << maxTrackTime / CLHEP::ns << " ns;"
0047       << " MaxZCentralCMS = " << maxZCentralCMS / CLHEP::m << " m"
0048       << " MaxTrackTimeForward = " << maxTrackTimeForward / CLHEP::ns << " ns"
0049       << " MaxNumberOfSteps = " << maxNumberOfSteps << "\n"
0050       << "                 Names of special volumes: " << cmseName_ << "  " << trackerName_ << "  " << caloName_ << "  "
0051       << btlName_;
0052
0053   numberTimes = maxTrackTimes.size();
0054   if (numberTimes > 0) {
0055     for (unsigned int i = 0; i < numberTimes; i++) {
0056       edm::LogVerbatim("SimG4CoreApplication")
0057           << "Phase2SteppingAction::MaxTrackTime for " << maxTimeNames[i] << " is " << maxTrackTimes[i] << " ns ";
0058       maxTrackTimes[i] *= CLHEP::ns;
0059     }
0060   }
0061
0062   ndeadRegions = deadRegionNames.size();
0063   if (ndeadRegions > 0) {
0064     edm::LogVerbatim("SimG4CoreApplication")
0065         << "Phase2SteppingAction: Number of DeadRegions where all trackes are killed " << ndeadRegions;
0066     for (unsigned int i = 0; i < ndeadRegions; ++i) {
0067       edm::LogVerbatim("SimG4CoreApplication")
0068           << "Phase2SteppingAction: DeadRegion " << i << ".  " << deadRegionNames[i];
0069     }
0070   }
0071   numberEkins = ekinNames.size();
0072   numberPart = ekinParticles.size();
0073   if (0 == numberPart) {
0074     numberEkins = 0;
0075   }
0076
0077   if (numberEkins > 0) {
0078     edm::LogVerbatim("SimG4CoreApplication")
0079         << "Phase2SteppingAction::Kill following " << numberPart << " particles in " << numberEkins << " volumes";
0080     for (unsigned int i = 0; i < numberPart; ++i) {
0081       edm::LogVerbatim("SimG4CoreApplication") << "Phase2SteppingAction::Particle " << i << " " << ekinParticles[i]
0082                                                << "  Threshold = " << ekinMins[i] << " MeV";
0083       ekinMins[i] *= CLHEP::MeV;
0084     }
0085     for (unsigned int i = 0; i < numberEkins; ++i) {
0086       edm::LogVerbatim("SimG4CoreApplication") << "Phase2SteppingAction::LogVolume[" << i << "] = " << ekinNames[i];
0087     }
0088   }
0089 }
0090
0091 void Phase2SteppingAction::UserSteppingAction(const G4Step* aStep) {
0092   if (!initialized) {
0093     initialized = initPointer();
0094   }
0095
0096   m_g4StepSignal(aStep);
0097
0098   G4Track* theTrack = aStep->GetTrack();
0099   TrackStatus tstat = (theTrack->GetTrackStatus() == fAlive) ? sAlive : sKilledByProcess;
0100
0101   if (theTrack->GetKineticEnergy() < 0.0) {
0102     if (nWarnings < 2) {
0103       ++nWarnings;
0104       edm::LogWarning("SimG4CoreApplication")
0105           << "Phase2SteppingAction::UserPhase2SteppingAction: Track #" << theTrack->GetTrackID() << " "
0106           << theTrack->GetDefinition()->GetParticleName()
0107           << " Ekin(MeV)= " << theTrack->GetKineticEnergy() / CLHEP::MeV;
0108     }
0109     theTrack->SetKineticEnergy(0.0);
0110   }
0111
0112   const G4StepPoint* preStep = aStep->GetPreStepPoint();
0113   const G4StepPoint* postStep = aStep->GetPostStepPoint();
0114
0115   // the track is killed by the process
0116   if (tstat == sKilledByProcess) {
0117     if (nullptr != steppingVerbose) {
0118       steppingVerbose->nextStep(aStep, fpSteppingManager, false);
0119     }
0120     return;
0121   }
0122
0123   if (sAlive == tstat && theTrack->GetCurrentStepNumber() > maxNumberOfSteps) {
0124     tstat = sNumberOfSteps;
0125     if (nWarnings < 5) {
0126       ++nWarnings;
0127       edm::LogWarning("SimG4CoreApplication")
0128           << "Track #" << theTrack->GetTrackID() << " " << theTrack->GetDefinition()->GetParticleName()
0129           << " E(MeV)= " << preStep->GetKineticEnergy() / CLHEP::MeV << " Nstep= " << theTrack->GetCurrentStepNumber()
0130           << " is killed due to limit on number of steps;/n  PV: " << preStep->GetPhysicalVolume()->GetName() << " at "
0131           << theTrack->GetPosition() << " StepLen(mm)= " << aStep->GetStepLength();
0132     }
0133   }
0134   const double time = theTrack->GetGlobalTime();
0135
0136   // check Z-coordinate
0137   if (sAlive == tstat && std::abs(theTrack->GetPosition().z()) >= maxZCentralCMS) {
0138     tstat = (time > maxTrackTimeForward) ? sOutOfTime : sVeryForward;
0139   }
0140
0141   // check G4Region
0142   if (sAlive == tstat) {
0143     // next logical volume and next region
0144     const G4LogicalVolume* lv = postStep->GetPhysicalVolume()->GetLogicalVolume();
0145     const G4Region* theRegion = lv->GetRegion();
0146
0147     // kill in dead regions
0148     if (isInsideDeadRegion(theRegion))
0149       tstat = sDeadRegion;
0150
0151     // kill out of time
0152     if (sAlive == tstat) {
0153       if (isOutOfTimeWindow(theRegion, time))
0154         tstat = sOutOfTime;
0155     }
0156
0157     // kill low-energy in volumes on demand
0158     if (sAlive == tstat && numberEkins > 0) {
0159       if (isLowEnergy(lv, theTrack))
0160         tstat = sLowEnergy;
0161     }
0162
0163     // kill low-energy in vacuum
0164     if (sAlive == tstat && killBeamPipe) {
0165       if (theTrack->GetKineticEnergy() < theCriticalEnergyForVacuum &&
0166           theTrack->GetDefinition()->GetPDGCharge() != 0.0 && lv->GetMaterial()->GetDensity() <= theCriticalDensity) {
0167         tstat = sLowEnergyInVacuum;
0168       }
0169     }
0170   }
0171   // check transition tracker/btl and tracker/calo
0172   bool isKilled = false;
0173   if (sAlive == tstat || sVeryForward == tstat) {
0174     // store TrackInformation about transition from one envelope to another
0175     if (preStep->GetPhysicalVolume() == tracker && postStep->GetPhysicalVolume() == btl) {
0176       // store transition tracker -> BTL only for tracks entering BTL for the first time
0177       TrackInformation* trkinfo = static_cast<TrackInformation*>(theTrack->GetUserInformation());
0178       if (!trkinfo->isFromTtoBTL() && !trkinfo->isFromBTLtoT()) {
0179         trkinfo->setFromTtoBTL();
0180 #ifdef EDM_ML_DEBUG
0181         LogDebug("SimG4CoreApplication") << "Setting flag for Tracker -> BTL " << trkinfo->isFromTtoBTL()
0182                                          << " IdAtBTLentrance = " << trkinfo->mcTruthID();
0183 #endif
0184       } else {
0185         trkinfo->setBTLlooper();
0186 #ifdef EDM_ML_DEBUG
0187         LogDebug("SimG4CoreApplication") << "Setting flag for BTL looper " << trkinfo->isBTLlooper();
0188 #endif
0189       }
0190     } else if (preStep->GetPhysicalVolume() == btl && postStep->GetPhysicalVolume() == tracker) {
0191       // store transition BTL -> tracker
0192       TrackInformation* trkinfo = static_cast<TrackInformation*>(theTrack->GetUserInformation());
0193       if (!trkinfo->isFromBTLtoT()) {
0194         trkinfo->setFromBTLtoT();
0195 #ifdef EDM_ML_DEBUG
0196         LogDebug("SimG4CoreApplication") << "Setting flag for BTL -> Tracker " << trkinfo->isFromBTLtoT();
0197 #endif
0198       }
0199     } else if (preStep->GetPhysicalVolume() == tracker && postStep->GetPhysicalVolume() == calo) {
0200       // store transition tracker -> calo
0201       TrackInformation* trkinfo = static_cast<TrackInformation*>(theTrack->GetUserInformation());
0202       if (!trkinfo->crossedBoundary()) {
0203         trkinfo->setCrossedBoundary(theTrack);
0204       }
0205     } else if ((preStep->GetPhysicalVolume() == calo && postStep->GetPhysicalVolume() == tracker) ||
0206                (preStep->GetPhysicalVolume() == cmse && postStep->GetPhysicalVolume() == tracker)) {
0207       // accounting for both geometries with direct CALO -> Tracker transitions or older versions with CMSE in the middle
0208       TrackInformation* trkinfo = static_cast<TrackInformation*>(theTrack->GetUserInformation());
0209       trkinfo->setInTrkFromBackscattering();
0210 #ifdef EDM_ML_DEBUG
0211       LogDebug("SimG4CoreApplication") << "Setting flag for backscattering from CALO "
0212                                        << trkinfo->isInTrkFromBackscattering();
0213 #endif
0214     }
0215   } else {
0216     theTrack->SetTrackStatus(fStopAndKill);
0217     isKilled = true;
0218 #ifdef EDM_ML_DEBUG
0219     PrintKilledTrack(theTrack, tstat);
0220 #endif
0221   }
0222   if (nullptr != steppingVerbose) {
0223     steppingVerbose->nextStep(aStep, fpSteppingManager, isKilled);
0224   }
0225 }
0226
0227 bool Phase2SteppingAction::isLowEnergy(const G4LogicalVolume* lv, const G4Track* theTrack) const {
0228   const double ekin = theTrack->GetKineticEnergy();
0229   int pCode = theTrack->GetDefinition()->GetPDGEncoding();
0230
0231   for (auto& vol : ekinVolumes) {
0232     if (lv == vol) {
0233       for (unsigned int i = 0; i < numberPart; ++i) {
0234         if (pCode == ekinPDG[i]) {
0235           return (ekin <= ekinMins[i]);
0236         }
0237       }
0238       break;
0239     }
0240   }
0241   return false;
0242 }
0243
0244 bool Phase2SteppingAction::initPointer() {
0245   const G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
0246   for (auto const& pvcite : *pvs) {
0247     const std::string& pvname = (std::string)(DD4hep2DDDName::namePV(pvcite->GetName(), dd4hep_));
0248     if (pvname == trackerName_) {
0249       tracker = pvcite;
0250     } else if (pvname == caloName_) {
0251       calo = pvcite;
0252     } else if (pvname == btlName_) {
0253       btl = pvcite;
0254     } else if (pvname == cmseName_) {
0255       cmse = pvcite;
0256     }
0257     if (tracker && calo && btl && cmse)
0258       break;
0259   }
0260   edm::LogVerbatim("SimG4CoreApplication")
0261       << "Phase2SteppingAction: pointer for Tracker " << tracker << " and for Calo " << calo << " and for BTL " << btl;
0262
0263   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0264   if (numberEkins > 0) {
0265     ekinVolumes.resize(numberEkins, nullptr);
0266     for (auto const& lvcite : *lvs) {
0267       std::string lvname = (std::string)(DD4hep2DDDName::nameMatterLV(lvcite->GetName(), dd4hep_));
0268       for (unsigned int i = 0; i < numberEkins; ++i) {
0269         if (lvname == ekinNames[i]) {
0270           ekinVolumes[i] = lvcite;
0271           break;
0272         }
0273       }
0274     }
0275     for (unsigned int i = 0; i < numberEkins; ++i) {
0276       edm::LogVerbatim("SimG4CoreApplication") << ekinVolumes[i]->GetName() << " with pointer " << ekinVolumes[i];
0277     }
0278   }
0279
0280   if (numberPart > 0) {
0281     G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
0282     ekinPDG.resize(numberPart, 0);
0283     for (unsigned int i = 0; i < numberPart; ++i) {
0284       const G4ParticleDefinition* part = theParticleTable->FindParticle(ekinParticles[i]);
0285       if (nullptr != part)
0286         ekinPDG[i] = part->GetPDGEncoding();
0287       edm::LogVerbatim("SimG4CoreApplication") << "Particle " << ekinParticles[i] << " with PDG code " << ekinPDG[i]
0288                                                << " and KE cut off " << ekinMins[i] / CLHEP::MeV << " MeV";
0289     }
0290   }
0291
0292   const G4RegionStore* rs = G4RegionStore::GetInstance();
0293   if (numberTimes > 0) {
0294     maxTimeRegions.resize(numberTimes, nullptr);
0295     for (auto const& rcite : *rs) {
0296       const G4String& rname = rcite->GetName();
0297       for (unsigned int i = 0; i < numberTimes; ++i) {
0298         if (rname == (G4String)(maxTimeNames[i])) {
0299           maxTimeRegions[i] = rcite;
0300           break;
0301         }
0302       }
0303     }
0304   }
0305   if (ndeadRegions > 0) {
0306     deadRegions.resize(ndeadRegions, nullptr);
0307     for (auto const& rcite : *rs) {
0308       const G4String& rname = rcite->GetName();
0309       for (unsigned int i = 0; i < ndeadRegions; ++i) {
0310         if (rname == (G4String)(deadRegionNames[i])) {
0311           deadRegions[i] = rcite;
0312           break;
0313         }
0314       }
0315     }
0316   }
0317   return true;
0318 }
0319
0320 void Phase2SteppingAction::PrintKilledTrack(const G4Track* aTrack, const TrackStatus& tst) const {
0321   std::string vname = "";
0322   std::string rname = "";
0323   std::string typ = " ";
0324   switch (tst) {
0325     case sDeadRegion:
0326       typ = " in dead region ";
0327       break;
0328     case sOutOfTime:
0329       typ = " out of time window ";
0330       break;
0331     case sLowEnergy:
0332       typ = " low energy limit ";
0333       break;
0334     case sLowEnergyInVacuum:
0335       typ = " low energy limit in vacuum ";
0336       break;
0337     case sEnergyDepNaN:
0338       typ = " energy deposition is NaN ";
0339       break;
0340     case sVeryForward:
0341       typ = " very forward track ";
0342       break;
0343     case sNumberOfSteps:
0344       typ = " too many steps ";
0345       break;
0346     default:
0347       break;
0348   }
0349   G4VPhysicalVolume* pv = aTrack->GetNextVolume();
0350   vname = pv->GetLogicalVolume()->GetName();
0351   rname = pv->GetLogicalVolume()->GetRegion()->GetName();
0352
0353   const double ekin = aTrack->GetKineticEnergy();
0354   if (ekin < 2 * CLHEP::MeV) {
0355     return;
0356   }
0357   edm::LogWarning("SimG4CoreApplication")
0358       << "Track #" << aTrack->GetTrackID() << " StepN= " << aTrack->GetCurrentStepNumber() << " "
0359       << aTrack->GetDefinition()->GetParticleName() << " E(MeV)=" << ekin / CLHEP::MeV
0360       << " T(ns)=" << aTrack->GetGlobalTime() / CLHEP::ns << " is killed due to " << typ << "\n  LV: " << vname << " ("
0361       << rname << ") at " << aTrack->GetPosition() << " step(cm)=" << aTrack->GetStep()->GetStepLength() / CLHEP::cm;
0362 }