#include "FillTrue.h"
#include "larcorealg/GeoAlgo/GeoAlgo.h"
#include "RecoUtils/RecoUtils.h"
#include <functional>
#include <algorithm>

Namespaces
	caf
	Common Analysis Files.

Macros
#define	MATCH_PROCESS(name) if (process_name == #name) {return caf::kG4 ## name;}

#define	MATCH_PROCESS_NAMED(strname, id) if (process_name == #strname) {return caf::kG4 ## id;}

Functions
caf::SRTrackTruth	MatchTrack2Truth (const detinfo::DetectorClocksData &clockData, const std::vector< caf::SRTrueParticle > &particles, const std::vector< art::Ptr< recob::Hit >> &hits, const std::map< int, caf::HitsEnergy > &all_hits_map)

caf::SRTruthMatch	MatchSlice2Truth (const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< art::Ptr< simb::MCTruth >> &neutrinos, const caf::SRTruthBranch &srtruth, const cheat::ParticleInventoryService &inventory_service, const detinfo::DetectorClocksData &clockData)

float	ContainedLength (const TVector3 &v0, const TVector3 &v1, const std::vector< geoalgo::AABox > &boxes)

bool	FRFillNumuCC (const simb::MCTruth &mctruth, const std::vector< art::Ptr< sim::MCTrack >> &mctracks, const std::vector< geo::BoxBoundedGeo > &volumes, TRandom &rand, caf::SRFakeReco &fakereco)

bool	FRFillNueCC (const simb::MCTruth &mctruth, const std::vector< caf::SRTrueParticle > &srparticle, const std::vector< geo::BoxBoundedGeo > &volumes, TRandom &rand, caf::SRFakeReco &fakereco)

bool	isFromNuVertex (const simb::MCTruth &mc, const sim::MCTrack &track, float distance=5.0)

double	PDGMass (int pdg)

double	SmearLepton (const caf::SRTrueParticle &lepton, TRandom &rand)

double	SmearHadron (const caf::SRTrueParticle &hadron, TRandom &rand)

void	CopyTMatrixDToVector (const TMatrixD &m, std::vector< float > &v)

void	caf::FillSRGlobal (const sbn::evwgh::EventWeightParameterSet &pset, caf::SRGlobal &srglobal, std::map< std::string, unsigned int > &weightPSetIndex)

void	caf::FillTrackTruth (const std::vector< art::Ptr< recob::Hit >> &hits, const std::map< int, caf::HitsEnergy > &id_hits_map, const std::vector< caf::SRTrueParticle > &particles, const detinfo::DetectorClocksData &clockData, caf::SRTrack &srtrack, bool allowEmpty)

void	caf::FillShowerTruth (const std::vector< art::Ptr< recob::Hit >> &hits, const std::map< int, caf::HitsEnergy > &id_hits_map, const std::vector< caf::SRTrueParticle > &particles, const detinfo::DetectorClocksData &clockData, caf::SRShower &srshower, bool allowEmpty)

void	caf::FillStubTruth (const std::vector< art::Ptr< recob::Hit >> &hits, const std::map< int, caf::HitsEnergy > &id_hits_map, const std::vector< caf::SRTrueParticle > &particles, const detinfo::DetectorClocksData &clockData, caf::SRStub &srstub, bool allowEmpty)

void	caf::FillSliceTruth (const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< art::Ptr< simb::MCTruth >> &neutrinos, const caf::SRTruthBranch &srmc, const cheat::ParticleInventoryService &inventory_service, const detinfo::DetectorClocksData &clockData, caf::SRSlice &srslice, bool allowEmpty)

void	caf::FillMeVPrtlTruth (const evgen::ldm::MeVPrtlTruth &truth, const std::vector< geo::BoxBoundedGeo > &active_volumes, caf::SRMeVPrtl &srtruth)

void	caf::FillSliceFakeReco (const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< art::Ptr< simb::MCTruth >> &neutrinos, const caf::SRTruthBranch &srmc, const cheat::ParticleInventoryService &inventory_service, const detinfo::DetectorClocksData &clockData, caf::SRSlice &srslice, const std::vector< caf::SRTrueParticle > &srparticles, const std::vector< art::Ptr< sim::MCTrack >> &mctracks, const std::vector< geo::BoxBoundedGeo > &volumes, TRandom &rand)

