Inheritance diagram for numu::NuMuEfficiencyStudy:

Public Member Functions
	NuMuEfficiencyStudy (fhicl::ParameterSet const &p)

	NuMuEfficiencyStudy (NuMuEfficiencyStudy const &)=delete

	NuMuEfficiencyStudy (NuMuEfficiencyStudy &&)=delete

NuMuEfficiencyStudy &	operator= (NuMuEfficiencyStudy const &)=delete

NuMuEfficiencyStudy &	operator= (NuMuEfficiencyStudy &&)=delete

void	analyze (art::Event const &e) override

Private Member Functions
bool	InFV (TVector3 pos)

void	InitHistos ()

void	InitHisto (Histos &h, const std::string &prefix)

void	FillTrue (Histos &h, const simb::MCTruth &interaction, const simb::MCParticle &G4lepton)

void	FillMatchingMu (Histos &h, float purity, float completeness)

void	FillMatchingSlc (Histos &h, float purity, float completeness)

float	ParticleLength (const simb::MCParticle &particle)

Private Attributes
Histos	fTruth

std::vector< Histos >	fReco

std::vector< std::vector < geo::BoxBoundedGeo > >	fTPCVolumes

std::vector< geo::BoxBoundedGeo >	fActiveVolumes

std::array< float, 6 >	fFiducialInset

std::vector< geo::BoxBoundedGeo >	fFiducialVolumes

Static Private Attributes
static const unsigned	n_reco_eff = 14

Detailed Description

Definition at line 86 of file NuMuEfficiencyStudy_module.cc.

Constructor & Destructor Documentation

numu::NuMuEfficiencyStudy::NuMuEfficiencyStudy ( fhicl::ParameterSet const & p )

explicit

Definition at line 188 of file NuMuEfficiencyStudy_module.cc.

189 : EDAnalyzer{p},

190 fFiducialInset({p.get<float>("xmin"), p.get<float>("xmax"), p.get<float>("ymin"), p.get<float>("ymax"), p.get<float>("zmin"), p.get<float>("zmax")})

p

pdgs p

Definition: selectors.fcl:22

numu::NuMuEfficiencyStudy::fFiducialInset

std::array< float, 6 > fFiducialInset

Definition: NuMuEfficiencyStudy_module.cc:110

numu::NuMuEfficiencyStudy::NuMuEfficiencyStudy ( NuMuEfficiencyStudy const & )

delete

numu::NuMuEfficiencyStudy::NuMuEfficiencyStudy ( NuMuEfficiencyStudy && )

delete

Member Function Documentation

void numu::NuMuEfficiencyStudy::analyze ( art::Event const & e )

override

