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