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