Definition at line 265 of file NuMuEfficiencyStudy_module.cc.

 {
   art::ServiceHandle<cheat::BackTrackerService> backtracker;
   art::ServiceHandle<cheat::ParticleInventoryService> inventory_service;
 
   auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(ev);
 
   art::ValidHandle<std::vector<simb::MCTruth>> mctruth_handle = ev.getValidHandle<std::vector<simb::MCTruth>>("generator");
   const std::vector<simb::MCTruth> &mctruth_list = *mctruth_handle;
   std::vector<art::Ptr<simb::MCTruth>> mctruth_ptrs;
   art::fill_ptr_vector(mctruth_ptrs, mctruth_handle);
 
   art::FindManyP<simb::MCParticle, sim::GeneratedParticleInfo> truth_to_particles(mctruth_handle, ev, "largeant");
   const std::vector<simb::MCParticle> &g4_particles = *ev.getValidHandle<std::vector<simb::MCParticle>>("largeant");
   std::map<int, const simb::MCParticle*> g4_particle_map;
   for (unsigned i_g4 = 0; i_g4 < g4_particles.size(); i_g4++) {
     g4_particle_map[g4_particles[i_g4].TrackId()] =  &g4_particles[i_g4];
   }
 
   std::vector<int> truth_to_numucc;
   std::vector<simb::MCTruth> numucc;
   std::vector<simb::MCParticle> G4muons;
   // iterate over the true interactions
   for (unsigned i = 0; i < mctruth_list.size(); i++) {
     const simb::MCTruth &truth = mctruth_list[i];
     TVector3 vertex(truth.GetNeutrino().Nu().Position().X(), truth.GetNeutrino().Nu().Position().Y(), truth.GetNeutrino().Nu().Position().Z());
     // fiducial numu CC interactions
     if (InFV(vertex) && truth.GetNeutrino().CCNC() == 0 && abs(truth.GetNeutrino().Nu().PdgCode()) == 14) {
       // find the muon in G4
       const simb::MCParticle *G4muon = Genie2G4MCParticle(truth.GetNeutrino().Lepton(), truth, truth_to_particles.at(i), truth_to_particles.data(i));
       FillTrue(fTruth, truth, *G4muon);
       numucc.push_back(truth);
       G4muons.push_back(*G4muon);
       truth_to_numucc.push_back(numucc.size()-1);
     }
     else truth_to_numucc.push_back(-1);
   }
 
   // no true numu CC? Return
   if (numucc.size() == 0) return;
 
   // get the deposited energy for each muon
   std::vector<float> muonVisE(G4muons.size(), 0.);
   for (unsigned i = 0; i < G4muons.size(); i++) {
     std::vector<const sim::IDE*> particle_ides(backtracker->TrackIdToSimIDEs_Ps(G4muons[i].TrackId()));
     for (const sim::IDE *ide: particle_ides) muonVisE[i] += ide->energy / 1000.;
   }
 
   // get the deposited energy from primary particles in the neutrino
   std::vector<float> numuVisE(numucc.size(), 0.);
   for (unsigned i_g4 = 0; i_g4 < g4_particles.size(); i_g4++) {
     const simb::MCParticle &particle = g4_particles[i_g4];
     if (particle.Process() != "primary") continue;
     art::Ptr<simb::MCTruth> truth = inventory_service->TrackIdToMCTruth_P(particle.TrackId());
     for (unsigned i_truth = 0; i_truth < mctruth_ptrs.size(); i_truth++) {
       if (truth == mctruth_ptrs[i_truth]) {
         if (truth_to_numucc[i_truth] != -1) {
           std::vector<const sim::IDE*> particle_ides(backtracker->TrackIdToSimIDEs_Ps(particle.TrackId()));
           for (const sim::IDE *ide: particle_ides) numuVisE[truth_to_numucc[i_truth]] += ide->energy / 1000.;
         }
         break;
       }
     }
   } 
 
   // now collect the reco information 
   //
   // does a muon track exist?
   std::vector<int> matching_track(G4muons.size(), -1);
   std::vector<float> track_purity(G4muons.size(), -1);
   std::vector<float> track_completeness(G4muons.size(), -1);
 
   art::ValidHandle<std::vector<recob::PFParticle>> particle_handle = ev.getValidHandle<std::vector<recob::PFParticle>>("pandora");
   art::FindManyP<recob::Track> tracks(particle_handle, ev, "pandoraTrack");
   art::FindManyP<larpandoraobj::PFParticleMetadata> metadatas(particle_handle, ev, "pandora");
   art::FindManyP<recob::Vertex> particles_to_vertex(particle_handle, ev, "pandora");
   for (unsigned i_part = 0; i_part < particle_handle->size(); i_part++) {
     if (!tracks.at(i_part).size()) continue;
     art::FindManyP<recob::Hit> fmHits(tracks.at(i_part), ev, "pandoraTrack");
     const std::vector<art::Ptr<recob::Hit>> &this_hits = fmHits.at(0);
 
     std::vector<std::pair<int, float>> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clock_data, this_hits);
 
     // sort by matching energy
     std::sort(matches.begin(), matches.end(), 
       [] (auto const &lhs, auto const &rhs) { return lhs.second > rhs.second; });
 
     for (unsigned i_cc = 0; i_cc < G4muons.size(); i_cc++) {
       if (matches.size() && matches[0].first == G4muons[i_cc].TrackId() && matches[0].second / muonVisE[i_cc] > 0.5) {
         matching_track[i_cc] = particle_handle->at(i_part).Self();
         track_purity[i_cc] = matches[0].second / 
           std::accumulate(matches.begin(), matches.end(), 0.f,
             [](auto &a, auto &b){return a + b.second;});
         track_completeness[i_cc] = matches[0].second / muonVisE[i_cc];
       }
     }
   }
 
   // also find a matching slice
   std::vector<int> matching_slice(numucc.size(), -1);
   std::vector<float> matching_purity(numucc.size(), -1);
   std::vector<float> matching_completeness(numucc.size(), -1);
 
   art::ValidHandle<std::vector<recob::Slice>> slice_handle = ev.getValidHandle<std::vector<recob::Slice>>("pandora"); 
   art::FindManyP<recob::Hit> slice_hits(slice_handle, ev, "pandora");
   art::FindManyP<recob::PFParticle> slice_particles(slice_handle, ev, "pandora");
 
   for (unsigned i_slice = 0; i_slice < slice_handle->size(); i_slice++) {
     std::vector<float> matchingE(numuVisE.size(), 0.);
     std::vector<float> matchingPrimaryE(numuVisE.size(), 0.);
     float totalE = 0.;
 
     std::vector<std::pair<int, float>> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clock_data, slice_hits.at(i_slice));
     for (unsigned i_match = 0; i_match < matches.size(); i_match++) {
       totalE += matches[i_match].second / 1000.;
     }
 
     for (unsigned i_match = 0; i_match < matches.size(); i_match++) {
       if (g4_particle_map.count(matches[i_match].first)) {
         const simb::MCParticle &part = *g4_particle_map.at(matches[i_match].first);
         art::Ptr<simb::MCTruth> truth = inventory_service->TrackIdToMCTruth_P(part.TrackId());
         for (unsigned i_truth = 0; i_truth < mctruth_ptrs.size(); i_truth++) {
           if (truth == mctruth_ptrs[i_truth]) {
             if (truth_to_numucc[i_truth] != -1) {
               if (part.Process() == "primary") matchingPrimaryE[truth_to_numucc[i_truth]] += matches[i_match].second / 1000.;
               matchingE[truth_to_numucc[i_truth]] += matches[i_match].second / 1000.;
             }
             break;
           }
         }
       }
     }
 
     unsigned best_match = std::distance(matchingE.begin(), std::max_element(matchingE.begin(), matchingE.end()));
     if (matchingPrimaryE[best_match] / numuVisE[best_match] > 0.5) {
       matching_slice[best_match] = i_slice;
       matching_purity[best_match] = matchingE[best_match] / totalE; 
       matching_completeness[best_match] = matchingPrimaryE[best_match] / numuVisE[best_match];
     }
   }
 
   // gather information for efficiency calculation
   for (unsigned i = 0; i < numucc.size(); i++) {
     bool has_muon_track = matching_track[i] != -1;
     bool has_slice = matching_slice[i] != -1;
 
     bool has_pure_slice = has_slice && matching_purity[i] > 0.5;
 
     bool muon_is_reco = false;
     bool muon_is_fid = false;
     if (has_muon_track) {
       const larpandoraobj::PFParticleMetadata &meta = *metadatas.at(matching_track[i]).at(0);
       auto const &properties = meta.GetPropertiesMap();
       std::cout << "Muon Track Properties:\n";
       for (auto pair: properties) std::cout << pair.first << std::endl;
       if (properties.count("IsClearCosmic")) {
         muon_is_reco = false;
       }
       else {
         muon_is_reco = true;
       }
 
       auto muon_vert = tracks.at(matching_track[i]).at(0)->Start();
       TVector3 muon_vert_v(muon_vert.X(), muon_vert.Y(), muon_vert.Z());
             
       muon_is_fid = InFV(muon_vert_v);
     }
 
     bool slice_is_neutrino = false;
     bool slice_is_reco = false;
     bool slice_is_fiducial = false;
     bool slice_has_muon = false;
     if (has_slice) {
       const std::vector<art::Ptr<recob::PFParticle>> &this_slice_particles = slice_particles.at(matching_slice[i]); 
       for (art::Ptr<recob::PFParticle> pfp: this_slice_particles) {
         if (pfp->IsPrimary()) {
           const larpandoraobj::PFParticleMetadata &meta = *metadatas.at(pfp.key()).at(0);
           auto const &properties = meta.GetPropertiesMap();
           std::cout << "Particle Properties:\n";
           for (auto pair: properties) std::cout << pair.first << " " << pair.second << std::endl;
 
           if (properties.count("IsNeutrino")) slice_is_neutrino = true;
           if (properties.count("IsNeutrino") && pfp->PdgCode() == 14) slice_is_reco = true;
 
           // check for vertex
           if (particles_to_vertex.at(pfp.key()).size()) {
             const recob::Vertex &vert = *particles_to_vertex.at(pfp.key()).at(0);
             TVector3 vect(vert.position().X(), vert.position().Y(), vert.position().Z());
             slice_is_fiducial = InFV(vect);
           }
         }
         if (has_muon_track && matching_track[i] == (int)pfp->Self()) {
           slice_has_muon = true;
         }
       }
     }
 
     bool has_pure_nu_slice = has_pure_slice && slice_is_neutrino;
     bool has_fid_nu_slice = slice_is_fiducial && slice_is_neutrino;
     bool has_fid_pure_slice  = slice_is_fiducial  && has_pure_slice;
     bool has_reco_fid_slice = slice_is_fiducial && slice_is_reco;
     bool slice_has_reco_muon = slice_has_muon && muon_is_reco;
 
     std::array<bool, n_reco_eff> pass = {has_muon_track, muon_is_reco, muon_is_fid, has_slice, has_pure_slice, slice_is_reco, slice_is_fiducial, slice_has_muon, has_fid_pure_slice, has_reco_fid_slice, slice_has_reco_muon, has_fid_nu_slice, has_pure_nu_slice, slice_is_neutrino};
     for (unsigned j = 0; j < n_reco_eff; j++) {
       for (unsigned k = 0; k < n_reco_eff; k++) {
         if (pass[j] && pass[k]) {
           FillTrue(fReco[j * n_reco_eff + k], numucc[i], G4muons[i]);
           FillMatchingMu(fReco[j * n_reco_eff + k], track_purity[i], track_completeness[i]);
           FillMatchingSlc(fReco[j * n_reco_eff + k], matching_purity[i], matching_completeness[i]);
         }
       }
     }
   }
   
 }

