#include <StackingAction.h>

Inheritance diagram for StackingAction:

Public Member Functions
virtual G4ClassificationOfNewTrack	ClassifyNewTrack (const G4Track *aTrack)

virtual void	NewStage ()

virtual void	PrepareNewEvent ()

	StackingAction (const edm::ParameterSet &ps)

virtual	~StackingAction ()

Private Member Functions
void	initPointer ()

int	isItFromPrimary (const G4Track &, int) const

bool	isItLongLived (const G4Track *) const

int	isItPrimaryDecayProductOrConversion (const G4Track *, const G4Track &) const

bool	isThisVolume (const G4VTouchable , std::vector< G4LogicalVolume > &) const

Private Attributes
std::vector< G4LogicalVolume * >	calo

bool	killDeltaRay

bool	killHeavy

double	kmaxIon

double	kmaxNeutron

double	kmaxProton

std::vector< std::string >	maxTimeNames

std::vector< G4Region * >	maxTimeRegions

double	maxTrackTime

std::vector< double >	maxTrackTimes

std::vector< G4LogicalVolume * >	muon

bool	saveFirstSecondary

bool	savePDandCinCalo

bool	savePDandCinMuon

bool	savePDandCinTracker

std::vector< G4LogicalVolume * >	tracker

bool	trackNeutrino

Detailed Description

Definition at line 14 of file StackingAction.h.

Constructor & Destructor Documentation

StackingAction::StackingAction ( const edm::ParameterSet & ps )

Definition at line 18 of file StackingAction.cc.

References edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), i, initPointer(), killDeltaRay, killHeavy, kmaxIon, kmaxNeutron, kmaxProton, maxTimeNames, maxTrackTime, maxTrackTimes, saveFirstSecondary, savePDandCinCalo, savePDandCinMuon, savePDandCinTracker, and trackNeutrino.

                                                         {
 
   trackNeutrino  = p.getParameter<bool>("TrackNeutrino");
   killHeavy      = p.getParameter<bool>("KillHeavy");
   kmaxIon        = p.getParameter<double>("IonThreshold")*MeV;
   kmaxProton     = p.getParameter<double>("ProtonThreshold")*MeV;
   kmaxNeutron    = p.getParameter<double>("NeutronThreshold")*MeV;
   killDeltaRay   = p.getParameter<bool>("KillDeltaRay");
   maxTrackTime   = p.getParameter<double>("MaxTrackTime")*ns;
   maxTrackTimes  = p.getParameter<std::vector<double> >("MaxTrackTimes");
   maxTimeNames   = p.getParameter<std::vector<std::string> >("MaxTimeNames");
   savePDandCinTracker = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInTracker",false);
   savePDandCinCalo    = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInCalo",false);
   savePDandCinMuon    = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInMuon",false);
   saveFirstSecondary  = p.getUntrackedParameter<bool>("SaveFirstLevelSecondary",false);
 
   edm::LogInfo("SimG4CoreApplication") << "StackingAction initiated with"
                        << " flag for saving decay products in "
                        << " Tracker: " << savePDandCinTracker
                                        << " in Calo: " << savePDandCinCalo
                                        << " in Muon: " << savePDandCinMuon
                        << "\n               saveFirstSecondary"
                        << ": " << saveFirstSecondary
                        << " Flag for tracking neutrino: "
                        << trackNeutrino << " Killing Flag "
                        << killHeavy << " protons below " 
                        << kmaxProton <<" MeV, neutrons below "
                        << kmaxNeutron << " MeV and ions"
                        << " below " << kmaxIon << " MeV\n"
                        << "               kill tracks with "
                        << "time larger than " << maxTrackTime
                        << " ns and kill Delta Ray flag set to "
                        << killDeltaRay;
   for (unsigned int i=0; i<maxTrackTimes.size(); i++) {
     maxTrackTimes[i] *= ns;
     edm::LogInfo("SimG4CoreApplication") << "StackingAction::MaxTrackTime for "
                      << maxTimeNames[i] << " is " 
                      << maxTrackTimes[i];
   }
   initPointer();
 }

