protoana::ProtoDUNETruthUtils Class Reference

#include <ProtoDUNETruthUtils.h>

Public Member Functions
	ProtoDUNETruthUtils ()
	~ProtoDUNETruthUtils ()
std::vector< const recob::Hit * >	FillSharedHits (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &mcpart, const std::vector< const recob::Hit * > &hitsVec, bool delta_ray, bool use_eve=true) const
std::vector< const recob::Hit * >	GetSharedHits (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &mcpart, const recob::PFParticle &pfpart, const art::Event &evt, std::string pfparticleModule, bool delta_ray=false) const
std::vector< const recob::Hit * >	GetSharedHits (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &mcpart, const recob::Track &track, const art::Event &evt, std::string trackModule, bool delta_ray=false) const
std::vector< const recob::Hit * >	GetSharedHits (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &mcpart, const recob::Shower &shower, const art::Event &evt, std::string showerModule, bool delta_ray=false) const
std::vector< const recob::Hit * >	GetMCParticleHits (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &mcpart, const art::Event &evt, std::string hitModule, bool use_eve=true) const
template<typename T >
double	GetCompleteness (detinfo::DetectorClocksData const &clockData, const T &recobj, const art::Event &evt, std::string recoModule, std::string hitModule) const
double	GetPurity (detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, const art::Event &evt, std::string pfparticleModule) const
double	GetPurity (detinfo::DetectorClocksData const &clockData, const recob::Track &track, const art::Event &evt, std::string trackModule) const
double	GetPurity (detinfo::DetectorClocksData const &clockData, const recob::Shower &shower, const art::Event &evt, std::string showerModule) const
std::vector< std::pair< const simb::MCParticle *, double > >	GetMCParticleListFromTrackHits (detinfo::DetectorClocksData const &clockData, const std::vector< const recob::Hit * > &hitVec, bool use_eve=true) const
std::vector< std::pair< const simb::MCParticle *, double > >	GetMCParticleListFromShowerHits (detinfo::DetectorClocksData const &clockData, const std::vector< const recob::Hit * > &hitVec) const
std::vector< std::pair< const simb::MCParticle *, double > >	GetMCParticleListFromPFParticle (detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, art::Event const &evt, std::string pfparticleModule) const
std::vector< std::pair< const recob::PFParticle *, double > >	GetPFParticleListFromMCParticle (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string pfparticleModule, bool use_eve=true) const
std::vector< std::pair< const simb::MCParticle *, double > >	GetMCParticleListFromRecoTrack (detinfo::DetectorClocksData const &clockData, const recob::Track &track, art::Event const &evt, std::string trackModule) const
std::vector< std::pair< const recob::Track *, double > >	GetRecoTrackListFromMCParticle (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string trackModule) const
std::vector< std::pair< const simb::MCParticle *, double > >	GetMCParticleListFromRecoShower (detinfo::DetectorClocksData const &clockData, const recob::Shower &shower, art::Event const &evt, std::string showerModule) const
std::vector< std::pair< const recob::Shower *, double > >	GetRecoShowerListFromMCParticle (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string showerModule) const
const simb::MCParticle *	GetMCParticleFromPFParticle (detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, art::Event const &evt, std::string pfparticleModule) const
const recob::PFParticle *	GetPFParticleFromMCParticle (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string pfparticleModule) const
const simb::MCParticle *	GetMCParticleFromRecoTrack (detinfo::DetectorClocksData const &clockData, const recob::Track &track, art::Event const &evt, std::string trackModule) const
const recob::Track *	GetRecoTrackFromMCParticle (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string trackModule) const
const simb::MCParticle *	GetMCParticleFromRecoShower (detinfo::DetectorClocksData const &clockData, const recob::Shower &shower, art::Event const &evt, std::string showerModule) const
const recob::Shower *	GetRecoShowerFromMCParticle (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string showerModule) const
const simb::MCParticle *	GetMCParticleFromReco (detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, art::Event const &evt, std::string pfparticleModule) const
const simb::MCParticle *	GetMCParticleFromReco (detinfo::DetectorClocksData const &clockData, const recob::Track &track, art::Event const &evt, std::string trackModule) const
const simb::MCParticle *	GetMCParticleFromReco (detinfo::DetectorClocksData const &clockData, const recob::Shower &shower, art::Event const &evt, std::string showerModule) const
std::vector< std::pair< const simb::MCParticle *, double > >	GetMCParticleListFromReco (detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, art::Event const &evt, std::string pfparticleModule) const
std::vector< std::pair< const simb::MCParticle *, double > >	GetMCParticleListFromReco (detinfo::DetectorClocksData const &clockData, const recob::Track &track, art::Event const &evt, std::string trackModule) const
std::vector< std::pair< const simb::MCParticle *, double > >	GetMCParticleListFromReco (detinfo::DetectorClocksData const &clockData, const recob::Shower &shower, art::Event const &evt, std::string showerModule) const
template<typename T >
std::vector< MCParticleSharedHits >	GetMCParticleListByHits (detinfo::DetectorClocksData const &clockData, const T &recobj, const art::Event &evt, std::string recoModule, std::string hitModule) const
template<typename T >
const MCParticleSharedHits	GetMCParticleByHits (detinfo::DetectorClocksData const &clockData, const T &recobj, const art::Event &evt, std::string recoModule, std::string hitModule) const
const simb::MCParticle *	MatchPduneMCtoG4 (const simb::MCParticle &pDunePart, const art::Event &evt)
const simb::MCParticle *	GetGeantGoodParticle (const simb::MCTruth &genTruth, const art::Event &evt) const
const float	ConvertTrueTimeToPandoraTimeNano (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part) const
const float	ConvertTrueTimeToPandoraTimeNano (detinfo::DetectorClocksData const &clockData, const float trueTime) const
const float	ConvertTrueTimeToPandoraTimeMicro (detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part) const
const float	ConvertTrueTimeToPandoraTimeMicro (detinfo::DetectorClocksData const &clockData, const float trueTime) const
int	GetProcessKey (std::string process)
double	GetDepEnergyMC (const art::Event &evt, geo::GeometryCore const *fGeom, int trackid, int whichview) const
int	GetFirstTrajectoryPointInTPCActiveVolume (const simb::MCParticle &mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh)
double	GetMCParticleLengthInTPCActiveVolume (const simb::MCParticle &mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh)
double	GetDepEnergyAtLastTrajPoint (const simb::MCParticle &mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh)
double	GetKinEnergyAtVertex (const simb::MCParticle &mcpart, double kinene_lastpoint=0.0)
std::map< size_t, std::vector< const sim::IDE * > >	GetSimIDEs (const simb::MCParticle &mcpart)
std::vector< const sim::IDE * >	GetSimIDEsBetweenPoints (const simb::MCParticle &mcpart, const TLorentzVector &p1, const TLorentzVector &p2)
std::map< int, std::vector< int > >	GetMapMCToPFPs_ByHits (detinfo::DetectorClocksData const &clockData, const art::Event &evt, std::string pfpTag, std::string hitTag)
template<typename T >
std::vector< protoana::MCParticleSharedHits >	GetMCParticleListByHits (detinfo::DetectorClocksData const &clockData, const T &recobj, const art::Event &evt, std::string recoModule, std::string hitModule) const
template<typename T >
const protoana::MCParticleSharedHits	GetMCParticleByHits (detinfo::DetectorClocksData const &clockData, const T &recobj, const art::Event &evt, std::string recoModule, std::string hitModule) const