void	caf::FillTrueNeutrino (const art::Ptr< simb::MCTruth > mctruth, const simb::MCFlux &mcflux, const simb::GTruth &gtruth, const std::vector< caf::SRTrueParticle > &srparticles, const std::map< int, std::vector< art::Ptr< recob::Hit >>> &id_to_truehit_map, caf::SRTrueInteraction &srneutrino, size_t i, const std::vector< geo::BoxBoundedGeo > &active_volumes)

void	caf::FillEventWeight (const sbn::evwgh::EventWeightMap &wgtmap, caf::SRTrueInteraction &srint, const std::map< std::string, unsigned int > &weightPSetIndex)

void	caf::FillTrueG4Particle (const simb::MCParticle &particle, const std::vector< geo::BoxBoundedGeo > &active_volumes, const std::vector< std::vector< geo::BoxBoundedGeo >> &tpc_volumes, const std::map< int, std::vector< std::pair< geo::WireID, const sim::IDE * >>> &id_to_ide_map, const std::map< int, std::vector< art::Ptr< recob::Hit >>> &id_to_truehit_map, const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, const std::vector< art::Ptr< simb::MCTruth >> &neutrinos, caf::SRTrueParticle &srparticle)

void	caf::FillFakeReco (const std::vector< art::Ptr< simb::MCTruth >> &mctruths, const std::vector< caf::SRTrueParticle > &srparticles, const std::vector< art::Ptr< sim::MCTrack >> &mctracks, const std::vector< geo::BoxBoundedGeo > &volumes, TRandom &rand, std::vector< caf::SRFakeReco > &srfakereco)

std::map< int, caf::HitsEnergy >	caf::SetupIDHitEnergyMap (const std::vector< art::Ptr< recob::Hit >> &allHits, const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker)

std::map< int, std::vector < art::Ptr< recob::Hit > > >	caf::PrepTrueHits (const std::vector< art::Ptr< recob::Hit >> &allHits, const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker)

std::map< int, std::vector < std::pair< geo::WireID, const sim::IDE * > > >	caf::PrepSimChannels (const std::vector< art::Ptr< sim::SimChannel >> &simchannels, const geo::GeometryCore &geo)

Macro Definition Documentation

#define MATCH_PROCESS ( name ) if (process_name == #name) {return caf::kG4 ## name;}

