#include <SteppingAction.h>

Inheritance diagram for SteppingAction:

Public Member Functions
	SteppingAction (EventAction ea, const edm::ParameterSet &ps, const CMSSteppingVerbose , bool hasW)

void	UserSteppingAction (const G4Step *aStep) final

	~SteppingAction () override

Public Attributes
SimActivityRegistry::G4StepSignal	m_g4StepSignal

Private Member Functions
bool	initPointer ()

bool	isInsideDeadRegion (const G4Region *reg) const

bool	isLowEnergy (const G4Step *aStep) const

bool	isOutOfTimeWindow (G4Track theTrack, const G4Region reg) const

bool	isThisVolume (const G4VTouchable touch, const G4VPhysicalVolume pv) const

void	PrintKilledTrack (const G4Track *, const TrackStatus &) const

Private Attributes
const G4VPhysicalVolume *	calo

std::vector< std::string >	deadRegionNames

std::vector< const G4Region * >	deadRegions

std::vector< double >	ekinMins

std::vector< std::string >	ekinNames

std::vector< std::string >	ekinParticles

std::vector< int >	ekinPDG

std::vector< G4LogicalVolume * >	ekinVolumes

EventAction *	eventAction_

bool	hasWatcher

bool	initialized

bool	killBeamPipe

std::vector< std::string >	maxTimeNames

std::vector< const G4Region * >	maxTimeRegions

double	maxTrackTime

std::vector< double >	maxTrackTimes

unsigned int	ndeadRegions

unsigned int	numberEkins

unsigned int	numberPart

unsigned int	numberTimes

unsigned int	nWarnings

const CMSSteppingVerbose *	steppingVerbose

double	theCriticalDensity

double	theCriticalEnergyForVacuum

const G4VPhysicalVolume *	tracker

Detailed Description

Definition at line 30 of file SteppingAction.h.

Constructor & Destructor Documentation

SteppingAction::SteppingAction	(	EventAction *	ea,
		const edm::ParameterSet &	ps,
		const CMSSteppingVerbose *	sv,
		bool	hasW
	)

explicit

Definition at line 18 of file SteppingAction.cc.

References deadRegionNames, ekinMins, ekinNames, ekinParticles, g, edm::ParameterSet::getParameter(), mps_fire::i, killBeamPipe, maxTimeNames, maxTrackTime, maxTrackTimes, MeV, ndeadRegions, numberEkins, numberPart, numberTimes, theCriticalDensity, and theCriticalEnergyForVacuum.

     : eventAction_(e),
       tracker(nullptr),
       calo(nullptr),
       steppingVerbose(sv),
       nWarnings(0),
       initialized(false),
       killBeamPipe(false),
       hasWatcher(hasW) {
   theCriticalEnergyForVacuum = (p.getParameter<double>("CriticalEnergyForVacuum") * CLHEP::MeV);
   if (0.0 < theCriticalEnergyForVacuum) {
     killBeamPipe = true;
   }
   theCriticalDensity = (p.getParameter<double>("CriticalDensity") * CLHEP::g / CLHEP::cm3);
   maxTrackTime = p.getParameter<double>("MaxTrackTime") * CLHEP::ns;
   maxTrackTimes = p.getParameter<std::vector<double> >("MaxTrackTimes");
   maxTimeNames = p.getParameter<std::vector<std::string> >("MaxTimeNames");
   deadRegionNames = p.getParameter<std::vector<std::string> >("DeadRegions");
   ekinMins = p.getParameter<std::vector<double> >("EkinThresholds");
   ekinNames = p.getParameter<std::vector<std::string> >("EkinNames");
   ekinParticles = p.getParameter<std::vector<std::string> >("EkinParticles");
 
   edm::LogVerbatim("SimG4CoreApplication")
       << "SteppingAction:: KillBeamPipe = " << killBeamPipe
       << " CriticalDensity = " << theCriticalDensity * CLHEP::cm3 / CLHEP::g << " g/cm3;"
       << " CriticalEnergyForVacuum = " << theCriticalEnergyForVacuum / CLHEP::MeV << " Mev;"
       << " MaxTrackTime = " << maxTrackTime / CLHEP::ns << " ns";
 
   numberTimes = maxTrackTimes.size();
   if (numberTimes > 0) {
     for (unsigned int i = 0; i < numberTimes; i++) {
       edm::LogVerbatim("SimG4CoreApplication")
           << "SteppingAction::MaxTrackTime for " << maxTimeNames[i] << " is " << maxTrackTimes[i] << " ns ";
       maxTrackTimes[i] *= ns;
     }
   }
 
   ndeadRegions = deadRegionNames.size();
   if (ndeadRegions > 0) {
     edm::LogVerbatim("SimG4CoreApplication")
         << "SteppingAction: Number of DeadRegions where all trackes are killed " << ndeadRegions;
     for (unsigned int i = 0; i < ndeadRegions; ++i) {
       edm::LogVerbatim("SimG4CoreApplication") << "SteppingAction: DeadRegion " << i << ".  " << deadRegionNames[i];
     }
   }
   numberEkins = ekinNames.size();
   numberPart = ekinParticles.size();
   if (0 == numberPart) {
     numberEkins = 0;
   }
 
   if (numberEkins > 0) {
     edm::LogVerbatim("SimG4CoreApplication")
         << "SteppingAction::Kill following " << numberPart << " particles in " << numberEkins << " volumes";
     for (unsigned int i = 0; i < numberPart; ++i) {
       edm::LogVerbatim("SimG4CoreApplication")
           << "SteppingAction::Particle " << i << " " << ekinParticles[i] << "  Threshold = " << ekinMins[i] << " MeV";
       ekinMins[i] *= CLHEP::MeV;
     }
     for (unsigned int i = 0; i < numberEkins; ++i) {
       edm::LogVerbatim("SimG4CoreApplication") << "SteppingAction::LogVolume[" << i << "] = " << ekinNames[i];
     }
   }
 }

