ParticleListAction_service.h

Go to the documentation of this file.

//
// accumulate a list of particles modeled by Geant4.
//

// Include guard
#ifndef PARTICLELISTACTION_SERVICE_H
#define PARTICLELISTACTION_SERVICE_H

#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceDeclarationMacros.h"

#include "lardataobj/Simulation/GeneratedParticleInfo.h"
#include "lardataobj/Simulation/ParticleAncestryMap.h"

#include "artg4tk/actionBase/EventActionBase.hh"
#include "artg4tk/actionBase/SteppingActionBase.hh"
#include "artg4tk/actionBase/TrackingActionBase.hh"

#include "nug4/ParticleNavigation/ParticleList.h"

#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/simb.h" // simb::GeneratedParticleIndex_t

#include "larcore/Geometry/Geometry.h"
#include "larcorealg/CoreUtils/ParticleFilters.h"

#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceDeclarationMacros.h"

#include "canvas/Persistency/Common/Assns.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Provenance/ProductID.h"

namespace art {
  class EDProductGetter;
}

#include "messagefacility/MessageLogger/MessageLogger.h"

#include "Geant4/globals.hh"

class G4Event;
class G4Track;
class G4Step;

class TLorentzVector;

#include <map>
#include <memory>
#include <set>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

namespace larg4 {