Detailed Description

Definition at line 39 of file ProtoDUNETruthUtils.h.

Constructor & Destructor Documentation

protoana::ProtoDUNETruthUtils::ProtoDUNETruthUtils

(

)

Definition at line 16 of file ProtoDUNETruthUtils.cxx.

                                                {

}

protoana::ProtoDUNETruthUtils::~ProtoDUNETruthUtils

(

)

Definition at line 20 of file ProtoDUNETruthUtils.cxx.

                                                 {

}

Member Function Documentation

const float protoana::ProtoDUNETruthUtils::ConvertTrueTimeToPandoraTimeMicro	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	part
	)		const

Definition at line 909 of file ProtoDUNETruthUtils.cxx.

                                                                                                              {
  return ConvertTrueTimeToPandoraTimeMicro(clockData, part.T());
}

const float protoana::ProtoDUNETruthUtils::ConvertTrueTimeToPandoraTimeMicro	(	detinfo::DetectorClocksData const &	clockData,
		const float	trueTime
	)		const

Definition at line 914 of file ProtoDUNETruthUtils.cxx.

                                                                                                        {

  return clockData.G4ToElecTime(trueTime);
}

const float protoana::ProtoDUNETruthUtils::ConvertTrueTimeToPandoraTimeNano	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	part
	)		const

Definition at line 898 of file ProtoDUNETruthUtils.cxx.

                                                                                                             {
  return ConvertTrueTimeToPandoraTimeNano(clockData, part.T());
}

const float protoana::ProtoDUNETruthUtils::ConvertTrueTimeToPandoraTimeNano	(	detinfo::DetectorClocksData const &	clockData,
		const float	trueTime
	)		const

Definition at line 903 of file ProtoDUNETruthUtils.cxx.

                                                                                                       {
  return 1000. * ConvertTrueTimeToPandoraTimeMicro(clockData, trueTime);
}

std::vector< const recob::Hit * > protoana::ProtoDUNETruthUtils::FillSharedHits	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	mcpart,
		const std::vector< const recob::Hit * > &	hitsVec,
		bool	delta_ray,
		bool	use_eve = true
	)		const

Definition at line 25 of file ProtoDUNETruthUtils.cxx.

                                         {

  std::vector<const recob::Hit*> outVec;

  // Push back all hits that are shared with this MCParticle
  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  for(const recob::Hit* hitp : hitsVec) {
    if (use_eve) {
      for(const sim::TrackIDE & ide : bt_serv->HitToEveTrackIDEs(clockData, *hitp)) {
        int trackId = ide.trackID;
        if( !delta_ray && ( trackId == mcpart.TrackId() ) ){
          outVec.push_back(hitp);
          break;
        }
        else if( delta_ray && ( trackId == (-1*mcpart.TrackId()) ) ) {
          outVec.push_back(hitp);
          break;
        }
      }
    }
    else {
      for(const int & trackId : bt_serv->HitToTrackIds(clockData, *hitp)) {
        if( !delta_ray && ( trackId == mcpart.TrackId() ) ){
          outVec.push_back(hitp);
          break;
        }
        else if( delta_ray && ( trackId == (-1*mcpart.TrackId()) ) ) {
          outVec.push_back(hitp);
          break;
        }
      }
    }
  }

  return outVec;
}

template<typename T >

double protoana::ProtoDUNETruthUtils::GetCompleteness	(	detinfo::DetectorClocksData const &	clockData,
		const T &	recobj,
		const art::Event &	evt,
		std::string	recoModule,
		std::string	hitModule
	)		const

Definition at line 150 of file ProtoDUNETruthUtils.cxx.

{
  // Find the hits of the corresponding MCParticle.
  const simb::MCParticle* mcmatch = GetMCParticleFromReco(clockData, recobj, evt, recoModule);
  if(mcmatch == 0x0) { return 0; }
  // std::cout << "MCmatch: " << mcmatch << '\n';
  const std::vector<const recob::Hit*> MCHits = GetMCParticleHits(clockData, *mcmatch, evt, hitModule);
  if(MCHits.size() == 0) { return 0; }
  // Get shared hits between reco object and MCParticle.
  const std::vector<const recob::Hit*> sharedHits = GetSharedHits(clockData, *mcmatch, recobj, evt, recoModule);

  return (double)sharedHits.size() / MCHits.size();
}

double protoana::ProtoDUNETruthUtils::GetDepEnergyAtLastTrajPoint	(	const simb::MCParticle &	mcpart,
		double	tpcactiveXlow,
		double	tpcactiveXhigh,
		double	tpcactiveYlow,
		double	tpcactiveYhigh,
		double	tpcactiveZlow,
		double	tpcactiveZhigh
	)

Definition at line 1029 of file ProtoDUNETruthUtils.cxx.

                                                                                                                                                                                                                                  {

  // First trajectory point inside TPC active volume
  int firstpoint = GetFirstTrajectoryPointInTPCActiveVolume(mcpart, tpcactiveXlow, tpcactiveXhigh, tpcactiveYlow, tpcactiveYhigh, tpcactiveZlow, tpcactiveZhigh);
  if(firstpoint < 0) return 0.0;

  // True track length in TPC active volume
  double length = GetMCParticleLengthInTPCActiveVolume(mcpart, tpcactiveXlow, tpcactiveXhigh, tpcactiveYlow, tpcactiveYhigh, tpcactiveZlow, tpcactiveZhigh);

  // Get the true track length up to the next to the last trajectory point
  int ntrajpoints = mcpart.NumberTrajectoryPoints();
  double plength = 0.0;
  TVector3 pos =  mcpart.Position(firstpoint).Vect();
  for(int i = firstpoint+1; i < ntrajpoints-1; ++i) {
    if(mcpart.Vx(i) >= tpcactiveXlow && mcpart.Vx(i) <= tpcactiveXhigh && mcpart.Vy(i) >= tpcactiveYlow && mcpart.Vy(i) <= tpcactiveYhigh && mcpart.Vz(i) >= tpcactiveZlow && mcpart.Vz(i) <= tpcactiveZhigh){

      pos -= mcpart.Position(i).Vect();
      plength += pos.Mag();
      pos = mcpart.Position(i).Vect();
    }
  }
  if(plength == 0.0) return 0.0;

  double length_laststep = length - plength;
  if(length_laststep <= 0.0) return 0.0;
  
  // Kinetic energy of the track
  double kinene_start = std::sqrt(mcpart.P(firstpoint)*mcpart.P(firstpoint) + mcpart.Mass()*mcpart.Mass()) - mcpart.Mass();
  double kinene_end = std::sqrt(mcpart.P(ntrajpoints-2)*mcpart.P(ntrajpoints-2) + mcpart.Mass()*mcpart.Mass()) - mcpart.Mass();
  double eneloss = (kinene_start - kinene_end);
  if(eneloss < 0.0 || kinene_start <= 0.0) return 0.0;

  // Average energy loss due to last trajectory point
  double aveneloss_laststep = length_laststep*eneloss/plength;

  return aveneloss_laststep;
  
}

double protoana::ProtoDUNETruthUtils::GetDepEnergyMC	(	const art::Event &	evt,
		geo::GeometryCore const *	fGeom,
		int	trackid,
		int	whichview
	)		const