SteppingAction::~SteppingAction ( )

override

Definition at line 83 of file SteppingAction.cc.

83 {}

Member Function Documentation

bool SteppingAction::initPointer ( )

private

Definition at line 187 of file SteppingAction.cc.

References calo, deadRegionNames, deadRegions, ekinMins, ekinNames, ekinParticles, ekinPDG, ekinVolumes, mps_fire::i, LogDebug, maxTimeNames, maxTimeRegions, MeV, ndeadRegions, numberEkins, numberPart, numberTimes, FSQDQM_cfi::pvs, and tracker.

Referenced by UserSteppingAction().

                                  {
   const G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
   if (pvs) {
     std::vector<G4VPhysicalVolume*>::const_iterator pvcite;
     for (pvcite = pvs->begin(); pvcite != pvs->end(); ++pvcite) {
       if ((*pvcite)->GetName() == "Tracker")
         tracker = (*pvcite);
       if ((*pvcite)->GetName() == "CALO")
         calo = (*pvcite);
       if (tracker && calo)
         break;
     }
     if (tracker || calo) {
       edm::LogVerbatim("SimG4CoreApplication") << "Pointer for Tracker " << tracker << " and for Calo " << calo;
       if (tracker)
         LogDebug("SimG4CoreApplication") << "Tracker vol name " << tracker->GetName();
       if (calo)
         LogDebug("SimG4CoreApplication") << "Calorimeter vol name " << calo->GetName();
     }
   }
 
   const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
   if (numberEkins > 0) {
     if (lvs) {
       ekinVolumes.resize(numberEkins, nullptr);
       std::vector<G4LogicalVolume*>::const_iterator lvcite;
       for (lvcite = lvs->begin(); lvcite != lvs->end(); ++lvcite) {
         for (unsigned int i = 0; i < numberEkins; ++i) {
           if ((*lvcite)->GetName() == (G4String)(ekinNames[i])) {
             ekinVolumes[i] = (*lvcite);
             break;
           }
         }
       }
     }
     for (unsigned int i = 0; i < numberEkins; ++i) {
       edm::LogVerbatim("SimG4CoreApplication") << ekinVolumes[i]->GetName() << " with pointer " << ekinVolumes[i];
     }
   }
 
   if (numberPart > 0) {
     G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
     G4String partName;
     ekinPDG.resize(numberPart, 0);
     for (unsigned int i = 0; i < numberPart; ++i) {
       ekinPDG[i] = theParticleTable->FindParticle(partName = ekinParticles[i])->GetPDGEncoding();
       edm::LogVerbatim("SimG4CoreApplication") << "Particle " << ekinParticles[i] << " with PDG code " << ekinPDG[i]
                                                << " and KE cut off " << ekinMins[i] / MeV << " MeV";
     }
   }
 
   const G4RegionStore* rs = G4RegionStore::GetInstance();
   if (numberTimes > 0) {
     maxTimeRegions.resize(numberTimes, nullptr);
     std::vector<G4Region*>::const_iterator rcite;
     for (rcite = rs->begin(); rcite != rs->end(); ++rcite) {
       for (unsigned int i = 0; i < numberTimes; ++i) {
         if ((*rcite)->GetName() == (G4String)(maxTimeNames[i])) {
           maxTimeRegions[i] = (*rcite);
           break;
         }
       }
     }
   }
   if (ndeadRegions > 0) {
     deadRegions.resize(ndeadRegions, nullptr);
     std::vector<G4Region*>::const_iterator rcite;
     for (rcite = rs->begin(); rcite != rs->end(); ++rcite) {
       for (unsigned int i = 0; i < ndeadRegions; ++i) {
         if ((*rcite)->GetName() == (G4String)(deadRegionNames[i])) {
           deadRegions[i] = (*rcite);
           break;
         }
       }
     }
   }
   return true;
 }

