NeutrinoShowerEff_module.cc
Go to the documentation of this file.
1 // LArSoft includes
2 #include "lardata/ArtDataHelper/MVAReader.h"
3 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
4 #include "lardataobj/RecoBase/Cluster.h"
5 #include "lardataobj/RecoBase/PFParticle.h"
6 #include "lardataobj/RecoBase/Shower.h"
7 #include "larsim/MCCheater/BackTrackerService.h"
8 #include "larsim/MCCheater/ParticleInventoryService.h"
9 #include "nusimdata/SimulationBase/MCParticle.h"
10 #include "nusimdata/SimulationBase/MCTruth.h"
11 
12 // Framework includes
13 #include "art/Framework/Core/EDAnalyzer.h"
14 #include "art/Framework/Core/ModuleMacros.h"
15 #include "art/Framework/Principal/Event.h"
16 #include "art/Framework/Principal/Handle.h"
17 #include "art/Framework/Services/Registry/ServiceHandle.h"
18 #include "art_root_io/TFileService.h"
19 #include "canvas/Persistency/Common/FindManyP.h"
20 #include "fhiclcpp/ParameterSet.h"
21 #include "messagefacility/MessageLogger/MessageLogger.h"
22 
23 // ROOT includes
24 #include "TEfficiency.h"
25 #include "TGraphAsymmErrors.h"
26 #include "TH1.h"
27 #include "TTree.h"
28 
29 #define MAX_SHOWERS 1000
30 using namespace std;
31 
32 //========================================================================
33 
34 namespace DUNE {
35 
36  class NeutrinoShowerEff : public art::EDAnalyzer {
37  public:
38  explicit NeutrinoShowerEff(fhicl::ParameterSet const& pset);
39 
40  private:
41  void beginJob() override;
42  void endJob() override;
43  void beginRun(const art::Run& run) override;
44  void analyze(const art::Event& evt) override;
45 
46  void processEff(detinfo::DetectorClocksData const& clockData,
47  const art::Event& evt,
48  bool& isFiducial);
49  void truthMatcher(detinfo::DetectorClocksData const& clockData,
50  std::vector<art::Ptr<recob::Hit>> all_hits,
51  std::vector<art::Ptr<recob::Hit>> shower_hits,
52  const simb::MCParticle*& MCparticle,
53  double& Efrac,
54  double& Ecomplet);
55  template <size_t N>
56  void checkCNNtrkshw(detinfo::DetectorClocksData const& clockData,
57  const art::Event& evt,
58  std::vector<art::Ptr<recob::Hit>> all_hits);
59  bool insideFV(double vertex[4]);
60  void doEfficiencies();
61  void reset();
62 
63  // the parameters we'll read from the .fcl
64 
65  art::InputTag fMCTruthModuleLabel;
66  art::InputTag fHitModuleLabel;
67  art::InputTag fShowerModuleLabel;
68  art::InputTag fCNNEMModuleLabel;
69 
70  int fNeutrinoPDGcode;
71  int fLeptonPDGcode;
72  double fMaxNeutrinoE;
73  bool fSaveMCTree;
74  double fMaxEfrac;
75  double fMinCompleteness;
76 
77  TTree* fEventTree;
78  // TTree *fHitsTree;
79 
80  TH1D* h_Ev_den;
81  TH1D* h_Ev_num;
82  TH1D* h_Ev_num_dEdx;
83 
84  TH1D* h_Ee_den;
85  TH1D* h_Ee_num;
86  TH1D* h_Ee_num_dEdx;
87 
88  TH1D* h_Pe_den;
89  TH1D* h_Pe_num;
90 
91  TH1D* h_theta_den;
92  TH1D* h_theta_num;
93 
94  TH1D* h_Efrac_shContamination;
95  TH1D* h_Efrac_shPurity;
96  TH1D* h_Ecomplet_lepton;
97 
98  TH1D* h_Efrac_NueCCPurity; //Signal
99  TH1D* h_Ecomplet_NueCC;
100  TH1D* h_HighestHitsProducedParticlePDG_NueCC;
101  TH1D* h_HighestHitsProducedParticlePDG_bkg;
102  TH1D* h_Efrac_bkgPurity; //Background
103  TH1D* h_Ecomplet_bkg;
104 
105  TEfficiency* h_Eff_Ev = 0;
106  TEfficiency* h_Eff_Ev_dEdx = 0;
107  TEfficiency* h_Eff_Ee = 0;
108  TEfficiency* h_Eff_Ee_dEdx = 0;
109  TEfficiency* h_Eff_Pe = 0;
110  TEfficiency* h_Eff_theta = 0;
111 
112  //Electron shower Best plane number
113  TH1D* h_esh_bestplane_NueCC;
114  TH1D* h_esh_bestplane_NC;
115  TH1D* h_dEdX_electronorpositron_NueCC;
116  TH1D* h_dEdX_electronorpositron_NC;
117  TH1D* h_dEdX_photon_NueCC;
118  TH1D* h_dEdX_photon_NC;
119  TH1D* h_dEdX_proton_NueCC;
120  TH1D* h_dEdX_proton_NC;
121  TH1D* h_dEdX_neutron_NueCC;
122  TH1D* h_dEdX_neutron_NC;
123  TH1D* h_dEdX_chargedpion_NueCC;
124  TH1D* h_dEdX_chargedpion_NC;
125  TH1D* h_dEdX_neutralpion_NueCC;
126  TH1D* h_dEdX_neutralpion_NC;
127  TH1D* h_dEdX_everythingelse_NueCC;
128  TH1D* h_dEdX_everythingelse_NC;
129 
130  TH1D* h_dEdXasymm_electronorpositron_NueCC;
131  TH1D* h_dEdXasymm_photon_NC;
132 
133  TH1D* h_mpi0_electronorpositron_NueCC;
134  TH1D* h_mpi0_photon_NC;
135 
136  //Study CNN track/shower id
137  TH1D* h_trklike_em;
138  TH1D* h_trklike_nonem;
139 
140  //Study the angle between the reconstructed shower direction w.r.t MC true particle direction
141  TH1D* h_CosThetaShDirwrtTrueparticle_electronorpositron_NueCC;
142  TH1D* h_CosThetaShDirwrtTrueparticle_electronorpositron_NC;
143  TH1D* h_CosThetaShDirwrtTrueparticle_photon_NueCC;
144  TH1D* h_CosThetaShDirwrtTrueparticle_photon_NC;
145  TH1D* h_CosThetaShDirwrtTrueparticle_proton_NueCC;
146  TH1D* h_CosThetaShDirwrtTrueparticle_proton_NC;
147  TH1D* h_CosThetaShDirwrtTrueparticle_chargedpion_NueCC;
148  TH1D* h_CosThetaShDirwrtTrueparticle_chargedpion_NC;
149 
150  //Study the reconstructed shower start position (x,y,z) w.r.t MC true particle start position
151  TH1D* h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NueCC;
152  TH1D* h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NueCC;
153  TH1D* h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NueCC;
154 
155  TH1D* h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NC;
156  TH1D* h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NC;
157  TH1D* h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NC;
158 
159  TH1D* h_ShStartXwrtTrueparticleStartXDiff_photon_NueCC;
160  TH1D* h_ShStartYwrtTrueparticleStartYDiff_photon_NueCC;
161  TH1D* h_ShStartZwrtTrueparticleStartZDiff_photon_NueCC;
162 
163  TH1D* h_ShStartXwrtTrueparticleStartXDiff_photon_NC;
164  TH1D* h_ShStartYwrtTrueparticleStartYDiff_photon_NC;
165  TH1D* h_ShStartZwrtTrueparticleStartZDiff_photon_NC;
166 
167  //True photon end position comparison with the reconstructed shower start position
168  TH1D* h_ShStartXwrtTrueparticleEndXDiff_photon_NueCC;
169  TH1D* h_ShStartYwrtTrueparticleEndYDiff_photon_NueCC;
170  TH1D* h_ShStartZwrtTrueparticleEndZDiff_photon_NueCC;
171 
172  TH1D* h_ShStartXwrtTrueparticleEndXDiff_photon_NC;
173  TH1D* h_ShStartYwrtTrueparticleEndYDiff_photon_NC;
174  TH1D* h_ShStartZwrtTrueparticleEndZDiff_photon_NC;
175 
176  TH1D* h_ShStartXwrtTrueparticleStartXDiff_proton_NueCC;
177  TH1D* h_ShStartYwrtTrueparticleStartYDiff_proton_NueCC;
178  TH1D* h_ShStartZwrtTrueparticleStartZDiff_proton_NueCC;
179 
180  TH1D* h_ShStartXwrtTrueparticleStartXDiff_proton_NC;
181  TH1D* h_ShStartYwrtTrueparticleStartYDiff_proton_NC;
182  TH1D* h_ShStartZwrtTrueparticleStartZDiff_proton_NC;
183 
184  TH1D* h_ShStartXwrtTrueparticleStartXDiff_chargedpion_NueCC;
185  TH1D* h_ShStartYwrtTrueparticleStartYDiff_chargedpion_NueCC;
186  TH1D* h_ShStartZwrtTrueparticleStartZDiff_chargedpion_NueCC;
187 
188  TH1D* h_ShStartXwrtTrueparticleStartXDiff_chargedpion_NC;
189  TH1D* h_ShStartYwrtTrueparticleStartYDiff_chargedpion_NC;
190  TH1D* h_ShStartZwrtTrueparticleStartZDiff_chargedpion_NC;
191 
192  // Event
193  int Event;
194  int Run;
195  int SubRun;
196 
197  //MC truth
198  int MC_incoming_PDG;
199  int MC_lepton_PDG;
200  int MC_isCC;
201  int MC_channel;
202  int MC_target;
203  double MC_Q2;
204  double MC_W;
205  double MC_vertex[4];
206  double MC_incoming_P[4];
207  double MC_lepton_startMomentum[4];
208  double MC_lepton_endMomentum[4];
209  double MC_lepton_startXYZT[4];
210  double MC_lepton_endXYZT[4];
211  double MC_lepton_theta;
212  int MC_leptonID;
213  int MC_LeptonTrack;
214 
215  double sh_direction_X[MAX_SHOWERS];
216  double sh_direction_Y[MAX_SHOWERS];
217  double sh_direction_Z[MAX_SHOWERS];
218  double sh_start_X[MAX_SHOWERS];
219  double sh_start_Y[MAX_SHOWERS];
220  double sh_start_Z[MAX_SHOWERS];
221  double sh_energy[MAX_SHOWERS][3];
222  double sh_MIPenergy[MAX_SHOWERS][3];
223  double sh_dEdx[MAX_SHOWERS][3];
224  int sh_bestplane[MAX_SHOWERS];
225  double sh_length[MAX_SHOWERS];
226  int sh_hasPrimary_e[MAX_SHOWERS];
227  double sh_Efrac_contamination[MAX_SHOWERS];
228  double sh_purity[MAX_SHOWERS];
229  double sh_completeness[MAX_SHOWERS];
230  int sh_nHits[MAX_SHOWERS];
231  int sh_largest;
232  int sh_pdg[MAX_SHOWERS];
233  double sh_dEdxasymm[MAX_SHOWERS];
234  double sh_mpi0;
235 
236  int n_recoShowers;
237  double sh_Efrac_best;
238 
239  float fFidVolCutX;
240  float fFidVolCutY;
241  float fFidVolCutZ;
242 
243  float fFidVolXmin;
244  float fFidVolXmax;
245  float fFidVolYmin;
246  float fFidVolYmax;
247  float fFidVolZmin;
248  float fFidVolZmax;
249 
250  art::ServiceHandle<geo::Geometry const> geom;
251 
252  }; // class NeutrinoShowerEff
253 
254  //========================================================================
255  NeutrinoShowerEff::NeutrinoShowerEff(fhicl::ParameterSet const& p) : EDAnalyzer(p)
256  {
257  fMCTruthModuleLabel = p.get<art::InputTag>("MCTruthModuleLabel");
258  fHitModuleLabel = p.get<art::InputTag>("HitModuleLabel");
259  fShowerModuleLabel = p.get<art::InputTag>("ShowerModuleLabel");
260  fCNNEMModuleLabel = p.get<art::InputTag>("CNNEMModuleLabel", "");
261  fLeptonPDGcode = p.get<int>("LeptonPDGcode");
262  fNeutrinoPDGcode = p.get<int>("NeutrinoPDGcode");
263  fMaxNeutrinoE = p.get<double>("MaxNeutrinoE");
264  fMaxEfrac = p.get<double>("MaxEfrac");
265  fMinCompleteness = p.get<double>("MinCompleteness");
266  fSaveMCTree = p.get<bool>("SaveMCTree");
267  fFidVolCutX = p.get<float>("FidVolCutX");
268  fFidVolCutY = p.get<float>("FidVolCutY");
269  fFidVolCutZ = p.get<float>("FidVolCutZ");
270  }
271  //========================================================================
272  //========================================================================
273  void
274  NeutrinoShowerEff::beginJob()
275  {
276  cout << "job begin..." << endl;
277 
278  // Get geometry.
279  auto const* geo = lar::providerFrom<geo::Geometry>();
280  // Define histogram boundaries (cm).
281  // For now only draw cryostat=0.
282  double minx = 1e9;
283  double maxx = -1e9;
284  double miny = 1e9;
285  double maxy = -1e9;
286  double minz = 1e9;
287  double maxz = -1e9;
288  for (size_t i = 0; i < geo->NTPC(); ++i) {
289  double local[3] = {0., 0., 0.};
290  double world[3] = {0., 0., 0.};
291  const geo::TPCGeo& tpc = geo->TPC(i);
292  tpc.LocalToWorld(local, world);
293  if (minx > world[0] - geo->DetHalfWidth(i)) minx = world[0] - geo->DetHalfWidth(i);
294  if (maxx < world[0] + geo->DetHalfWidth(i)) maxx = world[0] + geo->DetHalfWidth(i);
295  if (miny > world[1] - geo->DetHalfHeight(i)) miny = world[1] - geo->DetHalfHeight(i);
296  if (maxy < world[1] + geo->DetHalfHeight(i)) maxy = world[1] + geo->DetHalfHeight(i);
297  if (minz > world[2] - geo->DetLength(i) / 2.) minz = world[2] - geo->DetLength(i) / 2.;
298  if (maxz < world[2] + geo->DetLength(i) / 2.) maxz = world[2] + geo->DetLength(i) / 2.;
299  }
300 
301  fFidVolXmin = minx + fFidVolCutX;
302  fFidVolXmax = maxx - fFidVolCutX;
303  fFidVolYmin = miny + fFidVolCutY;
304  fFidVolYmax = maxy - fFidVolCutY;
305  fFidVolZmin = minz + fFidVolCutZ;
306  fFidVolZmax = maxz - fFidVolCutZ;
307 
308  std::cout << "Fiducial volume:"
309  << "\n"
310  << fFidVolXmin << "\t< x <\t" << fFidVolXmax << "\n"
311  << fFidVolYmin << "\t< y <\t" << fFidVolYmax << "\n"
312  << fFidVolZmin << "\t< z <\t" << fFidVolZmax << "\n";
313 
314  art::ServiceHandle<art::TFileService const> tfs;
315 
316  double E_bins[21] = {0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4, 4.5, 5.0,
317  5.5, 6.0, 7.0, 8.0, 10.0, 12.0, 14.0, 17.0, 20.0, 25.0};
318  double theta_bin[44] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.,
319  11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 22.,
320  24., 26., 28., 30., 32., 34., 36., 38., 40., 42., 44.,
321  46., 48., 50., 55., 60., 65., 70., 75., 80., 85., 90.};
322  // double Pbins[18] ={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};
323 
324  h_Ev_den =
325  tfs->make<TH1D>("h_Ev_den",
326  "Neutrino Energy; Neutrino Energy (GeV); Shower Reconstruction Efficiency",
327  20,
328  E_bins);
329  h_Ev_den->Sumw2();
330  h_Ev_num =
331  tfs->make<TH1D>("h_Ev_num",
332  "Neutrino Energy; Neutrino Energy (GeV); Shower Reconstruction Efficiency",
333  20,
334  E_bins);
335  h_Ev_num->Sumw2();
336  h_Ev_num_dEdx =
337  tfs->make<TH1D>("h_Ev_num_dEdx",
338  "Neutrino Energy; Neutrino Energy (GeV); Shower Reconstruction Efficiency",
339  20,
340  E_bins);
341  h_Ev_num_dEdx->Sumw2();
342 
343  h_Ee_den =
344  tfs->make<TH1D>("h_Ee_den",
345  "Electron Energy; Electron Energy (GeV); Shower Reconstruction Efficiency",
346  20,
347  E_bins);
348  h_Ee_den->Sumw2();
349  h_Ee_num =
350  tfs->make<TH1D>("h_Ee_num",
351  "Electron Energy; Electron Energy (GeV); Shower Reconstruction Efficiency",
352  20,
353  E_bins);
354  h_Ee_num->Sumw2();
355  h_Ee_num_dEdx =
356  tfs->make<TH1D>("h_Ee_num_dEdx",
357  "Electron Energy; Electron Energy (GeV); Shower Reconstruction Efficiency",
358  20,
359  E_bins);
360  h_Ee_num_dEdx->Sumw2();
361 
362  h_Pe_den = tfs->make<TH1D>(
363  "h_Pe_den",
364  "Electron Momentum; Electron Momentum (GeV); Shower reconstruction Efficiency",
365  20,
366  E_bins);
367  h_Pe_den->Sumw2();
368  h_Pe_num = tfs->make<TH1D>(
369  "h_Pe_num",
370  "Electron Momentum; Electron Momentum (GeV); Shower reconstruction Efficiency",
371  20,
372  E_bins);
373  h_Pe_num->Sumw2();
374 
375  h_theta_den = tfs->make<TH1D>(
376  "h_theta_den",
377  "CosTheta; CosTheta w.r.t beam direction (Degrees); Shower reconstruction Efficiency",
378  43,
379  theta_bin);
380  h_theta_den->Sumw2();
381  h_theta_num = tfs->make<TH1D>(
382  "h_theta_num",
383  "CosTheta; CosTheta w.r.t beam direction (Degrees); Shower reconstruction Efficiency",
384  43,
385  theta_bin);
386  h_theta_num->Sumw2();
387 
388  h_Efrac_shContamination = tfs->make<TH1D>(
389  "h_Efrac_shContamination", "Efrac Lepton; Energy fraction (contamination);", 60, 0, 1.2);
390  h_Efrac_shContamination->Sumw2();
391  h_Efrac_shPurity =
392  tfs->make<TH1D>("h_Efrac_shPurity", "Efrac Lepton; Energy fraction (Purity);", 60, 0, 1.2);
393  h_Efrac_shPurity->Sumw2();
394  h_Ecomplet_lepton =
395  tfs->make<TH1D>("h_Ecomplet_lepton", "Ecomplet Lepton; Shower Completeness;", 60, 0, 1.2);
396  h_Ecomplet_lepton->Sumw2();
397 
398  h_HighestHitsProducedParticlePDG_NueCC =
399  tfs->make<TH1D>("h_HighestHitsProducedParticlePDG_NueCC",
400  "PDG Code; PDG Code;",
401  4,
402  -0.5,
403  3.5); //0 for undefined, 1=electron, 2=photon, 3=anything else //Signal
404  h_HighestHitsProducedParticlePDG_NueCC->Sumw2();
405  h_HighestHitsProducedParticlePDG_bkg =
406  tfs->make<TH1D>("h_HighestHitsProducedParticlePDG_bkg",
407  "PDG Code; PDG Code;",
408  4,
409  -0.5,
410  3.5); //0 for undefined, 1=electron, 2=photon, 3=anything else //bkg
411  h_HighestHitsProducedParticlePDG_bkg->Sumw2();
412 
413  h_Efrac_NueCCPurity = tfs->make<TH1D>(
414  "h_Efrac_NueCCPurity", "Efrac NueCC; Energy fraction (Purity);", 60, 0, 1.2); //Signal
415  h_Efrac_NueCCPurity->Sumw2();
416  h_Ecomplet_NueCC =
417  tfs->make<TH1D>("h_Ecomplet_NueCC", "Ecomplet NueCC; Shower Completeness;", 60, 0, 1.2);
418  h_Ecomplet_NueCC->Sumw2();
419 
420  h_Efrac_bkgPurity = tfs->make<TH1D>(
421  "h_Efrac_bkgPurity", "Efrac bkg; Energy fraction (Purity);", 60, 0, 1.2); //Background
422  h_Efrac_bkgPurity->Sumw2();
423  h_Ecomplet_bkg =
424  tfs->make<TH1D>("h_Ecomplet_bkg", "Ecomplet bkg; Shower Completeness;", 60, 0, 1.2);
425  h_Ecomplet_bkg->Sumw2();
426 
427  h_esh_bestplane_NueCC =
428  tfs->make<TH1D>("h_esh_bestplane_NueCC", "Best plane; Best plane;", 4, -0.5, 3.5);
429  h_esh_bestplane_NC =
430  tfs->make<TH1D>("h_esh_bestplane_NC", "Best plane; Best plane;", 4, -0.5, 3.5);
431  h_dEdX_electronorpositron_NueCC = tfs->make<TH1D>("h_dEdX_electronorpositron_NueCC",
432  "dE/dX; Electron or Positron dE/dX (MeV/cm);",
433  100,
434  0.0,
435  15.0);
436  h_dEdX_electronorpositron_NC = tfs->make<TH1D>("h_dEdX_electronorpositron_NC",
437  "dE/dX; Electron or Positron dE/dX (MeV/cm);",
438  100,
439  0.0,
440  15.0);
441  h_dEdX_photon_NueCC =
442  tfs->make<TH1D>("h_dEdX_photon_NueCC", "dE/dX; photon dE/dX (MeV/cm);", 100, 0.0, 15.0);
443  h_dEdX_photon_NC =
444  tfs->make<TH1D>("h_dEdX_photon_NC", "dE/dX; photon dE/dX (MeV/cm);", 100, 0.0, 15.0);
445  h_dEdX_proton_NueCC =
446  tfs->make<TH1D>("h_dEdX_proton_NueCC", "dE/dX; proton dE/dX (MeV/cm);", 100, 0.0, 15.0);
447  h_dEdX_proton_NC =
448  tfs->make<TH1D>("h_dEdX_proton_NC", "dE/dX; proton dE/dX (MeV/cm);", 100, 0.0, 15.0);
449  h_dEdX_neutron_NueCC =
450  tfs->make<TH1D>("h_dEdX_neutron_NueCC", "dE/dX; neutron dE/dX (MeV/cm);", 100, 0.0, 15.0);
451  h_dEdX_neutron_NC =
452  tfs->make<TH1D>("h_dEdX_neutron_NC", "dE/dX; neutron dE/dX (MeV/cm);", 100, 0.0, 15.0);
453  h_dEdX_chargedpion_NueCC = tfs->make<TH1D>(
454  "h_dEdX_chargedpion_NueCC", "dE/dX; charged pion dE/dX (MeV/cm);", 100, 0.0, 15.0);
455  h_dEdX_chargedpion_NC = tfs->make<TH1D>(
456  "h_dEdX_chargedpion_NC", "dE/dX; charged pion dE/dX (MeV/cm);", 100, 0.0, 15.0);
457  h_dEdX_neutralpion_NueCC = tfs->make<TH1D>(
458  "h_dEdX_neutralpion_NueCC", "dE/dX; neutral pion dE/dX (MeV/cm);", 100, 0.0, 15.0);
459  h_dEdX_neutralpion_NC = tfs->make<TH1D>(
460  "h_dEdX_neutralpion_NC", "dE/dX; neutral pion dE/dX (MeV/cm);", 100, 0.0, 15.0);
461  h_dEdX_everythingelse_NueCC = tfs->make<TH1D>(
462  "h_dEdX_everythingelse_NueCC", "dE/dX; everythingelse dE/dX (MeV/cm);", 100, 0.0, 15.0);
463  h_dEdX_everythingelse_NC = tfs->make<TH1D>(
464  "h_dEdX_everythingelse_NC", "dE/dX; everythingelse dE/dX (MeV/cm);", 100, 0.0, 15.0);
465 
466  h_dEdXasymm_electronorpositron_NueCC =
467  tfs->make<TH1D>("h_dEdXasymm_electronorpositron_NueCC",
468  "dE/dX asymmetry; Electron or Positron dE/dX asymmetry;",
469  60,
470  0.0,
471  1.2);
472  h_dEdXasymm_photon_NC = tfs->make<TH1D>(
473  "h_dEdXasymm_photon_NC", "dE/dX asymmetry; photon dE/dx asymmetry;", 60, 0.0, 1.2);
474 
475  h_mpi0_electronorpositron_NueCC =
476  tfs->make<TH1D>("h_mpi0_electronorpositron_NueCC",
477  "m(#gamma#gamma); Electron or Positron dE/dX (MeV/cm);",
478  100,
479  0,
480  1);
481  h_mpi0_photon_NC = tfs->make<TH1D>(
482  "h_mpi0_photon_NC", "m(#gamma#gamma); Electron or Positron dE/dX (MeV/cm);", 100, 0, 1);
483 
484  h_esh_bestplane_NueCC->Sumw2();
485  h_esh_bestplane_NC->Sumw2();
486  h_dEdX_electronorpositron_NueCC->Sumw2();
487  h_dEdX_electronorpositron_NC->Sumw2();
488  h_dEdX_photon_NueCC->Sumw2();
489  h_dEdX_photon_NC->Sumw2();
490  h_dEdX_proton_NueCC->Sumw2();
491  h_dEdX_proton_NC->Sumw2();
492  h_dEdX_neutron_NueCC->Sumw2();
493  h_dEdX_neutron_NC->Sumw2();
494  h_dEdX_chargedpion_NueCC->Sumw2();
495  h_dEdX_chargedpion_NC->Sumw2();
496  h_dEdX_neutralpion_NueCC->Sumw2();
497  h_dEdX_neutralpion_NC->Sumw2();
498  h_dEdX_everythingelse_NueCC->Sumw2();
499  h_dEdX_everythingelse_NC->Sumw2();
500 
501  h_dEdXasymm_electronorpositron_NueCC->Sumw2();
502  h_dEdXasymm_photon_NC->Sumw2();
503 
504  h_mpi0_electronorpositron_NueCC->Sumw2();
505  h_mpi0_photon_NC->Sumw2();
506 
507  h_trklike_em = tfs->make<TH1D>("h_trklike_em", "EM hits; Track-like Score;", 100, 0, 1);
508  h_trklike_nonem =
509  tfs->make<TH1D>("h_trklike_nonem", "Non-EM hits; Track-like Score;", 100, 0, 1);
510 
511  //Study the constheta angle between the reconstructed shower direction w.r.t MC true particle direction
512  h_CosThetaShDirwrtTrueparticle_electronorpositron_NueCC =
513  tfs->make<TH1D>("h_CosThetaShDirwrtTrueparticle_electronorpositron_NueCC",
514  "CosTheta; cos#theta;",
515  110,
516  -1.1,
517  1.1);
518  h_CosThetaShDirwrtTrueparticle_electronorpositron_NC =
519  tfs->make<TH1D>("h_CosThetaShDirwrtTrueparticle_electronorpositron_NC",
520  "CosTheta;cos#theta;",
521  110,
522  -1.1,
523  1.1);
524  h_CosThetaShDirwrtTrueparticle_photon_NueCC = tfs->make<TH1D>(
525  "h_CosThetaShDirwrtTrueparticle_photon_NueCC", "CosTheta;cos#theta;", 110, -1.1, 1.1);
526  h_CosThetaShDirwrtTrueparticle_photon_NC = tfs->make<TH1D>(
527  "h_CosThetaShDirwrtTrueparticle_photon_NC", "CosTheta;cos#theta;", 110, -1.1, 1.1);
528  h_CosThetaShDirwrtTrueparticle_proton_NueCC = tfs->make<TH1D>(
529  "h_CosThetaShDirwrtTrueparticle_proton_NueCC", "CosTheta;cos#theta;", 110, -1.1, 1.1);
530  h_CosThetaShDirwrtTrueparticle_proton_NC = tfs->make<TH1D>(
531  "h_CosThetaShDirwrtTrueparticle_proton_NC", "CosTheta;cos#theta;", 110, -1.1, 1.1);
532  h_CosThetaShDirwrtTrueparticle_chargedpion_NueCC = tfs->make<TH1D>(
533  "h_CosThetaShDirwrtTrueparticle_chargedpion_NueCC", "CosTheta;cos#theta;", 110, -1.1, 1.1);
534  h_CosThetaShDirwrtTrueparticle_chargedpion_NC = tfs->make<TH1D>(
535  "h_CosThetaShDirwrtTrueparticle_chargedpion_NC", "CosTheta;cos#theta;", 110, -1.1, 1.1);
536 
537  //Study the reconstructed shower start position (x,y,z) w.r.t MC true particle start position
538  h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NueCC =
539  tfs->make<TH1D>("h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NueCC",
540  "ShVx-TrueParticleVx; ShVx-TrueParticleVx (cm);",
541  100,
542  -5.0,
543  5.0);
544  h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NueCC =
545  tfs->make<TH1D>("h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NueCC",
546  "ShVy-TrueParticleVy; ShVy-TrueParticleVy (cm);",
547  100,
548  -5.0,
549  5.0);
550  h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NueCC =
551  tfs->make<TH1D>("h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NueCC",
552  "ShVz-TrueParticleVz; ShVz-TrueParticleVz (cm);",
553  100,
554  -5.0,
555  5.0);
556 
557  h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NC =
558  tfs->make<TH1D>("h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NC",
559  "ShVx-TrueParticleVx; ShVx-TrueParticleVx (cm);",
560  100,
561  -5.0,
562  5.0);
563  h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NC =
564  tfs->make<TH1D>("h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NC",
565  "ShVy-TrueParticleVy; ShVy-TrueParticleVy (cm);",
566  100,
567  -5.0,
568  5.0);
569  h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NC =
570  tfs->make<TH1D>("h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NC",
571  "ShVz-TrueParticleVz; ShVz-TrueParticleVz (cm);",
572  100,
573  -5.0,
574  5.0);
575 
576  h_ShStartXwrtTrueparticleStartXDiff_photon_NueCC =
577  tfs->make<TH1D>("h_ShStartXwrtTrueparticleStartXDiff_photon_NueCC",
578  "ShVx-TrueParticleVx; ShVx-TrueParticleVx (cm);",
579  100,
580  -5.0,
581  5.0);
582  h_ShStartYwrtTrueparticleStartYDiff_photon_NueCC =
583  tfs->make<TH1D>("h_ShStartYwrtTrueparticleStartYDiff_photon_NueCC",
584  "ShVy-TrueParticleVy; ShVy-TrueParticleVy (cm);",
585  100,
586  -5.0,
587  5.0);
588  h_ShStartZwrtTrueparticleStartZDiff_photon_NueCC =
589  tfs->make<TH1D>("h_ShStartZwrtTrueparticleStartZDiff_photon_NueCC",
590  "ShVz-TrueParticleVz; ShVz-TrueParticleVz (cm);",
591  100,
592  -5.0,
593  5.0);
594 
595  h_ShStartXwrtTrueparticleStartXDiff_photon_NC =
596  tfs->make<TH1D>("h_ShStartXwrtTrueparticleStartXDiff_photon_NC",
597  "ShVx-TrueParticleVx; ShVx-TrueParticleVx (cm);",
598  100,
599  -5.0,
600  5.0);
601  h_ShStartYwrtTrueparticleStartYDiff_photon_NC =
602  tfs->make<TH1D>("h_ShStartYwrtTrueparticleStartYDiff_photon_NC",
603  "ShVy-TrueParticleVy; ShVy-TrueParticleVy (cm);",
604  100,
605  -5.0,
606  5.0);
607  h_ShStartZwrtTrueparticleStartZDiff_photon_NC =
608  tfs->make<TH1D>("h_ShStartZwrtTrueparticleStartZDiff_photon_NC",
609  "ShVz-TrueParticleVz; ShVz-TrueParticleVz (cm);",
610  100,
611  -5.0,
612  5.0);
613 
614  h_ShStartXwrtTrueparticleEndXDiff_photon_NueCC =
615  tfs->make<TH1D>("h_ShStartXwrtTrueparticleEndXDiff_photon_NueCC",
616  "ShVx-TrueParticleVx; ShVx-TrueParticleVx (cm);",
617  100,
618  -5.0,
619  5.0);
620  h_ShStartYwrtTrueparticleEndYDiff_photon_NueCC =
621  tfs->make<TH1D>("h_ShStartYwrtTrueparticleEndYDiff_photon_NueCC",
622  "ShVy-TrueParticleVy; ShVy-TrueParticleVy (cm);",
623  100,
624  -5.0,
625  5.0);
626  h_ShStartZwrtTrueparticleEndZDiff_photon_NueCC =
627  tfs->make<TH1D>("h_ShStartZwrtTrueparticleEndZDiff_photon_NueCC",
628  "ShVz-TrueParticleVz; ShVz-TrueParticleVz (cm);",
629  100,
630  -5.0,
631  5.0);
632 
633  h_ShStartXwrtTrueparticleEndXDiff_photon_NC =
634  tfs->make<TH1D>("h_ShStartXwrtTrueparticleEndXDiff_photon_NC",
635  "ShVx-TrueParticleVx; ShVx-TrueParticleVx (cm);",
636  100,
637  -5.0,
638  5.0);
639  h_ShStartYwrtTrueparticleEndYDiff_photon_NC =
640  tfs->make<TH1D>("h_ShStartYwrtTrueparticleEndYDiff_photon_NC",
641  "ShVy-TrueParticleVy; ShVy-TrueParticleVy (cm);",
642  100,
643  -5.0,
644  5.0);
645  h_ShStartZwrtTrueparticleEndZDiff_photon_NC =
646  tfs->make<TH1D>("h_ShStartZwrtTrueparticleEndZDiff_photon_NC",
647  "ShVz-TrueParticleVz; ShVz-TrueParticleVz (cm);",
648  100,
649  -5.0,
650  5.0);
651 
652  h_ShStartXwrtTrueparticleStartXDiff_proton_NueCC =
653  tfs->make<TH1D>("h_ShStartXwrtTrueparticleStartXDiff_proton_NueCC",
654  "ShVx-TrueParticleVx; ShVx-TrueParticleVx (cm);",
655  100,
656  -5.0,
657  5.0);
658  h_ShStartYwrtTrueparticleStartYDiff_proton_NueCC =
659  tfs->make<TH1D>("h_ShStartYwrtTrueparticleStartYDiff_proton_NueCC",
660  "ShVy-TrueParticleVy; ShVy-TrueParticleVy (cm);",
661  100,
662  -5.0,
663  5.0);
664  h_ShStartZwrtTrueparticleStartZDiff_proton_NueCC =
665  tfs->make<TH1D>("h_ShStartZwrtTrueparticleStartZDiff_proton_NueCC",
666  "ShVz-TrueParticleVz; ShVz-TrueParticleVz (cm);",
667  100,
668  -5.0,
669  5.0);
670 
671  h_ShStartXwrtTrueparticleStartXDiff_proton_NC =
672  tfs->make<TH1D>("h_ShStartXwrtTrueparticleStartXDiff_proton_NC",
673  "ShVx-TrueParticleVx; ShVx-TrueParticleVx (cm);",
674  100,
675  -5.0,
676  5.0);
677  h_ShStartYwrtTrueparticleStartYDiff_proton_NC =
678  tfs->make<TH1D>("h_ShStartYwrtTrueparticleStartYDiff_proton_NC",
679  "ShVy-TrueParticleVy; ShVy-TrueParticleVy (cm);",
680  100,
681  -5.0,
682  5.0);
683  h_ShStartZwrtTrueparticleStartZDiff_proton_NC =
684  tfs->make<TH1D>("h_ShStartZwrtTrueparticleStartZDiff_proton_NC",
685  "ShVz-TrueParticleVz; ShVz-TrueParticleVz (cm);",
686  100,
687  -5.0,
688  5.0);
689 
690  h_ShStartXwrtTrueparticleStartXDiff_chargedpion_NueCC =
691  tfs->make<TH1D>("h_ShStartXwrtTrueparticleStartXDiff_chargedpion_NueCC",
692  "ShVx-TrueParticleVx; ShVx-TrueParticleVx (cm);",
693  100,
694  -5.0,
695  5.0);
696  h_ShStartYwrtTrueparticleStartYDiff_chargedpion_NueCC =
697  tfs->make<TH1D>("h_ShStartYwrtTrueparticleStartYDiff_chargedpion_NueCC",
698  "ShVy-TrueParticleVy; ShVy-TrueParticleVy (cm);",
699  100,
700  -5.0,
701  5.0);
702  h_ShStartZwrtTrueparticleStartZDiff_chargedpion_NueCC =
703  tfs->make<TH1D>("h_ShStartZwrtTrueparticleStartZDiff_chargedpion_NueCC",
704  "ShVz-TrueParticleVz; ShVz-TrueParticleVz (cm);",
705  100,
706  -5.0,
707  5.0);
708 
709  h_ShStartXwrtTrueparticleStartXDiff_chargedpion_NC =
710  tfs->make<TH1D>("h_ShStartXwrtTrueparticleStartXDiff_chargedpion_NC",
711  "ShVx-TrueParticleVx; ShVx-TrueParticleVx (cm);",
712  100,
713  -5.0,
714  5.0);
715  h_ShStartYwrtTrueparticleStartYDiff_chargedpion_NC =
716  tfs->make<TH1D>("h_ShStartYwrtTrueparticleStartYDiff_chargedpion_NC",
717  "ShVy-TrueParticleVy; ShVy-TrueParticleVy (cm);",
718  100,
719  -5.0,
720  5.0);
721  h_ShStartZwrtTrueparticleStartZDiff_chargedpion_NC =
722  tfs->make<TH1D>("h_ShStartZwrtTrueparticleStartZDiff_chargedpion_NC",
723  "ShVz-TrueParticleVz; ShVz-TrueParticleVz (cm);",
724  100,
725  -5.0,
726  5.0);
727 
728  //Study the constheta angle between the reconstructed shower direction w.r.t MC true particle direction
729  h_CosThetaShDirwrtTrueparticle_electronorpositron_NueCC->Sumw2();
730  h_CosThetaShDirwrtTrueparticle_electronorpositron_NC->Sumw2();
731  h_CosThetaShDirwrtTrueparticle_photon_NueCC->Sumw2();
732  h_CosThetaShDirwrtTrueparticle_photon_NC->Sumw2();
733  h_CosThetaShDirwrtTrueparticle_proton_NueCC->Sumw2();
734  h_CosThetaShDirwrtTrueparticle_proton_NC->Sumw2();
735  h_CosThetaShDirwrtTrueparticle_chargedpion_NueCC->Sumw2();
736  h_CosThetaShDirwrtTrueparticle_chargedpion_NC->Sumw2();
737 
738  //Study the reconstructed shower start position (x,y,z) w.r.t MC true particle start position
739  h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NueCC->Sumw2();
740  h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NueCC->Sumw2();
741  h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NueCC->Sumw2();
742 
743  h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NC->Sumw2();
744  h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NC->Sumw2();
745  h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NC->Sumw2();
746 
747  h_ShStartXwrtTrueparticleStartXDiff_photon_NueCC->Sumw2();
748  h_ShStartYwrtTrueparticleStartYDiff_photon_NueCC->Sumw2();
749  h_ShStartZwrtTrueparticleStartZDiff_photon_NueCC->Sumw2();
750 
751  h_ShStartXwrtTrueparticleStartXDiff_photon_NC->Sumw2();
752  h_ShStartYwrtTrueparticleStartYDiff_photon_NC->Sumw2();
753  h_ShStartZwrtTrueparticleStartZDiff_photon_NC->Sumw2();
754 
755  h_ShStartXwrtTrueparticleEndXDiff_photon_NueCC->Sumw2();
756  h_ShStartYwrtTrueparticleEndYDiff_photon_NueCC->Sumw2();
757  h_ShStartZwrtTrueparticleEndZDiff_photon_NueCC->Sumw2();
758 
759  h_ShStartXwrtTrueparticleEndXDiff_photon_NC->Sumw2();
760  h_ShStartYwrtTrueparticleEndYDiff_photon_NC->Sumw2();
761  h_ShStartZwrtTrueparticleEndZDiff_photon_NC->Sumw2();
762 
763  h_ShStartXwrtTrueparticleStartXDiff_proton_NueCC->Sumw2();
764  h_ShStartYwrtTrueparticleStartYDiff_proton_NueCC->Sumw2();
765  h_ShStartZwrtTrueparticleStartZDiff_proton_NueCC->Sumw2();
766 
767  h_ShStartXwrtTrueparticleStartXDiff_proton_NC->Sumw2();
768  h_ShStartYwrtTrueparticleStartYDiff_proton_NC->Sumw2();
769  h_ShStartZwrtTrueparticleStartZDiff_proton_NC->Sumw2();
770 
771  h_ShStartXwrtTrueparticleStartXDiff_chargedpion_NueCC->Sumw2();
772  h_ShStartYwrtTrueparticleStartYDiff_chargedpion_NueCC->Sumw2();
773  h_ShStartZwrtTrueparticleStartZDiff_chargedpion_NueCC->Sumw2();
774 
775  h_ShStartXwrtTrueparticleStartXDiff_chargedpion_NC->Sumw2();
776  h_ShStartYwrtTrueparticleStartYDiff_chargedpion_NC->Sumw2();
777  h_ShStartZwrtTrueparticleStartZDiff_chargedpion_NC->Sumw2();
778 
779  if (fSaveMCTree) {
780  fEventTree = new TTree("Event", "Event Tree from Sim & Reco");
781  fEventTree->Branch("eventNo", &Event);
782  fEventTree->Branch("runNo", &Run);
783  fEventTree->Branch("subRunNo", &SubRun);
784  fEventTree->Branch("mc_incoming_PDG", &MC_incoming_PDG);
785  fEventTree->Branch("mc_lepton_PDG", &MC_lepton_PDG);
786  fEventTree->Branch("mc_isCC", &MC_isCC);
787  fEventTree->Branch("mc_target", &MC_target);
788  fEventTree->Branch("mc_channel", &MC_channel);
789  fEventTree->Branch("mc_Q2", &MC_Q2);
790  fEventTree->Branch("mc_W", &MC_W);
791  fEventTree->Branch("mc_vertex", &MC_vertex, "mc_vertex[4]/D");
792  fEventTree->Branch("mc_incoming_P", &MC_incoming_P, "mc_incoming_P[4]/D");
793  fEventTree->Branch(
794  "mc_lepton_startMomentum", &MC_lepton_startMomentum, "mc_lepton_startMomentum[4]/D");
795  fEventTree->Branch(
796  "mc_lepton_endMomentum", &MC_lepton_endMomentum, "mc_lepton_endMomentum[4]/D");
797  fEventTree->Branch("mc_lepton_startXYZT", &MC_lepton_startXYZT, "mc_lepton_startXYZT[4]/D");
798  fEventTree->Branch("mc_lepton_endXYZT", &MC_lepton_endXYZT, "mc_lepton_endXYZT[4]/D");
799  fEventTree->Branch("mc_lepton_theta", &MC_lepton_theta, "mc_lepton_theta/D");
800  fEventTree->Branch("mc_leptonID", &MC_leptonID, "mc_leptonID/I");
801 
802  fEventTree->Branch("n_showers", &n_recoShowers);
803  fEventTree->Branch("sh_direction_X", &sh_direction_X, "sh_direction_X[n_showers]/D");
804  fEventTree->Branch("sh_direction_Y", &sh_direction_Y, "sh_direction_Y[n_showers]/D");
805  fEventTree->Branch("sh_direction_Z", &sh_direction_Z, "sh_direction_Z[n_showers]/D");
806  fEventTree->Branch("sh_start_X", &sh_start_X, "sh_start_X[n_showers]/D");
807  fEventTree->Branch("sh_start_Y", &sh_start_Y, "sh_start_Y[n_showers]/D");
808  fEventTree->Branch("sh_start_Z", &sh_start_Z, "sh_start_Z[n_showers]/D");
809  fEventTree->Branch("sh_energy", &sh_energy, "sh_energy[n_showers][3]/D");
810  fEventTree->Branch("sh_MIPenergy", &sh_MIPenergy, "sh_MIPenergy[n_showers][3]/D");
811  fEventTree->Branch("sh_dEdx", &sh_dEdx, "sh_dEdx[n_showers][3]/D");
812  fEventTree->Branch("sh_bestplane", &sh_bestplane, "sh_bestplane[n_showers]/I");
813  fEventTree->Branch("sh_length", &sh_length, "sh_length[n_showers]/D");
814  fEventTree->Branch("sh_hasPrimary_e", &sh_hasPrimary_e, "sh_hasPrimary_e[n_showers]/I");
815  fEventTree->Branch(
816  "sh_Efrac_contamination", &sh_Efrac_contamination, "sh_Efrac_contamination[n_showers]/D");
817  fEventTree->Branch("sh_purity", &sh_purity, "sh_purity[n_showers]/D");
818  fEventTree->Branch("sh_completeness", &sh_completeness, "sh_completeness[n_showers]/D");
819  fEventTree->Branch("sh_Efrac_best", &sh_Efrac_best, "sh_Efrac_best/D");
820  fEventTree->Branch("sh_nHits", &sh_nHits, "sh_nHits[n_showers]/I");
821  fEventTree->Branch("sh_largest", &sh_largest, "sh_largest/I");
822  fEventTree->Branch("sh_pdg", &sh_pdg, "sh_pdg[n_showers]/I");
823  fEventTree->Branch("sh_dEdxasymm", &sh_dEdxasymm, "sh_dEdxasymm[n_showers]/D");
824  fEventTree->Branch("sh_mpi0", &sh_mpi0, "sh_mpi0/D");
825  }
826  }
827  //========================================================================
828  void
829  NeutrinoShowerEff::endJob()
830  {
831  doEfficiencies();
832  }
833  //========================================================================
834  void
835  NeutrinoShowerEff::beginRun(const art::Run& /*run*/)
836  {
837  mf::LogInfo("NeutrinoShowerEff") << "begin run..." << endl;
838  }
839  //========================================================================
840  void
841  NeutrinoShowerEff::analyze(const art::Event& event)
842  {
843 
844  reset();
845 
846  Event = event.id().event();
847  Run = event.run();
848  SubRun = event.subRun();
849  bool isFiducial = false;
850 
851  auto const clockData =
852  art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
853 
854  processEff(clockData, event, isFiducial);
855  if (fSaveMCTree) {
856  if (isFiducial) fEventTree->Fill();
857  }
858  }
859  //========================================================================
860  void
861  NeutrinoShowerEff::processEff(detinfo::DetectorClocksData const& clockData,
862  const art::Event& event,
863  bool& isFiducial)
864  {
865 
866  //!save neutrino's interaction info
867  art::Handle<std::vector<simb::MCTruth>> MCtruthHandle;
868  event.getByLabel(fMCTruthModuleLabel, MCtruthHandle);
869  std::vector<art::Ptr<simb::MCTruth>> MCtruthlist;
870  art::fill_ptr_vector(MCtruthlist, MCtruthHandle);
871  art::Ptr<simb::MCTruth> MCtruth;
872 
873  //For now assume that there is only one neutrino interaction...
874  int MCinteractions = MCtruthlist.size();
875  for (int i = 0; i < MCinteractions; i++) {
876  MCtruth = MCtruthlist[i];
877  if (MCtruth->NeutrinoSet()) {
878  simb::MCNeutrino nu = MCtruth->GetNeutrino();
879  if (nu.CCNC() == 0)
880  MC_isCC = 1;
881  else if (nu.CCNC() == 1)
882  MC_isCC = 0;
883  simb::MCParticle neutrino = nu.Nu();
884  MC_target = nu.Target();
885  MC_incoming_PDG = std::abs(nu.Nu().PdgCode());
886  MC_Q2 = nu.QSqr();
887  MC_channel = nu.InteractionType();
888  MC_W = nu.W();
889  const TLorentzVector& nu_momentum = nu.Nu().Momentum(0);
890  nu_momentum.GetXYZT(MC_incoming_P);
891  const TLorentzVector& vertex = neutrino.Position(0);
892  vertex.GetXYZT(MC_vertex);
893  simb::MCParticle lepton = nu.Lepton();
894  MC_lepton_PDG = lepton.PdgCode();
895  }
896  }
897 
898  //!save lepton
899  simb::MCParticle* MClepton = NULL;
900  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
901  const sim::ParticleList& plist = pi_serv->ParticleList();
902  simb::MCParticle* particle = 0;
903 
904  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar) {
905  particle = ipar->second;
906  if (std::abs(particle->PdgCode()) == fLeptonPDGcode &&
907  particle->Mother() == 0) { //primary lepton
908  const TLorentzVector& lepton_momentum = particle->Momentum(0);
909  const TLorentzVector& lepton_position = particle->Position(0);
910  const TLorentzVector& lepton_positionEnd = particle->EndPosition();
911  const TLorentzVector& lepton_momentumEnd = particle->EndMomentum();
912  lepton_momentum.GetXYZT(MC_lepton_startMomentum);
913  lepton_position.GetXYZT(MC_lepton_startXYZT);
914  lepton_positionEnd.GetXYZT(MC_lepton_endXYZT);
915  lepton_momentumEnd.GetXYZT(MC_lepton_endMomentum);
916  MC_leptonID = particle->TrackId();
917  MClepton = particle;
918  }
919  }
920  //Saving denominator histograms for lepton pions and protons
921  isFiducial = insideFV(MC_vertex);
922  if (!isFiducial) return;
923  double Pe = sqrt(pow(MC_lepton_startMomentum[0], 2) + pow(MC_lepton_startMomentum[1], 2) +
924  pow(MC_lepton_startMomentum[2], 2));
925  double Pv =
926  sqrt(pow(MC_incoming_P[0], 2) + pow(MC_incoming_P[1], 2) + pow(MC_incoming_P[2], 2));
927  double theta_e = acos((MC_incoming_P[0] * MC_lepton_startMomentum[0] +
928  MC_incoming_P[1] * MC_lepton_startMomentum[1] +
929  MC_incoming_P[2] * MC_lepton_startMomentum[2]) /
930  (Pv * Pe));
931  theta_e *= (180.0 / 3.14159);
932 
933  //save CC events within the fiducial volume with the favorite neutrino flavor
934  if (MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG)) {
935  if (MClepton) {
936  h_Ev_den->Fill(MC_incoming_P[3]);
937  h_Ee_den->Fill(MC_lepton_startMomentum[3]);
938  h_Pe_den->Fill(Pe);
939  h_theta_den->Fill(theta_e);
940  }
941  }
942 
943  //========================================================================
944  //========================================================================
945  // Reco stuff
946  //========================================================================
947  //========================================================================
948  art::Handle<std::vector<recob::Shower>> showerHandle;
949  if (!event.getByLabel(fShowerModuleLabel, showerHandle)) {
950  mf::LogError("NeutrinoShowerEff")
951  << "Could not find shower with label " << fShowerModuleLabel.encode();
952  return;
953  }
954  std::vector<art::Ptr<recob::Shower>> showerlist;
955  art::fill_ptr_vector(showerlist, showerHandle);
956 
957  art::Handle<std::vector<recob::Hit>> hitHandle;
958  std::vector<art::Ptr<recob::Hit>> all_hits;
959  if (event.getByLabel(fHitModuleLabel, hitHandle)) { art::fill_ptr_vector(all_hits, hitHandle); }
960 
961  n_recoShowers = showerlist.size();
962  //if ( n_recoShowers == 0 || n_recoShowers> MAX_SHOWERS ) return;
963  art::FindManyP<recob::Hit> sh_hitsAll(showerHandle, event, fShowerModuleLabel);
964  cout << "Found this many showers " << n_recoShowers << endl;
965  double Efrac_contamination = 999.0;
966  double Efrac_contaminationNueCC = 999.0;
967 
968  double Ecomplet_lepton = 0.0;
969  double Ecomplet_NueCC = 0.0;
970  int ParticlePDG_HighestShHits = 0; //undefined
971  int shower_bestplane = 0;
972  double Showerparticlededx_inbestplane = 0.0;
973  double dEdxasymm_largestshw = -1.;
974 
975  int showerPDGwithHighestHitsforFillingdEdX =
976  0; //0=undefined,1=electronorpositronshower,2=photonshower,3=protonshower,4=neutronshower,5=chargedpionshower,6=neutralpionshower,7=everythingelseshower
977 
978  double ShAngle = -9999.0, ShVxTrueParticleVxDiff = -9999.0, ShVyTrueParticleVyDiff = -9999.0,
979  ShVzTrueParticleVzDiff = -9999.0, ShStartVxTrueParticleEndVxDiff = -9999.0,
980  ShStartVyTrueParticleEndVyDiff = -9999.0, ShStartVzTrueParticleEndVzDiff = -9999.0;
981 
982  const simb::MCParticle* MClepton_reco = NULL;
983  int nHits = 0;
984 
985  TVector3 p1, p2;
986  double E1st = 0;
987  double E2nd = 0;
988 
989  for (int i = 0; i < n_recoShowers; i++) {
990 
991  art::Ptr<recob::Shower> shower = showerlist[i];
992  sh_direction_X[i] = shower->Direction().X();
993  sh_direction_Y[i] = shower->Direction().Y();
994  sh_direction_Z[i] = shower->Direction().Z();
995  sh_start_X[i] = shower->ShowerStart().X();
996  sh_start_Y[i] = shower->ShowerStart().Y();
997  sh_start_Z[i] = shower->ShowerStart().Z();
998  sh_bestplane[i] = shower->best_plane();
999  sh_length[i] = shower->Length();
1000  for (size_t j = 0; j < shower->Energy().size(); j++)
1001  sh_energy[i][j] = shower->Energy()[j];
1002  for (size_t j = 0; j < shower->MIPEnergy().size(); j++)
1003  sh_MIPenergy[i][j] = shower->MIPEnergy()[j];
1004  for (size_t j = 0; j < shower->dEdx().size(); j++)
1005  sh_dEdx[i][j] = shower->dEdx()[j];
1006 
1007  double dEdxasymm = -1;
1008  double dEdx0 = 0;
1009  if (shower->best_plane() >= 0 && shower->best_plane() < int(shower->dEdx().size())) {
1010  dEdx0 = shower->dEdx()[shower->best_plane()];
1011  }
1012  double dEdx1 = 0;
1013  double maxE = 0;
1014  for (int j = 0; j < 3; ++j) {
1015  if (j == shower->best_plane()) continue;
1016  if (j >= int(shower->Energy().size())) continue;
1017  if (shower->Energy()[j] > maxE) {
1018  maxE = shower->Energy()[j];
1019  dEdx1 = shower->dEdx()[j];
1020  }
1021  }
1022  if (dEdx0 || dEdx1) { dEdxasymm = std::abs(dEdx0 - dEdx1) / (dEdx0 + dEdx1); }
1023  sh_dEdxasymm[i] = dEdxasymm;
1024 
1025  if (shower->best_plane() >= 0 && shower->best_plane() < int(shower->Energy().size())) {
1026  if (shower->Energy()[shower->best_plane()] > E1st) {
1027  if (p1.Mag()) {
1028  E2nd = E1st;
1029  p2 = p1;
1030  }
1031  E1st = shower->Energy()[shower->best_plane()];
1032  p1[0] = E1st * shower->Direction().X();
1033  p1[1] = E1st * shower->Direction().Y();
1034  p1[2] = E1st * shower->Direction().Z();
1035  }
1036  else {
1037  if (shower->Energy()[shower->best_plane()] > E2nd) {
1038  E2nd = shower->Energy()[shower->best_plane()];
1039  p2[0] = E2nd * shower->Direction().X();
1040  p2[1] = E2nd * shower->Direction().Y();
1041  p2[2] = E2nd * shower->Direction().Z();
1042  }
1043  }
1044  }
1045 
1046  std::vector<art::Ptr<recob::Hit>> sh_hits = sh_hitsAll.at(i);
1047 
1048  if (!sh_hits.size()) {
1049  //no shower hits found, try pfparticle
1050  // PFParticles
1051  art::Handle<std::vector<recob::PFParticle>> pfpHandle;
1052  std::vector<art::Ptr<recob::PFParticle>> pfps;
1053  if (event.getByLabel(fShowerModuleLabel, pfpHandle)) art::fill_ptr_vector(pfps, pfpHandle);
1054  // Clusters
1055  art::Handle<std::vector<recob::Cluster>> clusterHandle;
1056  std::vector<art::Ptr<recob::Cluster>> clusters;
1057  if (event.getByLabel(fShowerModuleLabel, clusterHandle))
1058  art::fill_ptr_vector(clusters, clusterHandle);
1059  art::FindManyP<recob::PFParticle> fmps(showerHandle, event, fShowerModuleLabel);
1060  art::FindManyP<recob::Cluster> fmcp(pfpHandle, event, fShowerModuleLabel);
1061  art::FindManyP<recob::Hit> fmhc(clusterHandle, event, fShowerModuleLabel);
1062  if (fmps.isValid()) {
1063  std::vector<art::Ptr<recob::PFParticle>> pfs = fmps.at(i);
1064  for (size_t ipf = 0; ipf < pfs.size(); ++ipf) {
1065  if (fmcp.isValid()) {
1066  std::vector<art::Ptr<recob::Cluster>> clus = fmcp.at(pfs[ipf].key());
1067  for (size_t iclu = 0; iclu < clus.size(); ++iclu) {
1068  if (fmhc.isValid()) {
1069  std::vector<art::Ptr<recob::Hit>> hits = fmhc.at(clus[iclu].key());
1070  for (size_t ihit = 0; ihit < hits.size(); ++ihit) {
1071  sh_hits.push_back(hits[ihit]);
1072  }
1073  }
1074  }
1075  }
1076  }
1077  }
1078  }
1079 
1080  const simb::MCParticle* particle;
1081  double tmpEfrac_contamination =
1082  0.0; // fraction of non EM energy contatiminatio (see truthMatcher for
1083  // definition)
1084  double tmpEcomplet = 0;
1085 
1086  int tmp_nHits = sh_hits.size();
1087 
1088  truthMatcher(clockData, all_hits, sh_hits, particle, tmpEfrac_contamination, tmpEcomplet);
1089  if (!particle) continue;
1090 
1091  sh_Efrac_contamination[i] = tmpEfrac_contamination;
1092  sh_purity[i] = 1 - tmpEfrac_contamination;
1093  sh_completeness[i] = tmpEcomplet;
1094  sh_nHits[i] = tmp_nHits;
1095  sh_hasPrimary_e[i] = 0;
1096  sh_pdg[i] = particle->PdgCode();
1097 
1098  //Shower with highest hits
1099  if (tmp_nHits > nHits) {
1100  sh_largest = i;
1101  dEdxasymm_largestshw = dEdxasymm;
1102  nHits = tmp_nHits;
1103  Ecomplet_NueCC = tmpEcomplet;
1104  Efrac_contaminationNueCC = tmpEfrac_contamination;
1105  //Calculate Shower anagle w.r.t True particle
1106  double ShDirMag =
1107  sqrt(pow(sh_direction_X[i], 2) + pow(sh_direction_Y[i], 2) + pow(sh_direction_Z[i], 2));
1108  ShAngle = (sh_direction_X[i] * particle->Px() + sh_direction_Y[i] * particle->Py() +
1109  sh_direction_Z[i] * particle->Pz()) /
1110  (ShDirMag * particle->P());
1111 
1112  ShVxTrueParticleVxDiff = sh_start_X[i] - particle->Vx();
1113  ShVyTrueParticleVyDiff = sh_start_Y[i] - particle->Vy();
1114  ShVzTrueParticleVzDiff = sh_start_Z[i] - particle->Vz();
1115 
1116  //put overflow and underflow at top and bottom bins:
1117  if (ShVxTrueParticleVxDiff > 5)
1118  ShVxTrueParticleVxDiff = 4.99;
1119  else if (ShVxTrueParticleVxDiff < -5)
1120  ShVxTrueParticleVxDiff = -5;
1121  if (ShVyTrueParticleVyDiff > 5)
1122  ShVyTrueParticleVyDiff = 4.99;
1123  else if (ShVyTrueParticleVyDiff < -5)
1124  ShVyTrueParticleVyDiff = -5;
1125  if (ShVzTrueParticleVzDiff > 5)
1126  ShVzTrueParticleVzDiff = 4.99;
1127  else if (ShVzTrueParticleVzDiff < -5)
1128  ShVzTrueParticleVzDiff = -5;
1129 
1130  ShStartVxTrueParticleEndVxDiff = sh_start_X[i] - particle->EndX();
1131  ShStartVyTrueParticleEndVyDiff = sh_start_Y[i] - particle->EndY();
1132  ShStartVzTrueParticleEndVzDiff = sh_start_Z[i] - particle->EndZ();
1133 
1134  if (std::abs(particle->PdgCode()) == 11) { ParticlePDG_HighestShHits = 1; }
1135  else if (particle->PdgCode() == 22) {
1136  ParticlePDG_HighestShHits = 2;
1137  }
1138  else {
1139  ParticlePDG_HighestShHits = 3;
1140  }
1141 
1142  //dedx for different showers
1143  //Highest hits shower pdg for the dEdx study 0=undefined,1=electronorpositronshower,2=photonshower,3=protonshower,4=neutronshower,5=chargedpionshower,6=neutralpionshower,7=everythingelseshower
1144  shower_bestplane = shower->best_plane();
1145  if (shower_bestplane < 0 || shower_bestplane >= int(shower->dEdx().size())) {
1146  //bestplane is not set properly, just pick the first plane that has dEdx
1147  for (size_t i = 0; i < shower->dEdx().size(); ++i) {
1148  if (shower->dEdx()[i]) {
1149  shower_bestplane = i;
1150  break;
1151  }
1152  }
1153  }
1154  if (shower_bestplane < 0 || shower_bestplane >= int(shower->dEdx().size())) {
1155  //still a problem? just set it to 0
1156  shower_bestplane = 0;
1157  }
1158 
1159  if (shower_bestplane >= 0 and shower_bestplane < int(shower->dEdx().size()))
1160  Showerparticlededx_inbestplane = shower->dEdx()[shower_bestplane];
1161 
1162  if (std::abs(particle->PdgCode()) == 11) { //lepton shower
1163  showerPDGwithHighestHitsforFillingdEdX = 1;
1164  }
1165  else if (particle->PdgCode() == 22) { //photon shower
1166  showerPDGwithHighestHitsforFillingdEdX = 2;
1167  }
1168  else if (particle->PdgCode() == 2212) { //proton shower
1169  showerPDGwithHighestHitsforFillingdEdX = 3;
1170  }
1171  else if (particle->PdgCode() == 2112) { //neutron shower
1172  showerPDGwithHighestHitsforFillingdEdX = 4;
1173  }
1174  else if (std::abs(particle->PdgCode()) == 211) { //charged pion shower
1175  showerPDGwithHighestHitsforFillingdEdX = 5;
1176  }
1177  else if (particle->PdgCode() == 111) { //neutral pion shower
1178  showerPDGwithHighestHitsforFillingdEdX = 6;
1179  }
1180  else { //everythingelse shower
1181  showerPDGwithHighestHitsforFillingdEdX = 7;
1182  }
1183  }
1184 
1185  if (particle->PdgCode() == fLeptonPDGcode && particle->TrackId() == MC_leptonID)
1186  sh_hasPrimary_e[i] = 1;
1187  // save the best shower based on non EM and number of hits
1188 
1189  if (std::abs(particle->PdgCode()) == fLeptonPDGcode && particle->TrackId() == MC_leptonID) {
1190 
1191  if (tmpEcomplet > Ecomplet_lepton) {
1192 
1193  Ecomplet_lepton = tmpEcomplet;
1194 
1195  Efrac_contamination = tmpEfrac_contamination;
1196  MClepton_reco = particle;
1197  sh_Efrac_best = Efrac_contamination;
1198  }
1199  }
1200  } //end of looping all the showers
1201 
1202  if (p1.Mag() && p2.Mag()) { sh_mpi0 = sqrt(pow(p1.Mag() + p2.Mag(), 2) - (p1 + p2).Mag2()); }
1203  else
1204  sh_mpi0 = 0;
1205 
1206  if (MClepton_reco && MClepton) {
1207  if (MC_isCC && (fNeutrinoPDGcode == std::abs(MC_incoming_PDG))) {
1208  h_Efrac_shContamination->Fill(Efrac_contamination);
1209  h_Efrac_shPurity->Fill(1 - Efrac_contamination);
1210  h_Ecomplet_lepton->Fill(Ecomplet_lepton);
1211 
1212  // Selecting good showers requires completeness of gretaer than 70 % and Purity > 70 %
1213  if (Efrac_contamination < fMaxEfrac && Ecomplet_lepton > fMinCompleteness) {
1214 
1215  h_Pe_num->Fill(Pe);
1216  h_Ev_num->Fill(MC_incoming_P[3]);
1217  h_Ee_num->Fill(MC_lepton_startMomentum[3]);
1218  h_theta_num->Fill(theta_e);
1219 
1220  if (Showerparticlededx_inbestplane > 1 && Showerparticlededx_inbestplane < 3) {
1221  h_Ev_num_dEdx->Fill(MC_incoming_P[3]);
1222  h_Ee_num_dEdx->Fill(MC_lepton_startMomentum[3]);
1223  }
1224  }
1225  }
1226  }
1227 
1228  //NueCC SIgnal and background Completeness
1229  if (MC_isCC == 1 && (fNeutrinoPDGcode == std::abs(MC_incoming_PDG)) && isFiducial) {
1230  h_HighestHitsProducedParticlePDG_NueCC->Fill(ParticlePDG_HighestShHits);
1231 
1232  if (ParticlePDG_HighestShHits > 0) { // atleat one shower is reconstructed
1233  h_Ecomplet_NueCC->Fill(Ecomplet_NueCC);
1234  h_Efrac_NueCCPurity->Fill(1 - Efrac_contaminationNueCC);
1235 
1236  h_esh_bestplane_NueCC->Fill(shower_bestplane);
1237  if (showerPDGwithHighestHitsforFillingdEdX == 1) //electron or positron shower
1238  {
1239  h_dEdX_electronorpositron_NueCC->Fill(Showerparticlededx_inbestplane);
1240  //Study the angle between the reconstructed shower direction w.r.t MC true particle direction
1241  h_CosThetaShDirwrtTrueparticle_electronorpositron_NueCC->Fill(ShAngle);
1242 
1243  //Study the reconstructed shower start position (x,y,z) w.r.t MC true particle start position
1244  h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NueCC->Fill(
1245  ShVxTrueParticleVxDiff);
1246  h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NueCC->Fill(
1247  ShVyTrueParticleVyDiff);
1248  h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NueCC->Fill(
1249  ShVzTrueParticleVzDiff);
1250 
1251  h_dEdXasymm_electronorpositron_NueCC->Fill(dEdxasymm_largestshw);
1252 
1253  h_mpi0_electronorpositron_NueCC->Fill(sh_mpi0);
1254  }
1255  else if (showerPDGwithHighestHitsforFillingdEdX == 2) //photon shower
1256  {
1257  h_dEdX_photon_NueCC->Fill(Showerparticlededx_inbestplane);
1258  h_CosThetaShDirwrtTrueparticle_photon_NueCC->Fill(ShAngle);
1259  h_ShStartXwrtTrueparticleStartXDiff_photon_NueCC->Fill(ShVxTrueParticleVxDiff);
1260  h_ShStartYwrtTrueparticleStartYDiff_photon_NueCC->Fill(ShVyTrueParticleVyDiff);
1261  h_ShStartZwrtTrueparticleStartZDiff_photon_NueCC->Fill(ShVzTrueParticleVzDiff);
1262 
1263  h_ShStartXwrtTrueparticleEndXDiff_photon_NueCC->Fill(ShStartVxTrueParticleEndVxDiff);
1264  h_ShStartYwrtTrueparticleEndYDiff_photon_NueCC->Fill(ShStartVyTrueParticleEndVyDiff);
1265  h_ShStartZwrtTrueparticleEndZDiff_photon_NueCC->Fill(ShStartVzTrueParticleEndVzDiff);
1266  }
1267  else if (showerPDGwithHighestHitsforFillingdEdX == 3) //proton shower
1268  {
1269  h_dEdX_proton_NueCC->Fill(Showerparticlededx_inbestplane);
1270  h_CosThetaShDirwrtTrueparticle_proton_NueCC->Fill(ShAngle);
1271 
1272  h_ShStartXwrtTrueparticleStartXDiff_proton_NueCC->Fill(ShVxTrueParticleVxDiff);
1273  h_ShStartYwrtTrueparticleStartYDiff_proton_NueCC->Fill(ShVyTrueParticleVyDiff);
1274  h_ShStartZwrtTrueparticleStartZDiff_proton_NueCC->Fill(ShVzTrueParticleVzDiff);
1275  }
1276  else if (showerPDGwithHighestHitsforFillingdEdX == 4) //neutron shower
1277  {
1278  h_dEdX_neutron_NueCC->Fill(Showerparticlededx_inbestplane);
1279  }
1280  else if (showerPDGwithHighestHitsforFillingdEdX == 5) //charged pion shower
1281  {
1282  h_dEdX_chargedpion_NueCC->Fill(Showerparticlededx_inbestplane);
1283  h_CosThetaShDirwrtTrueparticle_chargedpion_NueCC->Fill(ShAngle);
1284  h_ShStartXwrtTrueparticleStartXDiff_chargedpion_NueCC->Fill(ShVxTrueParticleVxDiff);
1285  h_ShStartYwrtTrueparticleStartYDiff_chargedpion_NueCC->Fill(ShVyTrueParticleVyDiff);
1286  h_ShStartZwrtTrueparticleStartZDiff_chargedpion_NueCC->Fill(ShVzTrueParticleVzDiff);
1287  }
1288  else if (showerPDGwithHighestHitsforFillingdEdX == 6) //neutral pion shower
1289  {
1290  h_dEdX_neutralpion_NueCC->Fill(Showerparticlededx_inbestplane);
1291  }
1292  else if (showerPDGwithHighestHitsforFillingdEdX == 7) //everythingelse shower
1293  {
1294  h_dEdX_everythingelse_NueCC->Fill(Showerparticlededx_inbestplane);
1295  }
1296  }
1297  }
1298  else if (!MC_isCC && isFiducial) {
1299  h_HighestHitsProducedParticlePDG_bkg->Fill(ParticlePDG_HighestShHits);
1300 
1301  if (ParticlePDG_HighestShHits > 0) {
1302  h_Ecomplet_bkg->Fill(Ecomplet_NueCC);
1303  h_Efrac_bkgPurity->Fill(1 - Efrac_contaminationNueCC);
1304 
1305  h_esh_bestplane_NC->Fill(shower_bestplane);
1306  if (showerPDGwithHighestHitsforFillingdEdX == 1) //electron or positron shower
1307  {
1308  h_dEdX_electronorpositron_NC->Fill(Showerparticlededx_inbestplane);
1309  h_CosThetaShDirwrtTrueparticle_electronorpositron_NC->Fill(ShAngle);
1310 
1311  h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NC->Fill(ShVxTrueParticleVxDiff);
1312  h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NC->Fill(ShVyTrueParticleVyDiff);
1313  h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NC->Fill(ShVzTrueParticleVzDiff);
1314  }
1315  else if (showerPDGwithHighestHitsforFillingdEdX == 2) //photon shower
1316  {
1317  h_dEdX_photon_NC->Fill(Showerparticlededx_inbestplane);
1318  h_CosThetaShDirwrtTrueparticle_photon_NC->Fill(ShAngle);
1319 
1320  h_ShStartXwrtTrueparticleStartXDiff_photon_NC->Fill(ShVxTrueParticleVxDiff);
1321  h_ShStartYwrtTrueparticleStartYDiff_photon_NC->Fill(ShVyTrueParticleVyDiff);
1322  h_ShStartZwrtTrueparticleStartZDiff_photon_NC->Fill(ShVzTrueParticleVzDiff);
1323 
1324  h_ShStartXwrtTrueparticleEndXDiff_photon_NC->Fill(ShStartVxTrueParticleEndVxDiff);
1325  h_ShStartYwrtTrueparticleEndYDiff_photon_NC->Fill(ShStartVyTrueParticleEndVyDiff);
1326  h_ShStartZwrtTrueparticleEndZDiff_photon_NC->Fill(ShStartVzTrueParticleEndVzDiff);
1327 
1328  h_dEdXasymm_photon_NC->Fill(dEdxasymm_largestshw);
1329 
1330  h_mpi0_photon_NC->Fill(sh_mpi0);
1331  }
1332  else if (showerPDGwithHighestHitsforFillingdEdX == 3) //proton shower
1333  {
1334  h_dEdX_proton_NC->Fill(Showerparticlededx_inbestplane);
1335  h_CosThetaShDirwrtTrueparticle_proton_NC->Fill(ShAngle);
1336 
1337  h_ShStartXwrtTrueparticleStartXDiff_proton_NC->Fill(ShVxTrueParticleVxDiff);
1338  h_ShStartYwrtTrueparticleStartYDiff_proton_NC->Fill(ShVyTrueParticleVyDiff);
1339  h_ShStartZwrtTrueparticleStartZDiff_proton_NC->Fill(ShVzTrueParticleVzDiff);
1340  }
1341  else if (showerPDGwithHighestHitsforFillingdEdX == 4) //neutron shower
1342  {
1343  h_dEdX_neutron_NC->Fill(Showerparticlededx_inbestplane);
1344  }
1345  else if (showerPDGwithHighestHitsforFillingdEdX == 5) //charged pion shower
1346  {
1347  h_dEdX_chargedpion_NC->Fill(Showerparticlededx_inbestplane);
1348  h_CosThetaShDirwrtTrueparticle_chargedpion_NC->Fill(ShAngle);
1349  h_ShStartXwrtTrueparticleStartXDiff_chargedpion_NC->Fill(ShVxTrueParticleVxDiff);
1350  h_ShStartYwrtTrueparticleStartYDiff_chargedpion_NC->Fill(ShVyTrueParticleVyDiff);
1351  h_ShStartZwrtTrueparticleStartZDiff_chargedpion_NC->Fill(ShVzTrueParticleVzDiff);
1352  }
1353  else if (showerPDGwithHighestHitsforFillingdEdX == 6) //neutral pion shower
1354  {
1355  h_dEdX_neutralpion_NC->Fill(Showerparticlededx_inbestplane);
1356  }
1357  else if (showerPDGwithHighestHitsforFillingdEdX == 7) //everythingelse shower
1358  {
1359  h_dEdX_everythingelse_NC->Fill(Showerparticlededx_inbestplane);
1360  }
1361  } //if(ParticlePDG_HighestShHits>0)
1362  } //else if(!MC_isCC&&isFiducial)
1363 
1364  checkCNNtrkshw<4>(clockData, event, all_hits);
1365  }
1366 
1367  //========================================================================
1368  void
1369  NeutrinoShowerEff::truthMatcher(detinfo::DetectorClocksData const& clockData,
1370  std::vector<art::Ptr<recob::Hit>> all_hits,
1371  std::vector<art::Ptr<recob::Hit>> shower_hits,
1372  const simb::MCParticle*& MCparticle,
1373  double& Efrac,
1374  double& Ecomplet)
1375  {
1376 
1377  MCparticle = 0;
1378  Efrac = 1.0;
1379  Ecomplet = 0;
1380 
1381  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
1382  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
1383  std::map<int, double> trkID_E;
1384  for (size_t j = 0; j < shower_hits.size(); ++j) {
1385  art::Ptr<recob::Hit> hit = shower_hits[j];
1386  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToEveTrackIDEs(clockData, hit);
1387  for (size_t k = 0; k < TrackIDs.size(); k++) {
1388  if (trkID_E.find(std::abs(TrackIDs[k].trackID)) == trkID_E.end())
1389  trkID_E[std::abs(TrackIDs[k].trackID)] = 0;
1390  trkID_E[std::abs(TrackIDs[k].trackID)] += TrackIDs[k].energy;
1391  }
1392  }
1393  double max_E = -999.0;
1394  double total_E = 0.0;
1395  int TrackID = -999;
1396  double partial_E = 0.0;
1397 
1398  if (empty(trkID_E)) return; // Ghost shower???
1399 
1400  for (std::map<int, double>::iterator ii = trkID_E.begin(); ii != trkID_E.end(); ++ii) {
1401  total_E += ii->second;
1402  if ((ii->second) > max_E) {
1403  partial_E = ii->second;
1404  max_E = ii->second;
1405  TrackID = ii->first;
1406  }
1407  }
1408  MCparticle = pi_serv->TrackIdToParticle_P(TrackID);
1409 
1410  Efrac = 1 - (partial_E / total_E);
1411 
1412  //completeness
1413  double totenergy = 0;
1414  for (size_t k = 0; k < all_hits.size(); ++k) {
1415  art::Ptr<recob::Hit> hit = all_hits[k];
1416  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToEveTrackIDEs(clockData, hit);
1417  for (size_t l = 0; l < TrackIDs.size(); ++l) {
1418  if (std::abs(TrackIDs[l].trackID) == TrackID) { totenergy += TrackIDs[l].energy; }
1419  }
1420  }
1421  Ecomplet = partial_E / totenergy;
1422  }
1423 
1424  //========================================================================
1425  bool
1426  NeutrinoShowerEff::insideFV(double vertex[4])
1427  {
1428 
1429  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1430  //This is temporarily we should define a common FV
1431  double x = vertex[0];
1432  double y = vertex[1];
1433  double z = vertex[2];
1434 
1435  if (x > fFidVolXmin && x < fFidVolXmax && y > fFidVolYmin && y < fFidVolYmax &&
1436  z > fFidVolZmin && z < fFidVolZmax)
1437  return true;
1438  else
1439  return false;
1440  }
1441  //========================================================================
1442  void
1443  NeutrinoShowerEff::doEfficiencies()
1444  {
1445 
1446  art::ServiceHandle<art::TFileService const> tfs;
1447 
1448  if (TEfficiency::CheckConsistency(*h_Ev_num, *h_Ev_den)) {
1449  h_Eff_Ev = tfs->make<TEfficiency>(*h_Ev_num, *h_Ev_den);
1450  TGraphAsymmErrors* grEff_Ev = h_Eff_Ev->CreateGraph();
1451  grEff_Ev->Write("grEff_Ev");
1452  h_Eff_Ev->Write("h_Eff_Ev");
1453  }
1454 
1455  if (TEfficiency::CheckConsistency(*h_Ev_num_dEdx, *h_Ev_den)) {
1456  h_Eff_Ev_dEdx = tfs->make<TEfficiency>(*h_Ev_num_dEdx, *h_Ev_den);
1457  TGraphAsymmErrors* grEff_Ev_dEdx = h_Eff_Ev_dEdx->CreateGraph();
1458  grEff_Ev_dEdx->Write("grEff_Ev_dEdx");
1459  h_Eff_Ev_dEdx->Write("h_Eff_Ev_dEdx");
1460  }
1461 
1462  if (TEfficiency::CheckConsistency(*h_Ee_num, *h_Ee_den)) {
1463  h_Eff_Ee = tfs->make<TEfficiency>(*h_Ee_num, *h_Ee_den);
1464  TGraphAsymmErrors* grEff_Ee = h_Eff_Ee->CreateGraph();
1465  grEff_Ee->Write("grEff_Ee");
1466  h_Eff_Ee->Write("h_Eff_Ee");
1467  }
1468 
1469  if (TEfficiency::CheckConsistency(*h_Ee_num_dEdx, *h_Ee_den)) {
1470  h_Eff_Ee_dEdx = tfs->make<TEfficiency>(*h_Ee_num_dEdx, *h_Ee_den);
1471  TGraphAsymmErrors* grEff_Ee_dEdx = h_Eff_Ee_dEdx->CreateGraph();
1472  grEff_Ee_dEdx->Write("grEff_Ee_dEdx");
1473  h_Eff_Ee_dEdx->Write("h_Eff_Ee_dEdx");
1474  }
1475 
1476  if (TEfficiency::CheckConsistency(*h_Pe_num, *h_Pe_den)) {
1477  h_Eff_Pe = tfs->make<TEfficiency>(*h_Pe_num, *h_Pe_den);
1478  TGraphAsymmErrors* grEff_Pe = h_Eff_Pe->CreateGraph();
1479  grEff_Pe->Write("grEff_Pe");
1480  h_Eff_Pe->Write("h_Eff_Pe");
1481  }
1482  if (TEfficiency::CheckConsistency(*h_theta_num, *h_theta_den)) {
1483  h_Eff_theta = tfs->make<TEfficiency>(*h_theta_num, *h_theta_den);
1484  TGraphAsymmErrors* grEff_theta = h_Eff_theta->CreateGraph();
1485  grEff_theta->Write("grEff_theta");
1486  h_Eff_theta->Write("h_Eff_theta");
1487  }
1488  }
1489 
1490  //============================================
1491  // Check CNN track/shower ID
1492  //============================================
1493  template <size_t N>
1494  void
1495  NeutrinoShowerEff::checkCNNtrkshw(detinfo::DetectorClocksData const& clockData,
1496  const art::Event& evt,
1497  std::vector<art::Ptr<recob::Hit>> all_hits)
1498  {
1499  if (fCNNEMModuleLabel.empty()) return;
1500 
1501  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
1502  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
1503 
1504  auto hitResults = anab::MVAReader<recob::Hit, N>::create(evt, fCNNEMModuleLabel);
1505  if (hitResults) {
1506  int trkLikeIdx = hitResults->getIndex("track");
1507  int emLikeIdx = hitResults->getIndex("em");
1508  if ((trkLikeIdx < 0) || (emLikeIdx < 0)) {
1509  throw cet::exception("NeutrinoShowerEff")
1510  << "No em/track labeled columns in MVA data products." << std::endl;
1511  }
1512 
1513  for (size_t i = 0; i < all_hits.size(); ++i) {
1514  //find out if the hit was generated by an EM particle
1515  bool isEMparticle = false;
1516  int pdg = INT_MAX;
1517  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToEveTrackIDEs(clockData, all_hits[i]);
1518  if (!TrackIDs.size()) continue;
1519 
1520  int trkid = INT_MAX;
1521  double maxE = -1;
1522  for (size_t k = 0; k < TrackIDs.size(); k++) {
1523  if (TrackIDs[k].energy > maxE) {
1524  maxE = TrackIDs[k].energy;
1525  trkid = TrackIDs[k].trackID;
1526  }
1527  }
1528  if (trkid != INT_MAX) {
1529  auto* particle = pi_serv->TrackIdToParticle_P(trkid);
1530  if (particle) {
1531  pdg = particle->PdgCode();
1532  if (std::abs(pdg) == 11 || //electron/positron
1533  pdg == 22 || //photon
1534  pdg == 111) { //pi0
1535  isEMparticle = true;
1536  }
1537  }
1538  }
1539  auto vout = hitResults->getOutput(all_hits[i]);
1540  double trk_like = -1, trk_or_em = vout[trkLikeIdx] + vout[emLikeIdx];
1541  if (trk_or_em > 0) {
1542  trk_like = vout[trkLikeIdx] / trk_or_em;
1543  if (isEMparticle) { h_trklike_em->Fill(trk_like); }
1544  else {
1545  h_trklike_nonem->Fill(trk_like);
1546  }
1547  }
1548  }
1549  }
1550  else {
1551  std::cout << "Couldn't get hitResults." << std::endl;
1552  }
1553  }
1554 
1555  //========================================================================
1556  void
1557  NeutrinoShowerEff::reset()
1558  {
1559 
1560  MC_incoming_PDG = -999;
1561  MC_lepton_PDG = -999;
1562  MC_isCC = -999;
1563  MC_channel = -999;
1564  MC_target = -999;
1565  MC_Q2 = -999.0;
1566  MC_W = -999.0;
1567  MC_lepton_theta = -999.0;
1568  MC_leptonID = -999;
1569  MC_LeptonTrack = -999;
1570 
1571  for (int i = 0; i < MAX_SHOWERS; i++) {
1572  sh_direction_X[i] = -999.0;
1573  sh_direction_Y[i] = -999.0;
1574  sh_direction_Z[i] = -999.0;
1575  sh_start_X[i] = -999.0;
1576  sh_start_Y[i] = -999.0;
1577  sh_start_Z[i] = -999.0;
1578  sh_bestplane[i] = -999.0;
1579  sh_length[i] = -999.0;
1580  sh_hasPrimary_e[i] = -999.0;
1581  sh_Efrac_contamination[i] = -999.0;
1582  sh_purity[i] = -999.0;
1583  sh_completeness[i] = -999.0;
1584  sh_nHits[i] = -999.0;
1585  for (int j = 0; j < 3; j++) {
1586  sh_energy[i][j] = -999.0;
1587  sh_MIPenergy[i][j] = -999.0;
1588  sh_dEdx[i][j] = -999.0;
1589  }
1590  sh_pdg[i] = -999;
1591  sh_dEdxasymm[i] = -999;
1592  }
1593  sh_largest = -999;
1594  sh_mpi0 = -999;
1595  }
1596 
1597  //========================================================================
1598  DEFINE_ART_MODULE(NeutrinoShowerEff)
1599 
1600 }
rootstat.analyze
def analyze(root, level, gtrees, gbranches, doprint)
Definition: rootstat.py:69
recob::Shower::best_plane
int best_plane() const
Definition: Shower.h:200
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
MAX_SHOWERS
#define MAX_SHOWERS
Definition: NeutrinoShowerEff_module.cc:29
DUNE::NeutrinoShowerEff::h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NC
TH1D * h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NC
Definition: NeutrinoShowerEff_module.cc:155
DUNE::NeutrinoShowerEff::checkCNNtrkshw
void checkCNNtrkshw(detinfo::DetectorClocksData const &clockData, const art::Event &evt, std::vector< art::Ptr< recob::Hit >> all_hits)
Definition: NeutrinoShowerEff_module.cc:1495
recob::Shower::ShowerStart
const TVector3 & ShowerStart() const
Definition: Shower.h:192
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::MCNeutrino::QSqr
double QSqr() const
Definition: MCNeutrino.h:157
simb::MCTruth::GetNeutrino
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
generate_CCQE_events.energy
string energy
Definition: generate_CCQE_events.py:15
simb::MCParticle::Py
double Py(const int i=0) const
Definition: MCParticle.h:231
DUNE::NeutrinoShowerEff::sh_energy
double sh_energy[MAX_SHOWERS][3]
Definition: NeutrinoShowerEff_module.cc:221
simb::MCParticle::EndPosition
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:225
simb::MCParticle::EndZ
double EndZ() const
Definition: MCParticle.h:228
recob::Shower::Length
double Length() const
Definition: Shower.h:201
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
DUNE::NeutrinoShowerEff::h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NC
TH1D * h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NC
Definition: NeutrinoShowerEff_module.cc:156
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
recob::Shower::Energy
const std::vector< double > & Energy() const
Definition: Shower.h:195
simb::MCParticle::Px
double Px(const int i=0) const
Definition: MCParticle.h:230
DUNE::NeutrinoShowerEff::sh_dEdx
double sh_dEdx[MAX_SHOWERS][3]
Definition: NeutrinoShowerEff_module.cc:223
vector
struct vector vector
DUNE::NeutrinoShowerEff::analyze
void analyze(const art::Event &evt) override
Definition: NeutrinoShowerEff_module.cc:841
std
STL namespace.
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
DUNE::NeutrinoShowerEff::truthMatcher
void truthMatcher(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> all_hits, std::vector< art::Ptr< recob::Hit >> shower_hits, const simb::MCParticle *&MCparticle, double &Efrac, double &Ecomplet)
Definition: NeutrinoShowerEff_module.cc:1369
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:211
MCParticle.h
Particle class.
simb::MCParticle::EndY
double EndY() const
Definition: MCParticle.h:227
DUNE::NeutrinoShowerEff::h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NueCC
TH1D * h_ShStartXwrtTrueparticleStartXDiff_electronorpositron_NueCC
Definition: NeutrinoShowerEff_module.cc:151
l
static QStrList * l
Definition: config.cpp:1044
art::InputTag::encode
std::string encode() const
Definition: InputTag.cc:97
art::Run
Definition: Run.h:17
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:210
art::InputTag::empty
bool empty() const noexcept
Definition: InputTag.cc:73
DUNE::NeutrinoShowerEff::h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NueCC
TH1D * h_ShStartYwrtTrueparticleStartYDiff_electronorpositron_NueCC
Definition: NeutrinoShowerEff_module.cc:152
simb::MCNeutrino::InteractionType
int InteractionType() const
Definition: MCNeutrino.h:150
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
reset
QTextStream & reset(QTextStream &s)
Definition: qtextstream.cpp:2048
simb::MCNeutrino::Lepton
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
simb::MCNeutrino::W
double W() const
Definition: MCNeutrino.h:154
recob::Shower::dEdx
const std::vector< double > & dEdx() const
Definition: Shower.h:203
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
recob::Shower::MIPEnergy
const std::vector< double > & MIPEnergy() const
Definition: Shower.h:198
breakpoints::beginJob
void beginJob()
Definition: Breakpoints.cc:14
wirecell.dfp.graph.key
def key(type, name=None)
Definition: graph.py:13
simb::MCParticle::P
double P(const int i=0) const
Definition: MCParticle.h:234
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:271
test.p
p
Definition: test.py:223
recob::Shower::Direction
const TVector3 & Direction() const
Definition: Shower.h:189
TrackID
TrackID
Definition: ClusteringValidation_module.cc:57
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
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
DUNE::NeutrinoShowerEff::h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NueCC
TH1D * h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NueCC
Definition: NeutrinoShowerEff_module.cc:153
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:221
simb::MCNeutrino::Target
int Target() const
Definition: MCNeutrino.h:151
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:283
DUNE::NeutrinoShowerEff::processEff
void processEff(detinfo::DetectorClocksData const &clockData, const art::Event &evt, bool &isFiducial)
Definition: NeutrinoShowerEff_module.cc:861
cheat::BackTrackerService::HitToEveTrackIDEs
std::vector< sim::TrackIDE > HitToEveTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Definition: BackTrackerService_service.cc:203
simb::MCParticle::Momentum
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
simb::MCParticle::Pz
double Pz(const int i=0) const
Definition: MCParticle.h:232
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:223
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
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
DUNE::NeutrinoShowerEff::beginRun
void beginRun(const art::Run &run) override
Definition: NeutrinoShowerEff_module.cc:835
simb::MCParticle::EndX
double EndX() const
Definition: MCParticle.h:226
DUNE::NeutrinoShowerEff::sh_Efrac_contamination
double sh_Efrac_contamination[MAX_SHOWERS]
Definition: NeutrinoShowerEff_module.cc:227
DUNE::NeutrinoShowerEff::h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NC
TH1D * h_ShStartZwrtTrueparticleStartZDiff_electronorpositron_NC
Definition: NeutrinoShowerEff_module.cc:157
simb::MCNeutrino
Event generator information.
Definition: MCNeutrino.h:18
geo
LArSoft geometry interface.
Definition: ChannelGeo.h:16
DUNE::NeutrinoShowerEff::sh_MIPenergy
double sh_MIPenergy[MAX_SHOWERS][3]
Definition: NeutrinoShowerEff_module.cc:222
simb::MCParticle::EndMomentum
const TLorentzVector & EndMomentum() const
Definition: MCParticle.h:240
geo::TPCGeo::LocalToWorld
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
DUNE::NeutrinoShowerEff::geom
art::ServiceHandle< geo::Geometry const > geom
Definition: NeutrinoShowerEff_module.cc:250
anab::MVAReader::create
static std::unique_ptr< MVAReader > create(const art::Event &evt, const art::InputTag &tag)
Definition: MVAReader.h:110
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:222
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
util::empty
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:97
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
event
Event finding and building.
Definition: EventCheater_module.cc:32
DUNE::NeutrinoShowerEff::h_CosThetaShDirwrtTrueparticle_electronorpositron_NueCC
TH1D * h_CosThetaShDirwrtTrueparticle_electronorpositron_NueCC
Definition: NeutrinoShowerEff_module.cc:141
vertex
vertex reconstruction
Definition: AggregateVertexAna_module.cc:32
Generated by   doxygen 1.8.11