#include <ParticleListAction_service.h>

Inheritance diagram for larg4::ParticleListActionService:

Classes
struct	ParticleInfo_t

Public Member Functions
	ParticleListActionService (fhicl::ParameterSet const &)

void	preUserTrackingAction (const G4Track *) override

void	postUserTrackingAction (const G4Track *) override

void	userSteppingAction (const G4Step *) override

simb::GeneratedParticleIndex_t	GetPrimaryTruthIndex (int trackId) const
	Returns the index of primary truth (`sim::NoGeneratorIndex` if none). More...

void	beginOfEventAction (const G4Event *) override

void	endOfEventAction (const G4Event *) override

void	setInputCollections (std::vector< art::Handle< std::vector< simb::MCTruth >>> const &mclists)

void	setPtrInfo (art::ProductID pid, art::EDProductGetter const *productGetter)

std::unique_ptr< std::vector< simb::MCParticle > >	ParticleCollection ()

std::unique_ptr< art::Assns< simb::MCTruth, simb::MCParticle, sim::GeneratedParticleInfo > >	AssnsMCTruthToMCParticle ()

Private Member Functions
sim::ParticleList &&	YieldList ()

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::vector< std::string >	fkeepGenTrajectories

std::map< int, int >	fParentIDMap
	key is current track ID, value is parent ID More...

int	fCurrentTrackID

int	fTrackIDOffset

bool	fKeepEMShowerDaughters
	whether to keep EM shower secondaries, tertiaries, etc More...

std::vector< std::string >	fNotStoredPhysics
	Physics processes that will not be stored. More...

bool	fkeepOnlyPrimaryFullTraj

bool	fSparsifyTrajectories
	help reduce the number of trajectory points. More...

double	fSparsifyMargin
	set the sparsification margin More...

bool	fKeepTransportation
	tell whether or not to keep the transportation process More...

bool	fKeepSecondToLast
	tell whether or not to force keeping the second to last point More...

std::vector< art::Handle< std::vector< simb::MCTruth > > > const *	fMCLists
	MCTruthCollection input lists. More...

std::map< int, simb::GeneratedParticleIndex_t >	fPrimaryTruthMap
	Map: particle track ID -> index of primary information in MC truth. More...

std::map< int, size_t >	fMCTIndexMap
	Map: particle track ID -> index of primary parent in std::vector<simb::MCTruth> object. More...

std::map< int, bool >	fMCTPrimProcessKeepMap
	Map: particle trakc ID -> boolean decision to keep or not full trajectory points. More...

std::map< size_t, std::pair< std::string, G4bool > >	fMCTIndexToGeneratorMap
	Map: MCTruthIndex -> generator, input label of generator and keepGenerator decision. More...

std::unordered_map< std::string, int >	fNotStoredCounterUMap
	Map: not stored process and counter. More...

std::unique_ptr< std::vector< simb::MCParticle > >	partCol_

std::unique_ptr< art::Assns< simb::MCTruth, simb::MCParticle, sim::GeneratedParticleInfo > >	tpassn_

art::ProductID	pid_ {art::ProductID::invalid()}

art::EDProductGetter const *	productGetter_ {nullptr}

Detailed Description

Definition at line 44 of file ParticleListAction_service.h.

Constructor & Destructor Documentation

larg4::ParticleListActionService::ParticleListActionService ( fhicl::ParameterSet const & p )

explicit