Definition at line 964 of file ProtoDUNETruthUtils.cxx.

                                                                                                                                       {

  double edep = 0.0;

  auto simchannelHandle = evt.getHandle< std::vector<sim::SimChannel> >("largeant");
  if (simchannelHandle) {
    // Loop over sim channels
    for(auto const& simchannel : (*simchannelHandle)){
      // Only keep channels in the selected view
      if(fGeom->View(simchannel.Channel()) != whichview) continue;
      // Get all time slices
      auto const& alltimeslices = simchannel.TDCIDEMap();
      // Loop over time slices
      for(auto const& tslice : alltimeslices){
        auto const& simide = tslice.second;
        // Loop over energy deposits
        for(auto const& eDep : simide){
          if(eDep.trackID == trackid || eDep.trackID == -trackid)
            edep += eDep.energy;
        }
      }
    }
  }

  return edep;

}

int protoana::ProtoDUNETruthUtils::GetFirstTrajectoryPointInTPCActiveVolume	(	const simb::MCParticle &	mcpart,
		double	tpcactiveXlow,
		double	tpcactiveXhigh,
		double	tpcactiveYlow,
		double	tpcactiveYhigh,
		double	tpcactiveZlow,
		double	tpcactiveZhigh
	)

Definition at line 993 of file ProtoDUNETruthUtils.cxx.

                                                                                                                                                                                                                                            {

  int firstpoint = -999;
  for(unsigned int i = 0; i < mcpart.NumberTrajectoryPoints(); ++i) {
    if(mcpart.Vx(i) >= tpcactiveXlow && mcpart.Vx(i) <= tpcactiveXhigh && mcpart.Vy(i) >= tpcactiveYlow && mcpart.Vy(i) <= tpcactiveYhigh && mcpart.Vz(i) >= tpcactiveZlow && mcpart.Vz(i) <= tpcactiveZhigh){

      firstpoint = i;
      break;
    }
  }

  return firstpoint;
}

const simb::MCParticle * protoana::ProtoDUNETruthUtils::GetGeantGoodParticle	(	const simb::MCTruth &	genTruth,
		const art::Event &	evt
	)		const

Definition at line 860 of file ProtoDUNETruthUtils.cxx.

                                                                                                                             {

  // Get the good particle from the MCTruth
  simb::MCParticle goodPart;
  bool found = false;
  for(int t = 0; t < genTruth.NParticles(); ++t){
    simb::MCParticle part = genTruth.GetParticle(t);
    if(part.Process() == "primary"){
      goodPart = part;
      found = true;
      break;
    }
  }

  if(!found){
    std::cerr << "No good particle found, returning null pointer" << std::endl;
    return nullptr;
  }

  // Now loop over geant particles to find the one that matches
  // Get list of the g4 particles. plist should be a std::map< int, simb::MCParticle* >
  art::ServiceHandle< cheat::ParticleInventoryService > pi_serv;
  const sim::ParticleList & plist = pi_serv->ParticleList();

  for(auto const part : plist){
    if((goodPart.PdgCode() == part.second->PdgCode()) && fabs(part.second->E() - goodPart.E()) < 1e-5){
      return part.second;
    }
  }

  // If we get here something has gone wrong
  std::cerr << "No G4 version of the good particle was found, returning null pointer" << std::endl;
  return nullptr;

}

double protoana::ProtoDUNETruthUtils::GetKinEnergyAtVertex	(	const simb::MCParticle &	mcpart,
		double	kinene_lastpoint = 0.0
	)

Definition at line 1069 of file ProtoDUNETruthUtils.cxx.

                                                                                                             {

  int ntrajpoints = mcpart.NumberTrajectoryPoints();
  double kinene_end = std::sqrt(mcpart.P(ntrajpoints-2)*mcpart.P(ntrajpoints-2) + mcpart.Mass()*mcpart.Mass()) - mcpart.Mass();
  
  return (kinene_end - kinene_lastpoint);
}

std::map< int, std::vector< int > > protoana::ProtoDUNETruthUtils::GetMapMCToPFPs_ByHits	(	detinfo::DetectorClocksData const &	clockData,
		const art::Event &	evt,
		std::string	pfpTag,
		std::string	hitTag
	)

Definition at line 1140 of file ProtoDUNETruthUtils.cxx.

                                                                                                                                                   {

  std::map< int, std::vector< int > > results;

  auto allPFPs = evt.getValidHandle< std::vector< recob::PFParticle > >( pfpTag );

  for( size_t i = 0; i < allPFPs->size(); ++i ){
    auto thePFP = &(allPFPs->at(i));
    protoana::MCParticleSharedHits match = GetMCParticleByHits( clockData, *thePFP, evt, pfpTag, hitTag );
    if( match.particle ){
      results[ match.particle->TrackId() ].push_back( i );
    }
  }

  return results;
}

template<typename T >

const MCParticleSharedHits protoana::ProtoDUNETruthUtils::GetMCParticleByHits	(	detinfo::DetectorClocksData const &	clockData,
		const T &	recobj,
		const art::Event &	evt,
		std::string	recoModule,
		std::string	hitModule
	)		const

template<typename T >

const protoana::MCParticleSharedHits protoana::ProtoDUNETruthUtils::GetMCParticleByHits	(	detinfo::DetectorClocksData const &	clockData,
		const T &	recobj,
		const art::Event &	evt,
		std::string	recoModule,
		std::string	hitModule
	)		const

Definition at line 788 of file ProtoDUNETruthUtils.cxx.

                                                                                                                                                                                                        {

  auto list = GetMCParticleListByHits(clockData, recobj, evt, recoModule, hitModule );

  protoana::MCParticleSharedHits dummy;
  dummy.particle = 0x0;
  dummy.nSharedHits = 0;
  dummy.nSharedDeltaRayHits = 0;
  //if( !list.size() ) return 0x0;
  if( !list.size() ) return dummy;

  //return list[0].first;
  //return list[0].particle;
  return list[0];
}

const simb::MCParticle * protoana::ProtoDUNETruthUtils::GetMCParticleFromPFParticle	(	detinfo::DetectorClocksData const &	clockData,
		const recob::PFParticle &	pfpart,
		art::Event const &	evt,
		std::string	pfparticleModule
	)		const

Definition at line 565 of file ProtoDUNETruthUtils.cxx.

                                                                                         {

  const simb::MCParticle* outPart = 0x0;

  const std::vector<std::pair<const simb::MCParticle*, double>> inVec =
    GetMCParticleListFromPFParticle(clockData, pfpart, evt, pfparticleModule);

  if (inVec.size() != 0) outPart = inVec[0].first;

  return outPart;
}

const simb::MCParticle * protoana::ProtoDUNETruthUtils::GetMCParticleFromReco	(	detinfo::DetectorClocksData const &	clockData,
		const recob::PFParticle &	pfpart,
		art::Event const &	evt,
		std::string	pfparticleModule
	)		const

Definition at line 696 of file ProtoDUNETruthUtils.cxx.

                                                                                         {
  return GetMCParticleFromPFParticle(clockData, pfpart, evt, pfparticleModule);
}

const simb::MCParticle * protoana::ProtoDUNETruthUtils::GetMCParticleFromReco	(	detinfo::DetectorClocksData const &	clockData,
		const recob::Track &	track,
		art::Event const &	evt,
		std::string	trackModule
	)		const

Definition at line 701 of file ProtoDUNETruthUtils.cxx.

                                                                              {
  return GetMCParticleFromRecoTrack(clockData, track, evt, trackModule);
}