#define MATCH_PROCESS_NAMED	(	strname,
		id
	)	if (process_name == #strname) {return caf::kG4 ## id;}

Function Documentation

float ContainedLength	(	const TVector3 &	v0,
		const TVector3 &	v1,
		const std::vector< geoalgo::AABox > &	boxes
	)

Definition at line 1171 of file FillTrue.cxx.

                                                              {
   static const geoalgo::GeoAlgo algo;
 
   // if points are the same, return 0
   if ((v0 - v1).Mag() < 1e-6) return 0;
 
   // construct individual points
   geoalgo::Point_t p0(v0);
   geoalgo::Point_t p1(v1);
 
   // construct line segment
   geoalgo::LineSegment line(p0, p1);
 
   double length = 0;
 
   // total contained length is sum of lengths in all boxes
   // assuming they are non-overlapping
   for (auto const &box: boxes) {
     int n_contained = box.Contain(p0) + box.Contain(p1);
     // both points contained -- length is total length (also can break out of loop)
     if (n_contained == 2) {
       length = (v1 - v0).Mag();
       break;
     }
     // one contained -- have to find intersection point (which must exist)
     if (n_contained == 1) {
       auto intersections = algo.Intersection(line, box);
       // Because of floating point errors, it can sometimes happen
       // that there is 1 contained point but no "Intersections"
       // if one of the points is right on the edge
       if (intersections.size() == 0) {
         // determine which point is on the edge
         double tol = 1e-5;
         bool p0_edge = algo.SqDist(p0, box) < tol;
         bool p1_edge = algo.SqDist(p1, box) < tol;
         assert(p0_edge || p1_edge);
         // contained one is on edge -- can treat both as not contained
         //
         // In this case, no length
         if ((p0_edge && box.Contain(p0)) || (box.Contain(p1) && p1_edge))
           continue;
         // un-contaned one is on edge -- treat both as contained
         else if ((p0_edge && box.Contain(p1)) || (box.Contain(p0) && p1_edge)) {
     length = (v1 - v0).Mag();
     break;
         }
         else {
           assert(false); // bad
         }
       }
       // floating point errors can also falsely cause 2 intersection points
       //
       // in this case, one of the intersections must be very close to the
       // "contained" point, so the total contained length will be about
       // the same as the distance between the two intersection points
       else if (intersections.size() == 2) {
         length += (intersections.at(0).ToTLorentzVector().Vect() - intersections.at(1).ToTLorentzVector().Vect()).Mag();
         continue;
       }
       // "Correct"/ideal case -- 1 intersection point
       else if (intersections.size() == 1) {
         // get TVector at intersection point
         TVector3 int_tv(intersections.at(0).ToTLorentzVector().Vect());
         length += ( box.Contain(p0) ? (v0 - int_tv).Mag() : (v1 - int_tv).Mag() );
       }
       else assert(false); // bad
     }
     // none contained -- either must have zero or two intersections
     if (n_contained == 0) {
       auto intersections = algo.Intersection(line, box);
       if (!(intersections.size() == 0 || intersections.size() == 2)) {
         // more floating point error fixes...
         //
         // figure out which points are near the edge
         double tol = 1e-5;
         bool p0_edge = algo.SqDist(p0, box) < tol;
         bool p1_edge = algo.SqDist(p1, box) < tol;
         // and which points are near the intersection
         TVector3 vint = intersections.at(0).ToTLorentzVector().Vect();
 
         bool p0_int = (v0 - vint).Mag() < tol;
         bool p1_int = (v1 - vint).Mag() < tol;
         // exactly one of them should produce the intersection
         assert((p0_int && p0_edge) != (p1_int && p1_edge));
         // void variables when assert-ions are turned off
         (void) p0_int; (void) p1_int;
 
         // both close to edge -- full length is contained
         if (p0_edge && p1_edge) {
           length += (v0 - v1).Mag();
         }
         // otherwise -- one of them is not on an edge, no length is contained
         else {}
       }
       // assert(intersections.size() == 0 || intersections.size() == 2);
       else if (intersections.size() == 2) {
         TVector3 start(intersections.at(0).ToTLorentzVector().Vect());
         TVector3 end(intersections.at(1).ToTLorentzVector().Vect());
         length += (start - end).Mag();
       }
     }
   }
 
   return length;
 }//ContainedLength

void CopyTMatrixDToVector	(	const TMatrixD &	m,
		std::vector< float > &	v
	)

Definition at line 79 of file FillTrue.cxx.

 {
   const double& start = m(0, 0);
   v.insert(v.end(), &start, &start + m.GetNoElements());
 }

bool FRFillNueCC	(	const simb::MCTruth &	mctruth,
		const std::vector< caf::SRTrueParticle > &	srparticle,
		const std::vector< geo::BoxBoundedGeo > &	volumes,
		TRandom &	rand,
		caf::SRFakeReco &	fakereco
	)

Definition at line 873 of file FillTrue.cxx.

                                            {
   // Configuration -- TODO: make configurable?
 
   // Fiducial Volume
   float xmin = 10.;
   float xmax = 10.;
   float ymin = 10.;
   float ymax = 10.;
   float zmin = 10.;
   float zmax = 100.;
 
   // energy smearing
 /*  auto smear_lepton = [&rand](const caf::SRTrueParticle &lepton) -> float {
     const double smearing = 0.15;
     // Oscillation tech note says smear ionization deposition, but
     // code in old samples seems to use true energy.
     const double true_E = lepton.plane[0][2].visE + lepton.plane[1][2].visE;
     const double smeared_E = rand.Gaus(true_E, smearing * true_E / std::sqrt(true_E));
     return std::max(smeared_E, 0.0);
   };
   auto smear_hadron = [&rand](const caf::SRTrueParticle &hadron) -> float {
     const double smearing = 0.05;
     const double true_E = hadron.startE - PDGMass(hadron.pdg) / 1000;
     const double smeared_E = rand.Gaus(true_E, smearing * true_E);
     return std::max(smeared_E, 0.0);
   };
 */
   // visible energy threshold
   const float hadronic_energy_threshold = 0.021; // GeV
 
   fakereco.filled = false;
 
   // first check if neutrino exists
   if (!mctruth.NeutrinoSet()) return false;
 
   TVector3 nuVtx = mctruth.GetNeutrino().Nu().Position().Vect();
 
   // then check if fiducial
   int cryo_index = -1;
   for (unsigned i = 0; i < volumes.size(); i++) {
     const geo::BoxBoundedGeo &vol = volumes[i];
     geo::BoxBoundedGeo FV(vol.MinX() + xmin, vol.MaxX() - xmax, vol.MinY() + ymin, vol.MaxY() - ymax, vol.MinZ() + zmin, vol.MaxZ() - zmax);
     if (FV.ContainsPosition(nuVtx)) {
       cryo_index = i;
       break;
     }
   }
 
   if (cryo_index == -1) return false;
 
   std::vector<geoalgo::AABox> aa_volumes;
   const geo::BoxBoundedGeo &v = volumes.at(cryo_index);
   aa_volumes.emplace_back(v.MinX(), v.MinY(), v.MinZ(), v.MaxX(), v.MaxY(), v.MaxZ());
 
   //Showers arising from the vertex are identified and the reconstructed energy is found
   //from smearing the ionisation deposition. If more than one shower with energy above
   //100 MeV exists, the event is removed. This removes neutral pion events where the
   //pion decays into two photon showers.
 
   std::vector<const caf::SRTrueParticle*> lepton_candidates;
   for(const auto& particle: srparticles) {
     const int pdg = std::abs(particle.pdg);
     const float distance_from_vertex = std::hypot(nuVtx.X() - particle.start.x,
                                                   nuVtx.Y() - particle.start.y,
                                                   nuVtx.Z() - particle.start.z);
     if((pdg == 11 || pdg == 22) && distance_from_vertex < 5) {
       const auto parent = std::find_if(srparticles.begin(), srparticles.end(),
         [&particle](const caf::SRTrueParticle& parent_candidate) -> bool {
           return particle.parent == std::abs(parent_candidate.G4ID);
         });
       if((parent == srparticles.end()
            || !(std::abs((*parent).pdg) == 11 || std::abs((*parent).pdg) == 22))
          && SmearLepton(particle, rand) > 0.1) {
         lepton_candidates.push_back(&particle);
       }
     }
   }
   if(lepton_candidates.size() != 1) return false;
   
   const caf::SRTrueParticle* lepton = lepton_candidates[0];
   caf::SRFakeRecoParticle fake_lepton;
   fake_lepton.pid = 11;
   fake_lepton.ke = SmearLepton(*lepton, rand) - PDGMass(11) / 1000;
   // Should these be set for a shower?
   // fake_lepton.len = ???;
   // fake_lepton.costh = ???;
   // fake_lepton.contained = false;
 
   // Get hadrons
 
   std::vector<caf::SRFakeRecoParticle> fake_hadrons;
   float hadronic_E = 0.0;
   for(const auto& particle: srparticles) {
     const int pdg = std::abs(particle.pdg);
     const float ke = SmearHadron(particle, rand);
     const float distance_from_vertex = std::hypot(nuVtx.X() - particle.start.x,
                                                   nuVtx.Y() - particle.start.y,
                                                   nuVtx.Z() - particle.start.z);
 
     if((pdg == 2212 || pdg == 211 || pdg == 321) && distance_from_vertex < 5
        && particle.start_process == caf::kG4primary && ke > hadronic_energy_threshold) {
       caf::SRFakeRecoParticle fake_hadron;
       fake_hadron.pid = pdg;
       fake_hadron.len = ContainedLength(TVector3(particle.start), TVector3(particle.end), aa_volumes);
       fake_hadron.contained = false;
       for(const geo::BoxBoundedGeo &vol: volumes) {
         if(vol.ContainsPosition(TVector3(particle.start)) 
            && vol.ContainsPosition(TVector3(particle.end))) {
           fake_hadron.contained = true;
         }
       }
       fake_hadron.costh = TVector3(particle.start).CosTheta();
       fake_hadron.ke = ke;
       hadronic_E += ke;
       fake_hadrons.push_back(fake_hadron);
     }
   }
 
   // If there is only one photon candidate in the event, a conversion gap cut is applied.
   // If the vertex is deemed visible, due to the presence of at least 50 MeV of hadronic
   // kinetic energy, and the photon starts to shower further than 3 cm from the vertex the
   // event is removed.
   
   if(lepton->pdg == 22 && hadronic_E > 0.05 && (TVector3(lepton->end) - nuVtx).Mag() > 3)
     return false;
 
   // Remaining photons undergo a dE/dx cut resulting in a 94% background rejection.
 
   double weight = 1.0;
   if(lepton->pdg == 22) weight *= 0.06;
 
   // If a misidentified photon originates from a resonant numuCC interaction, the muon
   // lepton is identified. Events where a muon travels greater than 1 m are assumed to be
   // from numuCC interactions and the event is removed.
 
   if(std::abs(mctruth.GetNeutrino().Nu().PdgCode()) == 14 && mctruth.GetNeutrino().CCNC() == 0) {
     const caf::SRTrueParticle* muon = NULL;
     for(const auto& particle: srparticles) {
       const int pdg = std::abs(particle.pdg);
       const float distance_from_vertex = std::hypot(nuVtx.X() - particle.start.x,
                                                     nuVtx.Y() - particle.start.y,
                                                     nuVtx.Z() - particle.start.z);
       if(pdg == 13 && distance_from_vertex < 5 && particle.start_process == caf::kG4primary
          && (muon == NULL || muon->startE < particle.startE))
         muon = &particle;
     }
     if(muon != NULL && ContainedLength(TVector3(muon->start), TVector3(muon->end), aa_volumes) > 100)
       return false;
   }
 
   // Events where the shower has an energy less than 200 MeV are removed
 
   if((lepton->plane[0][2].visE + lepton->plane[1][2].visE) < 0.2) return false;
 
   // total up the energy to get the reco neutrino energy
   fakereco.nuE = fake_lepton.ke + PDGMass(11) / 1000 + hadronic_E;
 
   // save the particles
   fakereco.lepton = fake_lepton;
   fakereco.hadrons = fake_hadrons;
 
   // other info
   fakereco.vtx.x = nuVtx.X();
   fakereco.vtx.y = nuVtx.Y();
   fakereco.vtx.z = nuVtx.Z();
 
   fakereco.wgt = 0.8;
   fakereco.wgt *= weight;
 
   fakereco.nhad = fakereco.hadrons.size();
   fakereco.filled = true;
 
   return true;
 }

bool FRFillNumuCC	(	const simb::MCTruth &	mctruth,
		const std::vector< art::Ptr< sim::MCTrack >> &	mctracks,
		const std::vector< geo::BoxBoundedGeo > &	volumes,
		TRandom &	rand,
		caf::SRFakeReco &	fakereco
	)

Definition at line 710 of file FillTrue.cxx.

                                            {
   // Configuration -- TODO: make configurable?
 
   // Fiducial Volume
   float xmin = 10.;
   float xmax = 10.;
   float ymin = 10.;
   float ymax = 10.;
   float zmin = 10.;
   float zmax = 100.;
 
   // energy smearing
   float hadron_smearing = 0.05;
   float lepton_contained_smearing = 0.02;
   std::function<float (float)> lepton_exiting_smearing = [](float length) {
     float A = 0.102;
     float B = 0.000612;
     return -A * TMath::Log(B * length);
   };
 
   // visible energy threshold
   float hadronic_energy_threshold = 0.021; // GeV
 
   // length cut
   float contained_length_cut = 50.; // cm
   float exiting_length_cut = 100.; // cm
 
   fakereco.filled = false;
 
   // first check if neutrino exists
   if (!mctruth.NeutrinoSet()) return false;
 
   // then check if fiducial
   int cryo_index = -1;
   for (unsigned i = 0; i < volumes.size(); i++) {
     const geo::BoxBoundedGeo &vol = volumes[i];
     geo::BoxBoundedGeo FV(vol.MinX() + xmin, vol.MaxX() - xmax, vol.MinY() + ymin, vol.MaxY() - ymax, vol.MinZ() + zmin, vol.MaxZ() - zmax);
     if (FV.ContainsPosition(mctruth.GetNeutrino().Nu().Position().Vect())) {
       cryo_index = i;
       break;
     }
   }
 
   if (cryo_index == -1) return false;
 
   std::vector<geoalgo::AABox> aa_volumes;
   const geo::BoxBoundedGeo &v = volumes.at(cryo_index);
   aa_volumes.emplace_back(v.MinX(), v.MinY(), v.MinZ(), v.MaxX(), v.MaxY(), v.MaxZ());
 
   // look for CC lepton or a \pi^+/- which can "fake" a numu CC interaction
 
   int lepton_ind = -1;
   // CC lepton
   if (abs(mctruth.GetNeutrino().Nu().PdgCode()) == 14 && mctruth.GetNeutrino().CCNC() == 0) {
     for (int i = 0; i < (int)mctracks.size(); i++) {
       if (isFromNuVertex(mctruth, *mctracks[i]) && abs(mctracks[i]->PdgCode()) == 13 && mctracks[i]->Process() == "primary") {
         if (lepton_ind == -1 || mctracks[lepton_ind]->Start().E() < mctracks[i]->Start().E()) {
           lepton_ind = i;
         }
       }
     }
   }
   // NC pion
   else if (mctruth.GetNeutrino().CCNC() == 1) {
     for (int i = 0; i < (int)mctracks.size(); i++) {
       if (isFromNuVertex(mctruth, *mctracks[i]) && abs(mctracks[i]->PdgCode()) == 211 && mctracks[i]->Process() == "primary") {
         if (lepton_ind == -1 || mctracks[lepton_ind]->Start().E() < mctracks[i]->Start().E()) {
           lepton_ind = i;
         }
       }
     }
   }
 
   // no matching track -- no fake reco
   if (lepton_ind == -1) return false;
 
   // Now set the "lepton"
   caf::SRFakeRecoParticle fake_lepton;
 
   fake_lepton.pid = 13;
   fake_lepton.contained = false;
   for (const geo::BoxBoundedGeo &vol: volumes) {
     if (vol.ContainsPosition(mctracks[lepton_ind]->Start().Position().Vect()) && vol.ContainsPosition(mctracks[lepton_ind]->Start().Position().Vect())) {
       fake_lepton.contained = true;
     }
   }
   fake_lepton.len = ContainedLength(mctracks[lepton_ind]->Start().Position().Vect(), mctracks[lepton_ind]->End().Position().Vect(), aa_volumes);
   fake_lepton.costh = mctracks[lepton_ind]->Start().Position().Vect().CosTheta();
 
   // apply length cut
   if (fake_lepton.contained && fake_lepton.len < contained_length_cut) return false;
   if (!fake_lepton.contained && fake_lepton.len < exiting_length_cut) return false;
 
 
   // smear the lepton energy
   float smearing = fake_lepton.contained ? lepton_contained_smearing : lepton_exiting_smearing(fake_lepton.len);
   float ke = (mctracks[lepton_ind]->Start().E() - PDGMass(mctracks[lepton_ind]->PdgCode())) / 1000. /* MeV -> GeV*/;
   fake_lepton.ke = rand.Gaus(ke, smearing * ke);
   fake_lepton.ke = std::max(fake_lepton.ke, 0.f);
 
   // get the hadronic state
   std::vector<caf::SRFakeRecoParticle> hadrons;
 
   for (int i = 0; i < (int)mctracks.size(); i++) {
     if (isFromNuVertex(mctruth, *mctracks[i]) // from this interaction
      && (abs(mctracks[i]->PdgCode()) == 211 || abs(mctracks[i]->PdgCode()) == 321 || abs(mctracks[i]->PdgCode()) == 2212) // hadronic
      && mctracks[i]->Process() == "primary" // primary
      && i != lepton_ind // not the fake lepton
     ) {
       caf::SRFakeRecoParticle hadron;
       hadron.pid = abs(mctracks[i]->PdgCode());
       hadron.contained = false;
       for (const geo::BoxBoundedGeo &vol: volumes) {
         if (vol.ContainsPosition(mctracks[i]->Start().Position().Vect()) && vol.ContainsPosition(mctracks[i]->Start().Position().Vect())) {
           hadron.contained = true;
         }
       }
       hadron.len = ContainedLength(mctracks[i]->Start().Position().Vect(), mctracks[i]->Start().Position().Vect(), aa_volumes);
       hadron.costh = mctracks[i]->Start().Position().Vect().CosTheta();
 
       float ke = (mctracks[i]->Start().E() - PDGMass(mctracks[i]->PdgCode())) / 1000. /* MeV -> GeV*/;
       if (ke < hadronic_energy_threshold) continue;
 
       hadron.ke = rand.Gaus(ke, hadron_smearing * ke);
       hadron.ke = std::max(hadron.ke, 0.f);
 
       hadrons.push_back(hadron);
     }
   }
 
   // total up the energy to get the reco neutrino energy
   fakereco.nuE = fake_lepton.ke;
   for (const caf::SRFakeRecoParticle &had: hadrons) fakereco.nuE += had.ke;
   fakereco.nuE += PDGMass(13) / 1000.; // MeV -.> GeV
 
   // save the particles
   fakereco.lepton = fake_lepton;
   fakereco.hadrons = hadrons;
 
   // other info
   TVector3 vertex = mctruth.GetNeutrino().Nu().Position().Vect();
   fakereco.vtx.x = vertex.X();
   fakereco.vtx.y = vertex.Y();
   fakereco.vtx.z = vertex.Z();
 
   // signal
   if (abs(mctruth.GetNeutrino().Nu().PdgCode()) == 14 && mctruth.GetNeutrino().CCNC() == 0) {
     fakereco.wgt = 0.8;
   }
   // bkg
   else {
     fakereco.wgt = 1.;
   }
   fakereco.nhad = fakereco.hadrons.size();
   fakereco.filled = true;
 
   return true;
 }

bool isFromNuVertex	(	const simb::MCTruth &	mc,
		const sim::MCTrack &	track,
		float	distance = `5.0`
	)

Definition at line 36 of file FillTrue.cxx.

                                         {
   TVector3 nuVtx = mc.GetNeutrino().Nu().Trajectory().Position(0).Vect();
   TVector3 trkStart = track.Start().Position().Vect();
   return (trkStart - nuVtx).Mag() < distance;
 }

caf::SRTruthMatch MatchSlice2Truth	(	const std::vector< art::Ptr< recob::Hit >> &	hits,
		const std::vector< art::Ptr< simb::MCTruth >> &	neutrinos,
		const caf::SRTruthBranch &	srtruth,
		const cheat::ParticleInventoryService &	inventory_service,
		const detinfo::DetectorClocksData &	clockData
	)

Definition at line 1367 of file FillTrue.cxx.

                                                                                {
   caf::SRTruthMatch ret;
   float total_energy = CAFRecoUtils::TotalHitEnergy(clockData, hits);
   // speed optimization: if there are no truths, all the matching energy must be cosmic
   if (truths.empty()) {
     ret.visEinslc = total_energy / 1000. /* MeV -> GeV */;
     ret.visEcosmic = total_energy / 1000. /* MeV -> GeV */;
     ret.eff_cryo = -1.;
     ret.eff = -1;
     ret.pur = -1;
     ret.index = -1;
     return ret;
   }
   std::vector<std::pair<int, float>> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clockData, hits, true);
   std::vector<float> matching_energy(truths.size(), 0.);
   for (auto const &pair: matches) {
     art::Ptr<simb::MCTruth> truth;
     try {
       truth = inventory_service.TrackIdToMCTruth_P(pair.first);
     }
     // Ignore track ID's that cannot be looked up
     catch(...) {
       continue;
     }
     for (unsigned ind = 0; ind < truths.size(); ind++) {
       if (truth == truths[ind]) {
         matching_energy[ind] += pair.second;
         break;
       }
     }
   }
 
   float cosmic_energy = total_energy;
   for (float E: matching_energy) cosmic_energy -= E;
 
   float matching_frac = *std::max_element(matching_energy.begin(), matching_energy.end()) / total_energy;
   int index = (matching_frac > 0.5) ? std::distance(matching_energy.begin(), std::max_element(matching_energy.begin(), matching_energy.end())) : -1;
 
   ret.index = index;
   if (index >= 0) {
     ret.visEinslc = total_energy / 1000. /* MeV -> GeV */;
     ret.visEcosmic = cosmic_energy / 1000. /* MeV -> GeV */;
     ret.pur = matching_energy[index] / total_energy;
 
     double totVisE = 0.;
     for (int p = 0; p < 3; p++) {
       for (int c = 0; c < 2; c++) {
         totVisE += srmc.nu[index].plane[c][p].visE;
       }
     }
 
     ret.eff = (matching_energy[index] / 1000.) / totVisE;
 
     // Lookup the cryostat
     int icryo = -1;
     if (!hits.empty()) {
       icryo = hits[0]->WireID().Cryostat;
     }
 
     if (icryo >= 0) {
       ret.eff_cryo = (matching_energy[index] / 1000.) / 
           (srmc.nu[index].plane[icryo][0].visE + srmc.nu[index].plane[icryo][1].visE + srmc.nu[index].plane[icryo][2].visE);
     }
     else {
       ret.eff_cryo = -1;
     }
 
   }
   else {
     ret.pur = -1;
     ret.eff = -1;
     ret.eff_cryo = -1;
   }
   return ret;
 }//Slc2Truth

