larg4::ParticleListAction Class Reference

#include <ParticleListAction.h>

Inheritance diagram for larg4::ParticleListAction:

Classes
struct	ParticleInfo_t

Public Types
using	GeneratedParticleIndex_t = simb::GeneratedParticleIndex_t

Public Member Functions
	ParticleListAction (double energyCut, bool storeTrajectories=false, bool keepEMShowerDaughters=false, bool keepMCParticleList=true)
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< util::PositionInVolumeFilter > &&filter)
	Grabs a particle filter. More...
void	ResetTrackIDOffset ()
const sim::ParticleList *	GetList () const
std::map< int, GeneratedParticleIndex_t > const &	GetPrimaryTruthMap () const
bool	hasList () const
	Returns whether a particle list is being kept. More...
GeneratedParticleIndex_t	GetPrimaryTruthIndex (int trackId) const
	Returns the index of primary truth (sim::NoGeneratorIndex if none). More...
sim::ParticleList &&	YieldList ()

Static Public Member Functions
static int	GetCurrentTrackID ()
static int	GetCurrentPdgCode ()
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
std::unique_ptr< 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...
bool	fKeepEMShowerDaughters
	whether to keep EM shower secondaries, tertiaries, etc More...
std::unique_ptr< util::PositionInVolumeFilter >	fFilter
	filter for particles to be kept More...
std::map< int, GeneratedParticleIndex_t >	fPrimaryTruthMap
	Map: particle track ID -> index of primary information in MC truth. More...

Static Private Attributes
static int	fCurrentTrackID = sim::NoParticleId
static int	fCurrentPdgCode = 0
	pdg code of current particle More...
static int	fTrackIDOffset = 0

Detailed Description

Definition at line 45 of file ParticleListAction.h.

Member Typedef Documentation

using larg4::ParticleListAction::GeneratedParticleIndex_t = simb::GeneratedParticleIndex_t

Definition at line 47 of file ParticleListAction.h.

Constructor & Destructor Documentation

larg4::ParticleListAction::ParticleListAction	(	double	energyCut,
		bool	storeTrajectories = false,
		bool	keepEMShowerDaughters = false,
		bool	keepMCParticleList = true
	)

Definition at line 57 of file ParticleListAction.cxx.

    : fenergyCut(energyCut * CLHEP::GeV)
    , fparticleList(keepMCParticleList ? std::make_unique<sim::ParticleList>() : nullptr)
    , fstoreTrajectories(storeTrajectories)
    , fKeepEMShowerDaughters(keepEMShowerDaughters)
  {}

Member Function Documentation

void larg4::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 509 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 larg4::ParticleListAction::BeginOfEventAction ( const G4Event *  )

virtual

Definition at line 71 of file ParticleListAction.cxx.

  {
    // Clear any previous particle information.
    fCurrentParticle.clear();
    if (fparticleList) fparticleList->clear();
    fParentIDMap.clear();
    fCurrentTrackID = sim::NoParticleId;
    fCurrentPdgCode = 0;
  }

void larg4::ParticleListAction::EndOfEventAction ( const G4Event *  )

virtual

