Application/src/Phase2SteppingAction.cc

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