larg4::ParticleListActionService Class Reference

#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< std::map< int, std::set< int > > >	DroppedTracksCollection ()
std::unique_ptr< art::Assns < simb::MCTruth, simb::MCParticle, sim::GeneratedParticleInfo > >	AssnsMCTruthToMCParticle ()
std::map< int, int >	GetTargetIDMap ()

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...
std::map< int, int >	fTargetIDMap
	key is original track ID, value is ID to assign for downstream objs (e.g. SimEdeps) 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::map< int, std::set< int > >	fdroppedTracksMap
	map <ParentID, set: list of track ids for which no MCParticle was created> More...
std::unique_ptr< std::vector < simb::MCParticle > >	partCol_
std::unique_ptr< std::map< int, std::set< int > > >	droppedCol_
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 61 of file ParticleListAction_service.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 66 of file ParticleListAction.cc.

    : 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", {}))

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 674 of file ParticleListAction.cc.

  {
    // 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 111 of file ParticleListAction_service.h.

    {
      return std::move(tpassn_);
    }

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

override

Definition at line 130 of file ParticleListAction.cc.

  {
    // Clear any previous particle information.
    fCurrentParticle.clear();
    fParticleList.clear();
    fParentIDMap.clear();
    fTargetIDMap.clear();
    fMCTIndexMap.clear();
    fMCTPrimProcessKeepMap.clear();
    fCurrentTrackID = sim::NoParticleId;
    fTrackIDOffset = 0;
    fPrimaryTruthMap.clear();
    fMCTIndexToGeneratorMap.clear();
    fNotStoredCounterUMap.clear();
    fdroppedTracksMap.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)";
    }
  }

std::unique_ptr<std::map<int,std::set<int> > > larg4::ParticleListActionService::DroppedTracksCollection ( )

inline

Definition at line 105 of file ParticleListAction_service.h.

    {
      return std::move(droppedCol_);
    }

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

override

Definition at line 685 of file ParticleListAction.cc.

  {
    // -- 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>>();
    droppedCol_ = std::make_unique<std::map<int,std::set<int>>>();
    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;
        for(auto const& dropped :  fdroppedTracksMap )
          {
            droppedCol_->insert(dropped);
          }
       ++nMCTruths;
      }
    }
    fTrackIDOffset = 0;
  }

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

private