  class ParticleListActionService : public artg4tk::EventActionBase,
                                    public artg4tk::TrackingActionBase,
                                    public artg4tk::SteppingActionBase {
  public:
    // Standard constructors and destructors;
    explicit ParticleListActionService(fhicl::ParameterSet const&);

    // UserActions method that we'll override, to obtain access to
    // Geant4's particle tracks and trajectories.
    void preUserTrackingAction(const G4Track*) override;
    void postUserTrackingAction(const G4Track*) override;
    void userSteppingAction(const G4Step*) override;

    simb::GeneratedParticleIndex_t GetPrimaryTruthIndex(int trackId) const;

    // Called at the beginning of each event (note that this is after the
    // primaries have been generated and sent to the event manager)
    void beginOfEventAction(const G4Event*) override;

    // Called at the end of each event, right before GEANT's state switches
    // out of event processing and into closed geometry (last chance to access
    // the current event).
    void endOfEventAction(const G4Event*) override;

    // Set/get the current Art event
    void setInputCollections(std::vector<art::Handle<std::vector<simb::MCTruth>>> const& mclists)
    {
      fMCLists = &mclists;
    }

    void setPtrInfo(art::ProductID pid, art::EDProductGetter const* productGetter)
    {
      pid_ = pid;
      productGetter_ = productGetter;
    }

    std::unique_ptr<std::vector<simb::MCParticle>> ParticleCollection()
    {
      return std::move(partCol_);
    }
    std::unique_ptr<std::vector<simb::MCParticle>> DroppedParticleCollection()
    {
      return std::move(droppedPartCol_);
    }
    std::unique_ptr<sim::ParticleAncestryMap> DroppedTracksCollection()
    {
      return std::move(droppedCol_);
    }

    std::unique_ptr<art::Assns<simb::MCTruth, simb::MCParticle, sim::GeneratedParticleInfo>>
    AssnsMCTruthToMCParticle()
    {
      return std::move(tpassn_);
    }

    std::map<int, int> GetTargetIDMap() { return fTargetIDMap; }

    void CreateParticleFilter(std::vector<std::string> keepParticlesInVolumes,
                              std::unique_ptr<util::PositionInVolumeFilter>& filter)
    {
      // if we don't have favourite volumes, don't even bother creating a filter
      std::set<std::string> vol_names(keepParticlesInVolumes.begin(), keepParticlesInVolumes.end());

      if (empty(vol_names)) filter = {};

      auto const& geom = *art::ServiceHandle<geo::Geometry const>();

      std::vector<std::vector<TGeoNode const*>> node_paths = geom.FindAllVolumePaths(vol_names);

      // collection of interesting volumes
      util::PositionInVolumeFilter::AllVolumeInfo_t GeoVolumePairs;
      GeoVolumePairs.reserve(node_paths.size()); // because we are obsessed

      //for each interesting volume, follow the node path and collect
      //total rotations and translations
      for (size_t iVolume = 0; iVolume < node_paths.size(); ++iVolume) {
        std::vector<TGeoNode const*> path = node_paths[iVolume];

        auto pTransl = new TGeoTranslation(0., 0., 0.);
        auto pRot = new TGeoRotation();
        for (TGeoNode const* node : path) {
          TGeoTranslation thistranslate(*node->GetMatrix());
          TGeoRotation thisrotate(*node->GetMatrix());
          pTransl->Add(&thistranslate);
          *pRot = *pRot * thisrotate;
        }

        // for some reason, pRot and pTransl don't have tr and rot bits set
        // correctly make new translations and rotations so bits are set correctly
        auto pTransl2 = new TGeoTranslation(
          pTransl->GetTranslation()[0], pTransl->GetTranslation()[1], pTransl->GetTranslation()[2]);
        double phi = 0., theta = 0., psi = 0.;
        pRot->GetAngles(phi, theta, psi);
        auto pRot2 = new TGeoRotation();
        pRot2->SetAngles(phi, theta, psi);

        auto pTransf = new TGeoCombiTrans(*pTransl2, *pRot2);
        GeoVolumePairs.emplace_back(node_paths[iVolume].back()->GetVolume(), pTransf);
      }

      filter = std::make_unique<util::PositionInVolumeFilter>(std::move(GeoVolumePairs));
    }
    void ParticleFilter() { CreateParticleFilter(fKeepParticlesInVolumes, fFilter); }
    void DroppedParticleFilter()
    {
      CreateParticleFilter(fKeepDroppedParticlesInVolumes, fDroppedFilter);
    }
    bool storeDropped() const { return fStoreDroppedMCParticles; }

  private:
    struct ParticleInfo_t {
      simb::MCParticle* particle = nullptr; 
      bool keepFullTrajectory = false;      
      bool isInVolume = false;              
      bool isDropped = false;               

      simb::GeneratedParticleIndex_t truthIndex = simb::NoGeneratedParticleIndex;

      void clear()
      {
        particle = nullptr;
        keepFullTrajectory = false;
        isInVolume = false;
        isDropped = false;
        truthIndex = simb::NoGeneratedParticleIndex;
      }

      bool hasParticle() const { return particle; }

      bool isPrimary() const { return simb::isGeneratedParticleIndex(truthIndex); }

      simb::GeneratedParticleIndex_t truthInfoIndex() const { return truthIndex; }

      void print()
      {
        std::cout << "ParticleInfo_t: " << std::endl;
        std::cout << "  particle: " << particle << std::endl;
        std::cout << "  keepFullTrajectory: " << keepFullTrajectory << std::endl;
        std::cout << "  isInVolume: " << isInVolume << std::endl;
        std::cout << "  isDropped: " << isDropped << std::endl;
        std::cout << "  truthIndex: " << truthIndex << std::endl;
      }

    }; // ParticleInfo_t

    // Returns whether the particle was dropped
    bool isDropped(simb::MCParticle const* p);

    // Yields the ParticleList accumulated during the current event.
    sim::ParticleList&& YieldList();

    // Yields the (dropped) ParticleList accumulated during the current event.
    sim::ParticleList&& YieldDroppedList();

    // this method will loop over the fParentIDMap to get the
    // parentage of the provided trackid
    int GetParentage(int trackid) const;

    G4double fenergyCut;             
    ParticleInfo_t fCurrentParticle; 
    sim::ParticleList fParticleList; 
    G4bool fstoreTrajectories;       
    std::vector<std::string>
      fkeepGenTrajectories;          
    std::map<int, int> fParentIDMap; 
    std::map<int, int>
      fTargetIDMap; 
    int fCurrentTrackID;         
    mutable int fTrackIDOffset;  
    bool fKeepEMShowerDaughters; 
    std::vector<std::string> fNotStoredPhysics; 
    bool fkeepOnlyPrimaryFullTraj; 
    bool fSparsifyTrajectories;    
    double fSparsifyMargin;        
    bool fKeepTransportation;      
    bool fKeepSecondToLast; 

    std::vector<std::string>
      fKeepParticlesInVolumes; 

    std::vector<std::string> fKeepDroppedParticlesInVolumes; 
    bool
      fStoreDroppedMCParticles; 

    std::unique_ptr<sim::ParticleList>
      fdroppedParticleList; 
    std::vector<art::Handle<std::vector<simb::MCTruth>>> const*
      fMCLists; 

    std::map<int, simb::GeneratedParticleIndex_t> fPrimaryTruthMap;

    std::map<int, size_t> fMCTIndexMap;

    std::map<int, bool> fMCTPrimProcessKeepMap;

    std::map<size_t, std::pair<std::string, G4bool>> fMCTIndexToGeneratorMap;

    std::unordered_map<std::string, int> fNotStoredCounterUMap;

    std::map<int, std::set<int>> fdroppedTracksMap;

    std::unique_ptr<std::vector<simb::MCParticle>> partCol_;
    std::unique_ptr<std::vector<simb::MCParticle>> droppedPartCol_;
    std::unique_ptr<sim::ParticleAncestryMap> droppedCol_;
    std::unique_ptr<art::Assns<simb::MCTruth, simb::MCParticle, sim::GeneratedParticleInfo>>
      tpassn_;
    art::ProductID pid_{art::ProductID::invalid()};
    art::EDProductGetter const* productGetter_{nullptr};
    std::unique_ptr<util::PositionInVolumeFilter> fFilter; 
    std::unique_ptr<util::PositionInVolumeFilter>
      fDroppedFilter; 
    void AddPointToCurrentParticle(TLorentzVector const& pos,
                                   TLorentzVector const& mom,
                                   std::string const& process);
  };

} // namespace larg4

DECLARE_ART_SERVICE(larg4::ParticleListActionService, LEGACY)

#endif // PARTICLELISTACTION_SERVICE_H