Definition at line 440 of file ParticleListAction.cxx.

  {
    if (!fparticleList) return;

    // 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.get());

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

static int larg4::ParticleListAction::GetCurrentPdgCode ( )

inlinestatic

Definition at line 123 of file ParticleListAction.h.

    {
      return fCurrentPdgCode;
    }

static int larg4::ParticleListAction::GetCurrentTrackID ( )

inlinestatic

Definition at line 118 of file ParticleListAction.h.

    {
      return fCurrentTrackID;
    }

const sim::ParticleList * larg4::ParticleListAction::GetList

(

)

const

Definition at line 459 of file ParticleListAction.cxx.

  {
    if (!fparticleList) return nullptr;

    // 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;

    //Only change the fTrackIDOffset if there is in fact a particle to add to the event
    if ((fparticleList->size()) != 0) { fTrackIDOffset = highestID + 1; }

    return fparticleList.get();
  }

int larg4::ParticleListAction::GetParentage ( int  trackid ) const

private

Definition at line 87 of file ParticleListAction.cxx.

  {
    int parentid = sim::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;
  }

simb::GeneratedParticleIndex_t larg4::ParticleListAction::GetPrimaryTruthIndex

(

int

trackId

)

const

Returns the index of primary truth (sim::NoGeneratorIndex if none).

Definition at line 478 of file ParticleListAction.cxx.

  {
    auto const iInfo = GetPrimaryTruthMap().find(trackId);
    return (iInfo == GetPrimaryTruthMap().end()) ? simb::NoGeneratedParticleIndex : iInfo->second;
  } // ParticleListAction::GetPrimaryTruthIndex()

std::map<int, GeneratedParticleIndex_t> const& larg4::ParticleListAction::GetPrimaryTruthMap ( ) const

inline

Returns a map of truth record information index for each of the primary particles (by track ID).

Definition at line 140 of file ParticleListAction.h.

    {
      return fPrimaryTruthMap;
    }

bool larg4::ParticleListAction::hasList ( ) const

inline

Returns whether a particle list is being kept.

Definition at line 147 of file ParticleListAction.h.

    {
      return static_cast<bool>(fparticleList);
    }

bool larg4::ParticleListAction::isDropped ( simb::MCParticle const *  p )

static

returns whether the specified particle has been marked as dropped

Definition at line 50 of file ParticleListAction.cxx.

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

void larg4::ParticleListAction::ParticleFilter ( std::unique_ptr< util::PositionInVolumeFilter > &&  filter )

inline

Grabs a particle filter.

Definition at line 111 of file ParticleListAction.h.

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

void larg4::ParticleListAction::PostTrackingAction ( const G4Track *  aTrack )

virtual

Definition at line 260 of file ParticleListAction.cxx.

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

    // if we have found no reason to keep it, drop it!
    // (we might still need parentage information though)
    if (!fCurrentParticle.keep) {
      fparticleList->Archive(fCurrentParticle.particle.get());
      // after the particle is archived, it is deleted
      fCurrentParticle.clear();
      return;
    }

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

    // store truth record pointer, only if it is available
    if (fCurrentParticle.isPrimary()) {
      fPrimaryTruthMap[fCurrentParticle.particle->TrackId()] = fCurrentParticle.truthInfoIndex();
    }

    return;
  }

void larg4::ParticleListAction::PreTrackingAction ( const G4Track *  track )

virtual

Definition at line 110 of file ParticleListAction.cxx.

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

    // Get Geant4's ID number for this track.  This will be the same
    // ID number that we'll use in the ParticleList.
    // It is offset by the number of tracks accumulated from the previous Geant4
    // runs (if any)
    G4int trackID = track->GetTrackID() + fTrackIDOffset;
    fCurrentTrackID = trackID;
    fCurrentPdgCode = pdgCode;

    if (!fparticleList) {
      // the rest is about adding a new particle to the list: skip
      return; // note that fCurrentParticle is clear()'ed
    }

    // And the particle's parent (same offset as above):
    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();
    simb::GeneratedParticleIndex_t primaryIndex = simb::NoGeneratedParticleIndex;
    if (primaryParticle != 0) {
      const G4VUserPrimaryParticleInformation* gppi = primaryParticle->GetUserInformation();
      const g4b::PrimaryParticleInformation* ppi =
        dynamic_cast<const g4b::PrimaryParticleInformation*>(gppi);
      if (ppi != 0) {
        primaryIndex = ppi->MCParticleIndex();

        // 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;
      } // 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 - who wants to save
      // a buttload of electrons that arent from a CC interaction?
      process_name = track->GetCreatorProcess()->GetProcessName();
      if (!fKeepEMShowerDaughters && (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("Ion") != std::string::npos ||
                                      process_name.find("annihil") != std::string::npos)) {

        // 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 sim::IDE object for
        // the sim::SimChannel if we don't check it.
        if (!fparticleList->KnownParticle(fCurrentTrackID)) fCurrentTrackID = sim::NoParticleId;

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

      } // end if keeping EM shower daughters

      // 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 < 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;
      }

      // 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;
      }

    } // end if not a primary particle

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

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

    // 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.get());
  }

void larg4::ParticleListAction::ResetTrackIDOffset ( )

inline

Definition at line 129 of file ParticleListAction.h.

    {
      fTrackIDOffset = 0;
    }

void larg4::ParticleListAction::SteppingAction ( const G4Step *  step )

virtual

