#include <SteppingAction.h>

Public Member Functions
	SteppingAction (EventAction *ea, const edm::ParameterSet &ps)
void	UserSteppingAction (const G4Step *aStep)
	~SteppingAction ()
Public Attributes
SimActivityRegistry::G4StepSignal	m_g4StepSignal
Private Member Functions
bool	catchLongLived (const G4Step *aStep)
void	catchLowEnergyInVacuumHere (const G4Step *aStep)
void	catchLowEnergyInVacuumNext (const G4Step *aStep)
bool	initPointer ()
bool	isThisVolume (const G4VTouchable touch, G4VPhysicalVolume pv)
bool	killLowEnergy (const G4Step *aStep)
void	killTrack (const G4Step *aStep)
Private Attributes
G4VPhysicalVolume *	calo
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	initialized
bool	killBeamPipe
std::vector< std::string >	maxTimeNames
std::vector< G4Region * >	maxTimeRegions
double	maxTrackTime
std::vector< double >	maxTrackTimes
double	theCriticalDensity
double	theCriticalEnergyForVacuum
G4VPhysicalVolume *	tracker
int	verbose

Detailed Description

Definition at line 18 of file SteppingAction.h.

Constructor & Destructor Documentation

SteppingAction::SteppingAction	(	EventAction *	ea,
		const edm::ParameterSet &	ps
	)

Definition at line 13 of file SteppingAction.cc.

References ekinMins, ekinNames, ekinParticles, g, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), i, killBeamPipe, maxTimeNames, maxTrackTime, maxTrackTimes, theCriticalDensity, and theCriticalEnergyForVacuum.

  : eventAction_(e), initialized(false), tracker(0), calo(0) {

  killBeamPipe = (p.getParameter<bool>("KillBeamPipe"));
  theCriticalEnergyForVacuum = (p.getParameter<double>("CriticalEnergyForVacuum")*MeV);
  theCriticalDensity = (p.getParameter<double>("CriticalDensity")*g/cm3);
  maxTrackTime  = (p.getParameter<double>("MaxTrackTime")*ns);
  maxTrackTimes = (p.getParameter<std::vector<double> >("MaxTrackTimes"));
  maxTimeNames  = (p.getParameter<std::vector<std::string> >("MaxTimeNames"));
  ekinMins      = (p.getParameter<std::vector<double> >("EkinThresholds"));
  ekinNames     = (p.getParameter<std::vector<std::string> >("EkinNames"));
  ekinParticles = (p.getParameter<std::vector<std::string> >("EkinParticles"));
  verbose       = (p.getUntrackedParameter<int>("Verbosity",0));

  edm::LogInfo("SimG4CoreApplication") << "SteppingAction:: KillBeamPipe = "
                                       << killBeamPipe << " CriticalDensity = "
                                       << theCriticalDensity << " g/cc;"
                                       << " CriticalEnergyForVacuum = "
                                       << theCriticalEnergyForVacuum << " Mev;"
                                       << " MaxTrackTime = " << maxTrackTime 
                                       << " ns";
  for (unsigned int i=0; i<maxTrackTimes.size(); i++) {
    maxTrackTimes[i] *= ns;
    edm::LogInfo("SimG4CoreApplication") << "SteppingAction::MaxTrackTime for "
                                         << maxTimeNames[i] << " is " 
                                         << maxTrackTimes[i];
  }
  edm::LogInfo("SimG4CoreApplication") << "SteppingAction::Kill following "
                                       << ekinParticles.size() 
                                       << " particles in " << ekinNames.size()
                                       << " volumes";
  for (unsigned int i=0; i<ekinParticles.size(); i++) {
    ekinMins[i] *= GeV;
    edm::LogInfo("SimG4CoreApplication") << "SteppingAction::Particle " << i
                                         << " " << ekinParticles[i]
                                         << " (Threshold = " << ekinMins[i]
                                         << " MeV)";
  }
  for (unsigned int i=0; i<ekinNames.size(); i++) 
    edm::LogInfo("SimG4CoreApplication") << "SteppingAction::Volume[" << i
                                         << "] = " << ekinNames[i];
}