const simb::MCParticle * protoana::ProtoDUNETruthUtils::GetMCParticleFromReco	(	detinfo::DetectorClocksData const &	clockData,
		const recob::Shower &	shower,
		art::Event const &	evt,
		std::string	showerModule
	)		const

Definition at line 706 of file ProtoDUNETruthUtils.cxx.

                                                                                 {
  return GetMCParticleFromRecoShower(clockData, shower, evt, showerModule);
}

const simb::MCParticle * protoana::ProtoDUNETruthUtils::GetMCParticleFromRecoShower	(	detinfo::DetectorClocksData const &	clockData,
		const recob::Shower &	shower,
		art::Event const &	evt,
		std::string	showerModule
	)		const

Definition at line 653 of file ProtoDUNETruthUtils.cxx.

                                                                                 {

  const simb::MCParticle* outPart = 0x0;

  const std::vector<std::pair<const simb::MCParticle*, double>> inVec =
    GetMCParticleListFromRecoShower(clockData, shower, evt, showerModule);

  if (inVec.size() != 0) outPart = inVec[0].first;

  return outPart;
}

const simb::MCParticle * protoana::ProtoDUNETruthUtils::GetMCParticleFromRecoTrack	(	detinfo::DetectorClocksData const &	clockData,
		const recob::Track &	track,
		art::Event const &	evt,
		std::string	trackModule
	)		const

Definition at line 609 of file ProtoDUNETruthUtils.cxx.

                                                                              {

  const simb::MCParticle* outPart = 0x0;

  const std::vector<std::pair<const simb::MCParticle*, double>> inVec =
    GetMCParticleListFromRecoTrack(clockData, track, evt, trackModule);

  if (inVec.size() != 0) outPart = inVec[0].first;

  return outPart;
}

std::vector< const recob::Hit * > protoana::ProtoDUNETruthUtils::GetMCParticleHits	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	mcpart,
		const art::Event &	evt,
		std::string	hitModule,
		bool	use_eve = true
	)		const

Definition at line 103 of file ProtoDUNETruthUtils.cxx.

                        {
  std::vector<const recob::Hit*> outVec;
  
  auto hitHandle = evt.getHandle<std::vector<recob::Hit> >(hitModule);
  if (!hitHandle) {
    std::cout << "protoana::ProtoDUNETruthUtils::GetMCParticleHits: could not find hits in event.\n";
    return outVec;
  }

  // Backtrack all hits to verify whether they belong to the current MCParticle.
  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
  for(const recob::Hit& hit : *hitHandle) {
    if (use_eve) {
      for(const sim::TrackIDE & ide : bt_serv->HitToEveTrackIDEs(clockData, hit)) {
        int trackId = ide.trackID;
        if(pi_serv->TrackIdToParticle_P(trackId) == 
           pi_serv->TrackIdToParticle_P(mcpart.TrackId())) {
          outVec.push_back(&hit);
          break;
        }
      }
    }
    else {
      for(const int trackId : bt_serv->HitToTrackIds(clockData, hit)) {
        if(pi_serv->TrackIdToParticle_P(trackId) == 
           pi_serv->TrackIdToParticle_P(mcpart.TrackId())) {
          outVec.push_back(&hit);
          break;
        }
      }
    }
  }

  return outVec;
}

double protoana::ProtoDUNETruthUtils::GetMCParticleLengthInTPCActiveVolume	(	const simb::MCParticle &	mcpart,
		double	tpcactiveXlow,
		double	tpcactiveXhigh,
		double	tpcactiveYlow,
		double	tpcactiveYhigh,
		double	tpcactiveZlow,
		double	tpcactiveZhigh
	)

Definition at line 1008 of file ProtoDUNETruthUtils.cxx.

                                                                                                                                                                                                                                           {

  double length = 0.0;
  int firstpoint = GetFirstTrajectoryPointInTPCActiveVolume(mcpart, tpcactiveXlow, tpcactiveXhigh, tpcactiveYlow, tpcactiveYhigh, tpcactiveZlow, tpcactiveZhigh);

  if(firstpoint < 0) return length;

  TVector3 pos =  mcpart.Position(firstpoint).Vect();
  for(unsigned int i = firstpoint+1; i < mcpart.NumberTrajectoryPoints(); ++i) {
    if(mcpart.Vx(i) >= tpcactiveXlow && mcpart.Vx(i) <= tpcactiveXhigh && mcpart.Vy(i) >= tpcactiveYlow && mcpart.Vy(i) <= tpcactiveYhigh && mcpart.Vz(i) >= tpcactiveZlow && mcpart.Vz(i) <= tpcactiveZhigh){

      pos -= mcpart.Position(i).Vect();
      length += pos.Mag();
      pos = mcpart.Position(i).Vect();
    }
  }

  return length;
}

template<typename T >

std::vector< MCParticleSharedHits > protoana::ProtoDUNETruthUtils::GetMCParticleListByHits	(	detinfo::DetectorClocksData const &	clockData,
		const T &	recobj,
		const art::Event &	evt,
		std::string	recoModule,
		std::string	hitModule
	)		const

template<typename T >

std::vector< protoana::MCParticleSharedHits > protoana::ProtoDUNETruthUtils::GetMCParticleListByHits	(	detinfo::DetectorClocksData const &	clockData,
		const T &	recobj,
		const art::Event &	evt,
		std::string	recoModule,
		std::string	hitModule
	)		const

Definition at line 737 of file ProtoDUNETruthUtils.cxx.

                                                                                                                                                {
  
  //Get the particle list. This is sorted by energy contributed, but we'll turn this into nHits
  const std::vector< std::pair< const simb::MCParticle*, double > > mcparts = GetMCParticleListFromReco(clockData, recobj, evt, recoModule );

  //using weightedMCPair = std::pair<const simb::MCParticle*, size_t>;
  //std::vector< weightedMCPair > results;
  std::vector< protoana::MCParticleSharedHits > results;

  for( const std::pair< const simb::MCParticle*, double > part : mcparts ){
    //Get the non-delta ray hits
    const std::vector<const recob::Hit*> hits = GetSharedHits(clockData, *(part.first), recobj, evt, recoModule);

    //Get the delta ray hits
    const std::vector<const recob::Hit*> dr_hits = GetSharedHits(clockData, *(part.first), recobj, evt, recoModule, true);

    protoana::MCParticleSharedHits entry;
    entry.particle = part.first;
    entry.nSharedHits = hits.size();
    entry.nSharedDeltaRayHits = dr_hits.size();
    results.push_back( entry );
    //results.push_back( std::make_pair( part.first, hits.size() ) );
  }

  //std::sort(results.begin(), results.end(),
  //      [](weightedMCPair a, weightedMCPair b){ return a.second > b.second;});
  
  std::sort( results.begin(), results.end(),
        [](protoana::MCParticleSharedHits a, protoana::MCParticleSharedHits b)
          { return ( (a.nSharedHits + a.nSharedDeltaRayHits) > (b.nSharedHits + b.nSharedDeltaRayHits) ); });

  return results;
}

std::vector< std::pair< const simb::MCParticle *, double > > protoana::ProtoDUNETruthUtils::GetMCParticleListFromPFParticle	(	detinfo::DetectorClocksData const &	clockData,
		const recob::PFParticle &	pfpart,
		art::Event const &	evt,
		std::string	pfparticleModule
	)		const

