LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
NeutrinoTrackingEff_module.cc
Go to the documentation of this file.
1 //
2 //**Tracking Efficiency module***
3 //The basic idea is to loop over the hits from a given track and call BackTracker
4 //then look at std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackID(hit);
5 //then associete the hits to a G4 track ID (particle) that generate those hits(electrons)
6 //It was developed for CC neutrio interactions, it also can handle proton decay events p->k+nu_bar
7 //And protons, pion and muons from particle cannon by using the option isNeutrinoInt = false;
8 //All histograms are as a function of truth info (momentum, length)
9 //
10 // A. Higuera
11 // ahiguera@central.uh.edu
12 
13 // LArSoft includes
14 #include "larcore/CoreUtils/ServiceUtil.h"
15 #include "larcore/Geometry/Geometry.h"
16 #include "larcorealg/Geometry/TPCGeo.h"
17 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
18 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
19 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
20 #include "lardataobj/RecoBase/Track.h"
21 #include "larsim/MCCheater/BackTrackerService.h"
22 #include "larsim/MCCheater/ParticleInventoryService.h"
23 #include "nusimdata/SimulationBase/MCParticle.h"
24 #include "nusimdata/SimulationBase/MCTruth.h"
25 
26 // Framework includes
27 #include "art/Framework/Core/EDAnalyzer.h"
28 #include "art/Framework/Core/ModuleMacros.h"
29 #include "art/Framework/Principal/Event.h"
30 #include "art/Framework/Principal/Handle.h"
31 #include "art/Framework/Services/Registry/ServiceHandle.h"
32 #include "art_root_io/TFileService.h"
33 #include "canvas/Persistency/Common/FindManyP.h"
34 #include "fhiclcpp/ParameterSet.h"
35 #include "messagefacility/MessageLogger/MessageLogger.h"
36 
37 // ROOT includes
38 #include "TEfficiency.h"
39 #include "TGraphAsymmErrors.h"
40 #include "TH1.h"
41 
42 // C/C++ standard libraries
43 
44 #define MAX_TRACKS 1000
45 using namespace std;
46 
47 //========================================================================
48 
49 namespace DUNE {
50 
51  class NeutrinoTrackingEff : public art::EDAnalyzer {
52  public:
53  explicit NeutrinoTrackingEff(fhicl::ParameterSet const& pset);
54 
55  private:
56  void beginJob();
57  void endJob();
58  void beginRun(const art::Run& run);
59  void analyze(const art::Event& evt);
60 
61  void processEff(const art::Event& evt);
62  void truthMatcher(detinfo::DetectorClocksData const& clockData,
63  std::vector<art::Ptr<recob::Hit>> all_hits,
64  std::vector<art::Ptr<recob::Hit>> track_hits,
65  const simb::MCParticle*& MCparticle,
66  double& Efrac,
67  double& Ecomplet);
68  double truthLength(const detinfo::DetectorClocksData& clockData,
69  detinfo::DetectorPropertiesData const& detProp,
70  const simb::MCParticle* MCparticle);
71  bool insideFV(double vertex[4]);
72  void doEfficiencies();
73 
74  // the parameters we'll read from the .fcl
75  std::string fMCTruthModuleLabel;
76  std::string fTrackModuleLabel;
77  int fNeutrinoPDGcode;
78  int fLeptonPDGcode;
79  double fMaxNeutrinoE;
80  double fMaxLeptonP;
81  bool fisNeutrinoInt;
82 
83  int MC_isCC;
84  int MC_incoming_PDG;
85  double MC_incoming_P[4];
86  double MC_vertex[4];
87  double MC_lepton_startMomentum[4];
88 
89  int MC_leading_protonID;
90  int MC_leading_PionPlusID;
91  int MC_leading_PionMinusID;
92  int MC_leptonID;
93  int MC_kaonID;
94  int MC_michelID;
95 
96  double MC_leptonP;
97  double MC_leading_PionPlusP;
98  double MC_leading_ProtonP;
99  double MC_leading_PionMinusP;
100  double MC_kaonP;
101  double MC_michelP;
102 
103  TH1D* h_Ev_den;
104  TH1D* h_Ev_num;
105  TH1D* h_Pmu_den;
106  TH1D* h_Pmu_num;
107  TH1D* h_theta_den;
108  TH1D* h_theta_num;
109  TH1D* h_Pproton_den;
110  TH1D* h_Pproton_num;
111  TH1D* h_Ppion_plus_den;
112  TH1D* h_Ppion_plus_num;
113  TH1D* h_Ppion_minus_den;
114  TH1D* h_Ppion_minus_num;
115 
116  TH1D* h_Efrac_lepton;
117  TH1D* h_Ecomplet_lepton;
118  TH1D* h_Efrac_proton;
119  TH1D* h_Ecomplet_proton;
120  TH1D* h_Efrac_pion_plus;
121  TH1D* h_Ecomplet_pion_plus;
122  TH1D* h_Efrac_pion_minus;
123  TH1D* h_Ecomplet_pion_minus;
124  TH1D* h_trackRes_lepton;
125  TH1D* h_trackRes_proton;
126  TH1D* h_trackRes_pion_plus;
127  TH1D* h_trackRes_pion_minus;
128 
129  TH1D* h_muon_length;
130  TH1D* h_proton_length;
131  TH1D* h_pionp_length;
132  TH1D* h_pionm_length;
133  TH1D* h_muonwtrk_length;
134  TH1D* h_protonwtrk_length;
135  TH1D* h_pionpwtrk_length;
136  TH1D* h_pionmwtrk_length;
137 
138  TEfficiency* h_Eff_Ev = 0;
139  TEfficiency* h_Eff_Pmu = 0;
140  TEfficiency* h_Eff_theta = 0;
141  TEfficiency* h_Eff_Pproton = 0;
142  TEfficiency* h_Eff_Ppion_plus = 0;
143  TEfficiency* h_Eff_Ppion_minus = 0;
144 
145  TEfficiency* h_Eff_Lmuon = 0;
146  TEfficiency* h_Eff_Lproton = 0;
147  TEfficiency* h_Eff_Lpion_plus = 0;
148  TEfficiency* h_Eff_Lpion_minus = 0;
149 
150  //nucleon decay histograms
151  TH1D* h_Pkaon_den;
152  TH1D* h_Pkaon_num;
153  TH1D* h_Pmichel_e_den;
154  TH1D* h_Pmichel_e_num;
155  TH1D* h_Efrac_kaon;
156  TH1D* h_Ecomplet_kaon;
157  TH1D* h_trackRes_kaon;
158  TH1D* h_Efrac_michel;
159  TH1D* h_Ecomplet_michel;
160  TH1D* h_trackRes_michel;
161  TH1D* h_kaon_length;
162  TH1D* h_michel_length;
163  TH1D* h_kaonwtrk_length;
164  TH1D* h_michelwtrk_length;
165  TEfficiency* h_Eff_Pkaon = 0;
166  TEfficiency* h_Eff_Pmichel = 0;
167  TEfficiency* h_Eff_Lkaon = 0;
168  TEfficiency* h_Eff_Lmichel = 0;
169 
170  float fFidVolCutX;
171  float fFidVolCutY;
172  float fFidVolCutZ;
173 
174  float fFidVolXmin;
175  float fFidVolXmax;
176  float fFidVolYmin;
177  float fFidVolYmax;
178  float fFidVolZmin;
179  float fFidVolZmax;
180 
181  double fDriftVelocity; // in cm/ns
182  art::ServiceHandle<geo::Geometry const> geom;
183 
184  }; // class NeutrinoTrackingEff
185 
186  //========================================================================
187  NeutrinoTrackingEff::NeutrinoTrackingEff(fhicl::ParameterSet const& p) : EDAnalyzer(p)
188  {
189  fMCTruthModuleLabel = p.get<std::string>("MCTruthModuleLabel");
190  fTrackModuleLabel = p.get<std::string>("TrackModuleLabel");
191  fisNeutrinoInt = p.get<bool>("isNeutrinoInt");
192  fLeptonPDGcode = p.get<int>("LeptonPDGcode");
193  fNeutrinoPDGcode = p.get<int>("NeutrinoPDGcode");
194  fMaxNeutrinoE = p.get<double>("MaxNeutrinoE");
195  fMaxLeptonP = p.get<double>("MaxLeptonP");
196  fFidVolCutX = p.get<float>("FidVolCutX");
197  fFidVolCutY = p.get<float>("FidVolCutY");
198  fFidVolCutZ = p.get<float>("FidVolCutZ");
199 
200  auto const detProp =
201  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
202  fDriftVelocity = detProp.DriftVelocity() * 1e-3; // cm/ns
203  }
204  //========================================================================
205  void NeutrinoTrackingEff::beginJob()
206  {
207  std::cout << "job begin..." << std::endl;
208  auto const* geo = lar::providerFrom<geo::Geometry>();
209  // Define histogram boundaries (cm).
210  // For now only draw cryostat=0.
211  double minx = 1e9;
212  double maxx = -1e9;
213  double miny = 1e9;
214  double maxy = -1e9;
215  double minz = 1e9;
216  double maxz = -1e9;
217  for (auto const& tpc : geo->Iterate<geo::TPCGeo>(geo::CryostatID{0})) {
218  auto const world = tpc.GetCenter();
219  if (minx > world.X() - tpc.HalfWidth()) minx = world.X() - tpc.HalfWidth();
220  if (maxx < world.X() + tpc.HalfWidth()) maxx = world.X() + tpc.HalfWidth();
221  if (miny > world.Y() - tpc.HalfHeight()) miny = world.Y() - tpc.HalfHeight();
222  if (maxy < world.Y() + tpc.HalfHeight()) maxy = world.Y() + tpc.HalfHeight();
223  if (minz > world.Z() - tpc.Length() / 2.) minz = world.Z() - tpc.Length() / 2.;
224  if (maxz < world.Z() + tpc.Length() / 2.) maxz = world.Z() + tpc.Length() / 2.;
225  }
226 
227  fFidVolXmin = minx + fFidVolCutX;
228  fFidVolXmax = maxx - fFidVolCutX;
229  fFidVolYmin = miny + fFidVolCutY;
230  fFidVolYmax = maxy - fFidVolCutY;
231  fFidVolZmin = minz + fFidVolCutZ;
232  fFidVolZmax = maxz - fFidVolCutZ;
233 
234  std::cout << "Fiducial volume:"
235  << "\n"
236  << fFidVolXmin << "\t< x <\t" << fFidVolXmax << "\n"
237  << fFidVolYmin << "\t< y <\t" << fFidVolYmax << "\n"
238  << fFidVolZmin << "\t< z <\t" << fFidVolZmax << "\n";
239 
240  art::ServiceHandle<art::TFileService const> tfs;
241 
242  double E_bins[21] = {0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4, 4.5, 5.0,
243  5.5, 6.0, 7.0, 8.0, 10.0, 12.0, 14.0, 17.0, 20.0, 25.0};
244  double theta_bin[44] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.,
245  11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 22.,
246  24., 26., 28., 30., 32., 34., 36., 38., 40., 42., 44.,
247  46., 48., 50., 55., 60., 65., 70., 75., 80., 85., 90.};
248  double Pbins[18] = {
249  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0};
250 
251  for (int i = 0; i < 21; ++i)
252  E_bins[i] *= fMaxNeutrinoE / 25.;
253  for (int i = 0; i < 18; ++i)
254  Pbins[i] *= fMaxLeptonP / 3.0;
255 
256  h_Ev_den = tfs->make<TH1D>(
257  "h_Ev_den", "Neutrino Energy; Neutrino Energy (GeV); Tracking Efficiency", 20, E_bins);
258  h_Ev_num = tfs->make<TH1D>(
259  "h_Ev_num", "Neutrino Energy; Neutrino Energy (GeV); Tracking Efficiency", 20, E_bins);
260  h_Pmu_den = tfs->make<TH1D>(
261  "h_Pmu_den", "Muon Momentum; Muon Momentum (GeV); Tracking Efficiency", 20, E_bins);
262  h_Pmu_num = tfs->make<TH1D>(
263  "h_Pmu_num", "Muon Momentum; Muon Momentum (GeV); Tracking Efficiency", 20, E_bins);
264  h_theta_den =
265  tfs->make<TH1D>("h_theta_den",
266  "Theta; Theta w.r.t beam direction (Degrees); Tracking Efficiency",
267  43,
268  theta_bin);
269  h_theta_num =
270  tfs->make<TH1D>("h_theta_num",
271  "Theta; Theta w.r.t beam direction (Degrees); Tracking Efficiency",
272  43,
273  theta_bin);
274  h_Pproton_den = tfs->make<TH1D>(
275  "h_Pproton_den", "Protons; Proton Momentum (GeV); Tracking Efficiency", 17, Pbins);
276  h_Pproton_num = tfs->make<TH1D>(
277  "h_Pproton_num", "Protons; Proton Momentum (GeV); Tracking Efficiency", 17, Pbins);
278  h_Ppion_plus_den = tfs->make<TH1D>(
279  "h_Ppion_plus_den", "Pions Plus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
280  h_Ppion_plus_num = tfs->make<TH1D>(
281  "h_Ppion_plus_num", "Pions Plus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
282  h_Ppion_minus_den = tfs->make<TH1D>(
283  "h_Ppion_minus_den", "Pions Minus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
284  h_Ppion_minus_num = tfs->make<TH1D>(
285  "h_Ppion_minus_num", "Pions Minus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
286  h_Ev_den->Sumw2();
287  h_Ev_num->Sumw2();
288  h_Pmu_den->Sumw2();
289  h_Pmu_num->Sumw2();
290  h_theta_den->Sumw2();
291  h_theta_num->Sumw2();
292  h_Pproton_den->Sumw2();
293  h_Pproton_num->Sumw2();
294  h_Ppion_plus_den->Sumw2();
295  h_Ppion_plus_num->Sumw2();
296  h_Ppion_minus_den->Sumw2();
297  h_Ppion_minus_num->Sumw2();
298 
299  h_Efrac_lepton = tfs->make<TH1D>("h_Efrac_lepton", "Efrac Lepton; Track Purity;", 60, 0, 1.2);
300  h_Ecomplet_lepton =
301  tfs->make<TH1D>("h_Ecomplet_lepton", "Ecomplet Lepton; Track Completeness;", 60, 0, 1.2);
302  h_Efrac_proton = tfs->make<TH1D>("h_Efrac_proton", "Efrac Proton; Track Purity;", 60, 0, 1.2);
303  h_Ecomplet_proton =
304  tfs->make<TH1D>("h_Ecomplet_proton", "Ecomplet Proton; Track Completeness;", 60, 0, 1.2);
305  h_Efrac_pion_plus =
306  tfs->make<TH1D>("h_Efrac_pion_plus", "Efrac Pion +; Track Purity;", 60, 0, 1.2);
307  h_Ecomplet_pion_plus =
308  tfs->make<TH1D>("h_Ecomplet_pion_plus", "Ecomplet Pion +; Track Completeness;", 60, 0, 1.2);
309  h_Efrac_pion_minus =
310  tfs->make<TH1D>("h_Efrac_pion_minus", "Efrac Pion -; Track Purity;", 60, 0, 1.2);
311  h_Ecomplet_pion_minus =
312  tfs->make<TH1D>("h_Ecomplet_pion_minus", "Ecomplet Pion -; Track Completeness;", 60, 0, 1.2);
313  h_trackRes_lepton = tfs->make<TH1D>(
314  "h_trackRes_lepton", "Muon Residual; Truth length - Reco length (cm);", 200, -100, 100);
315  h_trackRes_proton = tfs->make<TH1D>(
316  "h_trackRes_proton", "Proton Residual; Truth length - Reco length (cm);", 200, -100, 100);
317  h_trackRes_pion_plus = tfs->make<TH1D>(
318  "h_trackRes_pion_plus", "Pion + Residual; Truth length - Reco length (cm);", 200, -100, 100);
319  h_trackRes_pion_minus = tfs->make<TH1D>(
320  "h_trackRes_pion_minus", "Pion - Residual; Truth length - Reco length (cm);", 200, -100, 100);
321  h_Efrac_lepton->Sumw2();
322  h_Ecomplet_lepton->Sumw2();
323  h_Efrac_proton->Sumw2();
324  h_Ecomplet_proton->Sumw2();
325  h_Efrac_pion_plus->Sumw2();
326  h_Ecomplet_pion_plus->Sumw2();
327  h_Efrac_pion_minus->Sumw2();
328  h_Ecomplet_pion_minus->Sumw2();
329  h_trackRes_lepton->Sumw2();
330  h_trackRes_proton->Sumw2();
331  h_trackRes_pion_plus->Sumw2();
332  h_trackRes_pion_minus->Sumw2();
333 
334  h_muon_length =
335  tfs->make<TH1D>("h_muon_length", "Muon Length; Muon Truth Length (cm)", 40, 0, 100);
336  h_proton_length =
337  tfs->make<TH1D>("h_proton_length", "Proton Length; Proton Truth Length (cm)", 40, 0, 100);
338  h_pionp_length =
339  tfs->make<TH1D>("h_pionp_length", "Pion + Length; Pion^{+} Truth Length (cm)", 40, 0, 100);
340  h_pionm_length =
341  tfs->make<TH1D>("h_pionm_length", "Pion - Length; Pion^{-} Truth Length (cm)", 40, 0, 100);
342 
343  h_muonwtrk_length =
344  tfs->make<TH1D>("h_muonwtrk_length", "Muon Length; Muon Truth Length (cm)", 40, 0, 100);
345  h_protonwtrk_length =
346  tfs->make<TH1D>("h_protonwtrk_length", "Proton Length; Proton Truth Length (cm)", 40, 0, 100);
347  h_pionpwtrk_length = tfs->make<TH1D>(
348  "h_pionpwtrk_length", "Pion + Length; Pion^{+} Truth Length (cm)", 40, 0, 100);
349  h_pionmwtrk_length = tfs->make<TH1D>(
350  "h_pionmwtrk_length", "Pion - Length; Pion^{-} Truth Length (cm)", 40, 0, 100);
351 
352  h_muon_length->Sumw2();
353  h_muonwtrk_length->Sumw2();
354  h_proton_length->Sumw2();
355  h_protonwtrk_length->Sumw2();
356  h_pionp_length->Sumw2();
357  h_pionpwtrk_length->Sumw2();
358  h_pionm_length->Sumw2();
359  h_pionmwtrk_length->Sumw2();
360 
361  h_Pkaon_den =
362  tfs->make<TH1D>("h_Pkaon_den", "Kaon; Kaon Momentum (GeV); Tracking Efficiency", 17, Pbins);
363  h_Pkaon_num =
364  tfs->make<TH1D>("h_Pkaon_num", "Kaon; Kaon Momentum (GeV); Tracking Efficiency", 17, Pbins);
365  h_Pmichel_e_den =
366  tfs->make<TH1D>("h_Pmichel_e_den",
367  "Michel Electron; Michele e Momentum (GeV); Tracking Efficiency",
368  17,
369  Pbins);
370  h_Pmichel_e_num =
371  tfs->make<TH1D>("h_Pmichel_e_num",
372  "Michel Electron; Michele e Momentum (GeV); Tracking Efficiency",
373  17,
374  Pbins);
375  h_Pkaon_den->Sumw2();
376  h_Pkaon_num->Sumw2();
377  h_Pmichel_e_den->Sumw2();
378  h_Pmichel_e_num->Sumw2();
379  h_Efrac_kaon = tfs->make<TH1D>("h_Efrac_kaon", "Efrac Kaon; Track Purity;", 60, 0, 1.2);
380  h_Ecomplet_kaon =
381  tfs->make<TH1D>("h_Ecomplet_kaon", "Ecomplet Kaon; Track Completeness;", 60, 0, 1.2);
382  h_trackRes_kaon = tfs->make<TH1D>(
383  "h_trackRes_kaon", "Kaon Residual; Truth length - Reco length (cm);", 200, -100, 100);
384  h_Efrac_michel =
385  tfs->make<TH1D>("h_Efrac_michel", "Efrac Michel; Track Energy fraction;", 60, 0, 1.2);
386  h_Ecomplet_michel =
387  tfs->make<TH1D>("h_Ecomplet_michel", "Ecomplet Michel; Track Completeness;", 60, 0, 1.2);
388  h_trackRes_michel = tfs->make<TH1D>(
389  "h_trackRes_michel", "Michel Residual; Truth length - Reco length (cm);", 200, -100, 100);
390  h_kaon_length =
391  tfs->make<TH1D>("h_kaon_length", "Kaon Length; Kaon Truth Length (cm)", 40, 0, 100);
392  h_kaonwtrk_length =
393  tfs->make<TH1D>("h_kaonwtrk_length", "Kaon Length; Kaon Truth Length (cm)", 40, 0, 100);
394  h_michel_length =
395  tfs->make<TH1D>("h_michel_length", "Michel Length; Michel e Truth Length (cm)", 40, 0, 100);
396  h_michelwtrk_length = tfs->make<TH1D>(
397  "h_michelwtrk_length", "Michel Length; Michel e Truth Length (cm)", 40, 0, 100);
398 
399  h_Efrac_kaon->Sumw2();
400  h_Ecomplet_kaon->Sumw2();
401  h_trackRes_kaon->Sumw2();
402  h_Efrac_michel->Sumw2();
403  h_Ecomplet_michel->Sumw2();
404  h_trackRes_michel->Sumw2();
405  h_kaon_length->Sumw2();
406  h_kaonwtrk_length->Sumw2();
407  h_michel_length->Sumw2();
408  h_michelwtrk_length->Sumw2();
409  }
410  //========================================================================
411  void NeutrinoTrackingEff::endJob()
412  {
413  doEfficiencies();
414  }
415  //========================================================================
416  void NeutrinoTrackingEff::beginRun(const art::Run& /*run*/)
417  {
418  mf::LogInfo("NeutrinoTrackingEff") << "begin run..." << std::endl;
419  }
420  //========================================================================
421  void NeutrinoTrackingEff::analyze(const art::Event& event)
422  {
423  if (event.isRealData()) return;
424 
425  processEff(event);
426  }
427  //========================================================================
428  void NeutrinoTrackingEff::processEff(const art::Event& event)
429  {
430  // Save neutrino's interaction info
431  art::Handle<std::vector<simb::MCTruth>> MCtruthHandle;
432  event.getByLabel(fMCTruthModuleLabel, MCtruthHandle);
433  std::vector<art::Ptr<simb::MCTruth>> MCtruthlist;
434  art::fill_ptr_vector(MCtruthlist, MCtruthHandle);
435  art::Ptr<simb::MCTruth> MCtruth;
436 
437  //For now assume that there is only one neutrino interaction...
438  MCtruth = MCtruthlist[0];
439  if (MCtruth->NeutrinoSet()) {
440  simb::MCNeutrino nu = MCtruth->GetNeutrino();
441  if (nu.CCNC() == 0)
442  MC_isCC = 1;
443  else if (nu.CCNC() == 1)
444  MC_isCC = 0;
445  simb::MCParticle neutrino = nu.Nu();
446  MC_incoming_PDG = nu.Nu().PdgCode();
447  const TLorentzVector& nu_momentum = nu.Nu().Momentum(0);
448  nu_momentum.GetXYZT(MC_incoming_P);
449  const TLorentzVector& vertex = neutrino.Position(0);
450  vertex.GetXYZT(MC_vertex);
451  }
452 
454 
455  double tmp_leadingPionPlusE = 0.0;
456  double tmp_leadingPionMinusE = 0.0;
457  double tmp_leadingProtonE = 0.0;
458 
459  simb::MCParticle* MClepton = nullptr;
460  simb::MCParticle* MCproton = nullptr;
461  simb::MCParticle* MCpion_plus = nullptr;
462  simb::MCParticle* MCpion_minus = nullptr;
463  simb::MCParticle* MCkaon = nullptr;
464  simb::MCParticle* MCmichel = nullptr;
465 
466  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
467  const sim::ParticleList& plist = pi_serv->ParticleList();
468  simb::MCParticle* particle = 0;
469  int i = 0; // particle index
470 
471  auto const clockData =
472  art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
473  auto const detProp =
474  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event);
475 
476  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar) {
477  particle = ipar->second;
478  if (particle->PdgCode() == fLeptonPDGcode && particle->Mother() == 0) { //primary lepton
479  const TLorentzVector& lepton_momentum = particle->Momentum(0);
480  lepton_momentum.GetXYZT(MC_lepton_startMomentum);
481  MC_leptonID = particle->TrackId();
482  MC_leptonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
483  pow(particle->Momentum().Pz(), 2));
484  MClepton = particle;
485  }
486  if (particle->Mother() == 0) { //save primary particle i.e. from the neutrino interaction
487  //save leading pion and proton
488  if (particle->PdgCode() == 2212) {
489  if (particle->Momentum().E() > tmp_leadingProtonE) {
490  tmp_leadingProtonE = particle->Momentum().E();
491  MC_leading_protonID = particle->TrackId();
492  MC_leading_ProtonP =
493  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
494  pow(particle->Momentum().Pz(), 2));
495  MCproton = particle;
496  }
497  }
498  else if (particle->PdgCode() == 211) {
499  if (particle->Momentum().E() > tmp_leadingPionPlusE) {
500  tmp_leadingPionPlusE = particle->Momentum().E();
501  MC_leading_PionPlusID = particle->TrackId();
502  MC_leading_PionPlusP =
503  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
504  pow(particle->Momentum().Pz(), 2));
505  MCpion_plus = particle;
506  }
507  }
508  else if (particle->PdgCode() == -211) {
509  if (particle->Momentum().E() > tmp_leadingPionMinusE) {
510  tmp_leadingPionMinusE = particle->Momentum().E();
511  MC_leading_PionMinusID = particle->TrackId();
512  MC_leading_PionMinusP =
513  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
514  pow(particle->Momentum().Pz(), 2));
515  MCpion_minus = particle;
516  }
517  }
518  i++; //paticle index
519  }
520 
521  //=======================================================================================
522  //add Nucleon decay stuff and particle cannon
523  //=======================================================================================
524  if (!fisNeutrinoInt) {
525  if (particle->Mother() == 0) {
526  const TLorentzVector& positionStart = particle->Position(0);
527  positionStart.GetXYZT(MC_vertex);
528  }
529  if (particle->PdgCode() == 321) { //save primary Kaon
530  MC_kaonID = particle->TrackId();
531  MC_kaonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
532  pow(particle->Momentum().Pz(), 2));
533  MCkaon = particle;
534  }
535  else if (particle->PdgCode() == fLeptonPDGcode) { // Particle cannon muon
536  MC_leptonID = particle->TrackId();
537  MC_leptonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
538  pow(particle->Momentum().Pz(), 2));
539  MClepton = particle;
540  }
541  else if (particle->PdgCode() == 2212) {
542  if (particle->Momentum().E() > tmp_leadingProtonE) {
543  tmp_leadingProtonE = particle->Momentum().E();
544  MC_leading_protonID = particle->TrackId();
545  MC_leading_ProtonP =
546  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
547  pow(particle->Momentum().Pz(), 2));
548  MCproton = particle;
549  }
550  }
551  else if (particle->PdgCode() == 211) {
552  if (particle->Momentum().E() > tmp_leadingPionPlusE) {
553  tmp_leadingPionPlusE = particle->Momentum().E();
554  MC_leading_PionPlusID = particle->TrackId();
555  MC_leading_PionPlusP =
556  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
557  pow(particle->Momentum().Pz(), 2));
558  MCpion_plus = particle;
559  }
560  }
561  else if (particle->PdgCode() == -211) {
562  if (particle->Momentum().E() > tmp_leadingPionMinusE) {
563  tmp_leadingPionMinusE = particle->Momentum().E();
564  MC_leading_PionMinusID = particle->TrackId();
565  MC_leading_PionMinusP =
566  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
567  pow(particle->Momentum().Pz(), 2));
568  MCpion_minus = particle;
569  }
570  }
571  else if (particle->Process() == "Decay" &&
572  particle->PdgCode() == -11) { // michel electron from muon decay
573  MC_michelID = particle->TrackId();
574  MC_michelP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
575  pow(particle->Momentum().Pz(), 2));
576  MCmichel = particle;
577  }
578  else if (TMath::Abs(particle->PdgCode() == 321)) { //save primary Kaon
579  MC_kaonID = particle->TrackId();
580  MC_kaonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
581  pow(particle->Momentum().Pz(), 2));
582  MCkaon = particle;
583  }
584  }
585  }
586  //===================================================================
587  //Saving denominator histograms
588  //===================================================================
589  if (not insideFV(MC_vertex)) return;
590  double Pv =
591  sqrt(pow(MC_incoming_P[0], 2) + pow(MC_incoming_P[1], 2) + pow(MC_incoming_P[2], 2));
592  double theta_mu = acos((MC_incoming_P[0] * MC_lepton_startMomentum[0] +
593  MC_incoming_P[1] * MC_lepton_startMomentum[1] +
594  MC_incoming_P[2] * MC_lepton_startMomentum[2]) /
595  (Pv * MC_leptonP));
596  theta_mu *= (180.0 / 3.14159);
597  double truth_lengthLepton = truthLength(clockData, detProp, MClepton);
598  double proton_length = truthLength(clockData, detProp, MCproton);
599  double pion_plus_length = truthLength(clockData, detProp, MCpion_plus);
600  double pion_minus_length = truthLength(clockData, detProp, MCpion_minus);
601  double kaonLength = truthLength(clockData, detProp, MCkaon);
602  double michelLength = truthLength(clockData, detProp, MCmichel);
603 
604  // save CC events within the fiducial volume with the favorite neutrino
605  // flavor
606  if (MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE)) {
607  if (MClepton) {
608  h_Ev_den->Fill(MC_incoming_P[3]);
609  h_Pmu_den->Fill(MC_leptonP);
610  h_theta_den->Fill(theta_mu);
611  h_muon_length->Fill(truth_lengthLepton);
612  }
613  if (MCproton) {
614  h_Pproton_den->Fill(MC_leading_ProtonP);
615  h_proton_length->Fill(proton_length);
616  }
617  if (MCpion_plus) {
618  h_Ppion_plus_den->Fill(MC_leading_PionPlusP);
619  h_pionp_length->Fill(pion_plus_length);
620  }
621  if (MCpion_minus) {
622  h_Ppion_minus_den->Fill(MC_leading_PionMinusP);
623  h_pionm_length->Fill(pion_minus_length);
624  }
625  if (MCkaon) {
626  h_Pkaon_den->Fill(MC_kaonP);
627  h_kaon_length->Fill(kaonLength);
628  }
629  }
630 
631  //save events for Nucleon decay and particle cannon
632  if (!fisNeutrinoInt) {
633  if (MClepton) {
634  h_Pmu_den->Fill(MC_leptonP);
635  h_muon_length->Fill(truth_lengthLepton);
636  }
637  if (MCkaon) {
638  h_Pkaon_den->Fill(MC_kaonP);
639  h_kaon_length->Fill(kaonLength);
640  }
641  if (MCproton) {
642  h_Pproton_den->Fill(MC_leading_ProtonP);
643  h_proton_length->Fill(proton_length);
644  }
645  if (MCpion_plus) {
646  h_Ppion_plus_den->Fill(MC_leading_PionPlusP);
647  h_pionp_length->Fill(pion_plus_length);
648  }
649  if (MCpion_minus) {
650  h_Ppion_minus_den->Fill(MC_leading_PionMinusP);
651  h_pionm_length->Fill(pion_minus_length);
652  }
653  if (MCmichel) {
654  h_Pmichel_e_den->Fill(MC_michelP);
655  h_michel_length->Fill(michelLength);
656  }
657  }
658 
659  //========================================================================
660  // Reco stuff, once we have selected a MC Particle let's find out if there is a track associated
661  //========================================================================
662  art::Handle<std::vector<recob::Track>> trackListHandle;
663  if (!event.getByLabel(fTrackModuleLabel, trackListHandle)) return;
664  std::vector<art::Ptr<recob::Track>> tracklist;
665  art::fill_ptr_vector(tracklist, trackListHandle);
666  int n_recoTrack = tracklist.size();
667 
668  art::FindManyP<recob::Hit> track_hits(trackListHandle, event, fTrackModuleLabel);
669  if (n_recoTrack == 0) {
670  MF_LOG_DEBUG("NeutrinoTrackingEff") << "There are no reco tracks... bye";
671  return;
672  }
673  MF_LOG_DEBUG("NeutrinoTrackingEff") << "Found this many reco tracks " << n_recoTrack;
674 
675  double Efrac_lepton = 0.0;
676  double Ecomplet_lepton = 0.0;
677  double Efrac_proton = 0.0;
678  double Ecomplet_proton = 0.0;
679  double Efrac_pionplus = 0.0;
680  double Ecomplet_pionplus = 0.0;
681  double Efrac_pionminus = 0.0;
682  double Ecomplet_pionminus = 0.0;
683  double Efrac_kaon = 0.0;
684  double Ecomplet_kaon = 0.0;
685  double Efrac_michel = 0.0;
686  double Ecomplet_michel = 0.0;
687  double trackLength_lepton = 0.0;
688  double trackLength_proton = 0.0;
689  double trackLength_pion_plus = 0.0;
690  double trackLength_pion_minus = 0.0;
691  double trackLength_kaon = 0.0;
692  double trackLength_michel = 0.0;
693  const simb::MCParticle* MClepton_reco = nullptr;
694  const simb::MCParticle* MCproton_reco = nullptr;
695  const simb::MCParticle* MCpion_plus_reco = nullptr;
696  const simb::MCParticle* MCpion_minus_reco = nullptr;
697  const simb::MCParticle* MCkaon_reco = nullptr;
698  const simb::MCParticle* MCmichel_reco = nullptr;
699 
700  std::vector<art::Ptr<recob::Hit>> tmp_all_trackHits = track_hits.at(0);
701  std::vector<art::Ptr<recob::Hit>> all_hits;
702  art::Handle<std::vector<recob::Hit>> hithandle;
703  auto const pd = event.getProductDescription(tmp_all_trackHits[0].id());
704  if (pd && event.getByLabel(pd->inputTag(), hithandle)) {
705  art::fill_ptr_vector(all_hits, hithandle);
706  }
707 
708  for (int i = 0; i < n_recoTrack; i++) {
709  art::Ptr<recob::Track> track = tracklist[i];
710  std::vector<art::Ptr<recob::Hit>> all_trackHits = track_hits.at(i);
711  double tmpEfrac = 0;
712  double tmpEcomplet = 0;
713  const simb::MCParticle* particle;
714  truthMatcher(clockData, all_hits, all_trackHits, particle, tmpEfrac, tmpEcomplet);
715  if (!particle) continue;
716  if ((particle->PdgCode() == fLeptonPDGcode) && (particle->TrackId() == MC_leptonID)) {
717  // save the best track ... based on completeness if there is more than
718  // one track if( tmpEfrac > Efrac_lepton ){ ///this was base on purity
719  if (tmpEcomplet > Ecomplet_lepton) {
720  Ecomplet_lepton = tmpEcomplet;
721  Efrac_lepton = tmpEfrac;
722  trackLength_lepton = track->Length();
723  MClepton_reco = particle;
724  }
725  }
726  else if ((particle->PdgCode() == 2212) && (particle->TrackId() == MC_leading_protonID)) {
727  //save the best track ... based on completeness if there is more than one track
728  if (tmpEcomplet > Ecomplet_proton) {
729  Ecomplet_proton = tmpEcomplet;
730  Efrac_proton = tmpEfrac;
731  trackLength_proton = track->Length();
732  MCproton_reco = particle;
733  }
734  }
735  else if ((particle->PdgCode() == 211) && (particle->TrackId() == MC_leading_PionPlusID)) {
736  //save the best track ... based on completeness if there is more than one track
737  if (tmpEcomplet > Ecomplet_pionplus) {
738  Ecomplet_pionplus = tmpEcomplet;
739  Efrac_pionplus = tmpEfrac;
740  trackLength_pion_plus = track->Length();
741  MCpion_plus_reco = particle;
742  }
743  }
744  else if ((particle->PdgCode() == -211) && (particle->TrackId() == MC_leading_PionMinusID)) {
745  //save the best track ... based on completeness if there is more than one track
746  if (tmpEcomplet > Ecomplet_pionminus) {
747  Ecomplet_pionminus = tmpEcomplet;
748  Efrac_pionminus = tmpEfrac;
749  trackLength_pion_minus = track->Length();
750  MCpion_minus_reco = particle;
751  }
752  }
753  //kaon from nucleon decay
754  else if ((TMath::Abs(particle->PdgCode()) == 321) && (particle->TrackId() == MC_kaonID)) {
755  //save the best track ... based on completeness if there is more than one track
756  if (tmpEcomplet > Ecomplet_kaon) {
757  Ecomplet_kaon = tmpEcomplet;
758  Efrac_kaon = tmpEfrac;
759  trackLength_kaon = track->Length();
760  MCkaon_reco = particle;
761  }
762  }
763  //michel from nucleon decay
764  else if ((particle->PdgCode() == -11) && (particle->TrackId() == MC_michelID)) {
765  //save the best track ... based on completeness if there is more than one track
766  if (tmpEcomplet > Ecomplet_michel) {
767  Ecomplet_michel = tmpEcomplet;
768  Efrac_michel = tmpEfrac;
769  trackLength_michel = track->Length();
770  MCmichel_reco = particle;
771  }
772  }
773  }
774 
775  double Reco_LengthRes = truth_lengthLepton - trackLength_lepton;
776  double Reco_LengthResProton = proton_length - trackLength_proton;
777  double Reco_LengthResPionPlus = pion_plus_length - trackLength_pion_plus;
778  double Reco_LengthResPionMinus = pion_minus_length - trackLength_pion_minus;
779 
780  if (MClepton_reco && MClepton) {
781  if (MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE)) {
782  h_Pmu_num->Fill(MC_leptonP);
783  h_Ev_num->Fill(MC_incoming_P[3]);
784  h_theta_num->Fill(theta_mu);
785  h_Efrac_lepton->Fill(Efrac_lepton);
786  h_Ecomplet_lepton->Fill(Ecomplet_lepton);
787  h_trackRes_lepton->Fill(Reco_LengthRes);
788  h_muonwtrk_length->Fill(truth_lengthLepton);
789  }
790  }
791  if (MCproton_reco && MCproton) {
792  if (MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE)) {
793  h_Pproton_num->Fill(MC_leading_ProtonP);
794  h_Efrac_proton->Fill(Efrac_proton);
795  h_Ecomplet_proton->Fill(Ecomplet_proton);
796  h_trackRes_proton->Fill(Reco_LengthResProton);
797  h_protonwtrk_length->Fill(proton_length);
798  }
799  }
800  if (MCpion_plus_reco && MCpion_plus) {
801  if (MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE)) {
802  h_Ppion_plus_num->Fill(MC_leading_PionPlusP);
803  h_Efrac_pion_plus->Fill(Efrac_pionplus);
804  h_Ecomplet_pion_plus->Fill(Ecomplet_pionplus);
805  h_trackRes_pion_plus->Fill(Reco_LengthResPionPlus);
806  h_pionpwtrk_length->Fill(pion_plus_length);
807  }
808  }
809  if (MCpion_minus_reco && MCpion_minus) {
810  if (MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE)) {
811  h_Ppion_minus_num->Fill(MC_leading_PionMinusP);
812  h_Efrac_pion_minus->Fill(Efrac_pionminus);
813  h_Ecomplet_pion_minus->Fill(Ecomplet_pionminus);
814  h_trackRes_pion_minus->Fill(Reco_LengthResPionMinus);
815  h_pionmwtrk_length->Fill(pion_minus_length);
816  }
817  }
818  if (MCkaon_reco && MCkaon) {
819  if (MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE)) {
820  h_Pkaon_num->Fill(MC_kaonP);
821  h_Efrac_kaon->Fill(Efrac_kaon);
822  h_Ecomplet_kaon->Fill(Ecomplet_kaon);
823  h_trackRes_kaon->Fill(kaonLength - trackLength_kaon);
824  h_kaonwtrk_length->Fill(kaonLength);
825  }
826  }
827  //Non neutrino events
828  //=========================================================
829  if (!fisNeutrinoInt) {
830  if (MClepton_reco && MClepton) {
831  h_Pmu_num->Fill(MC_leptonP);
832  h_Efrac_lepton->Fill(Efrac_lepton);
833  h_Ecomplet_lepton->Fill(Ecomplet_lepton);
834  h_trackRes_lepton->Fill(Reco_LengthRes);
835  h_muonwtrk_length->Fill(truth_lengthLepton);
836  }
837  if (MCkaon_reco && MCkaon) {
838  h_Pkaon_num->Fill(MC_kaonP);
839  h_Efrac_kaon->Fill(Efrac_kaon);
840  h_Ecomplet_kaon->Fill(Ecomplet_kaon);
841  h_trackRes_kaon->Fill(kaonLength - trackLength_kaon);
842  h_kaonwtrk_length->Fill(kaonLength);
843  }
844  if (MCproton_reco && MCproton) {
845  h_Pproton_num->Fill(MC_leading_ProtonP);
846  h_Efrac_proton->Fill(Efrac_proton);
847  h_Ecomplet_proton->Fill(Ecomplet_proton);
848  h_trackRes_proton->Fill(Reco_LengthResProton);
849  h_protonwtrk_length->Fill(proton_length);
850  }
851  if (MCpion_plus_reco && MCpion_plus) {
852  h_Ppion_plus_num->Fill(MC_leading_PionPlusP);
853  h_Efrac_pion_plus->Fill(Efrac_pionplus);
854  h_Ecomplet_pion_plus->Fill(Ecomplet_pionplus);
855  h_trackRes_pion_plus->Fill(Reco_LengthResPionPlus);
856  h_pionpwtrk_length->Fill(pion_plus_length);
857  }
858  if (MCpion_minus_reco && MCpion_minus) {
859  h_Ppion_minus_num->Fill(MC_leading_PionMinusP);
860  h_Efrac_pion_minus->Fill(Efrac_pionminus);
861  h_Ecomplet_pion_minus->Fill(Ecomplet_pionminus);
862  h_trackRes_pion_minus->Fill(Reco_LengthResPionMinus);
863  h_pionmwtrk_length->Fill(pion_minus_length);
864  }
865  if (MCmichel_reco && MCmichel) {
866  h_Pmichel_e_num->Fill(MC_michelP);
867  h_Efrac_michel->Fill(Efrac_michel);
868  h_Ecomplet_michel->Fill(Ecomplet_michel);
869  h_trackRes_michel->Fill(michelLength - trackLength_michel);
870  h_michelwtrk_length->Fill(michelLength);
871  }
872  }
873  }
874  //========================================================================
875  void NeutrinoTrackingEff::truthMatcher(detinfo::DetectorClocksData const& clockData,
876  std::vector<art::Ptr<recob::Hit>> all_hits,
877  std::vector<art::Ptr<recob::Hit>> track_hits,
878  const simb::MCParticle*& MCparticle,
879  double& Efrac,
880  double& Ecomplet)
881  {
882  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
883  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
884  std::map<int, double> trkID_E;
885  for (size_t j = 0; j < track_hits.size(); ++j) {
886  art::Ptr<recob::Hit> hit = track_hits[j];
887  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
888  for (size_t k = 0; k < TrackIDs.size(); k++) {
889  trkID_E[TrackIDs[k].trackID] += TrackIDs[k].energy;
890  }
891  }
892  double max_E = -999.0;
893  double total_E = 0.0;
894  int TrackID = -999;
895  double partial_E =
896  0.0; // amount of energy deposited by the particle that deposited more energy...
897 
898  // If the collection of hits have more than one particle associate
899  // save the particle w/ the highest energy deposition since we are
900  // looking for muons/pions/protons this should be enough
901  if (!trkID_E.size()) {
902  MCparticle = 0;
903  return; //Ghost track???
904  }
905  for (std::map<int, double>::iterator ii = trkID_E.begin(); ii != trkID_E.end(); ++ii) {
906  total_E += ii->second;
907  if ((ii->second) > max_E) {
908  partial_E = ii->second;
909  max_E = ii->second;
910  TrackID = ii->first;
911  }
912  }
913  MCparticle = pi_serv->TrackIdToParticle_P(TrackID);
914 
915  // In the current simulation, we do not save EM Shower daughters
916  // in GEANT. But we do save the energy deposition in TrackIDEs. If
917  // the energy deposition is from a particle that is the daughter
918  // of an EM particle, the negative of the parent track ID is saved
919  // in TrackIDE for the daughter particle we don't want to track
920  // gammas or any other EM activity
921  if (TrackID < 0) return;
922 
923  Efrac = (partial_E) / total_E;
924 
925  // Completeness
926  double totenergy = 0;
927  for (size_t k = 0; k < all_hits.size(); ++k) {
928  art::Ptr<recob::Hit> hit = all_hits[k];
929  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
930  for (size_t l = 0; l < TrackIDs.size(); ++l) {
931  if (TrackIDs[l].trackID == TrackID) totenergy += TrackIDs[l].energy;
932  }
933  }
934  Ecomplet = partial_E / totenergy;
935  }
936  //========================================================================
937  double NeutrinoTrackingEff::truthLength(const detinfo::DetectorClocksData& clockData,
938  const detinfo::DetectorPropertiesData& detProp,
939  const simb::MCParticle* MCparticle)
940  {
941  // Calculate the truth length considering only the part that is
942  // inside the TPC Base on a peace of code from
943  // dune/TrackingAna/TrackingEfficiency_module.cc
944 
945  if (!MCparticle) return -999.0;
946  int numberTrajectoryPoints = MCparticle->NumberTrajectoryPoints();
947  std::vector<double> TPCLengthHits(numberTrajectoryPoints, 0.0);
948  int FirstHit = 0, LastHit = 0;
949  double TPCLength = 0.0;
950  bool BeenInVolume = false;
951 
952  double const WindowSize = detProp.NumberTimeSamples() * clockData.TPCClock().TickPeriod() * 1e3;
953 
954  for (unsigned int MCHit = 0; MCHit < TPCLengthHits.size(); ++MCHit) {
955  if (MCHit != 0)
956  TPCLengthHits[MCHit] = sqrt(pow((MCparticle->Vx(MCHit - 1) - MCparticle->Vx(MCHit)), 2) +
957  pow((MCparticle->Vy(MCHit - 1) - MCparticle->Vy(MCHit)), 2) +
958  pow((MCparticle->Vz(MCHit - 1) - MCparticle->Vz(MCHit)), 2));
959  auto const position = geo::vect::toPoint(MCparticle->Position(MCHit).Vect());
960  geo::TPCID tpcid = geom->FindTPCAtPosition(position);
961  if (tpcid.isValid) {
962  // -- Check if hit is within drift window...
963  geo::TPCGeo const& tpc = geom->TPC(tpcid);
964  double XPlanePosition = tpc.Plane(0).GetCenter().X();
965  double DriftTimeCorrection = fabs(position.X() - XPlanePosition) / fDriftVelocity;
966  double TimeAtPlane = MCparticle->T() + DriftTimeCorrection;
967 
968  if (TimeAtPlane < trigger_offset(clockData) ||
969  TimeAtPlane > trigger_offset(clockData) + WindowSize)
970  continue;
971  LastHit = MCHit;
972  if (!BeenInVolume) {
973  BeenInVolume = true;
974  FirstHit = MCHit;
975  }
976  }
977  }
978  for (int Hit = FirstHit + 1; Hit <= LastHit; ++Hit)
979  TPCLength += TPCLengthHits[Hit];
980  return TPCLength;
981  }
982  //========================================================================
983  bool NeutrinoTrackingEff::insideFV(double vertex[4])
984  {
985  double const x = vertex[0];
986  double const y = vertex[1];
987  double const z = vertex[2];
988 
989  return x > fFidVolXmin && x < fFidVolXmax && y > fFidVolYmin && y < fFidVolYmax &&
990  z > fFidVolZmin && z < fFidVolZmax;
991  }
992  //========================================================================
993  void NeutrinoTrackingEff::doEfficiencies()
994  {
995  art::ServiceHandle<art::TFileService const> tfs;
996 
997  if (TEfficiency::CheckConsistency(*h_Ev_num, *h_Ev_den)) {
998  h_Eff_Ev = tfs->make<TEfficiency>(*h_Ev_num, *h_Ev_den);
999  TGraphAsymmErrors* grEff_Ev = h_Eff_Ev->CreateGraph();
1000  grEff_Ev->Write("grEff_Ev");
1001  h_Eff_Ev->Write("h_Eff_Ev");
1002  }
1003  if (TEfficiency::CheckConsistency(*h_Pmu_num, *h_Pmu_den)) {
1004  h_Eff_Pmu = tfs->make<TEfficiency>(*h_Pmu_num, *h_Pmu_den);
1005  TGraphAsymmErrors* grEff_Pmu = h_Eff_Pmu->CreateGraph();
1006  grEff_Pmu->Write("grEff_Pmu");
1007  h_Eff_Pmu->Write("h_Eff_Pmu");
1008  }
1009  if (TEfficiency::CheckConsistency(*h_theta_num, *h_theta_den)) {
1010  h_Eff_theta = tfs->make<TEfficiency>(*h_theta_num, *h_theta_den);
1011  TGraphAsymmErrors* grEff_theta = h_Eff_theta->CreateGraph();
1012  grEff_theta->Write("grEff_theta");
1013  h_Eff_theta->Write("h_Eff_theta");
1014  }
1015  if (TEfficiency::CheckConsistency(*h_Pproton_num, *h_Pproton_den)) {
1016  h_Eff_Pproton = tfs->make<TEfficiency>(*h_Pproton_num, *h_Pproton_den);
1017  TGraphAsymmErrors* grEff_Pproton = h_Eff_Pproton->CreateGraph();
1018  grEff_Pproton->Write("grEff_Pproton");
1019  h_Eff_Pproton->Write("h_Eff_Pproton");
1020  }
1021  if (TEfficiency::CheckConsistency(*h_Ppion_plus_num, *h_Ppion_plus_den)) {
1022  h_Eff_Ppion_plus = tfs->make<TEfficiency>(*h_Ppion_plus_num, *h_Ppion_plus_den);
1023  TGraphAsymmErrors* grEff_Ppion_plus = h_Eff_Ppion_plus->CreateGraph();
1024  grEff_Ppion_plus->Write("grEff_Ppion_plus");
1025  h_Eff_Ppion_plus->Write("h_Eff_Ppion_plus");
1026  }
1027  if (TEfficiency::CheckConsistency(*h_Ppion_minus_num, *h_Ppion_minus_den)) {
1028  h_Eff_Ppion_minus = tfs->make<TEfficiency>(*h_Ppion_minus_num, *h_Ppion_minus_den);
1029  TGraphAsymmErrors* grEff_Ppion_minus = h_Eff_Ppion_minus->CreateGraph();
1030  grEff_Ppion_minus->Write("grEff_Ppion_minus");
1031  h_Eff_Ppion_minus->Write("h_Eff_Ppion_minus");
1032  }
1033  if (TEfficiency::CheckConsistency(*h_muonwtrk_length, *h_muon_length)) {
1034  h_Eff_Lmuon = tfs->make<TEfficiency>(*h_muonwtrk_length, *h_muon_length);
1035  TGraphAsymmErrors* grEff_Lmuon = h_Eff_Lmuon->CreateGraph();
1036  grEff_Lmuon->Write("grEff_Lmuon");
1037  h_Eff_Lmuon->Write("h_Eff_Lmuon");
1038  }
1039  if (TEfficiency::CheckConsistency(*h_protonwtrk_length, *h_proton_length)) {
1040  h_Eff_Lproton = tfs->make<TEfficiency>(*h_protonwtrk_length, *h_proton_length);
1041  TGraphAsymmErrors* grEff_Lproton = h_Eff_Lproton->CreateGraph();
1042  grEff_Lproton->Write("grEff_Lproton");
1043  h_Eff_Lproton->Write("h_Eff_Lproton");
1044  }
1045  if (TEfficiency::CheckConsistency(*h_pionpwtrk_length, *h_pionp_length)) {
1046  h_Eff_Lpion_plus = tfs->make<TEfficiency>(*h_pionpwtrk_length, *h_pionp_length);
1047  TGraphAsymmErrors* grEff_Lpion_plus = h_Eff_Lpion_plus->CreateGraph();
1048  grEff_Lpion_plus->Write("grEff_Lpion_plus");
1049  h_Eff_Lpion_plus->Write("h_Eff_Lpion_plus");
1050  }
1051  if (TEfficiency::CheckConsistency(*h_pionpwtrk_length, *h_pionp_length)) {
1052  h_Eff_Lpion_minus = tfs->make<TEfficiency>(*h_pionmwtrk_length, *h_pionm_length);
1053  TGraphAsymmErrors* grEff_Lpion_minus = h_Eff_Lpion_minus->CreateGraph();
1054  grEff_Lpion_minus->Write("grEff_Lpion_minus");
1055  h_Eff_Lpion_minus->Write("h_Eff_Lpion_minus");
1056  }
1057  if (TEfficiency::CheckConsistency(*h_Pkaon_num, *h_Pkaon_den)) {
1058  h_Eff_Pkaon = tfs->make<TEfficiency>(*h_Pkaon_num, *h_Pkaon_den);
1059  TGraphAsymmErrors* grEff_Pkaon = h_Eff_Pkaon->CreateGraph();
1060  grEff_Pkaon->Write("grEff_Pkaon");
1061  h_Eff_Pkaon->Write("h_Eff_Pkaon");
1062  }
1063  if (TEfficiency::CheckConsistency(*h_kaonwtrk_length, *h_kaon_length)) {
1064  h_Eff_Lkaon = tfs->make<TEfficiency>(*h_kaonwtrk_length, *h_kaon_length);
1065  TGraphAsymmErrors* grEff_Lkaon = h_Eff_Lkaon->CreateGraph();
1066  grEff_Lkaon->Write("grEff_Lkaon");
1067  h_Eff_Lkaon->Write("h_Eff_Lkaon");
1068  }
1069  if (TEfficiency::CheckConsistency(*h_Pmichel_e_num, *h_Pmichel_e_den)) {
1070  h_Eff_Pmichel = tfs->make<TEfficiency>(*h_Pmichel_e_num, *h_Pmichel_e_den);
1071  TGraphAsymmErrors* grEff_Pmichel = h_Eff_Pmichel->CreateGraph();
1072  grEff_Pmichel->Write("grEff_Pmichel");
1073  h_Eff_Pmichel->Write("h_Eff_Pmichel");
1074  }
1075  if (TEfficiency::CheckConsistency(*h_michelwtrk_length, *h_michel_length)) {
1076  h_Eff_Lmichel = tfs->make<TEfficiency>(*h_michelwtrk_length, *h_michel_length);
1077  TGraphAsymmErrors* grEff_Lmichel = h_Eff_Lmichel->CreateGraph();
1078  grEff_Lmichel->Write("grEff_Lmichel");
1079  h_Eff_Lmichel->Write("h_Eff_Lmichel");
1080  }
1081  }
1082  //========================================================================
1083  DEFINE_ART_MODULE(NeutrinoTrackingEff)
1084 
1085 }
x
Float_t x
Definition: compare.C:6
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
DUNE::NeutrinoTrackingEff::geom
art::ServiceHandle< geo::Geometry const > geom
Definition: NeutrinoTrackingEff_module.cc:182
simb::MCParticle::NumberTrajectoryPoints
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:219
simb::MCParticle::Position
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:220
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:213
simb::MCNeutrino::CCNC
int CCNC() const
Definition: MCNeutrino.h:148
simb::MCTruth::GetNeutrino
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
ServiceUtil.h
Utilities related to art service access.
cheat::BackTrackerService::HitToTrackIDEs
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Definition: BackTrackerService.cc:163
geo::PlaneGeo::GetCenter
Point_t const & GetCenter() const
Returns the centre of the wire plane in world coordinates [cm].
Definition: PlaneGeo.h:435
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
TrackIDs
std::vector< TrackID > TrackIDs
Definition: ClusteringValidation_module.cc:60
cheat::ParticleInventoryService::TrackIdToParticle_P
const simb::MCParticle * TrackIdToParticle_P(int id) const
Definition: ParticleInventoryService.cc:154
simb::MCParticle::Mother
int Mother() const
Definition: MCParticle.h:214
y
Float_t y
Definition: compare.C:6
simb::MCNeutrino::Nu
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
sim::ParticleList::const_iterator
list_type::const_iterator const_iterator
Definition: ParticleList.h:132
z
Double_t z
Definition: plot.C:276
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:36
geo::vect::toPoint
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
Definition: geo_vectors_utils.h:1584
std
STL namespace.
art::Event::isRealData
bool isRealData() const
Definition: Event.cc:53
simb::MCParticle::Process
std::string Process() const
Definition: MCParticle.h:216
MCParticle.h
Particle class.
DUNE::NeutrinoTrackingEff::truthMatcher
void truthMatcher(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> all_hits, std::vector< art::Ptr< recob::Hit >> track_hits, const simb::MCParticle *&MCparticle, double &Efrac, double &Ecomplet)
Definition: NeutrinoTrackingEff_module.cc:875
art::Run
Definition: Run.h:37
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:211
DUNE::NeutrinoTrackingEff::processEff
void processEff(const art::Event &evt)
Definition: NeutrinoTrackingEff_module.cc:428
detinfo::ElecClock::TickPeriod
constexpr double TickPeriod() const noexcept
A single tick period in microseconds.
Definition: ElecClock.h:310
geo::GeometryCore::TPC
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
geo::GeometryCore::FindTPCAtPosition
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
Definition: GeometryCore.cxx:326
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:207
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
breakpoints::beginJob
void beginJob()
Definition: Breakpoints.cc:14
DUNE::NeutrinoTrackingEff::analyze
void analyze(const art::Event &evt)
Definition: NeutrinoTrackingEff_module.cc:421
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:314
Event.h
simb::MCParticle::T
double T(const int i=0) const
Definition: MCParticle.h:225
sim::ParticleList::begin
iterator begin()
Definition: ParticleList.h:305
detinfo::DetectorClocksData::TPCClock
ElecClock const & TPCClock() const noexcept
Borrow a const TPC clock with time set to Trigger time [us].
Definition: DetectorClocksData.h:375
Track.h
Provides recob::Track data product.
TrackID
TrackID
Definition: ClusteringValidation_module.cc:58
DUNE::NeutrinoTrackingEff::beginRun
void beginRun(const art::Run &run)
Definition: NeutrinoTrackingEff_module.cc:416
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
DUNE::NeutrinoTrackingEff::truthLength
double truthLength(const detinfo::DetectorClocksData &clockData, detinfo::DetectorPropertiesData const &detProp, const simb::MCParticle *MCparticle)
Definition: NeutrinoTrackingEff_module.cc:937
geo_types.h
Definition of data types for geometry description.
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
cheat::ParticleInventoryService::ParticleList
const sim::ParticleList & ParticleList() const
Definition: ParticleInventoryService.cc:138
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:222
art::ProductRetriever::getByLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: ProductRetriever.h:444
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:282
simb::MCParticle::Momentum
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:221
MF_LOG_DEBUG
#define MF_LOG_DEBUG(id)
Definition: MessageLogger.h:294
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:224
detinfo::trigger_offset
int trigger_offset(DetectorClocksData const &data)
Definition: DetectorClocksData.h:565
geo::TPCGeo::Plane
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:252
simb::MCTruth::NeutrinoSet
bool NeutrinoSet() const
Definition: MCTruth.h:78
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
e
Float_t e
Definition: plot.C:35
simb::MCNeutrino
Event generator information.
Definition: MCNeutrino.h:18
geo
Namespace collecting geometry-related classes utilities.
Definition: BeamFlashTrackMatchTaggerAlg.h:20
track
Float_t track
Definition: plot.C:35
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:223
Geometry.h
art framework interface to geometry description
event
Event finding and building.
Definition: EventCheater_module.cc:32
TPCGeo.h
Encapsulate the construction of a single detector plane.
geo::CryostatID
The data type to uniquely identify a cryostat.
Definition: geo_types.h:192
vertex
vertex reconstruction
Definition: AggregateVertexAna_module.cc:32