gar::garg4::ParticleListAction Class Reference

#include <ParticleListAction.h>

Inheritance diagram for gar::garg4::ParticleListAction:

Classes
struct	ParticleInfo_t

Public Member Functions
	ParticleListAction (double energyCut, bool storeTrajectories=false, bool keepEMShowerDaughters=false, std::string EMShowerDaughterMatRegex=".*")
virtual	~ParticleListAction ()
virtual void	BeginOfEventAction (const G4Event *)
virtual void	EndOfEventAction (const G4Event *)
virtual void	PreTrackingAction (const G4Track *)
virtual void	PostTrackingAction (const G4Track *)
virtual void	SteppingAction (const G4Step *)
void	ParticleFilter (std::unique_ptr< PositionInVolumeFilter > &&filter)
	Grabs a particle filter. More...
void	ResetTrackIDOffset ()
std::map< int, size_t >	TrackIDToMCTruthIndexMap () const
sim::ParticleList *	GetList () const
sim::ParticleList &&	YieldList ()

Static Public Member Functions
static int	GetCurrentTrackID ()
static int	GetTrackIDOffset ()
static bool	IsDropped (simb::MCParticle const *p)
	returns whether the specified particle has been marked as dropped More...

Private Member Functions
int	GetParentage (int trackid) const
void	AddPointToCurrentParticle (TLorentzVector const &pos, TLorentzVector const &mom, std::string const &process)
	Adds a trajectory point to the current particle, and runs the filter. More...

Private Attributes
G4double	fEnergyCut
ParticleInfo_t	fCurrentParticle
sim::ParticleList *	fParticleList
G4bool	fstoreTrajectories
	Whether to store particle trajectories with each particle. More...
std::map< int, int >	fParentIDMap
	key is current track ID, value is parent ID More...
std::map< int, size_t >	fTrackIDToMCTruthIndex
	map track ID to index of MCTruth in input list More...
bool	fKeepEMShowerDaughters
	whether to keep EM shower secondaries, tertiaries, etc More...
std::string	fEMShowerDaughterMatRegex
	if keeping EM shower daughters, save only in media matching this More...
std::unique_ptr< PositionInVolumeFilter >	fFilter
	filter for particles to be kept More...

Static Private Attributes
static int	fCurrentTrackID = sdp::NoParticleId
static int	fTrackIDOffset = 0

Detailed Description

Definition at line 39 of file ParticleListAction.h.

Constructor & Destructor Documentation

gar::garg4::ParticleListAction::ParticleListAction	(	double	energyCut,
		bool	storeTrajectories = false,
		bool	keepEMShowerDaughters = false,
		std::string	EMShowerDaughterMatRegex = ".*"
	)

Definition at line 55 of file ParticleListAction.cxx.

      : fEnergyCut            (energyCut)
      , fParticleList         (new sim::ParticleList())
      , fstoreTrajectories    (storeTrajectories)
      , fKeepEMShowerDaughters(keepEMShowerDaughters)
      , fEMShowerDaughterMatRegex(EMShowerDaughterMatRegex)
    {
      fParentIDMap.clear();
    }

gar::garg4::ParticleListAction::~ParticleListAction ( )

virtual

Definition at line 70 of file ParticleListAction.cxx.

    {
      // Delete anything that we created with "new'.
      delete fParticleList;
    }

Member Function Documentation

void gar::garg4::ParticleListAction::AddPointToCurrentParticle ( TLorentzVector const &  pos, TLorentzVector const &  mom, std::string const &  process  )

private

Adds a trajectory point to the current particle, and runs the filter.

Definition at line 541 of file ParticleListAction.cxx.

    {

      // Add the first point in the trajectory.
      fCurrentParticle.particle->AddTrajectoryPoint(pos, mom, process);

      // also see if we can decide to keep the particle
      if (!fCurrentParticle.keep)
        fCurrentParticle.keep = fFilter->mustKeep(pos);

    } // ParticleListAction::AddPointToCurrentParticle()

