File indexing completed on 2023-08-06 22:43:36
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
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
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
0128 if (sAlive == tstat && std::abs(theTrack->GetPosition().z()) >= maxZCentralCMS) {
0129 tstat = (time > maxTrackTimeForward) ? sOutOfTime : sVeryForward;
0130 }
0131
0132
0133 if (sAlive == tstat) {
0134
0135 const G4LogicalVolume* lv = postStep->GetPhysicalVolume()->GetLogicalVolume();
0136 const G4Region* theRegion = lv->GetRegion();
0137
0138
0139 if (isInsideDeadRegion(theRegion))
0140 tstat = sDeadRegion;
0141
0142
0143 if (sAlive == tstat) {
0144 if (isOutOfTimeWindow(theRegion, time))
0145 tstat = sOutOfTime;
0146 }
0147
0148
0149 if (sAlive == tstat && numberEkins > 0) {
0150 if (isLowEnergy(lv, theTrack))
0151 tstat = sLowEnergy;
0152 }
0153
0154
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
0163 bool isKilled = false;
0164 if (sAlive == tstat || sVeryForward == tstat) {
0165
0166 if (preStep->GetPhysicalVolume() == tracker && postStep->GetPhysicalVolume() == btl) {
0167
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
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
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
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 }