void numu::NuMuEfficiencyStudy::FillMatchingMu	(	Histos &	h,
		float	purity,
		float	completeness
	)

private

Definition at line 245 of file NuMuEfficiencyStudy_module.cc.

                                                                                         {
   h.muon_completeness->Fill(completeness);
   h.muon_purity->Fill(purity);
 }

void numu::NuMuEfficiencyStudy::FillMatchingSlc	(	Histos &	h,
		float	purity,
		float	completeness
	)

private

Definition at line 250 of file NuMuEfficiencyStudy_module.cc.

                                                                                          {
   h.slice_completeness->Fill(completeness);
   h.slice_purity->Fill(purity);
 }

void numu::NuMuEfficiencyStudy::FillTrue	(	Histos &	h,
		const simb::MCTruth &	interaction,
		const simb::MCParticle &	G4lepton
	)

private

Definition at line 255 of file NuMuEfficiencyStudy_module.cc.

                                                                                                           {
   h.neutrino_energy->Fill(truth.GetNeutrino().Nu().E());
   h.muon_momentum->Fill(truth.GetNeutrino().Lepton().Momentum().Vect().Mag());
   h.muon_length->Fill(ParticleLength(G4lepton));
 
   h.vertex_x->Fill(truth.GetNeutrino().Nu().Position().X());
   h.vertex_y->Fill(truth.GetNeutrino().Nu().Position().Y());
   h.vertex_z->Fill(truth.GetNeutrino().Nu().Position().Z());
 }

