#include <SteppingAction.h>

Inheritance diagram for SteppingAction:

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

virtual void	UserSteppingAction (const G4Step *aStep) final

virtual	~SteppingAction ()

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

const CMSSteppingVerbose *	steppingVerbose

double	theCriticalDensity

double	theCriticalEnergyForVacuum

const G4VPhysicalVolume *	tracker

Detailed Description

Definition at line 29 of file SteppingAction.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 17 of file SteppingAction.cc.

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

   : eventAction_(e), tracker(nullptr), calo(nullptr), steppingVerbose(sv),
     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::LogInfo("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::LogInfo("SimG4CoreApplication") << "SteppingAction::MaxTrackTime for "
                        << maxTimeNames[i] << " is " 
                        << maxTrackTimes[i] << " ns ";
       maxTrackTimes[i] *= ns;
     }
   }
 
   ndeadRegions =  deadRegionNames.size();
   if(ndeadRegions > 0) {
     edm::LogInfo("SimG4CoreApplication") 
       << "SteppingAction: Number of DeadRegions where all trackes are killed "
       << ndeadRegions;
     for(unsigned int i=0; i<ndeadRegions; ++i) {
       edm::LogInfo("SimG4CoreApplication") 
     << "SteppingAction: DeadRegion " << i << ".  " << deadRegionNames[i];
     }
   }
   numberEkins = ekinNames.size();
   numberPart  = ekinParticles.size();
   if(0 == numberPart) { numberEkins = 0; }
 
   if(numberEkins > 0) {
 
     edm::LogInfo("SimG4CoreApplication") << "SteppingAction::Kill following "
                      << numberPart 
                      << " particles in " << numberEkins
                      << " volumes";
     for (unsigned int i=0; i<numberPart; ++i) {
       edm::LogInfo("SimG4CoreApplication") << "SteppingAction::Particle " << i
                        << " " << ekinParticles[i]
                        << "  Threshold = " << ekinMins[i]
                        << " MeV";
       ekinMins[i] *= CLHEP::MeV;
     }
     for (unsigned int i=0; i<numberEkins; ++i) {
       edm::LogInfo("SimG4CoreApplication") << "SteppingAction::LogVolume[" << i
                        << "] = " << ekinNames[i];
     }
   }
 }

SteppingAction::~SteppingAction ( )

virtual

Definition at line 90 of file SteppingAction.cc.

90 {}

Member Function Documentation

bool SteppingAction::initPointer ( )

private

Definition at line 182 of file SteppingAction.cc.

References calo, deadRegionNames, deadRegions, ekinMins, ekinNames, ekinParticles, ekinPDG, ekinVolumes, 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::LogInfo("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, 0);
       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::LogInfo("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::LogInfo("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, 0);
     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, 0);
     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 74 of file SteppingAction.h.

References deadRegions, 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 159 of file SteppingAction.cc.

References ekinMins, ekinPDG, ekinVolumes, RemoveAddSevLevel::flag, 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 87 of file SteppingAction.h.

References 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 99 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 264 of file SteppingAction.cc.

References MeV, MetAnalyzer::pv(), sDeadRegion, 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;
   default:
     break;
   } 
   G4VPhysicalVolume* pv = aTrack->GetNextVolume();
   if(pv) { 
     vname = pv->GetLogicalVolume()->GetName(); 
     rname = pv->GetLogicalVolume()->GetRegion()->GetName(); 
   }
 
   edm::LogInfo("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 )

finalvirtual

Definition at line 92 of file SteppingAction.cc.

References EventAction::addTkCaloStateInfo(), calo, eventAction_, initialized, initPointer(), isInsideDeadRegion(), isLowEnergy(), isOutOfTimeWindow(), isThisVolume(), killBeamPipe, m_g4StepSignal, CMSSteppingVerbose::NextStep(), numberEkins, AlCaHLTBitMon_ParallelJobs::p, PrintKilledTrack(), sAlive, sDeadRegion, 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(0 == tstat && postStep->GetPhysicalVolume() != nullptr) {
 
     G4StepPoint* preStep = aStep->GetPreStepPoint();
     const G4Region* theRegion = 
       preStep->GetPhysicalVolume()->GetLogicalVolume()->GetRegion();
 
     // kill in dead regions
     if(isInsideDeadRegion(theRegion)) { tstat = sDeadRegion; }
 
     // kill out of time
     if(0 == tstat && isOutOfTimeWindow(theTrack, theRegion)) { tstat = sOutOfTime; }
 
     // kill low-energy in volumes on demand
     if(0 == tstat && numberEkins > 0 && isLowEnergy(aStep)) { tstat = sLowEnergy; }
 
     // kill low-energy in vacuum
     G4double kinEnergy = theTrack->GetKineticEnergy();
     if(0 == 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(0 == 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)); 
   }
 }