Definition at line 68 of file ParticleListAction_service.cxx.

     : artg4tk::EventActionBase("PLASEventActionBase")
     , artg4tk::TrackingActionBase("PLASTrackingActionBase")
     , artg4tk::SteppingActionBase("PLASSteppingActionBase")
     , fenergyCut(p.get<double>("EnergyCut", 0.0 * CLHEP::GeV))
     , fstoreTrajectories(p.get<bool>("storeTrajectories", true))
     , fkeepGenTrajectories(p.get<std::vector<std::string>>("keepGenTrajectories", {}))
     , fKeepEMShowerDaughters(p.get<bool>("keepEMShowerDaughters", true))
     , fNotStoredPhysics(p.get<std::vector<std::string>>("NotStoredPhysics", {}))
     , fkeepOnlyPrimaryFullTraj(p.get<bool>("keepOnlyPrimaryFullTrajectories", false))
     , fSparsifyTrajectories(p.get<bool>("SparsifyTrajectories", false))
     , fSparsifyMargin(p.get<double>("SparsifyMargin"))
     , fKeepTransportation(p.get<bool>("KeepTransportation", false))
     , fKeepSecondToLast(p.get<bool>("KeepSecondToLast", false))
   {
     // -- D.R. If a custom list of not storable physics is provided, use it, otherwise
     //    use the default list. This preserves the behavior of the keepEmShowerDaughters
     //    parameter
     bool customNotStored = not fNotStoredPhysics.empty();
     if (!fKeepEMShowerDaughters) {
       // -- Don't keep all processes
       if (!customNotStored) // -- Don't keep but haven't provided a list
       {                     // -- default list of not stored physics
         fNotStoredPhysics = {"conv",
                              "LowEnConversion",
                              "Pair",
                              "compt",
                              "Compt",
                              "Brem",
                              "phot",
                              "Photo",
                              "Ion",
                              "annihil"};
       }
 
       std::stringstream sstored;
       sstored << "The full tracking information will not be stored for particles"
               << " resulting from the following processes: \n{ ";
       for (auto const& i : fNotStoredPhysics) {
         sstored << "\"" << i << "\" ";
         fNotStoredCounterUMap.emplace(i, 0); // -- initialize counter
       }
       mf::LogInfo("ParticleListActionService") << sstored.str() << "}\n";
     }
     else { // -- Keep all processes
       mf::LogInfo("ParticleListActionService")
         << "Storing full tracking information for all processes. \n";
       if (customNotStored) // -- custom list will be ignored
       {
         mf::LogWarning("StoredPhysics")
           << "NotStoredPhysics provided, but will be ignored."
           << " To use NotStoredPhysics, set keepEMShowerDaughters to false";
       }
     }
 
     // -- sparsify info
     if (fSparsifyTrajectories)
       mf::LogInfo("ParticleListActionService")
         << "Trajectory sparsification enabled with SparsifyMargin : " << fSparsifyMargin << "\n";
   }

Member Function Documentation

void larg4::ParticleListActionService::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 691 of file ParticleListAction_service.cxx.

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

std::unique_ptr<art::Assns<simb::MCTruth, simb::MCParticle, sim::GeneratedParticleInfo> > larg4::ParticleListActionService::AssnsMCTruthToMCParticle ( )

inline

Definition at line 91 of file ParticleListAction_service.h.

     {
       return std::move(tpassn_);
     }

void larg4::ParticleListActionService::beginOfEventAction ( const G4Event * )

override

