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