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