Definition at line 132 of file ParticleListAction_service.cxx.

   {
     // Clear any previous particle information.
     fCurrentParticle.clear();
     fParticleList.clear();
     fParentIDMap.clear();
     fMCTIndexMap.clear();
     fMCTPrimProcessKeepMap.clear();
     fCurrentTrackID = sim::NoParticleId;
     fTrackIDOffset = 0;
 
     fPrimaryTruthMap.clear();
     fMCTIndexToGeneratorMap.clear();
     fNotStoredCounterUMap.clear();
 
     // -- D.R. If a custom list of keepGenTrajectories is provided, use it, otherwise
     //    keep or drop decision made based storeTrajectories parameter. This preserves
     //    the behavior of the storeTrajectories fhicl param
     bool customKeepTraj = not fkeepGenTrajectories.empty();
     if (!fstoreTrajectories) { // -- fstoreTrajectories : false
       mf::LogDebug("beginOfEventAction::Generator") << "Trajectory points will not be stored.";
     }
     else if (!customKeepTraj) { // -- fstoretrajectories : true and empty keepGenTrajectories list
       mf::LogDebug("beginOfEventAction::Generator")
         << "keepGenTrajectories list is empty. Will"
         << " store trajectory points for all generators";
     }
 
     // -- D.R. determine mapping between MCTruthIndex(s) and generator(s) for later reference
     size_t nKeep = 0;
     std::string generator_name = "unknown";
     for (size_t mcti = 0; mcti < fMCLists->size(); ++mcti) {
 
       std::stringstream sskeepgen;
       sskeepgen << "MCTruth object summary :";
       sskeepgen << "\n\tPrimary MCTIndex : " << mcti;
 
       // -- Obtain the generator (provenance) corresponding to the mctruth index:
       auto const& mclistHandle = (*fMCLists)[mcti];
       generator_name = mclistHandle.provenance()->inputTag().label();
       sskeepgen << "\n\tProvenance/Generator : " << generator_name;
 
       G4bool keepGen = false;
       if (fstoreTrajectories) // -- storeTrajectories set to true; check which
       {
         if (!customKeepTraj) { // -- no custom list, so keep all
           keepGen = true;
           nKeep++;
         }
         else { // -- custom list, so check the ones in the event against provided keep list
           for (auto keepableGen : fkeepGenTrajectories) {
             if (generator_name == keepableGen) { // -- exit upon finding match; false by default
               keepGen = true;
               nKeep++;
               break;
             }
           }
         }
       }
       fMCTIndexToGeneratorMap.emplace(mcti, std::make_pair(generator_name, keepGen));
       sskeepgen << "\n\tTrajectory points storable : " << (keepGen ? "true" : "false") << "\n";
       mf::LogDebug("beginOfEventAction::Generator") << sskeepgen.str();
     }
 
     if (nKeep == 0 && customKeepTraj && fstoreTrajectories) {
       mf::LogWarning("beginOfEventAction::keepableGenerators")
         << "storeTrajectories"
         << " set to true and a non-empty keepGenTrajectories list provided in configuration file, "
            "but"
         << " none of the generators in this list are present in the event! Double check list or "
            "don't"
         << " provide keepGenTrajectories in the configuration to keep all trajectories from all"
         << " generator labels. This may be expected for generators that have a nonzero probability "
            "of"
         << " producing no particles (e.g. some radiologicals)";
     }
   }

void larg4::ParticleListActionService::endOfEventAction ( const G4Event * )

override

Definition at line 702 of file ParticleListAction_service.cxx.

   {
     // -- End of Run Report
     if (!fNotStoredCounterUMap.empty()) { // -- Only if there is something to report
       std::stringstream sscounter;
       sscounter << "Not Stored Process summary:";
       for (auto const& [process, count] : fNotStoredCounterUMap) {
         sscounter << "\n\t" << process << " : " << count;
       }
       mf::LogInfo("ParticleListActionService") << sscounter.str();
     }
 
     partCol_ = std::make_unique<std::vector<simb::MCParticle>>();
     tpassn_ =
       std::make_unique<art::Assns<simb::MCTruth, simb::MCParticle, sim::GeneratedParticleInfo>>();
     // 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.
     std::for_each(
       fParticleList.begin(), fParticleList.end(), UpdateDaughterInformation{fParticleList});
 
     MF_LOG_INFO("endOfEventAction") << "MCTruth Handles Size: " << fMCLists->size();
 
     unsigned int nGeneratedParticles = 0;
     unsigned int nMCTruths = 0;
     sim::ParticleList particleList = YieldList();
 
     for (size_t mcl = 0; mcl < fMCLists->size(); ++mcl) {
       auto const& mclistHandle = (*fMCLists)[mcl];
       MF_LOG_INFO("endOfEventAction") << "mclistHandle Size: " << mclistHandle->size();
 
       for (size_t m = 0; m < mclistHandle->size(); ++m) {
         art::Ptr<simb::MCTruth> mct(mclistHandle, m);
 
         MF_LOG_INFO("endOfEventAction") << "Found " << mct->NParticles() << " particles";
 
         for (simb::MCParticle* p : particleList | ranges::views::values) {
           auto gen_index = fMCTIndexMap[p->TrackId()];
           if (gen_index != nMCTruths) continue;
           assert(p->NumberTrajectoryPoints() != 0ull);
 
           ++nGeneratedParticles;
           sim::GeneratedParticleInfo const truthInfo{GetPrimaryTruthIndex(p->TrackId())};
           if (!truthInfo.hasGeneratedParticleIndex() && p->Mother() == 0) {
             MF_LOG_WARNING("endOfEventAction") << "No GeneratedParticleIndex()!";
             // this means it's primary but with no information; logic error!!
             throw art::Exception(art::errors::LogicError)
               << "Failed to match primary particle:\n"
               << "\nwith particles from the truth record '" << mclistHandle.provenance()->inputTag()
               << "':\n\n";
           }
 
           partCol_->push_back(std::move(*p));
           art::Ptr<simb::MCParticle> mcp_ptr{pid_, partCol_->size() - 1, productGetter_};
           tpassn_->addSingle(mct, mcp_ptr, truthInfo);
         }
         mf::LogDebug("Offset") << "nGeneratedParticles = " << nGeneratedParticles;
         ++nMCTruths;
       }
     }
 
     fTrackIDOffset = 0;
   }

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