caf::SRTrackTruth MatchTrack2Truth	(	const detinfo::DetectorClocksData &	clockData,
		const std::vector< caf::SRTrueParticle > &	particles,
		const std::vector< art::Ptr< recob::Hit >> &	hits,
		const std::map< int, caf::HitsEnergy > &	all_hits_map
	)

Definition at line 1279 of file FillTrue.cxx.

                                                                                    {
 
   art::ServiceHandle<cheat::BackTrackerService> bt_serv;
 
   // this id is the same as the mcparticle ID as long as we got it from geant4
   std::vector<std::pair<int, float>> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clockData, hits, true);
   float total_energy = CAFRecoUtils::TotalHitEnergy(clockData, hits);
   std::map<int, caf::HitsEnergy> track_hits_map = caf::SetupIDHitEnergyMap(hits, clockData, *bt_serv.get());
 
   caf::SRTrackTruth ret;
 
   ret.visEintrk = total_energy / 1000. /* MeV -> GeV */;
 
   // setup the matches
   for (auto const &pair: matches) {
     caf::SRParticleMatch match;
     match.G4ID = pair.first;
     match.energy = pair.second / 1000. /* MeV -> GeV */;
 
     caf::HitsEnergy track_matched_hits = track_hits_map.find(match.G4ID)->second;
     caf::HitsEnergy all_matched_hits = all_hits_map.find(match.G4ID)->second;
 
     match.hit_purity = (hits.size() != 0) ? track_matched_hits.nHits / (float) hits.size() : 0.;
     match.energy_purity = (ret.visEintrk > 0) ? match.energy / ret.visEintrk : 0.;
     match.hit_completeness = (all_matched_hits.nHits != 0) ? track_matched_hits.nHits / (float) all_matched_hits.nHits : 0.;
     match.energy_completeness = (all_matched_hits.totE > 0) ? pair.second / all_matched_hits.totE : 0.;
     
     ret.matches.push_back(match);
   }
 
   // sort highest energy match to lowest
   std::sort(ret.matches.begin(), ret.matches.end(),
       [](const caf::SRParticleMatch &a, const caf::SRParticleMatch &b) {
         return a.energy > b.energy;
       }
       );
 
   bool found_bestmatch = false;
   if (ret.matches.size()) {
     ret.bestmatch = ret.matches.at(0);
     for (unsigned i_part = 0; i_part < particles.size(); i_part++) {
       if (particles[i_part].G4ID == ret.bestmatch.G4ID) {
         ret.p = particles[i_part];
         found_bestmatch = true;
         break;
       }
     }
   }
 
   // Calculate efficiency / purity
   if (found_bestmatch) {
     double match_total_energy = 0.;
     for (int p = 0; p < 3; p++) {
       for (int c = 0; c < 2; c++) {
         match_total_energy += ret.p.plane[c][p].visE;
       }
     }
 
     ret.eff = ret.matches[0].energy / match_total_energy; 
     ret.pur = ret.matches[0].energy / ret.visEintrk;
 
     int icryo = -1;
     if (!hits.empty()) {
       icryo = hits[0]->WireID().Cryostat;
     }
 
     assert(icryo < 2);
     if (icryo >= 0 && icryo < 2) {
       float match_cryo_energy = ret.p.plane[icryo][0].visE + ret.p.plane[icryo][1].visE + ret.p.plane[icryo][2].visE;
       ret.eff_cryo = ret.matches[0].energy / match_cryo_energy;
     }
     else {
       ret.eff_cryo = -1;
     }
 
   }
   else {
     ret.eff = -1.;
     ret.pur = -1.;
     ret.eff_cryo = -1.;
   }
 
   ret.nmatches = ret.matches.size();
 
   return ret;
 }//MatchTrack2Truth