Definition at line 292 of file ParticleListAction.cxx.

  {

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

    // Temporary fix for problem where  DeltaTime on the first step
    // of optical photon propagation is calculated incorrectly. -wforeman
    globalTime = step->GetTrack()->GetGlobalTime();
    velocity_G4 = step->GetTrack()->GetVelocity();
    velocity_step = step->GetStepLength() / step->GetDeltaTime();
    if ((step->GetTrack()->GetDefinition()->GetPDGEncoding() == 0) &&
        fabs(velocity_G4 - velocity_step) > 0.0001) {
      // Subtract the faulty step time from the global time,
      // and add the correct step time based on G4 velocity.
      step->GetPostStepPoint()->SetGlobalTime(globalTime - step->GetDeltaTime() +
                                              step->GetStepLength() / velocity_G4);
    }

    // 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 LArSoft uses cm, ns, GeV.
      TLorentzVector fourPos(position.x() / CLHEP::cm,
                             position.y() / CLHEP::cm,
                             position.z() / CLHEP::cm,
                             time / CLHEP::ns);

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

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

    } // end if this is the first step

    // At this point, the particle is being transported through the
    // simulation. This method is being called for every voxel that
    // the track passes through, but we don't want to update the
    // trajectory information if we're just updating voxels. To check
    // for this, look at the process name for the step, and compare it
    // against the voxelization process name (set in PhysicsList.cxx).
    G4String process = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
    G4bool ignoreProcess = process.contains("LArVoxel") || process.contains("OpDetReadout");

    MF_LOG_DEBUG("ParticleListAction::SteppingAction")
      << ": DEBUG - process='" << process << "'"
      << " ignoreProcess=" << ignoreProcess << " fstoreTrajectories=" << fstoreTrajectories;

    // 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 && !ignoreProcess) {
      // Get the post-step information from the G4Step.
      const G4StepPoint* postStepPoint = step->GetPostStepPoint();

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

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

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

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

sim::ParticleList && larg4::ParticleListAction::YieldList

(

)

Definition at line 487 of file ParticleListAction.cxx.

  {
    if (!fparticleList) {
      // check with `hasList()` before calling this method
      throw cet::exception("ParticleListAction")
        << "ParticleListAction::YieldList(): particle list not build by user request.\n";
    }
    // 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;

    //Only change the fTrackIDOffset if there is in fact a particle to add to the event
    if ((fparticleList->size()) != 0) { fTrackIDOffset = highestID + 1; }

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

Member Data Documentation

ParticleInfo_t larg4::ParticleListAction::fCurrentParticle

private

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

Definition at line 168 of file ParticleListAction.h.

int larg4::ParticleListAction::fCurrentPdgCode = 0

staticprivate

pdg code of current particle

Definition at line 176 of file ParticleListAction.h.

int larg4::ParticleListAction::fCurrentTrackID = sim::NoParticleId

staticprivate

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

Definition at line 174 of file ParticleListAction.h.

G4double larg4::ParticleListAction::fenergyCut

private

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

Definition at line 166 of file ParticleListAction.h.

std::unique_ptr<util::PositionInVolumeFilter> larg4::ParticleListAction::fFilter

private

filter for particles to be kept

Definition at line 181 of file ParticleListAction.h.

bool larg4::ParticleListAction::fKeepEMShowerDaughters

private

whether to keep EM shower secondaries, tertiaries, etc

Definition at line 179 of file ParticleListAction.h.

std::map<int, int> larg4::ParticleListAction::fParentIDMap

private

key is current track ID, value is parent ID

Definition at line 173 of file ParticleListAction.h.

std::unique_ptr<sim::ParticleList> larg4::ParticleListAction::fparticleList

private

The accumulated particle information for all particles in the event.

Definition at line 170 of file ParticleListAction.h.

std::map<int, GeneratedParticleIndex_t> larg4::ParticleListAction::fPrimaryTruthMap

private

Map: particle track ID -> index of primary information in MC truth.

Definition at line 184 of file ParticleListAction.h.

G4bool larg4::ParticleListAction::fstoreTrajectories

private

Whether to store particle trajectories with each particle.

Definition at line 172 of file ParticleListAction.h.

int larg4::ParticleListAction::fTrackIDOffset = 0

staticprivate

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

Definition at line 177 of file ParticleListAction.h.

---

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

• larsim/larsim/LegacyLArG4/ParticleListAction.h
• larsim/larsim/LegacyLArG4/ParticleListAction.cxx