void gar::garg4::ParticleListAction::BeginOfEventAction ( const G4Event *  )

virtual

Definition at line 78 of file ParticleListAction.cxx.

    {
      // Clear any previous particle information.
      fCurrentParticle.clear();
      fParticleList->clear();
      fParentIDMap.clear();
      fTrackIDToMCTruthIndex.clear();
      fCurrentTrackID = sdp::NoParticleId;

    }

void gar::garg4::ParticleListAction::EndOfEventAction ( const G4Event *  )

virtual

Definition at line 486 of file ParticleListAction.cxx.

    {
      // Set up the utility class for the "for_each" algorithm.  (We only
      // need a separate set-up for the utility class because we need to
      // give it the pointer to the particle list.  We're using the STL
      // "for_each" instead of the C++ "for loop" because it's supposed
      // to be faster.
      UpdateDaughterInformation updateDaughterInformation;
      updateDaughterInformation.SetParticleList( fParticleList );

      // Update the daughter information for each particle in the list.
      std::for_each(fParticleList->begin(),
                    fParticleList->end(),
                    updateDaughterInformation);
    }

static int gar::garg4::ParticleListAction::GetCurrentTrackID ( )

inlinestatic

Definition at line 79 of file ParticleListAction.h.

{ return fCurrentTrackID; }

sim::ParticleList * gar::garg4::ParticleListAction::GetList

(

)

const

Definition at line 504 of file ParticleListAction.cxx.

    {
      // check if the ParticleNavigator has entries, and if
      // so grab the highest track id value from it to
      // add to the fTrackIDOffset
      int highestID = 0;
      for(auto pn = fParticleList->begin(); pn != fParticleList->end(); ++pn)
        if( (*pn).first > highestID ) highestID = (*pn).first;

      fTrackIDOffset = highestID + 1;

      return fParticleList;
    }

int gar::garg4::ParticleListAction::GetParentage ( int  trackid ) const

private

Definition at line 93 of file ParticleListAction.cxx.

    {
      int parentid = sdp::NoParticleId;

      // search the fParentIDMap recursively until we have the parent id
      // of the first EM particle that led to this one
      std::map<int,int>::const_iterator itr = fParentIDMap.find(trackid);
      while( itr != fParentIDMap.end() ){
        MF_LOG_DEBUG("ParticleListAction")
          << "parentage for " << trackid
          << " " << (*itr).second;

        // set the parentid to the current parent ID, when the loop ends
        // this id will be the first EM particle
        parentid = (*itr).second;
        itr = fParentIDMap.find(parentid);
      }
      MF_LOG_DEBUG("ParticleListAction") << "final parent ID " << parentid;

      return parentid;
    }

static int gar::garg4::ParticleListAction::GetTrackIDOffset ( )

inlinestatic

Definition at line 80 of file ParticleListAction.h.

{ return fTrackIDOffset;  }

bool gar::garg4::ParticleListAction::IsDropped ( simb::MCParticle const *  p )

static

returns whether the specified particle has been marked as dropped

Definition at line 48 of file ParticleListAction.cxx.

    {
      return !p || p->Trajectory().empty();
    } // ParticleListAction::isDropped()

void gar::garg4::ParticleListAction::ParticleFilter ( std::unique_ptr< PositionInVolumeFilter > &&  filter )

inline

Grabs a particle filter.

Definition at line 74 of file ParticleListAction.h.

      { fFilter = std::move(filter); }

void gar::garg4::ParticleListAction::PostTrackingAction ( const G4Track *  aTrack )

virtual