Definition at line 325 of file ProtoDUNETruthUtils.cxx.

                                                                                          {

  std::vector<std::pair<const simb::MCParticle*, double>> outVec;

  // We must have MC for this module to make sense
  if(evt.isRealData()) return outVec;

  // Get the hits associated with this PFParticle
  protoana::ProtoDUNEPFParticleUtils pfpUtils;
  std::vector<const recob::Hit*> pfpHits = pfpUtils.GetPFParticleHits(pfpart, evt, pfparticleModule);

  // Call the appropriate function to determine MC matches.
  if(abs(pfpart.PdgCode()) == 11 || pfpart.PdgCode() == 22) {
    outVec = GetMCParticleListFromShowerHits(clockData, pfpHits);
  } else {
    outVec = GetMCParticleListFromTrackHits(clockData, pfpHits);
  }

  return outVec;
}

std::vector< std::pair< const simb::MCParticle *, double > > protoana::ProtoDUNETruthUtils::GetMCParticleListFromReco	(	detinfo::DetectorClocksData const &	clockData,
		const recob::PFParticle &	pfpart,
		art::Event const &	evt,
		std::string	pfparticleModule
	)		const

Definition at line 714 of file ProtoDUNETruthUtils.cxx.

                                                                                         {
  return GetMCParticleListFromPFParticle(clockData, pfpart, evt, pfparticleModule);
}

std::vector< std::pair< const simb::MCParticle *, double > > protoana::ProtoDUNETruthUtils::GetMCParticleListFromReco	(	detinfo::DetectorClocksData const &	clockData,
		const recob::Track &	track,
		art::Event const &	evt,
		std::string	trackModule
	)		const

Definition at line 720 of file ProtoDUNETruthUtils.cxx.

                                                                              {
  return GetMCParticleListFromRecoTrack(clockData, track, evt, trackModule);
}

std::vector< std::pair< const simb::MCParticle *, double > > protoana::ProtoDUNETruthUtils::GetMCParticleListFromReco	(	detinfo::DetectorClocksData const &	clockData,
		const recob::Shower &	shower,
		art::Event const &	evt,
		std::string	showerModule
	)		const

Definition at line 726 of file ProtoDUNETruthUtils.cxx.

                                                                                 {
  return GetMCParticleListFromRecoShower(clockData, shower, evt, showerModule);
}

std::vector< std::pair< const simb::MCParticle *, double > > protoana::ProtoDUNETruthUtils::GetMCParticleListFromRecoShower	(	detinfo::DetectorClocksData const &	clockData,
		const recob::Shower &	shower,
		art::Event const &	evt,
		std::string	showerModule
	)		const

Definition at line 487 of file ProtoDUNETruthUtils.cxx.

                                                                                 {
    
  // We must have MC for this module to make sense
  if(evt.isRealData()) return {};

  // Get the reconstructed shower hits
  protoana::ProtoDUNEShowerUtils showerUtils;
  std::vector<const recob::Hit*> showerHits = showerUtils.GetRecoShowerHits(shower, evt, showerModule);
  return GetMCParticleListFromShowerHits(clockData, showerHits);
}

std::vector< std::pair< const simb::MCParticle *, double > > protoana::ProtoDUNETruthUtils::GetMCParticleListFromRecoTrack	(	detinfo::DetectorClocksData const &	clockData,
		const recob::Track &	track,
		art::Event const &	evt,
		std::string	trackModule
	)		const

Definition at line 407 of file ProtoDUNETruthUtils.cxx.

                                                                               {

  // We must have MC for this module to make sense
  if(evt.isRealData()) return {};

  // Get the reconstructed track hits
  protoana::ProtoDUNETrackUtils trackUtils;
  std::vector<const recob::Hit*> trackHits = trackUtils.GetRecoTrackHits(track, evt, trackModule);
  return GetMCParticleListFromTrackHits(clockData, trackHits);
}

std::vector< std::pair< const simb::MCParticle *, double > > protoana::ProtoDUNETruthUtils::GetMCParticleListFromShowerHits	(	detinfo::DetectorClocksData const &	clockData,
		const std::vector< const recob::Hit * > &	hitVec
	)		const

Definition at line 286 of file ProtoDUNETruthUtils.cxx.

                                                   {

  using weightedMCPair = std::pair<const simb::MCParticle*, double>;
  std::vector<weightedMCPair> outVec;
  
  // Loop over all hits in the input vector and record the contributing MCParticles.
  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
  std::unordered_map<const simb::MCParticle*, double> mcEMap;
  double hitTotalE = 0;
  for(const recob::Hit* hitp : hitVec) {
    for(const sim::TrackIDE& ide : bt_serv->HitToEveTrackIDEs(clockData, *hitp)) {
      const simb::MCParticle* curr_part = pi_serv->TrackIdToParticle_P(ide.trackID);
      mcEMap[curr_part] += ide.energy;
      hitTotalE += ide.energy;
    }
  }

  // Fill and sort the output vector
  for (weightedMCPair const& p : mcEMap) {
    outVec.push_back(p);
  }
  std::sort(outVec.begin(), outVec.end(),
        [](weightedMCPair a, weightedMCPair b){ return a.second > b.second;});

  // Normalise the weights by the total track energy.
  if (hitTotalE < 1e-5) { hitTotalE = 1; } // Protect against zero division
  for (weightedMCPair& p : outVec) {
    p.second /= hitTotalE;
  }

  return outVec;
}

std::vector< std::pair< const simb::MCParticle *, double > > protoana::ProtoDUNETruthUtils::GetMCParticleListFromTrackHits	(	detinfo::DetectorClocksData const &	clockData,
		const std::vector< const recob::Hit * > &	hitVec,
		bool	use_eve = true
	)		const

Definition at line 238 of file ProtoDUNETruthUtils.cxx.

                                                                 {

  using weightedMCPair = std::pair<const simb::MCParticle*, double>;
  std::vector<weightedMCPair> outVec;
  
  // Loop over all hits in the input vector and record the contributing MCParticles.
  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
  std::unordered_map<const simb::MCParticle*, double> mcEMap;
  double hitTotalE = 0;
  for(const recob::Hit* hitp : hitVec) {
    if (use_eve) {
      for(const sim::TrackIDE& ide : bt_serv->HitToEveTrackIDEs(clockData, *hitp)) {
        const simb::MCParticle* curr_part = pi_serv->TrackIdToParticle_P(ide.trackID);
        mcEMap[curr_part] += ide.energy;
        hitTotalE += ide.energy;
      }
    }
    else {
      for(const sim::TrackIDE& ide : bt_serv->HitToTrackIDEs(clockData, *hitp)) {
        const simb::MCParticle* curr_part = pi_serv->TrackIdToParticle_P(ide.trackID);
        mcEMap[curr_part] += ide.energy;
        hitTotalE += ide.energy;
      }
    }
  }

  // Fill and sort the output vector
  for (weightedMCPair const& p : mcEMap) {
    outVec.push_back(p);
  }
  std::sort(outVec.begin(), outVec.end(),
        [](weightedMCPair a, weightedMCPair b){ return a.second > b.second;});

  // Normalise the weights by the total track energy.
  if (hitTotalE < 1e-5) { hitTotalE = 1; } // Protect against zero division
  for (weightedMCPair& p : outVec) {
    p.second /= hitTotalE;
  }

  return outVec;
}

const recob::PFParticle * protoana::ProtoDUNETruthUtils::GetPFParticleFromMCParticle	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	part,
		art::Event const &	evt,
		std::string	pfparticleModule
	)		const

