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