Definition at line 330 of file ParticleListAction.cxx.

    {
      if (!fCurrentParticle.hasParticle()) return;

      // if we have found no reason to keep it, drop it!
      // (we might still need parentage information though)
      if (!fCurrentParticle.keep) {

        MF_LOG_VERBATIM("ParticleListAction")
          << "dropping particle with track id "
          << fCurrentParticle.particle->TrackId();

        fParticleList->Archive(fCurrentParticle.particle);
        // after the particle is archived, it is deleted
        fCurrentParticle.clear();
        return;
      }

      if(aTrack){
        fCurrentParticle.particle->SetWeight(aTrack->GetWeight());
        if (aTrack->GetStep()->GetPostStepPoint()->GetProcessDefinedStep())
          {
            G4String process = aTrack->GetStep()->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
            fCurrentParticle.particle->SetEndProcess(process);
          }
      }

      return;
    }

void gar::garg4::ParticleListAction::PreTrackingAction ( const G4Track *  track )

virtual

Definition at line 117 of file ParticleListAction.cxx.

    {
      // Particle type.
      G4ParticleDefinition* particleDefinition = track->GetDefinition();
      G4int pdgCode = particleDefinition->GetPDGEncoding();

      art::ServiceHandle<geo::GeometryGAr const> geo;
      TGeoManager *geomanager = geo->ROOTGeoManager();

      // Get Geant4's ID number for this track.  This will be the same
      // ID number that we'll use in the ParticleList.
      G4int trackID       = track->GetTrackID() + fTrackIDOffset;
      fCurrentTrackID     = trackID;
      size_t mcTruthIndex = 0;

      // And the particle's parent:
      G4int parentID = track->GetParentID() + fTrackIDOffset;

      std::string process_name = "unknown";

      // Is there an MCTruth object associated with this G4Track?  We
      // have to go up a "chain" of information to find out:
      const G4DynamicParticle* dynamicParticle = track->GetDynamicParticle();
      const G4PrimaryParticle* primaryParticle = dynamicParticle->GetPrimaryParticle();
      if ( primaryParticle ){
        const G4VUserPrimaryParticleInformation* gppi = primaryParticle->GetUserInformation();
        const g4b::PrimaryParticleInformation* ppi = dynamic_cast<const g4b::PrimaryParticleInformation*>(gppi);
        if ( ppi != 0 ){
          // If we've made it this far, a PrimaryParticleInformation
          // object exists and we are using a primary particle, set the
          // process name accordingly
          process_name = "primary";

          // primary particles should have parentID = 0, even if there
          // are multiple MCTruths for this event
          parentID = 0;

          mcTruthIndex = ppi->MCTruthIndex();
        } // end else no primary particle information
      } // Is there a G4PrimaryParticle?
      // If this is not a primary particle...
      else{
        // check if this particle was made in an EM shower, don't put it in the particle
        // list as we don't care about secondaries, tertiaries, etc for these showers
        // figure out what process is making this track - skip it if it is
        // one of pair production, compton scattering, photoelectric effect
        // bremstrahlung, annihilation, any ionization - reduces memory and the
        // volume of the output files.

        // trj Dec. 4, 2019 -- add in a feature to optionally skip EM shower daughter MCParticle recording
        // in dense media.  This is so we can keep the shower daughters in TPC gas but not the calorimeters or rock

        process_name = track->GetCreatorProcess()->GetProcessName();
        bool is_em_shower_daughter = ( process_name.find("conv")            != std::string::npos  ||
                                       process_name.find("LowEnConversion") != std::string::npos  ||
                                       process_name.find("Pair")            != std::string::npos  ||
                                       process_name.find("compt")           != std::string::npos  ||
                                       process_name.find("Compt")           != std::string::npos  ||
                                       process_name.find("Brem")            != std::string::npos  ||
                                       process_name.find("phot")            != std::string::npos  ||
                                       process_name.find("Photo")           != std::string::npos  ||
                                       process_name.find("hIoni")           != std::string::npos  ||
                                       process_name.find("eIoni")           != std::string::npos  ||
                                       process_name.find("ionIoni")         != std::string::npos  ||
                                       (process_name.find("Ion")             != std::string::npos  &&
                                            process_name.find("mu")              != std::string::npos) ||
                                       process_name.find("annihil")         != std::string::npos );
        const G4ThreeVector& trackpos = track->GetPosition();
        std::string matname = "";
        auto geonode = geomanager->FindNode(trackpos[0]/CLHEP::cm,
                                            trackpos[1]/CLHEP::cm,
                                            trackpos[2]/CLHEP::cm);
        if (geonode)
          {
            matname = geonode->GetMedium()->GetMaterial()->GetName();
          }
        
        std::regex const re_material(fEMShowerDaughterMatRegex);

        if(  is_em_shower_daughter &&
             (! fKeepEMShowerDaughters || 
              ! std::regex_match(matname, re_material)))
          {       

            // figure out the ultimate parentage of this particle
            // first add this track id and its parent to the fParentIDMap
            fParentIDMap[trackID] = parentID;

            fCurrentTrackID = -1 * this->GetParentage(trackID);

            // check that fCurrentTrackID is in the particle list - it is possible
            // that this particle's parent is a particle that did not get tracked.
            // An example is a partent that was made due to muMinusCaptureAtRest
            // and the daughter was made by the phot process.  The parent likely
            // isn't saved in the particle list because it is below the energy cut
            // which will put a bogus track id value into the list if we don't check it.
            if(!fParticleList->KnownParticle(fCurrentTrackID))
              fCurrentTrackID = sdp::NoParticleId;

            // clear current particle as we are not stepping this particle and
            // adding trajectory points to it
            fCurrentParticle.clear();
            return;

          } // end if not keeping EM shower daughters or skipping some processes


        // Check the energy of the particle.  If it falls below the energy
        // cut, don't add it to our list.
        G4double energy = track->GetKineticEnergy();
        if( energy * CLHEP::MeV / CLHEP::GeV < fEnergyCut ){
          fCurrentParticle.clear();

          // do add the particle to the parent id map though
          // and set the current track id to be it's ultimate parent
          fParentIDMap[trackID] = parentID;

          fCurrentTrackID = -1 * this->GetParentage(trackID);

          return;
        }

        //std::cout << "In GArG4 particle list action: ->" << matname << "<- " << fEMShowerDaughterMatRegex << std::endl;
        //std::cout << "track position: " << trackpos[0]/CLHEP::cm << " " << trackpos[1]/CLHEP::cm << " " << trackpos[2]/CLHEP::cm << std::endl;
        //std::cout << "Keeping mcparticle " << process_name << " " << energy << std::endl;


        // check to see if the parent particle has been stored in the particle navigator
        // if not, then see if it is possible to walk up the fParentIDMap to find the
        // ultimate parent of this particle.  Use that ID as the parent ID for this
        // particle
        if( !fParticleList->KnownParticle(parentID) ){
          // do add the particle to the parent id map
          // just in case it makes a daughter that we have to track as well
          fParentIDMap[trackID] = parentID;
          int pid = this->GetParentage(parentID);

          // if we still can't find the parent in the particle navigator,
          // we have to give up
          if( !fParticleList->KnownParticle(pid) ){
            MF_LOG_WARNING("ParticleListAction")
              << "can't find parent id: "
              << parentID
              << " in the particle list, or fParentIDMap."
              << " Make "
              << parentID
              << " the mother ID for"
              << " track ID "
              << fCurrentTrackID
              << " in the hope that it will aid debugging.";
          }
          else
            parentID = pid;
        }

        // Attempt to find the MCTruth index corresponding to the
        // current particle.  If the fCurrentTrackID is not in the
        // map try the parent ID, if that is not there, throw an
        // exception
        if(fTrackIDToMCTruthIndex.count(fCurrentTrackID) > 0 )
          mcTruthIndex = fTrackIDToMCTruthIndex.at(fCurrentTrackID);
        else if(fTrackIDToMCTruthIndex.count(parentID) > 0 )
          mcTruthIndex = fTrackIDToMCTruthIndex.at(parentID);
        else
          throw cet::exception("ParticleListAction")
            << "Cannot find MCTruth index for track id "
            << fCurrentTrackID
            << " or "
            << parentID;

      }// end if not a primary particle

      // This is probably the PDG mass, but just in case:
      double mass = dynamicParticle->GetMass() * CLHEP::MeV / CLHEP::GeV;

      // Create the sim::Particle object.
      fCurrentParticle.clear();
      fCurrentParticle.particle = new simb::MCParticle(trackID,
                                                       pdgCode,
                                                       process_name,
                                                       parentID,
                                                       mass);

      // if we are not filtering, we have a decision already
      if (!fFilter) fCurrentParticle.keep = true;

      // Polarization.
      const G4ThreeVector& polarization = track->GetPolarization();
      fCurrentParticle.particle->SetPolarization( TVector3(polarization.x(),
                                                           polarization.y(),
                                                           polarization.z()) );

      // Save the particle in the ParticleList.
      fParticleList->Add( fCurrentParticle.particle );

      MF_LOG_DEBUG("ParticleListAction")
        << "There are now "
        << fParticleList->size()
        << " particles in the list";

      if(fTrackIDToMCTruthIndex.count(fCurrentTrackID) > 0)
        MF_LOG_WARNING("ParticleListAction")
          << "attempting to put "
          << fCurrentTrackID
          << " into fTrackIDToMCTruthIndex map "
          << " particle is\n"
          << *(fCurrentParticle.particle);

      fTrackIDToMCTruthIndex[fCurrentTrackID] = mcTruthIndex;

    }