private

Definition at line 216 of file ParticleListAction_service.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
     auto itr = fParentIDMap.find(trackid);
     while (itr != fParentIDMap.end()) {
 
       // 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);
     }
 
     return parentid;
   }

simb::GeneratedParticleIndex_t larg4::ParticleListActionService::GetPrimaryTruthIndex ( int trackId ) const

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

Definition at line 661 of file ParticleListAction_service.cxx.

   {
     auto const it = fPrimaryTruthMap.find(trackId);
     return it == fPrimaryTruthMap.end() ? simb::NoGeneratedParticleIndex : it->second;
   }

std::unique_ptr<std::vector<simb::MCParticle> > larg4::ParticleListActionService::ParticleCollection ( )

inline

Definition at line 85 of file ParticleListAction_service.h.

     {
       return std::move(partCol_);
     }

void larg4::ParticleListActionService::postUserTrackingAction ( const G4Track * aTrack )

override

Definition at line 446 of file ParticleListAction_service.cxx.

   {
     if (!fCurrentParticle.hasParticle()) return;
 
     if (aTrack) {
       fCurrentParticle.particle->SetWeight(aTrack->GetWeight());
 
       // Get the post-step information from the G4Step.
       const G4StepPoint* postStepPoint = aTrack->GetStep()->GetPostStepPoint();
       if (!postStepPoint->GetProcessDefinedStep()) {
         // Now we get to do some awkward cleanup because the
         // fparticleList was augmented during the
         // preUserTrackingAction.  We cannot call 'Archive' because
         // that only sets the mapped type of the entry to
         // nullptr...which is really bad whenever we iterate through
         // entries in the endOfEventAction and dereference the mapped
         // type.  We have to entirely erase the entry.
         auto key_to_erase = fParticleList.key(fCurrentParticle.particle);
         fParticleList.erase(key_to_erase);
         // after the particle is archived, it is deleted
         fCurrentParticle.clear();
         return;
       }
 
       G4String process = postStepPoint->GetProcessDefinedStep()->GetProcessName();
       fCurrentParticle.particle->SetEndProcess(process);
 
       // -- D.R. Store the final point only for particles that have not had intermediate trajectory
       //    points saved. This avoids double counting the final trajectory point for particles from
       //    generators with storable trajectory points.
 
       if (!fCurrentParticle.keepFullTrajectory) {
         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));
       }
       // -- particle has a full trajectory, apply SparsifyTrajectory method if enabled
       else if (fSparsifyTrajectories) {
         fCurrentParticle.particle->SparsifyTrajectory(fSparsifyMargin, fKeepSecondToLast);
       }
     }
 
     // store truth record pointer, only if it is available
     if (fCurrentParticle.isPrimary()) {
       fPrimaryTruthMap[fCurrentParticle.particle->TrackId()] = fCurrentParticle.truthInfoIndex();
     }
 
     return;
   }

void larg4::ParticleListActionService::preUserTrackingAction ( const G4Track * track )

override

