LArSoft  v08_35_01
Liquid Argon Software toolkit - http://larsoft.org/
NuShowerEff_module.cc
Go to the documentation of this file.
1 // Class: NuShowerEff
3 // Plugin Type: analyzer (art v2_11_03)
4 // File: NuShowerEff_module.cc
5 //
6 // Generated at Tue Sep 18 14:20:57 2018 by Wanwei Wu using cetskelgen
7 // from cetlib version v3_03_01.
9 
10 // Framework includes
15 #include "fhiclcpp/ParameterSet.h"
17 
20 #include "art_root_io/TFileService.h"
21 
22 // LArsoft includes
28 
33 
34 // ROOT includes
35 #include "TTree.h"
36 #include "TH1.h"
37 #include "TEfficiency.h"
38 #include "TGraphAsymmErrors.h"
39 
40 // user's defined constant
41 #define MAX_SHOWERS 1000
42 
43 using namespace std;
44 
45 class NuShowerEff;
46 
47 
48 class NuShowerEff : public art::EDAnalyzer {
49 public:
50  explicit NuShowerEff(fhicl::ParameterSet const & p);
51  // The compiler-generated destructor is fine for non-base
52  // classes without bare pointers or other resource use.
53 
54  // Plugins should not be copied or assigned.
55  NuShowerEff(NuShowerEff const &) = delete;
56  NuShowerEff(NuShowerEff &&) = delete;
57  NuShowerEff & operator = (NuShowerEff const &) = delete;
58  NuShowerEff & operator = (NuShowerEff &&) = delete;
59 
60 private:
61 
62  // Required functions.
63  void analyze(art::Event const & e) override;
64 
65  // user's defined functions
66  void beginJob() override;
67  void endJob() override;
68  void beginRun(art::Run const& run) override;
69  void doEfficiencies();
70  bool insideFV(double vertex[4]);
71  void reset();
72 
73  // Declare member data here.
74 
75  // read from fchicl: Input Tag
76  std::string fMCTruthModuleLabel;
77  std::string fHitModuleLabel;
78  std::string fShowerModuleLabel;
80 
81  // read from fchicl: user's defined parameters
85  bool fHitShowerThroughPFParticle; // an option for getting hits associated with shower
86  double fMinPurity; // for reference
87  double fMinCompleteness; // for reference
88  float fFidVolCutX;// Fiducail Volume cut [cm]
89  float fFidVolCutY;
90  float fFidVolCutZ;
91 
92  // Fiducial Volume parameters
93  float fFidVolXmin;
94  float fFidVolXmax;
95  float fFidVolYmin;
96  float fFidVolYmax;
97  float fFidVolZmin;
98  float fFidVolZmax;
99 
100  // Event
101  int Event;
102  int Run;
103  int SubRun;
104 
105  // MC Truth: Generator
108  int MC_isCC;
111  double MC_Q2;
112  double MC_W;
113  double MC_vertex[4];
114  double MC_incoming_P[4];
115  double MC_lepton_startMomentum[4];
116  double MC_lepton_endMomentum[40];
117  double MC_lepton_startXYZT[4];
118  double MC_lepton_endXYZT[4];
122 
123  double MC_incoming_E; // incoming neutrino energy
125 
126  // recob::Shower
128  double sh_direction_X[MAX_SHOWERS];
129  double sh_direction_Y[MAX_SHOWERS];
130  double sh_direction_Z[MAX_SHOWERS];
131  double sh_start_X[MAX_SHOWERS];
132  double sh_start_Y[MAX_SHOWERS];
133  double sh_start_Z[MAX_SHOWERS];
134  double sh_length[MAX_SHOWERS];
135  double sh_ehit_Q[MAX_SHOWERS];
136  int sh_TrackId[MAX_SHOWERS];
137  int sh_hasPrimary_e[MAX_SHOWERS];
138  double sh_allhit_Q[MAX_SHOWERS];
139  double sh_purity[MAX_SHOWERS];
140  double sh_completeness[MAX_SHOWERS];
141  double esh_1_purity; // largest shower in a CC event that contains electron contributions
143  double esh_each_purity[MAX_SHOWERS]; // each shower in a CC event that contains electron contributions
144  double esh_each_completeness[MAX_SHOWERS];
145  double esh_purity; // average over all electron showers in a CC event
147  int count_primary_e_in_Event;// number of showers containing electron contribution in a CC event
148 
149  // TTree
150  TTree *fEventTree;
151 
152  TH1D *h_Ev_den; // incoming neutrino energy from MC. den: denominator
153  TH1D *h_Ev_num; // recon. num: numerator
154 
155  TH1D *h_Ee_den; // primary electron energy from MC
156  TH1D *h_Ee_num;
157 
158  TEfficiency* h_Eff_Ev = 0;
159  TEfficiency* h_Eff_Ee = 0;
160 
161 
162 };
163 
164 // =====================================================================================
166  :
167  EDAnalyzer(p), // ,
168  // More initializers here.
169  //fMCTruthModuleLabel (p.get< std::string >("MCTruthModuleLabel", "generator")),
170  fMCTruthModuleLabel (p.get< std::string >("MCTruthModuleLabel")),// get parameter from fcl file
171  fHitModuleLabel (p.get< std::string >("HitModuleLabel")),
172  fShowerModuleLabel (p.get< std::string >("ShowerModuleLabel")),
173  fTruthMatchDataModuleLabel (p.get< std::string >("TruthMatchDataModuleLabel")),
174  fNeutrinoPDGcode (p.get<int>("NeutrinoPDGcode")),
175  fLeptonPDGcode (p.get<int>("LeptonPDGcode")),
176  fSaveMCTree (p.get<bool>("SaveMCTree")),
177  fHitShowerThroughPFParticle (p.get<bool>("HitShowerThroughPFParticle")),
178  fMinPurity (p.get<double>("MinPurity")),
179  fMinCompleteness (p.get<double>("MinCompleteness")),
180  fFidVolCutX (p.get<float>("FidVolCutX")),
181  fFidVolCutY (p.get<float>("FidVolCutY")),
182  fFidVolCutZ (p.get<float>("FidVolCutZ"))
183 {
184  //cout << "\n===== Please refer the fchicl for the values of preset parameters ====\n" << endl;
185 }
186 
187 //============================================================================
189  //cout << "\n===== function: beginJob() ====\n" << endl;
190 
191  // Get geometry: Fiducial Volume
192  auto const* geo = lar::providerFrom<geo::Geometry>();
193  double minx = 1e9; // [cm]
194  double maxx = -1e9;
195  double miny = 1e9;
196  double maxy = -1e9;
197  double minz = 1e9;
198  double maxz = -1e9;
199  //cout << "\nGeometry:\n\tgeo->NTPC(): " << geo->NTPC() << endl;
200  for (size_t i = 0; i<geo->NTPC(); ++i){
201  double local[3] = {0.,0.,0.};
202  double world[3] = {0.,0.,0.};
203  const geo::TPCGeo &tpc = geo->TPC(i);
204  tpc.LocalToWorld(local,world);
205  //cout << "\tlocal: " << local[0] << " ; " << local[1] << " ; " << local[2] << endl;
206  //cout << "\tworld: " << world[0] << " ; " << world[1] << " ; " << world[2] << endl;
207  //cout << "\tgeo->DetHalfWidth(" << i << "): " << geo->DetHalfWidth(i) << endl;
208  //cout << "\tgeo->DetHalfHeight(" << i << "): " << geo->DetHalfHeight(i) << endl;
209  //cout << "\tgeo->DetLength(" << i << "): " << geo->DetLength(i) << endl;
210 
211  if (minx > world[0] - geo->DetHalfWidth(i)) minx = world[0]-geo->DetHalfWidth(i);
212  if (maxx < world[0] + geo->DetHalfWidth(i)) maxx = world[0]+geo->DetHalfWidth(i);
213  if (miny > world[1] - geo->DetHalfHeight(i)) miny = world[1]-geo->DetHalfHeight(i);
214  if (maxy < world[1] + geo->DetHalfHeight(i)) maxy = world[1]+geo->DetHalfHeight(i);
215  if (minz > world[2] - geo->DetLength(i)/2.) minz = world[2]-geo->DetLength(i)/2.;
216  if (maxz < world[2] + geo->DetLength(i)/2.) maxz = world[2]+geo->DetLength(i)/2.;
217  }
218 
219  fFidVolXmin = minx + fFidVolCutX;
220  fFidVolXmax = maxx - fFidVolCutX;
221  fFidVolYmin = miny + fFidVolCutY;
222  fFidVolYmax = maxy - fFidVolCutY;
223  fFidVolZmin = minz + fFidVolCutZ;
224  fFidVolZmax = maxz - fFidVolCutZ;
225 
226  //cout << "\nFiducial Volume (length unit: cm):\n"
227  //<< "\t" << fFidVolXmin<<"\t< x <\t"<<fFidVolXmax<<"\n"
228  //<< "\t" << fFidVolYmin<<"\t< y <\t"<<fFidVolYmax<<"\n"
229  //<< "\t" << fFidVolZmin<<"\t< z <\t"<<fFidVolZmax<<"\n";
230 
232 
233  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};
234 // double theta_bins[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.};
235  // 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};
236 
237  h_Ev_den = tfs->make<TH1D>("h_Ev_den", "Neutrino Energy; Neutrino Energy (GeV); Shower Reconstruction Efficiency", 20, E_bins);
238  h_Ev_den->Sumw2();
239  h_Ev_num = tfs->make<TH1D>("h_Ev_num","Neutrino Energy; Neutrino Energy (GeV); Shower Reconstruction Efficiency",20,E_bins);
240  h_Ev_num->Sumw2();
241 
242  h_Ee_den = tfs->make<TH1D>("h_Ee_den","Electron Energy; Electron Energy (GeV); Shower Reconstruction Efficiency",20,E_bins);
243  h_Ee_den->Sumw2();
244  h_Ee_num = tfs->make<TH1D>("h_Ee_num","Electron Energy; Electron Energy (GeV); Shower Reconstruction Efficiency",20,E_bins);
245  h_Ee_num->Sumw2();
246 
247  if (fSaveMCTree) {
248  fEventTree = new TTree("Event", "Event Tree from Sim & Reco");
249  fEventTree->Branch("eventNo", &Event); // only select event in FV
250  fEventTree->Branch("runNo", &Run);
251  fEventTree->Branch("subRunNo", &SubRun);
252 
253  fEventTree->Branch("MC_incoming_E", &MC_incoming_E);
254  fEventTree->Branch("MC_lepton_startEnergy", &MC_lepton_startEnergy);
255 
256  fEventTree->Branch("n_showers", &n_recoShowers);
257  fEventTree->Branch("sh_hasPrimary_e", &sh_hasPrimary_e, "sh_hasPrimary_e[n_showers]/I");
258  //fEventTree->Branch("sh_purity", &sh_purity, "sh_purity[n_showers]/D");
259  //fEventTree->Branch("sh_completeness", &sh_completeness, "sh_completeness[n_showers]/D");
260  fEventTree->Branch("count_primary_e_in_Event", &count_primary_e_in_Event);
261  fEventTree->Branch("esh_1_purity", &esh_1_purity);
262  fEventTree->Branch("esh_1_completeness", &esh_1_completeness);
263  fEventTree->Branch("esh_each_purity", &esh_each_purity, "esh_each_purity[count_primary_e_in_Event]/D");
264  fEventTree->Branch("esh_each_completeness", &esh_each_completeness, "esh_each_completeness[count_primary_e_in_Event]/D");
265  fEventTree->Branch("esh_purity", &esh_purity);
266  fEventTree->Branch("esh_completeness", &esh_completeness);
267  }
268 
269 }
270 
271 //============================================================================
273  //cout << "\n===== function: endJob() =====\n" << endl;
274  doEfficiencies();
275 }
276 
277 //============================================================================
278 void NuShowerEff::beginRun(art::Run const & run){
279  //cout << "\n===== function: beginRun() =====\n" << endl;
280  mf::LogInfo("ShowerEff") << "==== begin run ... ====" << endl;
281 }
282 
283 //============================================================================
285 {
286  // Implementation of required member function here.
287  //cout << "\n===== function: analyze() =====\n" << endl;
288 
289  reset(); // for some variables
290 
291  Event = e.id().event();
292  Run = e.run();
293  SubRun = e.subRun();
294  //cout << "Event information:" << endl;
295  //cout << "\tEvent: " << Event << endl;
296  //cout << "\tRun: " << Run << endl;
297  //cout << "\tSubRun: " << SubRun << endl;
298 
299 
300  // -------- find Geant4 TrackId that corresponds to e+/e- from neutrino interaction ----------
301  // MCTruth: Generator
303  e.getByLabel(fMCTruthModuleLabel, MCtruthHandle);
304  std::vector<art::Ptr<simb::MCTruth>> MCtruthlist;
305  art::fill_ptr_vector(MCtruthlist, MCtruthHandle);
306  art::Ptr<simb::MCTruth> MCtruth;
307  // MC (neutrino) interaction
308  int MCinteractions = MCtruthlist.size();
309  //cout << "\nMCinteractions: " << MCinteractions << endl;
310  for (int i=0; i<MCinteractions; i++){
311  MCtruth = MCtruthlist[i];
312  if ( MCtruth->NeutrinoSet() ) { // NeutrinoSet(): whether the neutrino information has been set
313  simb::MCNeutrino nu = MCtruth->GetNeutrino();// GetNeutrino(): reference to neutrino info
314  if( nu.CCNC()==0 ) MC_isCC = 1;
315  else if (nu.CCNC()==1) MC_isCC = 0;
316  simb::MCParticle neutrino = nu.Nu(); // Nu(): the incoming neutrino
317  MC_target = nu.Target(); // Target(): nuclear target, as PDG code
318  MC_incoming_PDG = std::abs(nu.Nu().PdgCode());// here not use std::abs()
319  MC_Q2 = nu.QSqr();// QSqr(): momentum transfer Q^2, in GeV^2
320  MC_channel = nu.InteractionType();// 1001: CCQE
321  MC_W = nu.W(); // W(): hadronic invariant mass
322  const TLorentzVector& nu_momentum = nu.Nu().Momentum(0);
323  nu_momentum.GetXYZT(MC_incoming_P);
324  const TLorentzVector& vertex =neutrino.Position(0);
325  vertex.GetXYZT(MC_vertex);
326  simb::MCParticle lepton = nu.Lepton();// Lepton(): the outgoing lepton
327  MC_lepton_PDG = lepton.PdgCode();
328 
330 
331  //cout << "\tMCinteraction: " << i << "\n\t"
332  // << "neutrino PDG: " << MC_incoming_PDG << "\n\t"
333  // << "MC_lepton_PDG: " << MC_lepton_PDG << "\n\t"
334  // << "MC_channel: " << MC_channel << "\n\t"
335  // << "incoming E: " << MC_incoming_P[3] << "\n\t"
336  // << "MC_vertex: " << MC_vertex[0] << " , " << MC_vertex[1] << " , " <<MC_vertex[2] << " , " <<MC_vertex[3] << endl;
337  }
338  // MCTruth Generator Particles
339  int nParticles = MCtruthlist[0]->NParticles();
340  //cout << "\n\tNparticles: " << MCtruth->NParticles() <<endl;
341  for (int i=0; i<nParticles; i++){
342  simb::MCParticle pi = MCtruthlist[0]->GetParticle(i);
343  //cout << "\tparticle: " << i << "\n\t\t"
344  //<< "Pdgcode: " << pi.PdgCode() << "; Mother: " << pi.Mother() << "; TrackId: " << pi.TrackId() << endl;
345  // Mother(): mother = -1 means that this particle has no mother
346  // TrackId(): same as the index in the MCParticleList
347  }
348  }
349 
350  // Geant4: MCParticle -> lepton (e)
351  // Note: generator level MCPartilceList is different from Geant4 level MCParticleList. MCParticleList(Geant4) contains all particles in MCParticleList(generator) but their the indexes (TrackIds) are different.
352  simb::MCParticle *MClepton = NULL; //Geant4 level
354  const sim::ParticleList& plist = pi_serv->ParticleList();
355  simb::MCParticle *particle=0;
356 
357  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar!=plist.end();++ipar){
358  particle = ipar->second; // first is index(TrackId), second is value (point address)
359 
360  auto & truth = pi_serv->ParticleToMCTruth_P(particle);// beam neutrino only
361  if ( truth->Origin()==simb::kBeamNeutrino && std::abs(particle->PdgCode()) == fLeptonPDGcode && particle->Mother()==0){ // primary lepton; Mother() = 0 means e^{-} for v+n=e^{-}+p
362  const TLorentzVector& lepton_momentum =particle->Momentum(0);
363  const TLorentzVector& lepton_position =particle->Position(0);
364  const TLorentzVector& lepton_positionEnd = particle->EndPosition();
365  const TLorentzVector& lepton_momentumEnd = particle->EndMomentum();
366  lepton_momentum.GetXYZT(MC_lepton_startMomentum);
367  lepton_position.GetXYZT(MC_lepton_startXYZT);
368  lepton_positionEnd.GetXYZT(MC_lepton_endXYZT);
369  lepton_momentumEnd.GetXYZT(MC_lepton_endMomentum);
370  MC_leptonID = particle->TrackId();
371 
373  //cout << "\nGeant Track ID for electron/positron: " << endl;
374  //cout << "\tMClepton PDG: " << particle->PdgCode() << " ; MC_leptonID: " << MC_leptonID << endl;
375  MClepton = particle;
376  //cout << "\tMClepton PDG:" << MClepton->PdgCode() <<endl;
377  //cout << "\tipar->first (TrackId): " << ipar->first << endl;
378  }
379  }
380 
381  // check if the interaction is inside the Fiducial Volume
382  bool isFiducial = false;
383  isFiducial = insideFV( MC_vertex);
384  if (isFiducial) {
385  //cout <<"\nInfo: Interaction is inside the Fiducial Volume.\n" << endl;
386  }
387  else {
388  //cout << "\n********Interaction is NOT inside the Fiducial Volume. RETURN**********" << endl;
389  return;
390  }
391 
393  if (MClepton){
394  h_Ev_den->Fill(MC_incoming_P[3]);
396  }
397  }
398 
399  //if (MC_isCC && (fNeutrinoPDGcode == std::abs(MC_incoming_PDG)) && MC_incoming_E < 0.5) {
400  // cout << "\n------output CC event info for neutrino energy less than 0.5 GeV ------\n" << endl;
401  // cout << "\tEvent : " << Event << "\n"
402  // << "\tRun : " << Run << "\n"
403  // << "\tSubRun : " << SubRun << "\n"
404  // << "\tMC_incoming_E: " << MC_incoming_E << endl;
405  // }
406 
407  // recob::Hit
408  // ---- build a map for total hit charges corresponding to MC_leptonID (Geant4) ----
409 
411  std::vector<art::Ptr<recob::Hit>> all_hits;
412  if(e.getByLabel(fHitModuleLabel,hitHandle)){
413  art::fill_ptr_vector(all_hits, hitHandle);
414  }
415  //cout << "\nTotal hits:" << endl;
416  //cout << "\tall_hits.size(): " << all_hits.size() << endl;
417 
418  double ehit_total =0.0;
419 
420  art::FindManyP<simb::MCParticle,anab::BackTrackerHitMatchingData> fmhitmc(hitHandle, e, fTruthMatchDataModuleLabel);// need #include "lardataobj/AnalysisBase/BackTrackerMatchingData.h"
421 
422  std::map<int,double> all_hits_trk_Q;//Q for charge: Integral()
423  for (size_t i=0; i < all_hits.size(); ++i) {
424  art::Ptr<recob::Hit> hit = all_hits[i];
425  auto particles = fmhitmc.at(hit.key());// particles here is a pointer. A hit may come from multiple particles.
426  auto hitmatch = fmhitmc.data(hit.key());// metadata
427  double maxenergy = -1e9;
428  int hit_TrackId = 0;
429  for (size_t j = 0; j < particles.size(); ++j){
430  if (!particles[j]) continue;
431  if (!pi_serv->TrackIdToMotherParticle_P(particles[j]->TrackId())) continue;
432  size_t trkid = (pi_serv->TrackIdToMotherParticle_P(particles[j]->TrackId()))->TrackId();
433 
434  if ( (hitmatch[j]->energy) > maxenergy ){
435  maxenergy = hitmatch[j]->energy;
436  hit_TrackId = trkid;
437  }
438  }
439  if (hit_TrackId == MC_leptonID) ehit_total += hit->Integral();
440  all_hits_trk_Q[hit_TrackId] += hit->Integral(); // Integral(): the integral under the calibrated signal waveform of the hit, in tick x ADC units
441  }
442  //cout << "....ehit_total: " << ehit_total << endl;
443  //cout << "\tall_hits_trk_Q.size(): " << all_hits_trk_Q.size() << endl;
444 
445 
446  //--------- Loop over all showers: find the associated hits for each shower ------------
447  const simb::MCParticle *MClepton_reco = NULL;
448 
449  double temp_sh_ehit_Q = 0.0;
450  double temp_sh_allHit_Q = 0.0;
451  int temp_sh_TrackId = -999;
453 
455  std::vector<art::Ptr<recob::Shower>> showerlist;
456  if(e.getByLabel(fShowerModuleLabel,showerHandle)){
457  art::fill_ptr_vector(showerlist, showerHandle);
458  }
459  n_recoShowers= showerlist.size();
460  //cout << "\nRecon Showers: " << endl;
461  //cout << "\tn_recoShowers: " << n_recoShowers << endl;
462  art::FindManyP<recob::Hit> sh_hitsAll(showerHandle, e, fShowerModuleLabel);
463 
464  //std::vector<std::map<int,double>> showers_hits_trk_Q;
465  std::map<int,double> showers_trk_Q;
466  for (int i=0; i<n_recoShowers;i++){ // loop over showers
467  //const simb::MCParticle *particle;
468  sh_hasPrimary_e[i] = 0;
469 
470  std::map<int,double> sh_hits_trk_Q;//Q for charge: Integral()
471  sh_hits_trk_Q.clear();
472 
473  art::Ptr<recob::Shower> shower = showerlist[i];
474  sh_direction_X[i] = shower->Direction().X(); // Direction(): direction cosines at start of the shower
475  sh_direction_Y[i] = shower->Direction().Y();
476  sh_direction_Z[i] = shower->Direction().Z();
477  sh_start_X[i] = shower->ShowerStart().X();
478  sh_start_Y[i] = shower->ShowerStart().Y();
479  sh_start_Z[i] = shower->ShowerStart().Z();
480  sh_length[i] = shower->Length(); // shower length in [cm]
481  //cout << "\tInfo of shower " << i << "\n\t\t"
482  //<< "direction (cosines): " << sh_direction_X[i] << ", " << sh_direction_Y[i] << ", " << sh_start_Z[i] << "\n\t\t"
483  //<< "start position: " << sh_start_X[i] << ", " << sh_start_Y[i] << ", " << sh_start_Z[i] << "\n\t\t"
484  //<< "shower length: " << sh_length[i] << endl;
485 
486 
487  std::vector<art::Ptr<recob::Hit>> sh_hits;// associated hits for ith shower
488  //In mcc8, if we get hits associated with the shower through shower->hits association directly for pandora, the hit list is incomplete. The recommended way of getting hits is through association with pfparticle:
489  //shower->pfparticle->clusters->hits
490  //----------getting hits through PFParticle (an option here)-------------------
492  //cout << "\n==== Getting Hits associated with shower THROUGH PFParticle ====\n" << endl;
493  //cout << "\nHits in a shower through PFParticle:\n" << endl;
494 
495  // PFParticle
497  std::vector<art::Ptr<recob::PFParticle> > pfps;
498  if (e.getByLabel(fShowerModuleLabel, pfpHandle)) art::fill_ptr_vector(pfps, pfpHandle);
499  // Clusters
501  std::vector<art::Ptr<recob::Cluster> > clusters;
502  if (e.getByLabel(fShowerModuleLabel, clusterHandle)) art::fill_ptr_vector(clusters, clusterHandle);
503  art::FindManyP<recob::PFParticle> fmps(showerHandle, e, fShowerModuleLabel);// PFParticle in Shower
504  art::FindManyP<recob::Cluster> fmcp(pfpHandle, e, fShowerModuleLabel); // Cluster vs. PFParticle
505  art::FindManyP<recob::Hit> fmhc(clusterHandle, e, fShowerModuleLabel); // Hit in Shower
506  if (fmps.isValid()){
507  std::vector<art::Ptr<recob::PFParticle>> pfs = fmps.at(i);
508  for (size_t ipf = 0; ipf<pfs.size(); ++ipf){
509  if (fmcp.isValid()){
510  std::vector<art::Ptr<recob::Cluster>> clus = fmcp.at(pfs[ipf].key());
511  for (size_t iclu = 0; iclu<clus.size(); ++iclu){
512  if (fmhc.isValid()){
513  std::vector<art::Ptr<recob::Hit>> hits = fmhc.at(clus[iclu].key());
514  for (size_t ihit = 0; ihit<hits.size(); ++ihit){
515  sh_hits.push_back(hits[ihit]);
516  }
517  }
518  }
519  }
520  }
521  }
522 
523  // cout << "\tsh_hits.size(): " << sh_hits.size() << endl;
524 
525  // for (size_t k=0;k<sh_hits.size();k++){
526  // art::Ptr<recob::Hit> hit = sh_hits[k];
527  // cout << k << "\thit.key(): " << hit.key() << endl;
528  // cout << k << "\thit->Integral(): " << hit->Integral() << endl;
529  // }
530  } else {
531 
532  // ----- using shower->hit association for getting hits associated with shower-----
533  sh_hits = sh_hitsAll.at(i);// associated hits for ith shower (using association of hits and shower)
534  //cout << "\t\tsh_hits.size(): " << sh_hits.size() << endl;
535  } // we only choose one method for hits associated with shower here.
536 
537 
538  for (size_t k=0; k < sh_hits.size(); ++k) {
539  art::Ptr<recob::Hit> hit = sh_hits[k];
540  auto particles = fmhitmc.at(hit.key());
541  auto hitmatch = fmhitmc.data(hit.key());
542  double maxenergy = -1e9;
543  int hit_TrackId = 0;
544  for (size_t j = 0; j < particles.size(); ++j){
545  if (!particles[j]) continue;
546  if (!pi_serv->TrackIdToMotherParticle_P(particles[j]->TrackId())) continue;
547  size_t trkid = (pi_serv->TrackIdToMotherParticle_P(particles[j]->TrackId()))->TrackId();
548 
549  if ( (hitmatch[j]->energy) > maxenergy ){
550  maxenergy = hitmatch[j]->energy;
551  hit_TrackId = trkid;
552  }
553  }
554  if (hit_TrackId == MC_leptonID) {
555  sh_ehit_Q[i] += hit->Integral();
556  }
557  sh_allhit_Q[i] += hit->Integral();
558  sh_hits_trk_Q[hit_TrackId] += hit->Integral();// Q from the same hit_TrackId
559  }
560  //cout << "\tsh_hits_trk_Q.size(): " << sh_hits_trk_Q.size() << endl;
561  //showers_hits_trk_Q.push_back(sh_hits_trk_Q);
562 
563  // get TrackId for a shower
564  double maxshowerQ = -1.0e9;
565  //sh_TrackId[i] = 0;
566  for(std::map<int,double>::iterator k = sh_hits_trk_Q.begin(); k != sh_hits_trk_Q.end(); k++) {
567  //cout << k->first << "\t;\t" << k->second << endl;
568  if (k->second > maxshowerQ) {
569  maxshowerQ = k->second;
570  sh_TrackId[i] = k->first;//assign a sh_TrackId with TrackId from hit(particle) that contributing largest to the shower.
571  }
572  }
573 
574  //---------------------------------------------------------------------------------
575  //cout << "\nRecon Shower: " << i << endl;
576  //cout << "\t*****shower primary TrackId: " << sh_TrackId[i] << endl;
577 
578  if (sh_TrackId[i] == MC_leptonID && sh_ehit_Q[i] >0) {
579  temp_sh_TrackId = sh_TrackId[i];
580  sh_hasPrimary_e[i] = 1;
582  MClepton_reco = pi_serv->TrackIdToParticle_P(sh_TrackId[i]);
583 
584  if (sh_allhit_Q[i] >0 && sh_ehit_Q[i] <= sh_allhit_Q[i]){
585  sh_purity[i] = sh_ehit_Q[i]/sh_allhit_Q[i];
586  //cout << "\t*****shower purity: " << sh_purity[i] << endl;
587  } else {
588  sh_purity[i] = 0;
589  }
590  if(ehit_total >0 && sh_ehit_Q[i] <= sh_allhit_Q[i]){
591  sh_completeness[i] = sh_ehit_Q[i] / ehit_total;
592  //cout << "\t*****shower completeness: " << sh_completeness[i] << endl;
593  } else {
594  sh_completeness[i] = 0;
595  }
596  temp_sh_ehit_Q += sh_ehit_Q[i];
597  temp_sh_allHit_Q += sh_allhit_Q[i];
598  }
599 
600  showers_trk_Q[sh_TrackId[i]] += maxshowerQ;
601  //cout << "\tsh_TrackId: " << sh_TrackId[i] <<" ; maxshowerQ: " << maxshowerQ << endl;
602 
603  } // end: for loop over showers
604 
605  // ---------------------------------------------------------------
606  if (MClepton_reco && MClepton) {
607  if ((temp_sh_TrackId == MC_leptonID) && MC_isCC && (fNeutrinoPDGcode == std::abs(MC_incoming_PDG))) {
608  if ((temp_sh_allHit_Q >= temp_sh_ehit_Q) && (temp_sh_ehit_Q > 0.0)) {
609  esh_purity = temp_sh_ehit_Q/temp_sh_allHit_Q;
610  esh_completeness = temp_sh_ehit_Q/ehit_total;
611 
612  if (esh_purity > fMinPurity &&
614  //cout << "\nInfo: fill h_Ev_num ........\n" << endl;
615  h_Ev_num->Fill(MC_incoming_P[3]);
617  }
618  }
619  }
620  }
621  // --------------------------------------------------------------
622 
623  if ( (MClepton_reco && MClepton) && MC_isCC && (fNeutrinoPDGcode == std::abs(MC_incoming_PDG))){
624  //cout << "\n count_primary_e_in_Event: " << count_primary_e_in_Event << endl;
625  for (int j=0; j<count_primary_e_in_Event; j++){
626  esh_each_purity[j] = 0.0;
627  }
628 
629  double temp_esh_1_allhit = -1.0e9;
630  int temp_shower_index = -999;
631  int temp_esh_index = 0;
632  for (int i=0; i<n_recoShowers; i++) {
633  if (sh_TrackId[i] == MC_leptonID) {
634 
635  // for each electron shower
636  if (sh_ehit_Q[i] >0){
637  esh_each_purity[temp_esh_index] = sh_purity[i];
638  esh_each_completeness[temp_esh_index] = sh_completeness[i];
639  temp_esh_index += 1;
640  }
641 
642  // find largest shower
643  if ((sh_allhit_Q[i] > temp_esh_1_allhit) && (sh_ehit_Q[i] > 0.0) ) {
644  temp_esh_1_allhit = sh_allhit_Q[i];
645  temp_shower_index = i;
646  }
647  }
648  }
649  //if (temp_esh_index != count_primary_e_in_Event){
650  // cout << "wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww" << endl;
651  //}
652  // largest shower with electron contribution
653  esh_1_purity = sh_purity[temp_shower_index];
654  esh_1_completeness = sh_completeness[temp_shower_index];
655  }
656 
657  if (count_primary_e_in_Event>0 && MC_isCC && fSaveMCTree) { fEventTree->Fill();}// so far, only save CC event
658 }
659 
660 // ====================================================================================
662  //cout << "\n==== function: doEfficiencies() ====" << endl;
663 
665 
666  if(TEfficiency::CheckConsistency(*h_Ev_num,*h_Ev_den)){
667  h_Eff_Ev = tfs->make<TEfficiency>(*h_Ev_num,*h_Ev_den);
668  TGraphAsymmErrors *grEff_Ev = h_Eff_Ev->CreateGraph();
669  grEff_Ev->Write("grEff_Ev");
670  h_Eff_Ev->Write("h_Eff_Ev");
671  }
672 
673  if(TEfficiency::CheckConsistency(*h_Ee_num,*h_Ee_den)){
674  h_Eff_Ee = tfs->make<TEfficiency>(*h_Ee_num,*h_Ee_den);
675  TGraphAsymmErrors *grEff_Ee = h_Eff_Ee->CreateGraph();
676  grEff_Ee->Write("grEff_Ee");
677  h_Eff_Ee->Write("h_Eff_Ee");
678  }
679 
680 }
681 
682 // ====================================================================================
683 bool NuShowerEff::insideFV( double vertex[4] ){
684  //cout << "\n==== function: insideFV() ====" << endl;
685  double x = vertex[0];
686  double y = vertex[1];
687  double z = vertex[2];
688 
689  if (x>fFidVolXmin && x<fFidVolXmax&&
690  y>fFidVolYmin && y<fFidVolYmax&&
691  z>fFidVolZmin && z<fFidVolZmax)
692  {return true;}
693  else
694  {return false;}
695 }
696 
697 // ====================================================================================
699  //cout << "\n===== function: reset() =====\n" << endl;
700 
701  MC_incoming_PDG = -999;
702  MC_lepton_PDG =-999;
703  MC_isCC =-999;
704  MC_channel =-999;
705  MC_target =-999;
706  MC_Q2 =-999.0;
707  MC_W =-999.0;
708  MC_lepton_theta = -999.0;
709  MC_leptonID = -999;
710  MC_LeptonTrack = -999;
711 
712  for(int i=0; i<MAX_SHOWERS; i++){
713  sh_direction_X[i] = -999.0;
714  sh_direction_Y[i] = -999.0;
715  sh_direction_Z[i] = -999.0;
716  sh_start_X[i] = -999.0;
717  sh_start_Y[i] = -999.0;
718  sh_start_Z[i] = -999.0;
719  sh_length[i] = -999.0;
720  sh_hasPrimary_e[i] = -999;
721  sh_ehit_Q[i] = 0.0;
722  sh_allhit_Q[i] = 0.0;
723  sh_purity[i] = 0.0;
724  sh_completeness[i] = 0.0;
725  }
726 
727 }
Float_t x
Definition: compare.C:6
intermediate_table::iterator iterator
double esh_each_purity[MAX_SHOWERS]
double sh_direction_X[MAX_SHOWERS]
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 MC_lepton_startEnergy
TEfficiency * h_Eff_Ee
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:224
double Length() const
Definition: Shower.h:201
NuShowerEff(fhicl::ParameterSet const &p)
double sh_start_Y[MAX_SHOWERS]
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void endJob() override
const simb::MCParticle * TrackIdToParticle_P(int id) const
int Mother() const
Definition: MCParticle.h:212
Declaration of signal hit object.
Float_t y
Definition: compare.C:6
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
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
double sh_direction_Y[MAX_SHOWERS]
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::string fMCTruthModuleLabel
STL namespace.
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
#define MAX_SHOWERS
data_const_reference data(size_type i) const
Definition: FindManyP.h:448
double MC_lepton_startMomentum[4]
Particle class.
const art::Ptr< simb::MCTruth > & ParticleToMCTruth_P(const simb::MCParticle *p) const
Definition: Run.h:21
int TrackId() const
Definition: MCParticle.h:209
bool fHitShowerThroughPFParticle
int sh_TrackId[MAX_SHOWERS]
std::string fShowerModuleLabel
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
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:56
void beginJob()
Definition: Breakpoints.cc:14
std::string fHitModuleLabel
double esh_each_completeness[MAX_SHOWERS]
TEfficiency * h_Eff_Ev
int sh_hasPrimary_e[MAX_SHOWERS]
void analyze(art::Event const &e) override
bool insideFV(double vertex[4])
double energy
Definition: plottest35.C:25
key_type key() const
Definition: Ptr.h:238
iterator begin()
Definition: ParticleList.h:305
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:89
const TVector3 & Direction() const
Definition: Shower.h:189
double sh_start_Z[MAX_SHOWERS]
RunNumber_t run() const
Definition: DataViewImpl.cc:82
Declaration of cluster object.
double sh_purity[MAX_SHOWERS]
constexpr std::tuple_element_t< I, art::AssnsNode< L, R, D > > & get(art::AssnsNode< L, R, D > &t) noexcept
Detector simulation of raw signals on wires.
const sim::ParticleList & ParticleList() const
double sh_ehit_Q[MAX_SHOWERS]
double MC_lepton_startXYZT[4]
double MC_lepton_endXYZT[4]
double sh_direction_Z[MAX_SHOWERS]
int Target() const
Definition: MCNeutrino.h:151
double sh_completeness[MAX_SHOWERS]
double MC_lepton_endMomentum[40]
double sh_allhit_Q[MAX_SHOWERS]
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
void beginJob() override
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
EventNumber_t event() const
Definition: EventID.h:116
void beginRun(art::Run const &run) override
std::string fTruthMatchDataModuleLabel
bool NeutrinoSet() const
Definition: MCTruth.h:77
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:293
const simb::MCParticle * TrackIdToMotherParticle_P(int id) const
Float_t e
Definition: plot.C:34
Event generator information.
Definition: MCNeutrino.h:18
Namespace collecting geometry-related classes utilities.
double sh_start_X[MAX_SHOWERS]
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
EventID id() const
Definition: Event.cc:37
art framework interface to geometry description
double MC_incoming_P[4]
double sh_length[MAX_SHOWERS]
Beam neutrinos.
Definition: MCTruth.h:23
vertex reconstruction