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