Definition at line 582 of file ProtoDUNETruthUtils.cxx.

                                                                                       {

  const recob::PFParticle* outPFPart = 0x0;

  // We must have MC for this module to make sense
  if(evt.isRealData()) return outPFPart;

  // Get the list of contributing PFParticles for the MCParticle and see if any of
  // them have the MCParticle as a primary contributor.
  for (std::pair<const recob::PFParticle*, double> p :
         GetPFParticleListFromMCParticle(clockData, part, evt, pfparticleModule))
  {
    const recob::PFParticle* pfp = p.first;
    if (GetMCParticleFromPFParticle(clockData, *pfp, evt, pfparticleModule)->TrackId() == part.TrackId())
    {
      outPFPart = pfp;
      break;
    }
  }

  return outPFPart;
}

std::vector< std::pair< const recob::PFParticle *, double > > protoana::ProtoDUNETruthUtils::GetPFParticleListFromMCParticle	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	part,
		art::Event const &	evt,
		std::string	pfparticleModule,
		bool	use_eve = true
	)		const

Definition at line 349 of file ProtoDUNETruthUtils.cxx.

                                                                                                     {

  using weightedPFPair = std::pair<const recob::PFParticle*, double>;

  std::vector<weightedPFPair> outVec;

  // We must have MC for this module to make sense
  if(evt.isRealData()) return outVec;

  // Get the reconstructed PFParticles
  auto allPFParticles = evt.getValidHandle<std::vector<recob::PFParticle> >(pfparticleModule);

  // Loop over all PFParticles and save the energetic overlap with each MCParticle in a map
  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  std::unordered_map<const recob::PFParticle*, double> pfpEMap;
  double pfpTotalE = 0;
  for(const recob::PFParticle& pfpart : *allPFParticles) {
    for(const recob::Hit* hitp : GetSharedHits(clockData, part, pfpart, evt, pfparticleModule)) {
      if (use_eve) {
        for(const sim::TrackIDE& ide : bt_serv->HitToEveTrackIDEs(clockData, *hitp)) {
          if(ide.trackID == part.TrackId()) {
            pfpEMap[&pfpart] += ide.energy;
          }
          pfpTotalE += ide.energy;
        }
      }
      else {
        for(const sim::TrackIDE& ide : bt_serv->HitToTrackIDEs(clockData, *hitp)) {
          if(ide.trackID == part.TrackId()) {
            pfpEMap[&pfpart] += ide.energy;
          }
          pfpTotalE += ide.energy;
        }
      }
    }
  }

  // Fill and sort the output vector
  for (weightedPFPair const& p : pfpEMap) {
    outVec.push_back(p);
  }
  std::sort(outVec.begin(), outVec.end(),
        [](weightedPFPair a, weightedPFPair b){ return a.second > b.second;});

  // Normalise the weights by the total track energy.
  if (pfpTotalE < 1e-5) { pfpTotalE = 1; } // Protect against zero division
  for (weightedPFPair& p : outVec) {
    p.second /= pfpTotalE;
  }

  return outVec;
}

int protoana::ProtoDUNETruthUtils::GetProcessKey

(

std::string

process

)

Definition at line 921 of file ProtoDUNETruthUtils.cxx.

                                                               {

  if(process.compare("primary") == 0)                    return 0;
  else if(process.compare("hadElastic") == 0)            return 1;
  else if(process.compare("pi-Inelastic") == 0)          return 2;
  else if(process.compare("pi+Inelastic") == 0)          return 3;
  else if(process.compare("kaon-Inelastic") == 0)        return 4;
  else if(process.compare("kaon+Inelastic") == 0)        return 5;
  else if(process.compare("protonInelastic") == 0)       return 6;
  else if(process.compare("neutronInelastic") == 0)      return 7;
  else if(process.compare("kaon0SInelastic") == 0)       return 8;
  else if(process.compare("kaon0LInelastic") == 0)       return 9;
  else if(process.compare("lambdaInelastic") == 0)       return 10;
  else if(process.compare("omega-Inelastic") == 0)       return 11;
  else if(process.compare("sigma+Inelastic") == 0)       return 12;
  else if(process.compare("sigma-Inelastic") == 0)       return 13;
  else if(process.compare("sigma0Inelastic") == 0)       return 14;
  else if(process.compare("xi-Inelastic") == 0)          return 15;
  else if(process.compare("xi0Inelastic") == 0)          return 16;
  else if(process.compare("anti_protonInelastic") == 0)  return 20;
  else if(process.compare("anti_neutronInelastic") == 0) return 21;
  else if(process.compare("anti_lambdaInelastic") == 0)  return 22;
  else if(process.compare("anti_omega-Inelastic") == 0)  return 23;
  else if(process.compare("anti_sigma+Inelastic") == 0)  return 24;
  else if(process.compare("anti_sigma-Inelastic") == 0)  return 25;
  else if(process.compare("anti_xi-Inelastic") == 0)     return 26;
  else if(process.compare("anti_xi0Inelastic") == 0)     return 27;

  else if(process.compare("Decay") == 0)                 return 30;
  else if(process.compare("FastScintillation") == 0)     return 31;
  else if(process.compare("nKiller") == 0)               return 32; // Remove unwanted neutrons: neutron kinetic energy threshold (default 0) or time limit for neutron track
  else if(process.compare("nCapture") == 0)              return 33; // Neutron capture

  else if(process.compare("compt") == 0)                 return 40; // Compton Scattering
  else if(process.compare("rayleigh") == 0)              return 41; // Rayleigh Scattering
  else if(process.compare("phot") == 0)                  return 42; // Photoelectric Effect
  else if(process.compare("conv") == 0)                  return 43; // Pair production
  else if(process.compare("CoupledTransportation") == 0) return 44; //

  else return -1;
}

double protoana::ProtoDUNETruthUtils::GetPurity	(	detinfo::DetectorClocksData const &	clockData,
		const recob::PFParticle &	pfpart,
		const art::Event &	evt,
		std::string	pfparticleModule
	)		const

Definition at line 177 of file ProtoDUNETruthUtils.cxx.

                                                         {
  
  // Find the corresponding MCParticle and shared hits.
  const simb::MCParticle* mcmatch = GetMCParticleFromReco(clockData, pfpart, evt, pfparticleModule);
  if(mcmatch == 0x0) { return 0; }
  // Get shared hits between PFParticle and MCParticle.
  const std::vector<const recob::Hit*> sharedHits = GetSharedHits(clockData, *mcmatch, pfpart, evt, pfparticleModule);

  // Get the hits associated with this PFParticle
  protoana::ProtoDUNEPFParticleUtils pfpUtils;
  std::vector<const recob::Hit*> pfpHits = pfpUtils.GetPFParticleHits(pfpart, evt, pfparticleModule);

  return (double)sharedHits.size() / pfpHits.size();
}

double protoana::ProtoDUNETruthUtils::GetPurity	(	detinfo::DetectorClocksData const &	clockData,
		const recob::Track &	track,
		const art::Event &	evt,
		std::string	trackModule
	)		const

Definition at line 197 of file ProtoDUNETruthUtils.cxx.

                                                    {
  
  // Find the corresponding MCParticle and shared hits.
  const simb::MCParticle* mcmatch = GetMCParticleFromReco(clockData, track, evt, trackModule);
  if(mcmatch == 0x0) { return 0; }
  // Get shared hits between track and MCParticle.
  const std::vector<const recob::Hit*> sharedHits = GetSharedHits(clockData, *mcmatch, track, evt, trackModule);

  // Get the hits associated with this track
  protoana::ProtoDUNETrackUtils trackUtils;
  std::vector<const recob::Hit*> trackHits = trackUtils.GetRecoTrackHits(track, evt, trackModule);

  return (double)sharedHits.size() / trackHits.size();
}