double PDGMass ( int pdg )

Definition at line 44 of file FillTrue.cxx.

                         {
   const TDatabasePDG *PDGTable = TDatabasePDG::Instance();
   // regular particle
   if (pdg < 1000000000) {
     TParticlePDG* ple = PDGTable->GetParticle(pdg);
     if (ple == NULL) return -1;
     return ple->Mass() * 1000.0;
   }
   // ion
   else {
     int p = (pdg % 10000000) / 10000;
     int n = (pdg % 10000) / 10 - p;
     return (PDGTable->GetParticle(2212)->Mass() * p +
             PDGTable->GetParticle(2112)->Mass() * n) * 1000.0;
   }
 }

double SmearHadron	(	const caf::SRTrueParticle &	hadron,
		TRandom &	rand
	)

Definition at line 70 of file FillTrue.cxx.

                                                                    {
   const double smearing = 0.05;
   const double true_E = hadron.startE - PDGMass(hadron.pdg) / 1000;
   const double smeared_E = rand.Gaus(true_E, smearing * true_E);
   return std::max(smeared_E, 0.0);
 }

double SmearLepton	(	const caf::SRTrueParticle &	lepton,
		TRandom &	rand
	)

Definition at line 61 of file FillTrue.cxx.

                                                                    {
   const double smearing = 0.15;
   // Oscillation tech note says smear ionization deposition, but
   // code in old samples seems to use true energy.
   const double true_E = lepton.plane[0][2].visE + lepton.plane[1][2].visE;
   const double smeared_E = rand.Gaus(true_E, smearing * true_E / std::sqrt(true_E));
   return std::max(smeared_E, 0.0);
 }

Namespaces

Macros

Functions

Macro Definition Documentation

Function Documentation