KillSecondaries/src/KillSecondariesTrackAction.cc

0001 #include "SimG4Core/KillSecondaries/interface/KillSecondariesTrackAction.h"
0002
0003 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0004 #include "SimG4Core/Notification/interface/TrackInformation.h"
0005
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include "G4Track.hh"
0009
0010 using namespace CLHEP;
0011
0012 KillSecondariesTrackAction::KillSecondariesTrackAction(edm::ParameterSet const &p) {
0013   killHeavy = p.getParameter<bool>("KillHeavy");
0014   kmaxIon = p.getParameter<double>("IonThreshold") * MeV;
0015   kmaxProton = p.getParameter<double>("ProtonThreshold") * MeV;
0016   kmaxNeutron = p.getParameter<double>("NeutronThreshold") * MeV;
0017
0018   edm::LogInfo("KillSecondaries") << "KillSecondariesTrackAction:: Killing"
0019                                   << " Flag " << killHeavy << " protons below " << kmaxProton << " MeV, neutrons below "
0020                                   << kmaxNeutron << " MeV and ions below " << kmaxIon << " MeV\n";
0021 }
0022
0023 KillSecondariesTrackAction::~KillSecondariesTrackAction() {}
0024
0025 void KillSecondariesTrackAction::update(const BeginOfTrack *trk) {
0026   if (killHeavy) {
0027     G4Track *theTrack = (G4Track *)((*trk)());
0028     TrackInformation *trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
0029     if (trkInfo) {
0030       int pdg = theTrack->GetDefinition()->GetPDGEncoding();
0031       if (!(trkInfo->isPrimary())) {  // Only secondary particles
0032         double ke = theTrack->GetKineticEnergy() / MeV;
0033         if ((((pdg / 1000000000 == 1 && ((pdg / 10000) % 100) > 0 && ((pdg / 10) % 100) > 0)) && (ke < kmaxIon)) ||
0034             ((pdg == 2212) && (ke < kmaxProton)) || ((pdg == 2112) && (ke < kmaxNeutron))) {
0035           theTrack->SetTrackStatus(fStopAndKill);
0036           edm::LogInfo("KillSecondaries")
0037               << "Kill Track " << theTrack->GetTrackID() << " Type " << theTrack->GetDefinition()->GetParticleName()
0038               << " Kinetic Energy " << ke << " MeV";
0039         }
0040       }
0041     }
0042   }
0043 }