double protoana::ProtoDUNETruthUtils::GetPurity	(	detinfo::DetectorClocksData const &	clockData,
		const recob::Shower &	shower,
		const art::Event &	evt,
		std::string	showerModule
	)		const

Definition at line 217 of file ProtoDUNETruthUtils.cxx.

                                                     {
  
  // Find the corresponding MCParticle and shared hits.
  const simb::MCParticle* mcmatch = GetMCParticleFromReco(clockData, shower, evt, showerModule);
  if(mcmatch == 0x0) { return 0; }
  // Get shared hits between shower and MCParticle.
  const std::vector<const recob::Hit*> sharedHits = GetSharedHits(clockData, *mcmatch, shower, evt, showerModule);

  // Get the hits associated with this shower
  protoana::ProtoDUNEShowerUtils showerUtils;
  std::vector<const recob::Hit*> showerHits = showerUtils.GetRecoShowerHits(shower, evt, showerModule);

  return (double)sharedHits.size() / showerHits.size();
}

const recob::Shower * protoana::ProtoDUNETruthUtils::GetRecoShowerFromMCParticle	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	part,
		art::Event const &	evt,
		std::string	showerModule
	)		const

Definition at line 670 of file ProtoDUNETruthUtils.cxx.

                                                                                   {

  const recob::Shower* outShower = 0x0;

  // We must have MC for this module to make sense
  if(evt.isRealData()) return outShower;

  // Get the list of contributing showers for the MCParticle and see if any of
  // them have the MCParticle as a primary contributor.
  for (std::pair<const recob::Shower*, double> p :
         GetRecoShowerListFromMCParticle(clockData, part, evt, showerModule))
  {
    const recob::Shower* tr = p.first;
    if (GetMCParticleFromRecoShower(clockData, *tr, evt, showerModule)->TrackId() == part.TrackId())
    {
      outShower = tr;
      break;
    }
  }

  return outShower;
}

std::vector< std::pair< const recob::Shower *, double > > protoana::ProtoDUNETruthUtils::GetRecoShowerListFromMCParticle	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	part,
		art::Event const &	evt,
		std::string	showerModule
	)		const

Definition at line 504 of file ProtoDUNETruthUtils.cxx.

                                                                                  {

  using weightedShowerPair = std::pair<const recob::Shower*, double>;

  std::vector<weightedShowerPair> outVec;

  // We must have MC for this module to make sense
  if(evt.isRealData()) return outVec;

  // Get the reconstructed showers
  auto allRecoShowers = evt.getValidHandle<std::vector<recob::Shower> >(showerModule);

  // We need the association between the showers and the hits
  const art::FindManyP<recob::Hit> findShowerHits(allRecoShowers, evt, showerModule);

  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  ProtoDUNEShowerUtils shUtils;
  std::unordered_map<int, double> recoShower;

  // Record the energy contribution to the MCParticle of every relevant reco shower
  double part_E = 0;
  for (recob::Shower const & shower : *allRecoShowers)
  {
    const unsigned showerIndex = shUtils.GetShowerIndex(shower, evt, showerModule);
    // Since MCParticles do not have shower hits associated with them, loop over all shower hits
    for (art::Ptr<recob::Hit> hit : findShowerHits.at(showerIndex))
    {
      // Loop over hit IDEs to find our particle's part in it
      for (const sim::TrackIDE& ide : bt_serv->HitToEveTrackIDEs(clockData, hit))
      {
        if (ide.trackID == part.TrackId())
        {
          recoShower[showerIndex] += ide.energy;
        }
        part_E += ide.energy;
      } // sim::IDE*
    } // art::Ptr<recob::Hit>
  } // const recob::Shower&

  // Fill and sort the output vector
  for (std::pair<int, double> const& p : recoShower)
  {
    auto const showerIt = std::find_if(allRecoShowers->begin(), allRecoShowers->end(),
      [&](recob::Shower sh){ return shUtils.GetShowerIndex(sh, evt, showerModule) == p.first; });
    outVec.push_back(std::make_pair(&*showerIt, p.second));
  }
  std::sort(outVec.begin(), outVec.end(),
    [](weightedShowerPair a, weightedShowerPair b){ return a.second > b.second;});

  // Normalise the vector weights
  if (part_E < 1e-5) { part_E = 1; } // Protect against zero division
  std::for_each(outVec.begin(), outVec.end(),
          [&](weightedShowerPair& p){ p.second /= part_E; });

  return outVec;
}

const recob::Track * protoana::ProtoDUNETruthUtils::GetRecoTrackFromMCParticle	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	part,
		art::Event const &	evt,
		std::string	trackModule
	)		const

Definition at line 626 of file ProtoDUNETruthUtils.cxx.

                                                                                  {

  const recob::Track* outTrack = 0x0;

  // We must have MC for this module to make sense
  if(evt.isRealData()) return outTrack;

  // Get the list of contributing tracks for the MCParticle and see if any of
  // them have the MCParticle as a primary contributor.
  for (std::pair<const recob::Track*, double> p :
         GetRecoTrackListFromMCParticle(clockData, part, evt, trackModule))
  {
    const recob::Track* tr = p.first;
    if (GetMCParticleFromRecoTrack(clockData, *tr, evt, trackModule)->TrackId() == part.TrackId())
    {
      outTrack = tr;
      break;
    }
  }

  return outTrack;
}

std::vector< std::pair< const recob::Track *, double > > protoana::ProtoDUNETruthUtils::GetRecoTrackListFromMCParticle	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	part,
		art::Event const &	evt,
		std::string	trackModule
	)		const

Definition at line 424 of file ProtoDUNETruthUtils.cxx.

                                                                                 {

  using weightedTrackPair = std::pair<const recob::Track*, double>;

  std::vector<weightedTrackPair> outVec;

  // We must have MC for this module to make sense
  if(evt.isRealData()) return outVec;

  // Get the reconstructed tracks
  auto allRecoTracks = evt.getValidHandle<std::vector<recob::Track> >(trackModule);

  // We need the association between the tracks and the hits
  const art::FindManyP<recob::Hit> findTrackHits(allRecoTracks, evt, trackModule);

  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  std::unordered_map<int, double> recoTrack;

  // Record the energy contribution to the MCParticle of every relevant reco track
  double part_E = 0;
  for (recob::Track const & track : *allRecoTracks)
  {
    // This finds all hits that are shared between a reconstructed track and MCParticle
    for (art::Ptr<recob::Hit> const & hptr :
          bt_serv->TrackIdToHits_Ps(clockData, part.TrackId(), findTrackHits.at(track.ID())))
    {
      // Loop over hit IDEs to find our particle's part in it
      for (const sim::IDE* ide : bt_serv->HitToSimIDEs_Ps(clockData, hptr))
      {
        if (ide->trackID == part.TrackId())
        {
          recoTrack[track.ID()] += ide->energy;
        }
        part_E += ide->energy;
      } // sim::IDE*
    } // art::Ptr<recob::Hit>
  } // const recob::Track&

  // Fill and sort the output vector
  for (std::pair<int, double> const& p : recoTrack)
  {
    auto const trackIt = std::find_if(allRecoTracks->begin(), allRecoTracks->end(),
                               [&](recob::Track tr){ return tr.ID() == p.first; });
    outVec.push_back(std::make_pair(&*trackIt, p.second));
  }
  std::sort(outVec.begin(), outVec.end(),
    [](weightedTrackPair a, weightedTrackPair b){ return a.second > b.second;});

  // Normalise the vector weights
  if (part_E < 1e-5) { part_E = 1; } // Protect against zero division
  for (weightedTrackPair& p : outVec)
  {
    p.second /= part_E;
  }

  return outVec;
}

