Application/src/SteppingAction.cc

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