StackingAction::~StackingAction ( )

virtual

Definition at line 60 of file StackingAction.cc.

60 {}

Member Function Documentation

G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack ( const G4Track * aTrack )

virtual

Definition at line 62 of file StackingAction.cc.

References abs, calo, isItFromPrimary(), isItLongLived(), isItPrimaryDecayProductOrConversion(), isThisVolume(), ke, killDeltaRay, killHeavy, kmaxIon, kmaxNeutron, kmaxProton, LogDebug, muon, NewTrackAction::primary(), saveFirstSecondary, savePDandCinCalo, savePDandCinMuon, savePDandCinTracker, NewTrackAction::secondary(), CurrentG4Track::track(), tracker, and trackNeutrino.

                                                                                   {
 
   // G4 interface part
   G4ClassificationOfNewTrack classification = fUrgent;
   int flag = 0;
 
   NewTrackAction newTA;
   if (aTrack->GetCreatorProcess()==0 || aTrack->GetParentID()==0) {
     newTA.primary(aTrack);
   } else if (aTrack->GetTouchable() == 0) {
     edm::LogError("SimG4CoreApplication")
       << "StackingAction: no touchable for track " << aTrack->GetTrackID()
       << " from " << aTrack->GetParentID()
       << " with PDG code " << aTrack->GetDefinition()->GetParticleName();
     classification = fKill;
   } else {
     const G4Track * mother = CurrentG4Track::track();
     if ((savePDandCinTracker && isThisVolume(aTrack->GetTouchable(),tracker))||
     (savePDandCinCalo && isThisVolume(aTrack->GetTouchable(),calo)) ||
     (savePDandCinMuon && isThisVolume(aTrack->GetTouchable(),muon)))
       flag = isItPrimaryDecayProductOrConversion(aTrack, *mother);
     if (saveFirstSecondary) flag = isItFromPrimary(*mother, flag);
     newTA.secondary(aTrack, *mother, flag);
 
     if (aTrack->GetTrackStatus() == fStopAndKill) classification = fKill;
     if (killHeavy) {
       int    pdg = aTrack->GetDefinition()->GetPDGEncoding();
       double ke  = aTrack->GetKineticEnergy()/MeV;
       if (((pdg/1000000000 == 1) && (((pdg/10000)%100) > 0) && 
        (((pdg/10)%100) > 0) && (ke<kmaxIon)) || 
       ((pdg == 2212) && (ke < kmaxProton)) ||
       ((pdg == 2112) && (ke < kmaxNeutron))) classification = fKill;
     }
     if (!trackNeutrino) {
       int    pdg = std::abs(aTrack->GetDefinition()->GetPDGEncoding());
       if (pdg == 12 || pdg == 14 || pdg == 16 || pdg == 18) 
     classification = fKill;
     }
     if (isItLongLived(aTrack)) classification = fKill;
     if (killDeltaRay) {
       if (aTrack->GetCreatorProcess()->GetProcessType() == fElectromagnetic &&
       aTrack->GetCreatorProcess()->GetProcessSubType() == fIonisation)
     classification = fKill;
     }
 #ifdef DebugLog
     LogDebug("SimG4CoreApplication") << "StackingAction:Classify Track "
                      << aTrack->GetTrackID() << " Parent " 
                      << aTrack->GetParentID() << " Type "
                      << aTrack->GetDefinition()->GetParticleName() 
                      << " K.E. " << aTrack->GetKineticEnergy()/MeV
                      << " MeV from process/subprocess " 
                      << aTrack->GetCreatorProcess()->GetProcessType() << "|"
                      << aTrack->GetCreatorProcess()->GetProcessSubType()
                      << " as " << classification << " Flag " << flag;
 #endif
   }
   return classification;
 }

void StackingAction::initPointer ( )

private

Definition at line 125 of file StackingAction.cc.

References calo, i, maxTimeNames, maxTimeRegions, maxTrackTimes, muon, mergeVDriftHistosByStation::name, savePDandCinCalo, savePDandCinMuon, savePDandCinTracker, and tracker.