SteppingAction::~SteppingAction ( )

Definition at line 56 of file SteppingAction.cc.

{}

Member Function Documentation

bool SteppingAction::catchLongLived ( const G4Step * aStep ) [private]

Definition at line 133 of file SteppingAction.cc.

References i, killTrack(), maxTimeRegions, maxTrackTime, maxTrackTimes, and cond::rpcobgas::time.

Referenced by UserSteppingAction().

                                                        {

  bool flag   = true;
  double time = (aStep->GetPostStepPoint()->GetGlobalTime())/nanosecond;
  double tofM = maxTrackTime;
  if (maxTimeRegions.size() > 0) {
    G4Region* reg = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetRegion();
    for (unsigned int i=0; i<maxTimeRegions.size(); i++) {
      if (reg == maxTimeRegions[i]) {
        tofM = maxTrackTimes[i];
        break;
      }
    }
  }
  if (time > tofM) {
    killTrack(aStep);
    flag = false;
  }
  return flag;
}

void SteppingAction::catchLowEnergyInVacuumHere ( const G4Step * aStep ) [private]

Definition at line 92 of file SteppingAction.cc.

References g.

Referenced by UserSteppingAction().

                                                                    {
  G4Track * theTrack = aStep->GetTrack();
  double theKenergy = theTrack->GetKineticEnergy();
  if (theTrack->GetVolume()!=0) {
    double density = theTrack->GetVolume()->GetLogicalVolume()->GetMaterial()->GetDensity();
    if (theKenergy <= theCriticalEnergyForVacuum && theKenergy > 0.0 &&
        density <= theCriticalDensity && theTrack->GetDefinition()->GetPDGCharge() != 0 &&
        theTrack->GetTrackStatus() != fStopAndKill) {
      if (verbose>1)
        edm::LogInfo("SimG4CoreApplication") 
          <<   " SteppingAction: LoopCatchSteppingAction:catchLowEnergyInVacuumHere: "
          << " Track from " << theTrack->GetDefinition()->GetParticleName()
          << " of kinetic energy " << theKenergy/MeV << " MeV "
          << " killed in " << theTrack->GetVolume()->GetLogicalVolume()->GetName()
          << " of density " << density/(g/cm3) << " g/cm3" ;
      theTrack->SetTrackStatus(fStopAndKill);
    }
  }
}

void SteppingAction::catchLowEnergyInVacuumNext ( const G4Step * aStep ) [private]

Definition at line 112 of file SteppingAction.cc.

References g.

Referenced by UserSteppingAction().

                                                                    {
  G4Track * theTrack = aStep->GetTrack();
  double theKenergy = theTrack->GetKineticEnergy();
  if (theTrack->GetNextVolume()) {
    double density = theTrack->GetNextVolume()->GetLogicalVolume()->GetMaterial()->GetDensity();
    if (theKenergy <=  theCriticalEnergyForVacuum && theKenergy > 0.0 &&
        density <= theCriticalDensity && theTrack->GetDefinition()->GetPDGCharge() != 0 &&
        theTrack->GetTrackStatus() != fStopAndKill) {
      if (verbose>1)
        edm::LogInfo("SimG4CoreApplication") 
          << " SteppingAction: LoopCatchSteppingAction::catchLowEnergyInVacuumNext: "
          << " Track from " << theTrack->GetDefinition()->GetParticleName()
          << " of kinetic energy " << theKenergy/MeV << " MeV "
          << " stopped in " << theTrack->GetVolume()->GetLogicalVolume()->GetName()
          << " before going into "<< theTrack->GetNextVolume()->GetLogicalVolume()->GetName()
          << " of density " << density/(g/cm3) << " g/cm3" ;
      theTrack->SetTrackStatus(fStopButAlive);
    }
  }
}