Definition at line 237 of file ParticleListAction_service.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)
     int const trackID = track->GetTrackID() + fTrackIDOffset;
     fCurrentTrackID = trackID;
 
     // And the particle's parent (same offset as above):
     int parentID = track->GetParentID() + fTrackIDOffset;
 
     std::string process_name = "unknown";
     std::string mct_primary_process = "unknown";
     bool isFromMCTProcessPrimary = false;
 
     // 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;
     size_t primarymctIndex = 0;
     if (primaryParticle != nullptr) {
       const G4VUserPrimaryParticleInformation* gppi = primaryParticle->GetUserInformation();
       const g4b::PrimaryParticleInformation* ppi =
         dynamic_cast<const g4b::PrimaryParticleInformation*>(gppi);
       if (ppi != nullptr) {
         primaryIndex = ppi->MCParticleIndex();
         primarymctIndex = ppi->MCTruthIndex();
 
         mct_primary_process = ppi->GetMCParticle()->Process();
 
         // If we've made it this far, a PrimaryParticleInformation
         // object exists and we are using a primary particle, set the
         // process name accordingly
 
         // -- D.R. : if process == "primary" exactly (this will most likely be the case), mark it as
         //    a primary with keepable trajectory for itself and its descendants.
         // -- elsif it simply starts with "primary", accept it but mark it as non-keep
         // -- else force it to be primary because we have determined that it is a primary particle
         //    This was the original behavior of the code i.e. process was _set_ to "primary"
         //    regardless of what the process name was in the generator MCTruth object.
         //
         // -- NOTE: This enforces a convention for process names assigned in the gen stage.
         if (mct_primary_process.compare("primary") == 0) {
           process_name = "primary";
           isFromMCTProcessPrimary = true;
         }
         else if (mct_primary_process.find("primary") == 0) {
           process_name = mct_primary_process;
           isFromMCTProcessPrimary = false;
           mf::LogDebug("PrimaryParticle") << "MCTruth primary process name contains \"primary\" "
                                           << " but is not solely \"primary\" : " << process_name
                                           << ".\nWill not store full set of trajectory points.";
         }
         else { // -- override it
           process_name = "primary";
           isFromMCTProcessPrimary = true;
           mf::LogWarning("PrimaryParticle")
             << "MCTruth primary process does not beging with string"
             << " literal \"primary\" : " << process_name << "\nOVERRIDING it to \"primary\"";
         }
         // -- The process_name check is simply to allow an additional way to reduce memory usage,
         //    namely, by creating MCTruth input particles with multiple process labels (e.g. "primary"
         //    and "primaryBackground") that can be used to veto storage of trajectory points.
 
         // 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 undesirable process. For example:
       // if one is not interested in EM shower particles, don't put it in the particle
       // list as one wouldn'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, or ionization
       process_name = track->GetCreatorProcess()->GetProcessName();
       if (!fKeepEMShowerDaughters) {
         bool notstore = false;
         for (auto const& p : fNotStoredPhysics) {
           if (process_name.find(p) != std::string::npos) {
             notstore = true;
             mf::LogDebug("NotStoredPhysics") << "Found process : " << process_name;
 
             int old = 0;
             auto search = fNotStoredCounterUMap.find(p);
             if (search != fNotStoredCounterUMap.end()) { old = search->second; }
             fNotStoredCounterUMap.insert_or_assign(p, (old + 1));
             break;
           }
         }
 
         if (notstore) {
 
           // 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 parent 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 process matches an undesired process
       }   // 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("ParticleListActionService")
             << "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;
       }
 
       // Once the parentID is secured, inherit the MCTruth Index
       // which should have been set already
       if (auto it = fMCTIndexMap.find(parentID); it != cend(fMCTIndexMap)) {
         primarymctIndex = it->second;
       }
       else {
         throw art::Exception(art::errors::LogicError)
           << "Could not locate MCT Index for parent trackID of " << parentID;
       }
 
       // Inherit whether the parent is from a primary with MCTruth process_name == "primary"
       isFromMCTProcessPrimary = fMCTPrimProcessKeepMap[parentID];
     } // 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;
 
     fMCTIndexMap[trackID] = primarymctIndex;
 
     fMCTPrimProcessKeepMap[trackID] = isFromMCTProcessPrimary;
 
     // -- determine whether full set of trajectorie points should be stored or only the start and end points
     fCurrentParticle.keepFullTrajectory =
       (!fstoreTrajectories) ?
         false : /*don't want trajectory points at all, bail*/
         (!(fMCTIndexToGeneratorMap[primarymctIndex].second)) ?
         false : /*particle is not from a storable generator*/
           (!fkeepOnlyPrimaryFullTraj) ?
         true : /*want all primaries tracked for a storable generator*/
             (isFromMCTProcessPrimary) ?
         true :       /*only descendants from primaries with MCTruth process == "primary"*/
               false; /*not from MCTruth process "primary"*/
 
     // 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);
   }