bool SteppingAction::isInsideDeadRegion ( const G4Region * reg ) const

inlineprivate

Definition at line 73 of file SteppingAction.h.

References deadRegions, mps_fire::i, and ndeadRegions.

Referenced by UserSteppingAction().

                                                                         {
   bool res = false;
   for (unsigned int i = 0; i < ndeadRegions; ++i) {
     if (reg == deadRegions[i]) {
       res = true;
       break;
     }
   }
   return res;
 }

bool SteppingAction::isLowEnergy ( const G4Step * aStep ) const

private

Definition at line 165 of file SteppingAction.cc.

References ekinMins, ekinPDG, ekinVolumes, RemoveAddSevLevel::flag, mps_fire::i, numberEkins, and numberPart.

Referenced by UserSteppingAction().

                                                           {
   bool flag = false;
   const G4StepPoint* sp = aStep->GetPostStepPoint();
   G4LogicalVolume* lv = sp->GetPhysicalVolume()->GetLogicalVolume();
   for (unsigned int i = 0; i < numberEkins; ++i) {
     if (lv == ekinVolumes[i]) {
       flag = true;
       break;
     }
   }
   if (flag) {
     double ekin = sp->GetKineticEnergy();
     int pCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
     for (unsigned int i = 0; i < numberPart; ++i) {
       if (pCode == ekinPDG[i]) {
         return (ekin <= ekinMins[i]) ? true : false;
       }
     }
   }
   return false;
 }

bool SteppingAction::isOutOfTimeWindow	(	G4Track *	theTrack,
		const G4Region *	reg
	)		const

inlineprivate

Definition at line 84 of file SteppingAction.h.

References mps_fire::i, maxTimeRegions, maxTrackTime, maxTrackTimes, and numberTimes.

Referenced by UserSteppingAction().

                                                                                           {
   double tofM = maxTrackTime;
   for (unsigned int i = 0; i < numberTimes; ++i) {
     if (reg == maxTimeRegions[i]) {
       tofM = maxTrackTimes[i];
       break;
     }
   }
   return (theTrack->GetGlobalTime() > tofM) ? true : false;
 }

bool SteppingAction::isThisVolume	(	const G4VTouchable *	touch,
		const G4VPhysicalVolume *	pv
	)		const

inlineprivate

Definition at line 95 of file SteppingAction.h.

References personalPlayback::level, and MetAnalyzer::pv().

Referenced by UserSteppingAction().

                                                                                                      {
   int level = (touch->GetHistoryDepth()) + 1;
   return (level >= 3) ? (touch->GetVolume(level - 3) == pv) : false;
 }

void SteppingAction::PrintKilledTrack	(	const G4Track *	aTrack,
		const TrackStatus &	tst
	)		const

private

Definition at line 266 of file SteppingAction.cc.

References MeV, MetAnalyzer::pv(), rname, sDeadRegion, sEnergyDepNaN, sKilledByProcess, sLowEnergy, sLowEnergyInVacuum, sOutOfTime, and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by UserSteppingAction().

                                                                                          {
   std::string vname = "";
   std::string rname = "";
   std::string typ = " ";
   switch (tst) {
     case sKilledByProcess:
       typ = " G4Process ";
       break;
     case sDeadRegion:
       typ = " in dead region ";
       break;
     case sOutOfTime:
       typ = " out of time window ";
       break;
     case sLowEnergy:
       typ = " low energy limit ";
       break;
     case sLowEnergyInVacuum:
       typ = " low energy limit in vacuum ";
       break;
     case sEnergyDepNaN:
       typ = " energy deposition is NaN ";
       break;
     default:
       break;
   }
   G4VPhysicalVolume* pv = aTrack->GetNextVolume();
   if (pv) {
     vname = pv->GetLogicalVolume()->GetName();
     rname = pv->GetLogicalVolume()->GetRegion()->GetName();
   }
 
   edm::LogVerbatim("SimG4CoreApplication")
       << "Track #" << aTrack->GetTrackID() << " " << aTrack->GetDefinition()->GetParticleName()
       << " E(MeV)= " << aTrack->GetKineticEnergy() / MeV << " T(ns)= " << aTrack->GetGlobalTime() / ns
       << " is killed due to " << typ << " inside LV: " << vname << " (" << rname << ") at " << aTrack->GetPosition();
 }

