MRCCAna_module.cc
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////////
2 // MRCCAna_module.cc
3 // Mike Wallbank, University of Cincinnati <wallbank@fnal.gov>, July 2020
4 //
5 // Class: MRCCAna
6 // Type: analyzer
7 // Description: MRCC validation module to analyze the hits which were
8 // removed by the MuonRemoval process, in particular with
9 // regards to NC interactions.
10 /////////////////////////////////////////////////////////////////////////////
11 
12 // framework
20 #include "fhiclcpp/ParameterSet.h"
21 
22 // nova
25 #include "MCCheater/BackTracker.h"
26 #include "RecoBase/CellHit.h"
27 #include "RecoBase/Cluster.h"
28 #include "RecoBase/Vertex.h"
29 #include "RecoBase/FitSum.h"
30 #include "RecoBase/FilterList.h"
31 
32 // stl
33 #include <string>
34 #include <vector>
35 #include <map>
36 
37 // ROOT
38 #include "TH1I.h"
39 #include "TH1F.h"
40 #include "TProfile.h"
41 
42 // --------------------------------------------------------------------------
43 namespace murem {
44  class MRCCNeutrino;
45  class MRCCParticle;
46  class MRCCAna;
47 }
48 
49 // --------------------------------------------------------------------------
51 public:
52 
54  const std::vector<const sim::Particle*>& particles);
55 
56  void AddAssociatedCluster(unsigned int clusterID);
57  void AddAssociatedMRCCCluster(unsigned int clusterID);
58  void AddParticleHit(unsigned int trackID, art::Ptr<rb::CellHit> hit);
59  void AddParticleMRCCHit(unsigned int trackID, art::Ptr<rb::CellHit> hit);
60 
61  float Energy() const;
62  bool CC() const;
63  float ClusterCharge() const;
64  float ClusterMRCCCharge() const;
65  float ClusterHadronicCharge() const;
66  float ClusterMRCCHadronicCharge() const;
67  float ClusterLeptonicCharge() const;
68  float ClusterMRCCLeptonicCharge() const;
69 
70  unsigned int NumAssociatedClusters() const;
71  unsigned int NumAssociatedMRCCClusters() const;
72 
74  std::map<unsigned int, MRCCParticle*> Particles() const;
75  int Muon() const;
76 
77 private:
78 
80 
81  std::map<unsigned int, MRCCParticle*> fParticles;
82  unsigned int* fMuon;
83 
84  std::vector<unsigned int> fAssociatedClusters;
85  std::vector<unsigned int> fAssociatedMRCCClusters;
86 
87 };
88 
89 // --------------------------------------------------------------------------
91 public:
92 
93  MRCCParticle(const simb::MCParticle* particle);
94 
95  void AddHit(art::Ptr<rb::CellHit> hit);
96  void AddMRCCHit(art::Ptr<rb::CellHit> hit);
97 
98  float Energy() const;
99  float DepositedCharge() const;
100  float DepositedMRCCCharge() const;
101 
102  const simb::MCParticle* TrueParticle() const;
103 
104 private:
105 
108 
110  std::vector<art::Ptr<rb::CellHit> > fHits;
111  std::vector<art::Ptr<rb::CellHit> > fMRCCHits;
112 
113 };
114 
115 // --------------------------------------------------------------------------
117  const std::vector<const sim::Particle*>& particles) {
118  fTruth = neutrino_int;
119  fMuon = nullptr;
120  for (std::vector<const sim::Particle*>::const_iterator particleIt = particles.begin();
121  particleIt != particles.end(); ++particleIt) {
122  fParticles[(*particleIt)->TrackId()] = new MRCCParticle(*particleIt);
123  if (!neutrino_int->GetNeutrino().CCNC() and !fMuon and (*particleIt)->PdgCode() == 13)
124  fMuon = new unsigned int((*particleIt)->TrackId());
125  }
126 }
127 
128 // --------------------------------------------------------------------------
129 void murem::MRCCNeutrino::AddAssociatedCluster(unsigned int clusterID) {
130  fAssociatedClusters.push_back(clusterID);
131  return;
132 }
133 
134 // --------------------------------------------------------------------------
135 void murem::MRCCNeutrino::AddAssociatedMRCCCluster(unsigned int clusterID) {
136  fAssociatedMRCCClusters.push_back(clusterID);
137  return;
138 }
139 
140 // --------------------------------------------------------------------------
142  fParticles[trackID]->AddHit(hit);
143  return;
144 }
145 
146 // --------------------------------------------------------------------------
148  fParticles[trackID]->AddMRCCHit(hit);
149  return;
150 }
151 
152 // --------------------------------------------------------------------------
154  return fAssociatedClusters.size();
155 }
156 
157 // --------------------------------------------------------------------------
159  return fAssociatedMRCCClusters.size();
160 }
161 
162 // --------------------------------------------------------------------------
164  return fTruth->GetNeutrino().Nu().E();
165 }
166 
167 // --------------------------------------------------------------------------
169  return !fTruth->GetNeutrino().CCNC();
170 }
171 
172 // --------------------------------------------------------------------------
174  float charge = 0.;
175  for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt = fParticles.begin();
176  particleIt != fParticles.end(); ++particleIt)
177  charge += particleIt->second->DepositedCharge();
178  return charge;
179 }
180 
181 // --------------------------------------------------------------------------
183  float charge = 0.;
184  for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt = fParticles.begin();
185  particleIt != fParticles.end(); ++particleIt)
186  charge += particleIt->second->DepositedMRCCCharge();
187  return charge;
188 }
189 
190 // --------------------------------------------------------------------------
192  float charge = 0.;
193  for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt = fParticles.begin();
194  particleIt != fParticles.end(); ++particleIt) {
195  if (fMuon and particleIt->first == *fMuon)
196  continue;
197  charge += particleIt->second->DepositedCharge();
198  }
199  return charge;
200 }
201 
202 // --------------------------------------------------------------------------
204  float charge = 0.;
205  for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt = fParticles.begin();
206  particleIt != fParticles.end(); ++particleIt) {
207  if (fMuon and particleIt->first == *fMuon)
208  continue;
209  charge += particleIt->second->DepositedMRCCCharge();
210  }
211  return charge;
212 }
213 
214 // --------------------------------------------------------------------------
216  float charge = 0.;
217  if (fMuon)
218  charge = fParticles.at(*fMuon)->DepositedCharge();
219  return charge;
220 }
221 
222 // --------------------------------------------------------------------------
224  float charge = 0.;
225  if (fMuon)
226  charge = fParticles.at(*fMuon)->DepositedMRCCCharge();
227  return charge;
228 }
229 
230 // --------------------------------------------------------------------------
232  return fTruth;
233 }
234 
235 // --------------------------------------------------------------------------
236 std::map<unsigned int, murem::MRCCParticle*> murem::MRCCNeutrino::Particles() const {
237  return fParticles;
238 }
239 
240 // --------------------------------------------------------------------------
242  if (fMuon) return *fMuon;
243  return -1;
244 }
245 
246 // --------------------------------------------------------------------------
248  fTotalCharge(0.),
249  fTotalMRCCCharge(0.),
250  fTrueParticle(particle) {
251 }
252 
253 // --------------------------------------------------------------------------
255  fHits.push_back(hit);
256  fTotalCharge += hit->ADC();
257  return;
258 }
259 
260 // --------------------------------------------------------------------------
262  fMRCCHits.push_back(hit);
263  fTotalMRCCCharge += hit->ADC();
264  return;
265 }
266 
267 // --------------------------------------------------------------------------
269  return fTrueParticle->E();
270 }
271 
272 // --------------------------------------------------------------------------
274  return fTotalCharge;
275 }
276 
277 // --------------------------------------------------------------------------
279  return fTotalMRCCCharge;
280 }
281 
282 // --------------------------------------------------------------------------
284  return fTrueParticle;
285 }
286 
287 // --------------------------------------------------------------------------
289 public:
290 
291  explicit MRCCAna(const fhicl::ParameterSet& pset);
292  virtual ~MRCCAna();
293 
294  void analyze(const art::Event& evt);
295  void beginJob();
296 
297 private:
298 
299  // methods
300  bool SliceSelection(const art::Event& evt, const art::Ptr<rb::Cluster>& slice);
301 
302  // labels
311 
312  // services
315 
316  // data products
336  TH1D* hHadYMRCC;
337  TH1D* hHadYNoMRCC;
341 
342 };
343 
345 
346 // --------------------------------------------------------------------------
347 murem::MRCCAna::MRCCAna(const fhicl::ParameterSet& pset) :
348  EDAnalyzer(pset),
349  fTruthModuleLabel(pset.get<std::string>("TruthModuleLabel")),
350  fHitsModuleLabel(pset.get<std::string>("HitsModuleLabel")),
351  fClusterModuleLabel(pset.get<std::string>("ClusterModuleLabel")),
352  fElasticModuleLabel(pset.get<std::string>("ElasticModuleLabel")),
353  fFuzzy3DModuleLabel(pset.get<std::string>("Fuzzy3DModuleLabel")),
354  fBPFTrackModuleLabel(pset.get<std::string>("BPFTrackModuleLabel")),
355  fMRCCHitsModuleLabel(pset.get<std::string>("MRCCHitsModuleLabel")),
356  fMRCCClusterModuleLabel(pset.get<std::string>("MRCCClusterModuleLabel")) {
357 }
358 
359 // --------------------------------------------------------------------------
361 }
362 
363 // --------------------------------------------------------------------------
365 
366  hNumNuClusters
367  = tfs->make<TH1I>("NumNuClusters", ";Number of Clusters Associated with True Neutrino;", 5,0,5);
368  hNumNuClusters->GetXaxis()->CenterLabels();
369  hNumNuClustersMRCC
370  = tfs->make<TH1I>("NumNuClustersMRCC", ";Number of Clusters Associated with True Neutrino;", 5,0,5);
371  hNumNuClustersMRCC->GetXaxis()->CenterLabels();
372  hSliceHadronicFrac
373  = tfs->make<TProfile>("SliceHadronicFrac", ";Neutrino Energy (GeV);Hadronic Fraction of Total Slice Charge;", 100,0,40);
374  hSliceHadronicFrac->SetMinimum(0);
375  hSliceHadronicFrac->SetMaximum(1);
376  hSliceHadronicFracMRCC
377  = tfs->make<TProfile>("SliceHadronicFracMRCC", ";Neutrino Energy (GeV);Hadronic Fraction of Total Slice Charge;", 100,0,40);
378  hSliceHadronicFracMRCC->SetMinimum(0);
379  hSliceHadronicFracMRCC->SetMaximum(1);
380  hAnyMRCCNeutrinoEnergy
381  = tfs->make<TProfile>("AnyMRCCNeutrinoEnergy", ";Neutrino Energy (GeV);Fraction of Neutrino Interactions with Charge Removed;", 100,0,40);
382  hAnyMRCCNeutrinoEnergy->SetMinimum(0);
383  hAnyMRCCNeutrinoEnergy->SetMaximum(1);
384  hAnyMRCCMuonEnergy
385  = tfs->make<TProfile>("AnyMRCCMuonEnergy", ";Muon Energy (GeV);Fraction of CC Muons with Charge Removed;", 100,0,40);
386  hAnyMRCCMuonEnergy->SetMinimum(0);
387  hAnyMRCCMuonEnergy->SetMaximum(1);
388  hSliceEFracRemoved
389  = tfs->make<TProfile>("SliceEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Charge Removed;", 100,0,40);
390  hSliceEFracRemoved->SetMinimum(0);
391  hSliceEFracRemoved->SetMaximum(1);
392  hSliceHadEFracRemoved
393  = tfs->make<TProfile>("SliceHadEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Hadronic Charge Removed;", 100,0,40);
394  hSliceHadEFracRemoved->SetMinimum(0);
395  hSliceHadEFracRemoved->SetMaximum(1);
396  hSliceLepEFracRemoved
397  = tfs->make<TProfile>("SliceLepEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Leptonic Charge Removed;", 100,0,40);
398  hSliceLepEFracRemoved->SetMinimum(0);
399  hSliceLepEFracRemoved->SetMaximum(1);
400  hSliceDeltaE
401  = tfs->make<TH1F>("SliceDeltaE", ";Fraction of Slice Charge Removed (GeV);", 100,0,1);
402  hSliceDeltaLepE
403  = tfs->make<TH1F>("SliceDeltaLepE", ";Fraction of Slice Leptonic Charge Removed (GeV);", 100,0,1);
404  hSliceDeltaHadE
405  = tfs->make<TH1F>("SliceDeltaHadE", ";Fraction of Slice Hadronic Charge Removed (GeV);", 100,0,1);
406  hLeptonEFracRemoved
407  = tfs->make<TProfile>("LeptonEFracRemoved", ";Lepton Energy (GeV);Fraction of Charge Removed;", 100,0,40);
408  hLeptonEFracRemoved->SetMinimum(0);
409  hLeptonEFracRemoved->SetMaximum(1);
410  hHadronEFracRemoved
411  = tfs->make<TProfile>("HadronEFracRemoved", ";Hadron Energy (GeV);Fraction of Charge Removed;", 100,0,40);
412  hHadronEFracRemoved->SetMinimum(0);
413  hHadronEFracRemoved->SetMaximum(1);
414  hMRCCSliceEFracRemoved
415  = tfs->make<TProfile>("MRCCSliceEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Charge Removed;", 100,0,40);
416  hMRCCSliceEFracRemoved->SetMinimum(0);
417  hMRCCSliceEFracRemoved->SetMaximum(1);
418  hMRCCSliceHadEFracRemoved
419  = tfs->make<TProfile>("MRCCSliceHadEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Hadronic Charge Removed;", 100,0,40);
420  hMRCCSliceHadEFracRemoved->SetMinimum(0);
421  hMRCCSliceHadEFracRemoved->SetMaximum(1);
422  hMRCCSliceLepEFracRemoved
423  = tfs->make<TProfile>("MRCCSliceLepEFracRemoved", ";Neutrino Energy (GeV);Fraction of Slice Leptonic Charge Removed;", 100,0,40);
424  hMRCCSliceLepEFracRemoved->SetMinimum(0);
425  hMRCCSliceLepEFracRemoved->SetMaximum(1);
426  hMRCCLeptonEFracRemoved
427  = tfs->make<TProfile>("MRCCLeptonEFracRemoved", ";Lepton Energy (GeV);Fraction of Charge Removed;", 100,0,40);
428  hMRCCLeptonEFracRemoved->SetMinimum(0);
429  hMRCCLeptonEFracRemoved->SetMaximum(1);
430  hMRCCHadronEFracRemoved
431  = tfs->make<TProfile>("MRCCHadronEFracRemoved", ";Hadron Energy (GeV);Fraction of Charge Removed;", 100,0,40);
432  hMRCCHadronEFracRemoved->SetMinimum(0);
433  hMRCCHadronEFracRemoved->SetMaximum(1);
434  hHadYMRCC = tfs->make<TH1D>("HadYMRCC", ";Hadronic Y;", 100, 0, 1);
435  hHadYNoMRCC = tfs->make<TH1D>("HadYNoMRCC", ";Hadronic Y;", 100, 0, 1);
436  hSliceEFracRemovedQE
437  = tfs->make<TProfile>("SliceEFracRemovedQE", ";Fraction of Slice Charge Removed (QE);", 100,0,40);
438  hSliceEFracRemovedDIS
439  = tfs->make<TProfile>("SliceEFracRemovedDIS", ";Fraction of Slice Charge Removed (DIS);", 100,0,40);
440  hSliceEFracRemovedRES
441  = tfs->make<TProfile>("SliceEFracRemovedRES", ";Fraction of Slice Charge Removed (RES);", 100,0,40);
442 
443  return;
444 
445 }
446 
447 // --------------------------------------------------------------------------
449 
450  // Get the truth information
452  std::vector<art::Ptr<simb::MCTruth> > truths;
453  if (evt.getByLabel(fTruthModuleLabel, truthHandle))
454  art::fill_ptr_vector(truths, truthHandle);
455 
456  // Get all the hits before muon removal
458  std::vector<art::Ptr<rb::CellHit> > hits;
459  if (evt.getByLabel(fHitsModuleLabel, hitsHandle))
460  art::fill_ptr_vector(hits, hitsHandle);
461 
462  // Get the clusters before muon removal
463  art::Handle<std::vector<rb::Cluster> > clusterHandle;
464  std::vector<art::Ptr<rb::Cluster> > clusters;
465  if (evt.getByLabel(fClusterModuleLabel, clusterHandle))
466  art::fill_ptr_vector(clusters, clusterHandle);
467 
468  // Get the hits following muon removal
469  art::Handle<std::vector<rb::CellHit> > mrccHitsHandle;
470  std::vector<art::Ptr<rb::CellHit> > mrccHits;
471  if (evt.getByLabel(fMRCCHitsModuleLabel, mrccHitsHandle))
472  art::fill_ptr_vector(mrccHits, mrccHitsHandle);
473 
474  // Get the clusters following muon removal
475  art::Handle<std::vector<rb::Cluster> > mrccClusterHandle;
476  std::vector<art::Ptr<rb::Cluster> > mrccClusters;
477  if (evt.getByLabel(fMRCCClusterModuleLabel, mrccClusterHandle))
478  art::fill_ptr_vector(mrccClusters, mrccClusterHandle);
479 
480  // Container to hold true neutrino information
481  std::map<unsigned int, MRCCNeutrino*> true_neutrinos;
482 
483  // Grab and save the neutrino interaction truth information
484  for (std::vector<art::Ptr<simb::MCTruth> >::const_iterator truthIt = truths.begin(); truthIt != truths.end(); ++truthIt) {
485  std::vector<const sim::Particle*> sim_particles = bt->MCTruthToParticles(*truthIt);
486  true_neutrinos[truthIt->key()] = new MRCCNeutrino(*truthIt, sim_particles);
487  }
488 
489  // Add all hits to the relevant true particle for each neutrino interaction, before muon removal
490  for (std::vector<art::Ptr<rb::Cluster> >::const_iterator clusterIt = clusters.begin();
491  clusterIt != clusters.end(); ++clusterIt) {
492  if ((*clusterIt)->IsNoise() or !SliceSelection(evt, *clusterIt))
493  continue;
494  std::vector<cheat::NeutrinoEffPur> nus = bt->SliceToNeutrinoInteractions(*clusterIt, hits);
495  if (nus.size() != 1)
496  continue;
497  true_neutrinos[nus[0].neutrinoInt.key()]->AddAssociatedCluster(clusterIt->key());
499  for (art::PtrVector<rb::CellHit>::const_iterator hitIt = cluster_hits.begin();
500  hitIt != cluster_hits.end(); ++hitIt) {
501  const sim::Particle* true_particle = bt->HitToParticle(*(*hitIt));
502  if (true_particle == nullptr)
503  continue;
504  true_neutrinos[nus[0].neutrinoInt.key()]->AddParticleHit(true_particle->TrackId(), *hitIt);
505  }
506  }
507 
508  // Add all hits to the relevant true particle for each neutrino interaction, after muon removal
509  for (std::vector<art::Ptr<rb::Cluster> >::const_iterator clusterIt = mrccClusters.begin();
510  clusterIt != mrccClusters.end(); ++clusterIt) {
511  if ((*clusterIt)->IsNoise())// or !SliceSelection(evt, *clusterIt))
512  continue;
513  std::vector<cheat::NeutrinoEffPur> nus = bt->SliceToNeutrinoInteractions(*clusterIt, mrccHits);
514  if (nus.size() != 1)
515  continue;
516  true_neutrinos[nus[0].neutrinoInt.key()]->AddAssociatedMRCCCluster(clusterIt->key());
518  for (art::PtrVector<rb::CellHit>::const_iterator hitIt = cluster_hits.begin();
519  hitIt != cluster_hits.end(); ++hitIt) {
520  const sim::Particle* true_particle = bt->HitToParticle(*(*hitIt));
521  if (true_particle == nullptr)
522  continue;
523  true_neutrinos[nus[0].neutrinoInt.key()]->AddParticleMRCCHit(true_particle->TrackId(), *hitIt);
524  }
525  }
526 
527  // Plot things!
528  for (std::map<unsigned int, MRCCNeutrino*>::const_iterator neutrinoIt = true_neutrinos.begin();
529  neutrinoIt != true_neutrinos.end(); ++neutrinoIt) {
530 
531  const MRCCNeutrino* true_nu = neutrinoIt->second;
532 
533  // Neutrino-level plots
534  hNumNuClusters->Fill(true_nu->NumAssociatedClusters());
535  hNumNuClustersMRCC->Fill(true_nu->NumAssociatedMRCCClusters());
536 
537  // For the slice plots, only consider unique neutrino slices which weren't affected by the muon removal
538  if (true_nu->NumAssociatedClusters() != 1 or true_nu->NumAssociatedClusters() != true_nu->NumAssociatedMRCCClusters())
539  continue;
540 
541  // Slice-level plots
542  hSliceHadronicFrac->Fill(true_nu->Energy(), true_nu->ClusterHadronicCharge()/true_nu->ClusterCharge());
543  hSliceHadronicFracMRCC->Fill(true_nu->Energy(), true_nu->ClusterMRCCHadronicCharge()/true_nu->ClusterMRCCCharge());
544 
545  float charge_before = true_nu->ClusterCharge(), charge_after = true_nu->ClusterMRCCCharge();
546  bool mrcc = charge_before != charge_after;
547 
548  if (true_nu->CC())
549  hAnyMRCCNeutrinoEnergy->Fill(true_nu->Energy(), mrcc);
550 
551  if (charge_before) {
552  hSliceEFracRemoved->Fill(true_nu->Energy(), (charge_before-charge_after)/charge_before);
553  hSliceDeltaE->Fill((charge_before-charge_after)/charge_before);
554  if (true_nu->Truth()->GetNeutrino().InteractionType() == 1001)
555  hSliceEFracRemovedQE->Fill(true_nu->Energy(), (charge_before-charge_after)/charge_before);
556  if (true_nu->Truth()->GetNeutrino().InteractionType() >= 1003 and true_nu->Truth()->GetNeutrino().InteractionType() <= 1090)
557  hSliceEFracRemovedRES->Fill(true_nu->Energy(), (charge_before-charge_after)/charge_before);
558  if (true_nu->Truth()->GetNeutrino().InteractionType() == 1091)
559  hSliceEFracRemovedDIS->Fill(true_nu->Energy(), (charge_before-charge_after)/charge_before);
560  if (mrcc)
561  hMRCCSliceEFracRemoved->Fill(true_nu->Energy(), (charge_before-charge_after)/charge_before);
562  }
563 
564  float had_charge_before = true_nu->ClusterHadronicCharge(), had_charge_after = true_nu->ClusterMRCCHadronicCharge();
565  if (had_charge_before) {
566  hSliceHadEFracRemoved->Fill(true_nu->Energy(), (had_charge_before-had_charge_after)/had_charge_before);
567  hSliceDeltaHadE->Fill((had_charge_before-had_charge_after)/had_charge_before);
568  if (mrcc)
569  hMRCCSliceHadEFracRemoved->Fill(true_nu->Energy(), (had_charge_before-had_charge_after)/had_charge_before);
570  }
571 
572  float lep_charge_before = true_nu->ClusterLeptonicCharge(), lep_charge_after = true_nu->ClusterMRCCLeptonicCharge();
573  if (lep_charge_before) {
574  hSliceLepEFracRemoved->Fill(true_nu->Energy(), (lep_charge_before-lep_charge_after)/lep_charge_before);
575  hSliceDeltaLepE->Fill((lep_charge_before-lep_charge_after)/lep_charge_before);
576  if (mrcc)
577  hMRCCSliceLepEFracRemoved->Fill(true_nu->Energy(), (lep_charge_before-lep_charge_after)/lep_charge_before);
578  }
579 
580  if (mrcc)
581  hHadYMRCC->Fill(true_nu->Truth()->GetNeutrino().Y());
582  else
583  hHadYNoMRCC->Fill(true_nu->Truth()->GetNeutrino().Y());
584 
585  // Particle-level plots
586  std::map<unsigned int, MRCCParticle*> particles = true_nu->Particles();
587  int muon = true_nu->Muon();
588  for (std::map<unsigned int, MRCCParticle*>::const_iterator particleIt = particles.begin();
589  particleIt != particles.end(); ++particleIt) {
590  float p_charge_before = particleIt->second->DepositedCharge();
591  float p_charge_after = particleIt->second->DepositedMRCCCharge();
592  if (!p_charge_before)
593  continue;
594  if ((int)particleIt->first == muon) {
595  hLeptonEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
596  hAnyMRCCMuonEnergy->Fill(particleIt->second->Energy(), mrcc);
597  if (mrcc)
598  hMRCCLeptonEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
599  } else {
600  hHadronEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
601  if (mrcc)
602  hMRCCHadronEFracRemoved->Fill(particleIt->second->Energy(), (p_charge_before-p_charge_after)/p_charge_before);
603  }
604  }
605 
606  }
607 
608  return;
609 
610 }
611 
612 // --------------------------------------------------------------------------
614 
615  // Copied from MuonRemove_module
616 
617  // rock filtered?
618  if (rb::IsFiltered(evt, slice, "rockpresel"))
619  return false;
620 
621  // find vertex
622  art::FindManyP<rb::Vertex> fmElastic({slice}, evt, fElasticModuleLabel);
623  if (!fmElastic.isValid())
624  return false;
625 
626  std::vector<art::Ptr<rb::Vertex> > elastics = fmElastic.at(0);
627  if (elastics.empty())
628  return false;
629 
630  // find prongs
631  art::FindManyP<rb::Prong> fmFuzzyProng3D(elastics, evt, fFuzzy3DModuleLabel);
632  std::vector<art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(0);
633  if (prongs.empty())
634  return false;
635 
636  // find energy of parent muon
637  double muonEmax = -1;
638  art::FindManyP<rb::Track> fmt({slice}, evt, fBPFTrackModuleLabel);
639  std::vector<art::Ptr<rb::Track> > tracks = fmt.at(0);
640  for (unsigned int trackIt = 0; trackIt < tracks.size(); ++trackIt) {
641  art::FindManyP<rb::FitSum> fmBPFFitSums({tracks[trackIt]}, evt, fBPFTrackModuleLabel);
642  if (fmBPFFitSums.isValid()) {
643  std::vector<art::Ptr<rb::FitSum> > bpfFitSums = fmBPFFitSums.at(0);
644  if (abs(bpfFitSums[0]->PDG()) == 13) {
645  double muonE = bpfFitSums[0]->FourMom().E();
646  if (muonE > muonEmax)
647  muonEmax = muonE;
648  }
649  }
650  }
651  if (muonEmax <= 0)
652  return false;
653 
654  return true;
655 
656 }
double E(const int i=0) const
Definition: MCParticle.h:232
float ClusterCharge() const
virtual ~MRCCAna()
back track the reconstruction to the simulation
int CCNC() const
Definition: MCNeutrino.h:148
float ClusterMRCCCharge() const
MRCCNeutrino(art::Ptr< simb::MCTruth > neutrino_int, const std::vector< const sim::Particle * > &particles)
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
TProfile * hSliceEFracRemovedRES
TProfile * hSliceHadEFracRemoved
art::Ptr< simb::MCTruth > fTruth
rb::Cluster cluster_hits(std::vector< art::Ptr< rb::CellHit > > const &hits, int32_t t_start, int32_t t_end, int16_t adc_cut=0)
Definition: DDTHelpers.cxx:13
std::vector< art::Ptr< rb::CellHit > > fMRCCHits
iterator begin()
Definition: PtrVector.h:223
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
std::string fTruthModuleLabel
void AddHit(art::Ptr< rb::CellHit > hit)
float ClusterMRCCLeptonicCharge() const
void abs(TH1 *hist)
TProfile * hSliceEFracRemoved
unsigned int NumAssociatedClusters() const
DEFINE_ART_MODULE(TestTMapFile)
TProfile * hLeptonEFracRemoved
const simb::MCParticle * TrueParticle() const
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
std::string fFuzzy3DModuleLabel
int TrackId() const
Definition: MCParticle.h:209
TProfile * hMRCCSliceHadEFracRemoved
TProfile * hSliceLepEFracRemoved
int InteractionType() const
Definition: MCNeutrino.h:150
TProfile * hAnyMRCCNeutrinoEnergy
void AddAssociatedCluster(unsigned int clusterID)
TH1I * hNumNuClustersMRCC
void hits()
Definition: readHits.C:15
TProfile * hMRCCSliceEFracRemoved
void beginJob()
double Y() const
Definition: MCNeutrino.h:156
std::string fBPFTrackModuleLabel
bool SliceSelection(const art::Event &evt, const art::Ptr< rb::Cluster > &slice)
float DepositedMRCCCharge() const
std::string fHitsModuleLabel
int evt
iterator end()
Definition: PtrVector.h:237
float DepositedCharge() const
std::string fElasticModuleLabel
art::Ptr< simb::MCTruth > Truth() const
TProfile * hMRCCSliceLepEFracRemoved
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
float ClusterMRCCHadronicCharge() const
TProfile * hSliceHadronicFracMRCC
TProfile * hHadronEFracRemoved
float Energy() const
std::string fMRCCClusterModuleLabel
Vertex location in position and time.
TProfile * hMRCCLeptonEFracRemoved
TProfile * hAnyMRCCMuonEnergy
TProfile * hMRCCHadronEFracRemoved
art::ServiceHandle< cheat::BackTracker > bt
std::vector< art::Ptr< rb::CellHit > > fHits
data_t::const_iterator const_iterator
Definition: PtrVector.h:61
std::string fClusterModuleLabel
float ClusterHadronicCharge() const
TProfile * hSliceHadronicFrac
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void AddAssociatedMRCCCluster(unsigned int clusterID)
std::string fMRCCHitsModuleLabel
void AddMRCCHit(art::Ptr< rb::CellHit > hit)
std::vector< unsigned int > fAssociatedClusters
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
float Energy() const
unsigned int * fMuon
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::map< unsigned int, MRCCParticle * > Particles() const
std::map< unsigned int, MRCCParticle * > fParticles
Definition: structs.h:12
void AddParticleMRCCHit(unsigned int trackID, art::Ptr< rb::CellHit > hit)
TProfile * hSliceEFracRemovedQE
MRCCParticle(const simb::MCParticle *particle)
std::vector< unsigned int > fAssociatedMRCCClusters
void analyze(const art::Event &evt)
art::ServiceHandle< art::TFileService > tfs
TProfile * hSliceEFracRemovedDIS
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Int_t trackID
Definition: plot.C:84
void AddParticleHit(unsigned int trackID, art::Ptr< rb::CellHit > hit)
unsigned int NumAssociatedMRCCClusters() const
const simb::MCParticle * fTrueParticle
float ClusterLeptonicCharge() const