#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 84 of file SteppingAction.cc.

84 {}

Member Function Documentation

bool SteppingAction::initPointer ( )

private

Definition at line 193 of file SteppingAction.cc.

References calo, deadRegionNames, deadRegions, ekinMins, ekinNames, ekinParticles, ekinPDG, ekinVolumes, mps_fire::i, LogDebug, maxTimeNames, maxTimeRegions, MeV, ndeadRegions, numberEkins, numberPart, numberTimes, 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 76 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 170 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 89 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 101 of file SteppingAction.h.

References hcalDigis_cfi::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 274 of file SteppingAction.cc.

References MeV, MetAnalyzer::pv(), 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 86 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)); 
   }
 }