bool SteppingAction::initPointer ( ) [private]

Definition at line 186 of file SteppingAction.cc.

References calo, ekinMins, ekinNames, ekinParticles, ekinPDG, ekinVolumes, i, LogDebug, maxTimeNames, maxTimeRegions, maxTrackTimes, mergeVDriftHistosByStation::name, and tracker.

Referenced by UserSteppingAction().

                                 {

  bool flag = true;
  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;
    }
    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();
  } else {
    flag = false;
  }

  const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
  unsigned int num = ekinNames.size();
  if (num > 0) {
    if (lvs) {
      std::vector<G4LogicalVolume*>::const_iterator lvcite;
      for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
        for (unsigned int i=0; i<num; i++) {
          if ((*lvcite)->GetName() == (G4String)(ekinNames[i])) {
            ekinVolumes.push_back(*lvcite);
            break;
          }
        }
        if (ekinVolumes.size() == num) break;
      }
    }
    if (ekinVolumes.size() != num) flag = false;
    for (unsigned int i=0; i<ekinVolumes.size(); i++) {
      edm::LogInfo("SimG4CoreApplication") << ekinVolumes[i]->GetName()
                                           <<" with pointer " <<ekinVolumes[i];
    }
  }

  G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
  G4String particleName;
  for (unsigned int i=0; i<ekinParticles.size(); i++) {
    int pdg = theParticleTable->FindParticle(particleName=ekinParticles[i])->GetPDGEncoding();
    ekinPDG.push_back(pdg);
    edm::LogInfo("SimG4CoreApplication") << "Particle " << ekinParticles[i]
                                         << " with code " << ekinPDG[i]
                                         << " and KE cut off " << ekinMins[i];
  }
  if (!flag) edm::LogInfo("SimG4CoreApplication") << "SteppingAction fails to"
                                                  << " initialize some the "
                                                  << "LV pointers correctly";

  const G4RegionStore * rs = G4RegionStore::GetInstance();
  num = maxTimeNames.size();
  if (num > 0) {
    std::vector<double> tofs;
    if (rs) {
      std::vector<G4Region*>::const_iterator rcite;
      for (rcite = rs->begin(); rcite != rs->end(); rcite++) {
        for (unsigned int i=0; i<num; i++) {
          if ((*rcite)->GetName() == (G4String)(maxTimeNames[i])) {
            maxTimeRegions.push_back(*rcite);
            tofs.push_back(maxTrackTimes[i]);
            break;
          }
        }
        if (tofs.size() == num) break;
      }
    }
    for (unsigned int i=0; i<tofs.size(); i++) {
      maxTrackTimes[i] = tofs[i];
      G4String name = "Unknown";
      if (maxTimeRegions[i]) name = maxTimeRegions[i]->GetName();
      edm::LogInfo("SimG4CoreApplication") << name << " with pointer " 
                                           << maxTimeRegions[i]<<" KE cut off "
                                           << maxTrackTimes[i];
    }
    if (tofs.size() != num) 
      edm::LogInfo("SimG4CoreApplication") << "SteppingAction fails to "
                                           << "initialize some the region "
                                           << "pointers correctly";
  }
  return true;
}

bool SteppingAction::isThisVolume	(	const G4VTouchable *	touch,
		G4VPhysicalVolume *	pv
	)		`[private]`

Definition at line 275 of file SteppingAction.cc.

References testEve_cfg::level.

Referenced by UserSteppingAction().

                                                         {

  int level = ((touch->GetHistoryDepth())+1);
  if (level > 0 && level >= 3) {
    unsigned int ii = (unsigned int)(level - 3);
    return (touch->GetVolume(ii) == pv);
  }
  return false;
}

bool SteppingAction::killLowEnergy ( const G4Step * aStep ) [private]

Definition at line 154 of file SteppingAction.cc.