std::vector< const recob::Hit * > protoana::ProtoDUNETruthUtils::GetSharedHits	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	mcpart,
		const recob::PFParticle &	pfpart,
		const art::Event &	evt,
		std::string	pfparticleModule,
		bool	delta_ray = false
	)		const

Definition at line 67 of file ProtoDUNETruthUtils.cxx.

                                                                                                                                       {

  // Get the hits associated with this PFParticle
  protoana::ProtoDUNEPFParticleUtils pfpUtils;
  std::vector<const recob::Hit*> pfpHits = pfpUtils.GetPFParticleHits(pfpart, evt, pfparticleModule);

  return FillSharedHits( clockData, mcpart, pfpHits, delta_ray );
}

std::vector< const recob::Hit * > protoana::ProtoDUNETruthUtils::GetSharedHits	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	mcpart,
		const recob::Track &	track,
		const art::Event &	evt,
		std::string	trackModule,
		bool	delta_ray = false
	)		const

Definition at line 80 of file ProtoDUNETruthUtils.cxx.

                                                                                                                            {

  // Get the hits associated with this track
  protoana::ProtoDUNETrackUtils trackUtils;
  std::vector<const recob::Hit*> trackHits = trackUtils.GetRecoTrackHits(track, evt, trackModule);
  return FillSharedHits( clockData, mcpart, trackHits, delta_ray );
}

std::vector< const recob::Hit * > protoana::ProtoDUNETruthUtils::GetSharedHits	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle &	mcpart,
		const recob::Shower &	shower,
		const art::Event &	evt,
		std::string	showerModule,
		bool	delta_ray = false
	)		const

Definition at line 92 of file ProtoDUNETruthUtils.cxx.

                                                                                                                               {

  // Get the hits associated with this shower
  protoana::ProtoDUNEShowerUtils shUtils;
  std::vector<const recob::Hit*> shHits = shUtils.GetRecoShowerHits(shower, evt, showerModule);
  return FillSharedHits( clockData, mcpart, shHits, delta_ray );
}

std::map< size_t, std::vector< const sim::IDE * > > protoana::ProtoDUNETruthUtils::GetSimIDEs

(

const simb::MCParticle &

mcpart

)

Definition at line 1079 of file ProtoDUNETruthUtils.cxx.

{
  art::ServiceHandle< cheat::BackTrackerService >       bt_serv;
  //art::ServiceHandle< geo::Geometry >                   geom;

  const simb::MCTrajectory & mctraj = mcpart.Trajectory();

  //Get all ides from the MCParticle
  std::vector< const sim::IDE * > ides = bt_serv->TrackIdToSimIDEs_Ps( mcpart.TrackId()/*, geom->View(2)*/ );
  //Sort by z position, assume traveling in beam direction.
  std::sort( ides.begin(), ides.end(), sort_IDEs );

  std::map< size_t, std::vector< const sim::IDE * > > results;
  
  size_t ide_start = 0;

  for( size_t i = 0; i < mctraj.size(); ++i ){
    results[i] = std::vector< const sim::IDE * >();  

    for( size_t j = ide_start; j < ides.size(); ++j ){

      if( ( ides[j]->z > mctraj.Z(i) ) && ( ides[j]->z < mctraj.Z(i+1) ) ){
        results[i].push_back( ides[j] );
      }

      else if( ( ides[j]->z < mctraj.Z(i) ) && ( ides[j]->z < mctraj.Z(i+1) ) ){
        continue; 
      }

      else{
        ide_start = j;
        break;
      }
    }
  }

  return results;
}

std::vector< const sim::IDE * > protoana::ProtoDUNETruthUtils::GetSimIDEsBetweenPoints	(	const simb::MCParticle &	mcpart,
		const TLorentzVector &	p1,
		const TLorentzVector &	p2
	)

Definition at line 1119 of file ProtoDUNETruthUtils.cxx.

{
  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  art::ServiceHandle<geo::Geometry> geom;
  std::vector<const sim::IDE *> ides =
      bt_serv->TrackIdToSimIDEs_Ps(mcpart.TrackId(), geom->View(2));

  std::vector<const sim::IDE *> results;
  for (const sim::IDE * theIDE : ides) {
    if (theIDE->z > p1.Z() && theIDE->z < p2.Z() /*&&
        theIDE->x > p1.X() && theIDE->x < p2.X() &&
        theIDE->y > p1.Y() && theIDE->y < p2.Y()*/) {
      results.push_back(theIDE);
    }
  }

  return results;
}

const simb::MCParticle * protoana::ProtoDUNETruthUtils::MatchPduneMCtoG4	(	const simb::MCParticle &	pDunePart,
		const art::Event &	evt
	)

Definition at line 813 of file ProtoDUNETruthUtils.cxx.

{  // Function that will match the protoDUNE MC particle to the Geant 4 MC particle, and return the matched particle (or a null pointer).

  // Find the energy of the procided MC Particle
  double pDuneEnergy = pDunePart.E();

  // Get list of the g4 particles. plist should be a std::map< int, simb::MCParticle* >
  art::ServiceHandle< cheat::ParticleInventoryService > pi_serv;
  const sim::ParticleList & plist = pi_serv->ParticleList();

  // Check if plist is empty
  if ( !plist.size() ) {
    std::cerr << "\n\n#####################################\n"
              << "\nEvent " << evt.id().event() << "\n"
              << "sim::ParticleList from cheat::ParticleInventoryService is empty\n"
              << "A null pointer will be returned\n"
              << "#####################################\n\n";
    return nullptr;
  }

  // Go through the list of G4 particles
  for ( auto partIt = plist.begin() ; partIt != plist.end() ; partIt++ ) {

    const simb::MCParticle* pPart = partIt->second;
    if (!pPart) {
      std::cerr << "\n\n#####################################\n"
                << "\nEvent " << evt.id().event() << "\n"
                << "GEANT particle #" << partIt->first << " returned a null pointer\n"
                << "This is not necessarily bad. It just means at least one\n"
                << "of the G4 particles returned a null pointer. It may well\n"
                << "have still matched a PD particle and a G4 particle.\n"
                << "#####################################\n\n";
      continue;
    }

    // If the initial energy of the g4 particle is very close to the energy of the protoDUNE particle, call it a day and have a cuppa.
    if ( (pDunePart.PdgCode() == pPart->PdgCode()) && fabs(pPart->E() - pDuneEnergy) < 0.00001 ) {
      return pPart;
    }

  }  // G4 particle list loop end.

  std::cout << "No G4 particle was matched for Event " << evt.id().event() << ". Null pointer returned\n";
  return nullptr;

}  // End MatchPduneMCtoG4

---

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

• protoduneana/protoduneana/Utilities/ProtoDUNETruthUtils.h
• protoduneana/protoduneana/Utilities/ProtoDUNETruthUtils.cxx