bool numu::NuMuEfficiencyStudy::InFV ( TVector3 pos )

private

Definition at line 128 of file NuMuEfficiencyStudy_module.cc.

                                                {
   for (auto const &FV: fFiducialVolumes) {
     if (FV.ContainsPosition(pos)) return true;
   }
   return false;
 }

void numu::NuMuEfficiencyStudy::InitHisto	(	Histos &	h,
		const std::string &	prefix
	)

private

Definition at line 163 of file NuMuEfficiencyStudy_module.cc.

                                                                           {
   art::ServiceHandle<art::TFileService> tfs; 
 
   h.neutrino_energy = tfs->make<TH1D>((prefix + "_nuE").c_str(), "nuE", 100, 0, 4);
   h.muon_momentum = tfs->make<TH1D>((prefix + "_muP").c_str(), "muP", 100, 0, 4);
   h.muon_length = tfs->make<TH1D>((prefix + "_muL").c_str(), "muL", 100, 0, 600);
 
   float XMin = std::min_element(fFiducialVolumes.begin(), fFiducialVolumes.end(), [](auto &lhs, auto &rhs) { return lhs.MinX() < rhs.MinX(); })->MinX();
   float XMax = std::max_element(fFiducialVolumes.begin(), fFiducialVolumes.end(), [](auto &lhs, auto &rhs) { return lhs.MaxX() < rhs.MaxX(); })->MaxX();
   float YMin = std::min_element(fFiducialVolumes.begin(), fFiducialVolumes.end(), [](auto &lhs, auto &rhs) { return lhs.MinY() < rhs.MinY(); })->MinY();
   float YMax = std::max_element(fFiducialVolumes.begin(), fFiducialVolumes.end(), [](auto &lhs, auto &rhs) { return lhs.MaxY() < rhs.MaxY(); })->MaxY();
   float ZMin = std::min_element(fFiducialVolumes.begin(), fFiducialVolumes.end(), [](auto &lhs, auto &rhs) { return lhs.MinZ() < rhs.MinZ(); })->MinZ();
   float ZMax = std::max_element(fFiducialVolumes.begin(), fFiducialVolumes.end(), [](auto &lhs, auto &rhs) { return lhs.MaxZ() < rhs.MaxZ(); })->MaxZ();
                  
   h.vertex_x = tfs->make<TH1D>((prefix + "_v_x").c_str(), "v_x", 200, XMin, XMax); 
   h.vertex_y = tfs->make<TH1D>((prefix + "_v_y").c_str(), "v_y", 200, YMin, YMax); 
   h.vertex_z = tfs->make<TH1D>((prefix + "_v_z").c_str(), "v_z", 200, ZMin, ZMax); 
 
   h.muon_completeness = tfs->make<TH1D>((prefix + "_muC").c_str(), "muC", 100, 0., 1.);
   h.muon_purity = tfs->make<TH1D>((prefix + "_mupur").c_str(), "mupur", 100, 0., 1.);
 
   h.slice_completeness = tfs->make<TH1D>((prefix + "_nuC").c_str(), "nuC", 100, 0., 1.);
   h.slice_purity = tfs->make<TH1D>((prefix + "_nupur").c_str(), "nupur", 100, 0., 1.);
 }