void SteppingAction::UserSteppingAction ( const G4Step * aStep )

final

Definition at line 85 of file SteppingAction.cc.

References EventAction::addTkCaloStateInfo(), calo, eventAction_, initialized, initPointer(), isInsideDeadRegion(), isLowEnergy(), edm::isNotFinite(), isOutOfTimeWindow(), isThisVolume(), killBeamPipe, m_g4StepSignal, MeV, CMSSteppingVerbose::NextStep(), numberEkins, nWarnings, AlCaHLTBitMon_ParallelJobs::p, PrintKilledTrack(), sAlive, sDeadRegion, sEnergyDepNaN, sKilledByProcess, sLowEnergy, sLowEnergyInVacuum, sOutOfTime, steppingVerbose, theCriticalDensity, theCriticalEnergyForVacuum, tracker, x, y, and z.

                                                            {
   if (!initialized) {
     initialized = initPointer();
   }
 
   //if(hasWatcher) { m_g4StepSignal(aStep); }
   m_g4StepSignal(aStep);
 
   G4Track* theTrack = aStep->GetTrack();
   TrackStatus tstat = (theTrack->GetTrackStatus() == fAlive) ? sAlive : sKilledByProcess;
 
   G4StepPoint* postStep = aStep->GetPostStepPoint();
 
   if (sAlive == tstat && postStep->GetPhysicalVolume() != nullptr) {
     G4StepPoint* preStep = aStep->GetPreStepPoint();
     const G4Region* theRegion = preStep->GetPhysicalVolume()->GetLogicalVolume()->GetRegion();
 
     // kill in dead regions
     if (isInsideDeadRegion(theRegion)) {
       tstat = sDeadRegion;
     }
 
     // NaN energy deposit
     if (sAlive == tstat && edm::isNotFinite(aStep->GetTotalEnergyDeposit())) {
       tstat = sEnergyDepNaN;
       if (nWarnings < 20) {
         ++nWarnings;
         edm::LogWarning("SimG4CoreApplication")
             << "Track #" << theTrack->GetTrackID() << " " << theTrack->GetDefinition()->GetParticleName()
             << " E(MeV)= " << preStep->GetKineticEnergy() / MeV
             << " is killed due to edep=NaN inside PV: " << preStep->GetPhysicalVolume()->GetName() << " at "
             << theTrack->GetPosition() << " StepLen(mm)= " << aStep->GetStepLength();
       }
     }
 
     // kill out of time
     if (sAlive == tstat && isOutOfTimeWindow(theTrack, theRegion)) {
       tstat = sOutOfTime;
     }
 
     // kill low-energy in volumes on demand
     if (sAlive == tstat && numberEkins > 0 && isLowEnergy(aStep)) {
       tstat = sLowEnergy;
     }
 
     // kill low-energy in vacuum
     G4double kinEnergy = theTrack->GetKineticEnergy();
     if (sAlive == tstat && killBeamPipe && kinEnergy < theCriticalEnergyForVacuum &&
         theTrack->GetDefinition()->GetPDGCharge() != 0.0 && kinEnergy > 0.0 &&
         theTrack->GetNextVolume()->GetLogicalVolume()->GetMaterial()->GetDensity() <= theCriticalDensity) {
       tstat = sLowEnergyInVacuum;
     }
 
     // check transition tracker/calo
     if (sAlive == tstat) {
       if (isThisVolume(preStep->GetTouchable(), tracker) && isThisVolume(postStep->GetTouchable(), calo)) {
         math::XYZVectorD pos((preStep->GetPosition()).x(), (preStep->GetPosition()).y(), (preStep->GetPosition()).z());
 
         math::XYZTLorentzVectorD mom((preStep->GetMomentum()).x(),
                                      (preStep->GetMomentum()).y(),
                                      (preStep->GetMomentum()).z(),
                                      preStep->GetTotalEnergy());
 
         uint32_t id = theTrack->GetTrackID();
 
         std::pair<math::XYZVectorD, math::XYZTLorentzVectorD> p(pos, mom);
         eventAction_->addTkCaloStateInfo(id, p);
       }
     } else {
       theTrack->SetTrackStatus(fStopAndKill);
 #ifdef DebugLog
       PrintKilledTrack(theTrack, tstat);
 #endif
     }
   }
   if (nullptr != steppingVerbose) {
     steppingVerbose->NextStep(aStep, fpSteppingManager, (1 < tstat));
   }
 }