void gar::garg4::ParticleListAction::ResetTrackIDOffset ( )

inline

Definition at line 82 of file ParticleListAction.h.

{ fTrackIDOffset = 0;     }

void gar::garg4::ParticleListAction::SteppingAction ( const G4Step *  step )

virtual

Definition at line 362 of file ParticleListAction.cxx.

    {

      if ( !fCurrentParticle.hasParticle() ) return;

      // For the most part, we just want to add the post-step
      // information to the particle's trajectory.  There's one
      // exception: In PreTrackingAction, the correct time information
      // is not available.  So add the correct vertex information here.

      if ( fCurrentParticle.particle->NumberTrajectoryPoints() == 0 ){

        // Get the pre/along-step information from the G4Step.
        const G4StepPoint* preStepPoint = step->GetPreStepPoint();

        const G4ThreeVector position = preStepPoint->GetPosition();
        G4double time = preStepPoint->GetGlobalTime();

        // Remember that GArSoft uses cm, ns, GeV.
        TLorentzVector fourPos(position.x() / CLHEP::cm,
                               position.y() / CLHEP::cm,
                               position.z() / CLHEP::cm,
                               time);

        const G4ThreeVector momentum = preStepPoint->GetMomentum();
        const G4double energy = preStepPoint->GetTotalEnergy();
        TLorentzVector fourMom(momentum.x() * CLHEP::MeV / CLHEP::GeV,
                               momentum.y() * CLHEP::MeV / CLHEP::GeV,
                               momentum.z() * CLHEP::MeV / CLHEP::GeV,
                               energy * CLHEP::MeV / CLHEP::GeV);

        // Add the first point in the trajectory.
        AddPointToCurrentParticle( fourPos, fourMom, "Start" );

      } // end if this is the first step

      if (step->GetPostStepPoint()->GetProcessDefinedStep())
        {
          G4String process = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();

          // We store the initial creation point of the particle
          // and its final position (ie where it has no more energy, or at least < 1 eV) no matter
          // what, but whether we store the rest of the trajectory depends
          // on the process, and on a user switch.
          if ( fstoreTrajectories ){
            // Get the post-step information from the G4Step.
            const G4StepPoint* postStepPoint = step->GetPostStepPoint();

            const G4ThreeVector position = postStepPoint->GetPosition();
            G4double time = postStepPoint->GetGlobalTime();

            // Remember that GArSoft uses cm, ns, GeV.
            TLorentzVector fourPos(position.x() / CLHEP::cm,
                                   position.y() / CLHEP::cm,
                                   position.z() / CLHEP::cm,
                                   time);

            const G4ThreeVector momentum = postStepPoint->GetMomentum();
            const G4double energy = postStepPoint->GetTotalEnergy();
            TLorentzVector fourMom( momentum.x() * CLHEP::MeV / CLHEP::GeV,
                                    momentum.y() * CLHEP::MeV / CLHEP::GeV,
                                    momentum.z() * CLHEP::MeV / CLHEP::GeV,
                                    energy * CLHEP::MeV / CLHEP::GeV );

            // Add another point in the trajectory.
            AddPointToCurrentParticle( fourPos, fourMom, std::string(process) );

          }
        }
    }