Referenced by StackingAction().

                                  {
 
   const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
   if (lvs) {
     std::vector<G4LogicalVolume*>::const_iterator lvcite;
     for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
       if (savePDandCinTracker) {
         if (strcmp("Tracker",(*lvcite)->GetName().c_str()) == 0) tracker.push_back(*lvcite);
         if (strcmp("BEAM",(*lvcite)->GetName().substr(0,4).c_str()) == 0) tracker.push_back(*lvcite);
       }
       if (savePDandCinCalo) {
         if (strcmp("CALO",(*lvcite)->GetName().c_str()) == 0) calo.push_back(*lvcite);
         if (strcmp("VCAL",(*lvcite)->GetName().c_str()) == 0) calo.push_back(*lvcite);
       }
       if (savePDandCinMuon) {
         if (strcmp("MUON",(*lvcite)->GetName().c_str()) == 0) muon.push_back(*lvcite);
       }
     }
     edm::LogInfo("SimG4CoreApplication") << "# of LV for Tracker " 
                      << tracker.size() << " for Calo "
                                          << calo.size() << " for Muon "
                      << muon.size();
     for (unsigned int i=0; i<tracker.size(); ++i)
       edm::LogInfo("SimG4CoreApplication") << "Tracker vol " << i << " name "
                        << tracker[i]->GetName();
     for (unsigned int i=0; i<calo.size(); ++i)
       edm::LogInfo("SimG4CoreApplication") << "Calorimeter vol " <<i <<" name "
                        << calo[i]->GetName();
     for (unsigned int i=0; i<muon.size(); ++i)
       edm::LogInfo("SimG4CoreApplication") << "Muon vol " << i << " name "
                        << muon[i]->GetName();
   }
 
   const G4RegionStore * rs = G4RegionStore::GetInstance();
   unsigned int 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];
     }
   }
 
 }

int StackingAction::isItFromPrimary	(	const G4Track &	mother,
		int	flagIn
	)		const

private

Definition at line 216 of file StackingAction.cc.

References TrackInformation::isPrimary().

Referenced by ClassifyNewTrack().

                                                                             {
 
   int flag = flagIn;
   if (flag != 1) {
     TrackInformationExtractor extractor;
     const TrackInformation & motherInfo(extractor(mother));
     if (motherInfo.isPrimary()) flag = 3;
   }
   return flag;
 }

bool StackingAction::isItLongLived ( const G4Track * aTrack ) const

private

Definition at line 227 of file StackingAction.cc.

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

Referenced by ClassifyNewTrack().

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

int StackingAction::isItPrimaryDecayProductOrConversion	(	const G4Track *	aTrack,
		const G4Track &	mother
	)		const

private

Definition at line 201 of file StackingAction.cc.

References TrackInformation::isPrimary().

Referenced by ClassifyNewTrack().

                                                           {
 
   int flag = 0;
   TrackInformationExtractor extractor;
   const TrackInformation & motherInfo(extractor(mother));
   // Check whether mother is a primary
   if (motherInfo.isPrimary()) {
     if (aTrack->GetCreatorProcess()->GetProcessType() == fDecay) flag = 1;
     else if (aTrack->GetCreatorProcess()->GetProcessType() == fElectromagnetic &&
          aTrack->GetCreatorProcess()->GetProcessSubType() == fGammaConversion) flag = 2;
   }
   return flag;
 }

bool StackingAction::isThisVolume	(	const G4VTouchable *	touch,
		std::vector< G4LogicalVolume * > &	lvs
	)		const

private

Definition at line 187 of file StackingAction.cc.

References prof2calltree::count, and testEve_cfg::level.

Referenced by ClassifyNewTrack().

                                                            {
 
   bool flag = false;
   if (lvs.size() > 0 && touch !=0) {
     int level = ((touch->GetHistoryDepth())+1);
     if (level >= 3) {
       unsigned int  ii = (unsigned int)(level - 3);
       flag    = (std::count(lvs.begin(),lvs.end(),(touch->GetVolume(ii)->GetLogicalVolume())) != 0);
     }
   }
   return flag;
 }