Definition at line 215 of file ParticleListAction.cc.

  {
    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 644 of file ParticleListAction.cc.

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

std::map<int, int> larg4::ParticleListActionService::GetTargetIDMap ( )

inline

Definition at line 116 of file ParticleListAction_service.h.

{return fTargetIDMap;}

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

inline

Definition at line 101 of file ParticleListAction_service.h.

    {
      return std::move(partCol_);
    }

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

override

Definition at line 438 of file ParticleListAction.cc.

  {
    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 235 of file ParticleListAction.cc.

  {
    // 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;
    fTargetIDMap[trackID] = fCurrentTrackID;
    // 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;
          fTargetIDMap[trackID] = fCurrentTrackID;
          // clear current particle as we are not stepping this particle and
          // adding trajectory points to it
          fdroppedTracksMap[this->GetParentage(trackID)].insert(trackID);
          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) {
        fdroppedTracksMap[this->GetParentage(trackID)].insert(trackID);
        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);
        fTargetIDMap[trackID] = fCurrentTrackID;
        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 88 of file ParticleListAction_service.h.

    {
      fMCLists = &mclists;
    }

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

inline

Definition at line 94 of file ParticleListAction_service.h.

    {
      pid_ = pid;
      productGetter_ = productGetter;
    }

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

override

Definition at line 506 of file ParticleListAction.cc.

  {
    // 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.
    // change below:
    // This method is being called for every step that
    // the track passes through, but we don't want to update the
    // trajectory information if the step  was defined by the StepLimiter. 
    G4String process = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
    G4bool ignoreProcess = process.contains("StepLimiter");
    
    // 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 653 of file ParticleListAction.cc.

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

std::unique_ptr<std::map<int,std::set<int> > > larg4::ParticleListActionService::droppedCol_

private

Definition at line 214 of file ParticleListAction_service.h.

ParticleInfo_t larg4::ParticleListActionService::fCurrentParticle

private

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

Definition at line 166 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 179 of file ParticleListAction_service.h.

std::map<int,std::set<int> > larg4::ParticleListActionService::fdroppedTracksMap

private

map <ParentID, set: list of track ids for which no MCParticle was created>

Definition at line 211 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 164 of file ParticleListAction_service.h.

bool larg4::ParticleListActionService::fKeepEMShowerDaughters

private

whether to keep EM shower secondaries, tertiaries, etc

Definition at line 183 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 172 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 185 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 190 of file ParticleListAction_service.h.

bool larg4::ParticleListActionService::fKeepTransportation

private

tell whether or not to keep the transportation process

Definition at line 189 of file ParticleListAction_service.h.

std::vector<art::Handle<std::vector<simb::MCTruth> > > const* larg4::ParticleListActionService::fMCLists

private

MCTruthCollection input lists.

Definition at line 193 of file ParticleListAction_service.h.

std::map<int, size_t> larg4::ParticleListActionService::fMCTIndexMap

private

Map: particle track ID -> index of primary parent in std::vector<simb::MCTruth> object.

Definition at line 199 of file ParticleListAction_service.h.

std::map<size_t, std::pair<std::string, G4bool> > larg4::ParticleListActionService::fMCTIndexToGeneratorMap

private

Map: MCTruthIndex -> generator, input label of generator and keepGenerator decision.

Definition at line 205 of file ParticleListAction_service.h.

std::map<int, bool> larg4::ParticleListActionService::fMCTPrimProcessKeepMap

private

Map: particle trakc ID -> boolean decision to keep or not full trajectory points.

Definition at line 202 of file ParticleListAction_service.h.

std::unordered_map<std::string, int> larg4::ParticleListActionService::fNotStoredCounterUMap

private

Map: not stored process and counter.

Definition at line 208 of file ParticleListAction_service.h.

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

private

Physics processes that will not be stored.

Definition at line 184 of file ParticleListAction_service.h.

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

private

key is current track ID, value is parent ID

Definition at line 177 of file ParticleListAction_service.h.

sim::ParticleList larg4::ParticleListActionService::fParticleList

private

The accumulated particle information for all particles in the event.

Definition at line 168 of file ParticleListAction_service.h.

std::map<int, simb::GeneratedParticleIndex_t> larg4::ParticleListActionService::fPrimaryTruthMap

private

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

Definition at line 196 of file ParticleListAction_service.h.

double larg4::ParticleListActionService::fSparsifyMargin

private

set the sparsification margin

Definition at line 188 of file ParticleListAction_service.h.

bool larg4::ParticleListActionService::fSparsifyTrajectories

private

help reduce the number of trajectory points.

Definition at line 187 of file ParticleListAction_service.h.

G4bool larg4::ParticleListActionService::fstoreTrajectories

private

Whether to store particle trajectories with each particle.

Definition at line 170 of file ParticleListAction_service.h.

std::map<int, int> larg4::ParticleListActionService::fTargetIDMap

private

key is original track ID, value is ID to assign for downstream objs (e.g. SimEdeps)

Definition at line 178 of file ParticleListAction_service.h.

int larg4::ParticleListActionService::fTrackIDOffset

mutableprivate

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

Definition at line 181 of file ParticleListAction_service.h.

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

private

Definition at line 213 of file ParticleListAction_service.h.

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

private

Definition at line 217 of file ParticleListAction_service.h.

art::EDProductGetter const* larg4::ParticleListActionService::productGetter_ {nullptr}

private

Definition at line 218 of file ParticleListAction_service.h.

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

private

Definition at line 216 of file ParticleListAction_service.h.

---

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

• ParticleListAction_service.h
• ParticleListAction.cc