std::map< int, size_t > gar::garg4::ParticleListAction::TrackIDToMCTruthIndexMap

(

)

const

Definition at line 535 of file ParticleListAction.cxx.

    {
      return fTrackIDToMCTruthIndex;
    }

sim::ParticleList && gar::garg4::ParticleListAction::YieldList

(

)

Definition at line 520 of file ParticleListAction.cxx.

    {
      // check if the ParticleNavigator has entries, and if
      // so grab the highest track id value from it to
      // add to the fTrackIDOffset
      int highestID = 0;
      for(auto pn = fParticleList->begin(); pn != fParticleList->end(); ++pn)
        if( (*pn).first > highestID ) highestID = (*pn).first;

      fTrackIDOffset = highestID + 1;

      return std::move(*fParticleList);
    } // ParticleList&& ParticleListAction::YieldList()

Member Data Documentation

ParticleInfo_t gar::garg4::ParticleListAction::fCurrentParticle

private

information about the particle currently being simulated for a single particle.

Definition at line 103 of file ParticleListAction.h.

int gar::garg4::ParticleListAction::fCurrentTrackID = sdp::NoParticleId

staticprivate

track ID of the current particle, set to eve ID for EM shower particles

Definition at line 110 of file ParticleListAction.h.

std::string gar::garg4::ParticleListAction::fEMShowerDaughterMatRegex