void numu::NuMuEfficiencyStudy::InitHistos ( )

private

Definition at line 135 of file NuMuEfficiencyStudy_module.cc.

                                          {
   InitHisto(fTruth, "truth");
 
   std::array<std::string, n_reco_eff> effnames {"has_muon",
                                                 "has_reco_muon", 
                                                 "muon_is_fid", 
                                                 "has_slice", 
                                                 "has_pure_slice", 
                                                 "has_reco_slice", 
                                                 "has_fid_slice", 
                                                 "slice_has_muon", 
                                                 "has_fid_pure_slice", 
                                                 "has_reco_fid_slice", 
                                                 "slice_has_reco_muon",
                                                 "has_fid_nu_slice",
                                                 "has_pure_nu_slice",
                                                 "has_neutrino_slice"};
   for (unsigned i = 0; i < n_reco_eff; i++) {
     for (unsigned j = 0; j < n_reco_eff; j++) {
       std::string name;
       if (i == j) name = effnames[i];
       else name = effnames[i] + "_" + effnames[j];
       fReco.emplace_back();
       InitHisto(fReco.back(), name);
     }
   }
 }

NuMuEfficiencyStudy& numu::NuMuEfficiencyStudy::operator= ( NuMuEfficiencyStudy const & )

delete

NuMuEfficiencyStudy& numu::NuMuEfficiencyStudy::operator= ( NuMuEfficiencyStudy && )

delete

float numu::NuMuEfficiencyStudy::ParticleLength ( const simb::MCParticle & particle )

private

Definition at line 228 of file NuMuEfficiencyStudy_module.cc.

                                                                             {
   std::vector<geoalgo::AABox> aa_volumes;
   for (const geo::BoxBoundedGeo &AV: fActiveVolumes) {
     if (AV.ContainsPosition(particle.Position().Vect())) {
       aa_volumes.emplace_back(AV.MinX(), AV.MinY(), AV.MinZ(), AV.MaxX(), AV.MaxY(), AV.MaxZ());
       break;
     }
   }
 
   float length = 0.;
   for (unsigned i = 1; i < particle.NumberTrajectoryPoints(); i++) {
     length += ContainedLength(particle.Position(i-1).Vect(), particle.Position(i).Vect(), aa_volumes);
   }
 
   return length;
 }