File indexing completed on 2024-05-20 22:40:11
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
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
0137 if (sAlive == tstat && std::abs(theTrack->GetPosition().z()) >= maxZCentralCMS) {
0138 tstat = (time > maxTrackTimeForward) ? sOutOfTime : sVeryForward;
0139 }
0140
0141
0142 if (sAlive == tstat) {
0143
0144 const G4LogicalVolume* lv = postStep->GetPhysicalVolume()->GetLogicalVolume();
0145 const G4Region* theRegion = lv->GetRegion();
0146
0147
0148 if (isInsideDeadRegion(theRegion))
0149 tstat = sDeadRegion;
0150
0151
0152 if (sAlive == tstat) {
0153 if (isOutOfTimeWindow(theRegion, time))
0154 tstat = sOutOfTime;
0155 }
0156
0157
0158 if (sAlive == tstat && numberEkins > 0) {
0159 if (isLowEnergy(lv, theTrack))
0160 tstat = sLowEnergy;
0161 }
0162
0163
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
0172 bool isKilled = false;
0173 if (sAlive == tstat || sVeryForward == tstat) {
0174
0175 if (preStep->GetPhysicalVolume() == tracker && postStep->GetPhysicalVolume() == btl) {
0176
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
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
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
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 }