private

if keeping EM shower daughters, save only in media matching this

Definition at line 115 of file ParticleListAction.h.

G4double gar::garg4::ParticleListAction::fEnergyCut

private

The minimum energy for a particle to be included in the list.

Definition at line 101 of file ParticleListAction.h.

std::unique_ptr<PositionInVolumeFilter> gar::garg4::ParticleListAction::fFilter

private

filter for particles to be kept

Definition at line 117 of file ParticleListAction.h.

bool gar::garg4::ParticleListAction::fKeepEMShowerDaughters

private

whether to keep EM shower secondaries, tertiaries, etc

Definition at line 114 of file ParticleListAction.h.

std::map<int, int> gar::garg4::ParticleListAction::fParentIDMap

private

key is current track ID, value is parent ID

Definition at line 108 of file ParticleListAction.h.

sim::ParticleList* gar::garg4::ParticleListAction::fParticleList

private

The accumulated particle information for all particles in the event.

Definition at line 105 of file ParticleListAction.h.

G4bool gar::garg4::ParticleListAction::fstoreTrajectories

private

Whether to store particle trajectories with each particle.

Definition at line 107 of file ParticleListAction.h.

int gar::garg4::ParticleListAction::fTrackIDOffset = 0

staticprivate

offset added to track ids when running over multiple MCTruth objects.

Definition at line 112 of file ParticleListAction.h.

std::map<int, size_t> gar::garg4::ParticleListAction::fTrackIDToMCTruthIndex

private

map track ID to index of MCTruth in input list

Definition at line 109 of file ParticleListAction.h.

---

The documentation for this class was generated from the following files:

• garsoft/GArG4/ParticleListAction.h
• garsoft/GArG4/ParticleListAction.cxx