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