References ekinMins, ekinPDG, ekinVolumes, i, killTrack(), and convertSQLiteXML::ok.

                                                       {

  bool ok = true;
  if (ekinVolumes.size() > 0) {
    bool flag = false;
    G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
    for (unsigned int i=0; i<ekinVolumes.size(); i++) {
      if (lv == ekinVolumes[i]) {
        flag = true;
        break;
      }
    }
    if (flag) {
      G4Track * track = aStep->GetTrack();
      double    ekin  = track->GetKineticEnergy();
      double    ekinM = 0;
      int       pCode = track->GetDefinition()->GetPDGEncoding();
      for (unsigned int i=0; i<ekinPDG.size(); i++) {
        if (pCode == ekinPDG[i]) {
          ekinM = ekinMins[i];
          break;
        }
      }
      if (ekin > ekinM) {
        killTrack(aStep);
        ok = false;
      }
    }
  }
  return ok;
}

void SteppingAction::killTrack ( const G4Step * aStep ) [private]

Definition at line 286 of file SteppingAction.cc.

References GetVolume(), GetRecoTauVFromDQM_MC_cff::kk, and LogDebug.

Referenced by catchLongLived(), and killLowEnergy().

                                                   {
  
  aStep->GetTrack()->SetTrackStatus(fStopAndKill);
  G4TrackVector tv = *(aStep->GetSecondary());
  for (unsigned int kk=0; kk<tv.size(); kk++) {
    if (tv[kk]->GetVolume() == aStep->GetPreStepPoint()->GetPhysicalVolume())
      tv[kk]->SetTrackStatus(fStopAndKill);
  }
  LogDebug("SimG4CoreApplication") << "SteppingAction: Kills track " 
                                   << aStep->GetTrack()->GetTrackID() << " ("
                                   << aStep->GetTrack()->GetDefinition()->GetParticleName()
                                   << ") at " << aStep->GetPostStepPoint()->GetGlobalTime()/nanosecond 
                                   << " ns in " << aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName()
                                   << " from region " << aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetRegion()->GetName();
}

void SteppingAction::UserSteppingAction ( const G4Step * aStep )

Definition at line 58 of file SteppingAction.cc.

References EventAction::addTkCaloStateInfo(), calo, catchLongLived(), catchLowEnergyInVacuumHere(), catchLowEnergyInVacuumNext(), eventAction_, initialized, initPointer(), isThisVolume(), killBeamPipe, m_g4StepSignal, convertSQLiteXML::ok, AlCaHLTBitMon_ParallelJobs::p, pos, tracker, x, detailsBasic3DVector::y, and z.

                                                            {
  if (!initialized) initialized = initPointer();
  m_g4StepSignal(aStep);

  if (killBeamPipe) {
    catchLowEnergyInVacuumHere(aStep);
    catchLowEnergyInVacuumNext(aStep);
  }

  if (aStep->GetPostStepPoint()->GetPhysicalVolume() != 0) {
    bool ok = catchLongLived(aStep);
    if (ok) ok = catchLongLived(aStep);

    ok = (isThisVolume(aStep->GetPreStepPoint()->GetTouchable(),tracker)&&
          isThisVolume(aStep->GetPostStepPoint()->GetTouchable(),calo));
    if (ok) {

      math::XYZVectorD pos((aStep->GetPreStepPoint()->GetPosition()).x(),
                           (aStep->GetPreStepPoint()->GetPosition()).y(),
                           (aStep->GetPreStepPoint()->GetPosition()).z());
      
      math::XYZTLorentzVectorD mom((aStep->GetPreStepPoint()->GetMomentum()).x(),
                                   (aStep->GetPreStepPoint()->GetMomentum()).y(),
                                   (aStep->GetPreStepPoint()->GetMomentum()).z(),
                                   aStep->GetPreStepPoint()->GetTotalEnergy());
      
      uint32_t id = aStep->GetTrack()->GetTrackID();
      
      std::pair<math::XYZVectorD,math::XYZTLorentzVectorD> p(pos,mom);
      eventAction_->addTkCaloStateInfo(id,p);
    }
  }
}