void larg4::ParticleListActionService::setInputCollections ( std::vector< art::Handle< std::vector< simb::MCTruth >>> const & mclists )

inline

Definition at line 71 of file ParticleListAction_service.h.

     {
       fMCLists = &mclists;
     }

void larg4::ParticleListActionService::setPtrInfo	(	art::ProductID	pid,
		art::EDProductGetter const *	productGetter
	)

inline

Definition at line 77 of file ParticleListAction_service.h.

     {
       pid_ = pid;
       productGetter_ = productGetter;
     }

void larg4::ParticleListActionService::userSteppingAction ( const G4Step * step )

override

Definition at line 514 of file ParticleListAction_service.cxx.

   {
     // N.B. G4 guarantees that following are non-null:
     //  - step
     //  - step->GetPostStepPoint()
     if (!fCurrentParticle.hasParticle() || !step->GetPostStepPoint()->GetProcessDefinedStep()) {
       return;
     }
     // Temporary fix for problem where  DeltaTime on the first step
     // of optical photon propagation is calculated incorrectly. -wforeman
     double const globalTime = step->GetTrack()->GetGlobalTime();
     double const velocity_G4 = step->GetTrack()->GetVelocity();
     double const 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::LogDebug("ParticleListActionService::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.
     // -- D.R. Store additional trajectory points only for desired generators and processes
     if (!ignoreProcess && fCurrentParticle.keepFullTrajectory) {
 
       // 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::ParticleListActionService::YieldList ( )

private

Definition at line 670 of file ParticleListAction_service.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;
 
     //Only change the fTrackIDOffset if there is in fact a particle to add to the event
     if ((fParticleList.size()) != 0) {
       fTrackIDOffset = highestID + 1;
       mf::LogDebug("YieldList:fTrackIDOffset")
         << "highestID = " << highestID << "\nfTrackIDOffset= " << fTrackIDOffset;
     }
 
     return std::move(fParticleList);
   } // ParticleList&& ParticleListActionService::YieldList()

Member Data Documentation

ParticleInfo_t larg4::ParticleListActionService::fCurrentParticle

private

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

Definition at line 145 of file ParticleListAction_service.h.

int larg4::ParticleListActionService::fCurrentTrackID

private

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

Definition at line 157 of file ParticleListAction_service.h.

G4double larg4::ParticleListActionService::fenergyCut

private

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

Definition at line 143 of file ParticleListAction_service.h.

bool larg4::ParticleListActionService::fKeepEMShowerDaughters

private

whether to keep EM shower secondaries, tertiaries, etc

Definition at line 161 of file ParticleListAction_service.h.

std::vector<std::string> larg4::ParticleListActionService::fkeepGenTrajectories

private

List of generators for which fstoreTrajectories applies. if not provided and storeTrajectories is true, then all trajectories for all generators will be stored. If storeTrajectories is set to false, this list is ignored and all additional trajectory points are not stored.

Definition at line 151 of file ParticleListAction_service.h.

bool larg4::ParticleListActionService::fkeepOnlyPrimaryFullTraj

private

Whether to store trajectories only for primaries and their descendants with MCTruth process = "primary"

Definition at line 163 of file ParticleListAction_service.h.

bool larg4::ParticleListActionService::fKeepSecondToLast

private

tell whether or not to force keeping the second to last point

Definition at line 168 of file ParticleListAction_service.h.

bool larg4::ParticleListActionService::fKeepTransportation

private

tell whether or not to keep the transportation process