LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
PFPAna_module.cc
Go to the documentation of this file.
1 //
3 // PFPAna class
4 //
5 // Bruce Baller
6 //
8 
9 #include <TH1F.h>
10 #include <TProfile.h>
11 #include <array>
12 #include <iomanip>
13 #include <string>
14 
15 // Framework includes
21 #include "art_root_io/TFileService.h"
25 #include "fhiclcpp/ParameterSet.h"
27 
28 // LArSoft includes
41 
42 namespace pfpf {
43 
44  class PFPAna : public art::EDAnalyzer {
45  public:
46  explicit PFPAna(fhicl::ParameterSet const& pset);
47 
48  private:
49  void analyze(const art::Event& evt) override;
50  void beginJob() override;
51 
52  TH1F* fNClusters;
54  // Cosmic Rays
55  TH1F* fCREP2;
56  TH1F* fCRE;
57  TH1F* fCRP;
58  // Neutrino interactions
59  TH1F* fNuKE_elec;
60  TH1F* fNuKE_muon;
61  TH1F* fNuKE_pion;
62  TH1F* fNuKE_kaon;
63  TH1F* fNuKE_prot;
64  TH1F* fNuEP2_elec;
65  TH1F* fNuEP2_muon;
66  TH1F* fNuEP2_pion;
67  TH1F* fNuEP2_kaon;
68  TH1F* fNuEP2_prot;
69  TH1F* fNuE_elec;
70  TH1F* fNuE_muon;
71  TH1F* fNuE_pion;
72  TH1F* fNuE_kaon;
73  TH1F* fNuE_prot;
74  TH1F* fNuP_elec;
75  TH1F* fNuP_muon;
76  TH1F* fNuP_pion;
77  TH1F* fNuP_kaon;
78  TH1F* fNuP_prot;
79 
80  TH1F* fNuVtx_dx;
81  TH1F* fNuVtx_dy;
82  TH1F* fNuVtx_dz;
83 
84  TProfile* fNuEP2_KE_elec;
85  TProfile* fNuEP2_KE_muon;
86  TProfile* fNuEP2_KE_pion;
87  TProfile* fNuEP2_KE_kaon;
88  TProfile* fNuEP2_KE_prot;
89 
90  std::string fHitsModuleLabel;
91  std::string fClusterModuleLabel;
92  std::string fTrackModuleLabel;
94  std::string fVertexModuleLabel;
95  std::vector<float> fElecKERange;
96  std::vector<float> fMuonKERange;
97  std::vector<float> fPionKERange;
98  std::vector<float> fKaonKERange;
99  std::vector<float> fProtKERange;
103  short fPrintLevel;
104 
105  }; // class PFPAna
106 
107  //--------------------------------------------------------------------
109  : EDAnalyzer(pset)
110  , fHitsModuleLabel(pset.get<std::string>("HitsModuleLabel"))
111  , fClusterModuleLabel(pset.get<std::string>("ClusterModuleLabel"))
112  , fTrackModuleLabel(pset.get<std::string>("TrackModuleLabel"))
113  , fPFParticleModuleLabel(pset.get<std::string>("PFParticleModuleLabel"))
114  , fVertexModuleLabel(pset.get<std::string>("VertexModuleLabel"))
115  , fElecKERange(pset.get<std::vector<float>>("ElecKERange"))
116  , fMuonKERange(pset.get<std::vector<float>>("MuonKERange"))
117  , fPionKERange(pset.get<std::vector<float>>("PionKERange"))
118  , fKaonKERange(pset.get<std::vector<float>>("KaonKERange"))
119  , fProtKERange(pset.get<std::vector<float>>("ProtKERange"))
120  , fTrackWeightOption(pset.get<short>("TrackWeightOption"))
121  , fMergeDaughters(pset.get<bool>("MergeDaughters"))
122  , fSkipCosmics(pset.get<bool>("SkipCosmics"))
123  , fPrintLevel(pset.get<short>("PrintLevel"))
124  {}
125 
126  //------------------------------------------------------------------
128  {
130 
131  fNClusters = tfs->make<TH1F>("fNoClustersInEvent", "Number of Clusters", 40, 0, 400);
132  fNHitInCluster = tfs->make<TH1F>("fNHitInCluster", "NHitInCluster", 100, 0, 100);
133 
134  if (!fSkipCosmics) {
135  // Cosmic ray histos
136  fCREP2 = tfs->make<TH1F>("CREP2", "CREP2", 50, 0, 1);
137  fCRE = tfs->make<TH1F>("CRE", "CR Efficiency", 50, 0, 1);
138  fCRP = tfs->make<TH1F>("CRP", "CR Purity", 50, 0, 1);
139  }
140  // Neutrino Int histos
141  fNuKE_elec = tfs->make<TH1F>("NuKE_elec", "NuKE electron", 100, 0, 4000);
142  fNuKE_muon = tfs->make<TH1F>("NuKE_muon", "NuKE muon", 100, 0, 4000);
143  fNuKE_pion = tfs->make<TH1F>("NuKE_pion", "NuKE pion", 100, 0, 4000);
144  fNuKE_kaon = tfs->make<TH1F>("NuKE_kaon", "NuKE kaon", 100, 0, 4000);
145  fNuKE_prot = tfs->make<TH1F>("NuKE_prot", "NuKE proton", 100, 0, 4000);
146 
147  fNuEP2_elec = tfs->make<TH1F>("NuEP2_elec", "NuEP2 electron", 50, 0, 1);
148  fNuEP2_muon = tfs->make<TH1F>("NuEP2_muon", "NuEP2 muon", 50, 0, 1);
149  fNuEP2_pion = tfs->make<TH1F>("NuEP2_pion", "NuEP2 pion", 50, 0, 1);
150  fNuEP2_kaon = tfs->make<TH1F>("NuEP2_kaon", "NuEP2 kaon", 50, 0, 1);
151  fNuEP2_prot = tfs->make<TH1F>("NuEP2_prot", "NuEP2 proton", 50, 0, 1);
152 
153  fNuE_elec = tfs->make<TH1F>("NuE_elec", "Nu Efficiency electron", 50, 0, 1);
154  fNuE_muon = tfs->make<TH1F>("NuE_muon", "Nu Efficiency muon", 50, 0, 1);
155  fNuE_pion = tfs->make<TH1F>("NuE_pion", "Nu Efficiency pion", 50, 0, 1);
156  fNuE_kaon = tfs->make<TH1F>("NuE_kaon", "Nu Efficiency kaon", 50, 0, 1);
157  fNuE_prot = tfs->make<TH1F>("NuE_prot", "Nu Efficiency proton", 50, 0, 1);
158 
159  fNuP_elec = tfs->make<TH1F>("NuP_elec", "Nu Purity electron", 50, 0, 1);
160  fNuP_muon = tfs->make<TH1F>("NuP_muon", "Nu Purity muon", 50, 0, 1);
161  fNuP_pion = tfs->make<TH1F>("NuP_pion", "Nu Purity pion", 50, 0, 1);
162  fNuP_kaon = tfs->make<TH1F>("NuP_kaon", "Nu Purity kaon", 50, 0, 1);
163  fNuP_prot = tfs->make<TH1F>("NuP_prot", "Nu Purity proton", 50, 0, 1);
164 
165  // True - Reco vertex difference
166  fNuVtx_dx = tfs->make<TH1F>("Vtx dx", "Vtx dx", 80, -10, 10);
167  fNuVtx_dy = tfs->make<TH1F>("Vtx dy", "Vtx dy", 80, -10, 10);
168  fNuVtx_dz = tfs->make<TH1F>("Vtx dz", "Vtx dz", 80, -10, 10);
169 
170  fNuEP2_KE_elec = tfs->make<TProfile>("NuEP2_KE_elec", "NuEP2 electron vs KE", 20, 0, 2000);
171  fNuEP2_KE_muon = tfs->make<TProfile>("NuEP2_KE_muon", "NuEP2 muon vs KE", 20, 0, 2000);
172  fNuEP2_KE_pion = tfs->make<TProfile>("NuEP2_KE_pion", "NuEP2 pion vs KE", 20, 0, 2000);
173  fNuEP2_KE_kaon = tfs->make<TProfile>("NuEP2_KE_kaon", "NuEP2 kaon vs KE", 20, 0, 2000);
174  fNuEP2_KE_prot = tfs->make<TProfile>("NuEP2_KE_prot", "NuEP2 proton vs KE", 20, 0, 2000);
175  }
176 
178  {
179  // code stolen from TrackAna_module.cc
180  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();
181  if (wireReadoutGeom.Nplanes() > 3) return;
182 
183  // get all hits in the event
185  evt.getByLabel(fHitsModuleLabel, hitListHandle);
186  std::vector<art::Ptr<recob::Hit>> allhits;
187  art::fill_ptr_vector(allhits, hitListHandle);
188  if (empty(allhits)) return;
189 
190  // get clusters and cluster-hit associations
191  art::Handle<std::vector<recob::Cluster>> clusterListHandle;
192  evt.getByLabel(fClusterModuleLabel, clusterListHandle);
193  art::FindManyP<recob::Hit> fmh(clusterListHandle, evt, fClusterModuleLabel);
194  if (clusterListHandle->size() == 0) return;
195 
196  // get 3D vertices
197  art::Handle<std::vector<recob::Vertex>> vertexListHandle;
198  evt.getByLabel(fVertexModuleLabel, vertexListHandle);
199  art::PtrVector<recob::Vertex> recoVtxList;
200  double xyz[3] = {0, 0, 0};
201  for (unsigned int ii = 0; ii < vertexListHandle->size(); ++ii) {
202  art::Ptr<recob::Vertex> vertex(vertexListHandle, ii);
203  recoVtxList.push_back(vertex);
204  vertex->XYZ(xyz);
205  }
206 
207  // get PFParticles
209  evt.getByLabel(fPFParticleModuleLabel, PFPListHandle);
211  for (unsigned int ii = 0; ii < PFPListHandle->size(); ++ii) {
212  art::Ptr<recob::PFParticle> pfp(PFPListHandle, ii);
213  recoPFPList.push_back(pfp);
214  mf::LogVerbatim("PFPAna") << "PFParticle PDG " << pfp->PdgCode();
215  }
216 
217  // list of all true particles
220  sim::ParticleList const& plist = pi_serv->ParticleList();
221  // list of all true particles that will be considered
222  std::vector<const simb::MCParticle*> plist2;
223  // true (reconstructed) hits for each particle in plist2
224  std::vector<std::vector<art::Ptr<recob::Hit>>> hlist2;
225  // index of cluster matched to each particle in plist2 in each plane
226  std::vector<std::vector<short>> truToCl;
227  // number of true hits in each plane and cluster
228  std::vector<std::vector<unsigned short>> nTruHitInCl;
229  //number of reconstructed hits in all clusters
230  std::vector<unsigned short> nRecHitInCl;
231 
232  // calculate average EP2 for every event to facilitate code development
233  // Beam Neutrinos - muons and not-muons
234 
235  float aveNuEP2mu = 0.;
236  float numNuEP2mu = 0.;
237  float aveNuEP2nm = 0.;
238  float numNuEP2nm = 0.;
239  // Cosmic Rays
240  float aveCREP2 = 0.;
241  float numCREP2 = 0.;
242 
243  // track ID of the neutrino
244  int neutTrackID = -1;
245  std::vector<int> tidlist;
246  float neutEnergy = -1.;
247  int neutIntType = -1;
248  int neutCCNC = -1;
249 
250  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
251 
252  for (sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
253  const simb::MCParticle* part = (*ipart).second;
254  assert(part != 0);
255  int pdg = abs(part->PdgCode());
256  int trackID = part->TrackId();
257  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(trackID);
258  if (fSkipCosmics && theTruth->Origin() == simb::kCosmicRay) continue;
259 
260  if (fPrintLevel > 3)
261  mf::LogVerbatim("PFPAna") << "Pre-Cuts origin " << theTruth->Origin() << " trackID "
262  << trackID << " PDG " << part->PdgCode() << " E " << part->E()
263  << " mass " << part->Mass() << " Mother "
264  << part->Mother() + neutTrackID << " Proc " << part->Process();
265 
266  // Get the neutrino track ID. Assume that there is only one neutrino
267  // interaction and it is first in the list of BeamNeutrino particles
268  if (theTruth->Origin() == simb::kBeamNeutrino && neutTrackID < 0) {
269  neutTrackID = trackID;
270  simb::MCNeutrino theNeutrino = theTruth->GetNeutrino();
271  neutEnergy = 1000. * theNeutrino.Nu().E();
272  neutIntType = theNeutrino.InteractionType();
273  neutCCNC = theNeutrino.CCNC();
274  for (unsigned short iv = 0; iv < recoVtxList.size(); ++iv) {
275  recoVtxList[iv]->XYZ(xyz);
276  fNuVtx_dx->Fill(part->Vx() - xyz[0]);
277  fNuVtx_dy->Fill(part->Vy() - xyz[1]);
278  fNuVtx_dz->Fill(part->Vz() - xyz[2]);
279  } // iv
280  } // theTruth->Origin() == simb::kBeamNeutrino && neutTrackID <
281 
282  bool isCharged = (pdg == 11) || (pdg == 13) || (pdg == 211) || (pdg == 321) || (pdg == 2212);
283 
284  if (!isCharged) continue;
285 
286  float KE = 1000 * (part->E() - part->Mass());
287  // KE (MeV) cuts
288  if (pdg == 11) {
289  if (fElecKERange[0] < 0) continue;
290  // only allow primary electrons
291  if (part->Process() != "primary") continue;
292  if (KE < fElecKERange[0] || KE > fElecKERange[1]) continue;
293  }
294  if (pdg == 13) {
295  if (fMuonKERange[0] < 0) continue;
296  if (KE < fMuonKERange[0] || KE > fMuonKERange[1]) continue;
297  }
298  if (pdg == 211) {
299  if (fPionKERange[0] < 0) continue;
300  if (KE < fPionKERange[0] || KE > fPionKERange[1]) continue;
301  }
302  if (pdg == 321) {
303  if (fKaonKERange[0] < 0) continue;
304  if (KE < fKaonKERange[0] || KE > fKaonKERange[1]) continue;
305  }
306  if (pdg == 2212) {
307  if (fProtKERange[0] < 0) continue;
308  if (KE < fProtKERange[0] || KE > fProtKERange[1]) continue;
309  }
310  // ignore secondaries from neutron interactions
311  if (part->Process() == "NeutronInelastic") continue;
312  plist2.push_back(part);
313  tidlist.push_back(trackID);
314  // initialize the true->(cluster,plane) association
315  std::vector<short> temp{-1, -1, -1};
316  truToCl.push_back(temp);
317  // initialize the true hit count
318  std::vector<unsigned short> temp2(3);
319  nTruHitInCl.push_back(temp2);
320 
321  if (fPrintLevel > 2)
322  mf::LogVerbatim("PFPAna") << plist2.size() - 1 << " Origin " << theTruth->Origin()
323  << " trackID " << trackID << " PDG " << part->PdgCode() << " KE "
324  << (int)KE << " Mother " << part->Mother() + neutTrackID
325  << " Proc " << part->Process();
326  }
327 
328  if (empty(plist2)) return;
329 
330  // get the hits (in all planes) that are matched to the true tracks
331  hlist2 = bt_serv->TrackIdsToHits_Ps(clockData, tidlist, allhits);
332  if (hlist2.size() != plist2.size()) {
333  mf::LogError("PFPAna") << "MC particle list size " << plist2.size()
334  << " != size of MC particle true hits lists " << hlist2.size();
335  return;
336  }
337  tidlist.clear();
338 
339  // vector of (mother, daughter) pairs
340  std::vector<std::pair<unsigned short, unsigned short>> moda;
341  // Deal with mother-daughter tracks
342  if (fMergeDaughters && neutTrackID >= 0) {
343  // Assume that daughters appear later in the list. Step backwards
344  // to accumulate all generations of daughters
345  for (unsigned short dpl = plist2.size() - 1; dpl > 0; --dpl) {
346  // no mother
347  if (plist2[dpl]->Mother() == 0) continue;
348  // electron
349  if (abs(plist2[dpl]->PdgCode()) == 11) continue;
350  // the actual mother trackID is offset from the neutrino trackID
351  int motherID = neutTrackID + plist2[dpl]->Mother() - 1;
352  // ensure that we are only looking at BeamNeutrino daughters
353  if (motherID < 0) continue;
354  // count the number of daughters
355  int ndtr = 0;
356  for (unsigned short kpl = 0; kpl < plist2.size(); ++kpl) {
357  if (plist2[kpl]->Mother() == motherID) ++ndtr;
358  }
359  // require only one daughter
360  if (ndtr > 1) continue;
361  // find the mother in the list
362  int mpl = -1;
363  for (unsigned short jpl = dpl - 1; jpl > 0; --jpl) {
364  if (plist2[jpl]->TrackId() == motherID) {
365  mpl = jpl;
366  break;
367  }
368  } // jpl
369  // mother not found for some reason
370  if (mpl < 0) continue;
371  // ensure that PDG code for mother and daughter are the same
372  if (plist2[dpl]->PdgCode() != plist2[mpl]->PdgCode()) continue;
373  moda.push_back(std::make_pair(mpl, dpl));
374  } // dpl
375  } // MergeDaughters
376 
377  // Now match reconstructed clusters to true particles.
379  for (unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
380  art::Ptr<recob::Cluster> clusterHolder(clusterListHandle, ii);
381  clusters.push_back(clusterHolder);
382  }
383 
384  fNClusters->Fill(clusterListHandle->size());
385  nRecHitInCl.resize(clusters.size());
386 
387  // get the plane from the view. Perhaps there is a method that does
388  // this somewhere...
389  std::map<geo::View_t, unsigned int> ViewToPlane;
390  for (auto const& plane : wireReadoutGeom.Iterate<geo::PlaneGeo>(geo::TPCID{0, 0})) {
391  geo::View_t view = plane.View();
392  ViewToPlane[view] = plane.ID().Plane;
393  }
394  for (size_t icl = 0; icl < clusters.size(); ++icl) {
395  unsigned int plane = ViewToPlane[clusters[icl]->View()];
396  std::vector<art::Ptr<recob::Hit>> cluhits = fmh.at(icl);
397  fNHitInCluster->Fill(cluhits.size());
398  nRecHitInCl[icl] = cluhits.size();
399  // count the number of hits matched to each true particle in plist2
400  std::vector<unsigned short> nHitInPl2(plist2.size());
401  for (size_t iht = 0; iht < cluhits.size(); ++iht) {
402  // look for this hit in all of the truth hit lists
403  short hitInPl2 = -1;
404  for (unsigned short ipl = 0; ipl < plist2.size(); ++ipl) {
405  unsigned short imat = 0;
406  for (imat = 0; imat < hlist2[ipl].size(); ++imat) {
407  if (cluhits[iht] == hlist2[ipl][imat]) break;
408  } // imat
409  if (imat < hlist2[ipl].size()) {
410  hitInPl2 = ipl;
411  break;
412  }
413  } // ipl
414  if (hitInPl2 < 0) continue;
415  // Assign the hit count to the mother if this is a daughter.
416  // Mother-daughter pairs are entered in the moda vector in reverse
417  // order, so assign daughter hits to the highest generation mother.
418  for (unsigned short imd = 0; imd < moda.size(); ++imd) {
419  if (moda[imd].second == hitInPl2) hitInPl2 = moda[imd].first;
420  }
421  // count
422  ++nHitInPl2[hitInPl2];
423  } // iht
424  // Associate the cluster with the truth particle that has the highest
425  // number of cluster hits
426  unsigned short nhit = 0;
427  short imtru = -1;
428  for (unsigned int ipl = 0; ipl < nHitInPl2.size(); ++ipl) {
429  if (nHitInPl2[ipl] > nhit) {
430  nhit = nHitInPl2[ipl];
431  imtru = ipl;
432  }
433  } // ipl
434  // make the cluster->(true,plane) association and save the
435  // number of true hits in the cluster
436  if (imtru != -1) {
437  // clobber a previously made association?
438  if (nhit > nTruHitInCl[imtru][plane]) {
439  truToCl[imtru][plane] = icl;
440  nTruHitInCl[imtru][plane] = nhit;
441  }
442  } // imtru != 1
443  } // icl
444 
445  // ready to calculate Efficiency, Purity in each plane and EP2
446  for (unsigned short ipl = 0; ipl < plist2.size(); ++ipl) {
447  // ignore daughters
448  bool skipit = false;
449  for (unsigned short ii = 0; ii < moda.size(); ++ii) {
450  if (moda[ii].second == ipl) {
451  skipit = true;
452  break;
453  }
454  } // ii
455  if (skipit) continue;
456  // ignore true particles with few true hits. Outside the detector
457  // or not reconstructable?
458  if (hlist2[ipl].size() < 3) continue;
459 
460  int trackID = plist2[ipl]->TrackId();
461  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(trackID);
462  bool isCosmic = (theTruth->Origin() == simb::kCosmicRay);
463  float KE = 1000 * (plist2[ipl]->E() - plist2[ipl]->Mass());
464  int PDG = abs(plist2[ipl]->PdgCode());
465 
466  std::vector<short> nTru(wireReadoutGeom.Nplanes());
467  std::vector<short> nRec(wireReadoutGeom.Nplanes());
468  std::vector<short> nTruRec(wireReadoutGeom.Nplanes());
469  std::vector<float> eff(wireReadoutGeom.Nplanes());
470  std::vector<float> pur(wireReadoutGeom.Nplanes());
471  std::vector<float> ep(wireReadoutGeom.Nplanes());
472  for (unsigned int plane = 0; plane < wireReadoutGeom.Nplanes(); ++plane) {
473  // count the number of true hits in this plane for the true particle.
474  // First count the mother hits
475  for (unsigned short ii = 0; ii < hlist2[ipl].size(); ++ii) {
476  if (ViewToPlane[hlist2[ipl][ii]->View()] == plane) ++nTru[plane];
477  } // ii
478  // next look for daughters and count those hits in all generations
479  unsigned short mom = ipl;
480  std::vector<std::pair<unsigned short, unsigned short>>::reverse_iterator rit =
481  moda.rbegin();
482  while (rit != moda.rend()) {
483  if ((*rit).first == mom) {
484  unsigned short dau = (*rit).second;
485  for (unsigned short jj = 0; jj < hlist2[dau].size(); ++jj) {
486  if (ViewToPlane[hlist2[dau][jj]->View()] == plane) ++nTru[plane];
487  } // jj
488  // It is likely that one hit appears in the mother list
489  // as well as the daughter list, so subtract one from the count
490  --nTru[plane];
491  mom = (*rit).second;
492  } // (*rit).first == mom
493  ++rit;
494  } // rit
495  if (nTru[plane] == 0) { continue; }
496  short icl = truToCl[ipl][plane];
497  nRec[plane] = nRecHitInCl[icl];
498  nTruRec[plane] = nTruHitInCl[ipl][plane];
499  if (nTru[plane] > 0) eff[plane] = (float)nTruRec[plane] / (float)nTru[plane];
500  if (nRec[plane] > 0) pur[plane] = (float)nTruRec[plane] / (float)nRec[plane];
501  ep[plane] = eff[plane] * pur[plane];
502  } // plane
503  // sort the ep values in ascending order
504  std::vector<float> temp;
505  temp = ep;
506  std::sort(temp.begin(), temp.end());
507  // EP2 is the second highest value
508  unsigned short ii = temp.size() - 2;
509  float ep2 = temp[ii];
510  // find the plane that defined EP2
511  short ep2Plane = 0;
512  short ep2Cluster = 0;
513  for (unsigned short jj = 0; jj < temp.size(); ++jj) {
514  if (ep[jj] == ep2) {
515  ep2Plane = jj;
516  ep2Cluster = truToCl[ipl][ep2Plane];
517  break;
518  }
519  } // jj
520  // find the US and DS ends of the cluster for printing
521  std::array<double, 2> clBeg, clEnd;
522  if (ep2Cluster >= 0) {
523  clBeg[0] = clusters[ep2Cluster]->StartWire();
524  clBeg[1] = clusters[ep2Cluster]->StartTick();
525  clEnd[0] = clusters[ep2Cluster]->EndWire();
526  clEnd[1] = clusters[ep2Cluster]->EndTick();
527  }
528  else {
529  clBeg.fill(0.);
530  clEnd.fill(0.);
531  }
532  // fill histograms
533  if (isCosmic) {
534  fCREP2->Fill(ep2);
535  fCRE->Fill(eff[ep2Plane]);
536  fCRP->Fill(pur[ep2Plane]);
537  aveCREP2 += ep2;
538  numCREP2 += 1.;
539  if (fPrintLevel > 1)
540  mf::LogVerbatim("PFPAna")
541  << ">>>CREP2 " << std::fixed << std::setprecision(2) << ep2 << " E " << eff[ep2Plane]
542  << std::setprecision(2) << " P " << pur[ep2Plane] << " P:W:T " << ep2Plane << ":"
543  << (int)clBeg[0] << ":" << (int)clBeg[1] << "-" << ep2Plane << ":" << (int)clEnd[0]
544  << ":" << (int)clEnd[1] << " PDG " << PDG << " KE " << (int)KE << " MeV";
545  } // isCosmic
546  else {
547  float wght = 1.;
548  if (fTrackWeightOption == 1) wght = KE;
549  // accumulate statistics for muons and not-muons
550  if (PDG == 13) {
551  aveNuEP2mu += ep2 * wght;
552  numNuEP2mu += wght;
553  }
554  else {
555  aveNuEP2nm += ep2 * wght;
556  numNuEP2nm += wght;
557  }
558  if (PDG == 11) {
559  fNuKE_elec->Fill(KE, wght);
560  fNuE_elec->Fill(eff[ep2Plane], wght);
561  fNuP_elec->Fill(pur[ep2Plane], wght);
562  fNuEP2_elec->Fill(ep2, wght);
563  fNuEP2_KE_elec->Fill(KE, ep2, wght);
564  }
565  else if (PDG == 13) {
566  fNuKE_muon->Fill(KE, wght);
567  fNuE_muon->Fill(eff[ep2Plane], wght);
568  fNuP_muon->Fill(pur[ep2Plane], wght);
569  fNuEP2_muon->Fill(ep2, wght);
570  fNuEP2_KE_muon->Fill(KE, ep2, wght);
571  }
572  else if (PDG == 211) {
573  fNuKE_pion->Fill(KE, wght);
574  fNuE_pion->Fill(eff[ep2Plane], wght);
575  fNuP_pion->Fill(pur[ep2Plane], wght);
576  fNuEP2_pion->Fill(ep2, wght);
577  fNuEP2_KE_pion->Fill(KE, ep2, wght);
578  }
579  else if (PDG == 321) {
580  fNuKE_kaon->Fill(KE, wght);
581  fNuE_kaon->Fill(eff[ep2Plane], wght);
582  fNuP_kaon->Fill(pur[ep2Plane], wght);
583  fNuEP2_kaon->Fill(ep2, wght);
584  fNuEP2_KE_kaon->Fill(KE, ep2, wght);
585  }
586  else if (PDG == 2212) {
587  fNuKE_prot->Fill(KE, wght);
588  fNuE_prot->Fill(eff[ep2Plane], wght);
589  fNuP_prot->Fill(pur[ep2Plane], wght);
590  fNuEP2_prot->Fill(ep2, wght);
591  fNuEP2_KE_prot->Fill(KE, ep2, wght);
592  }
593  if (fPrintLevel > 1)
594  mf::LogVerbatim("PFPAna")
595  << ">>>NuEP2 " << std::fixed << std::setprecision(2) << ep2 << " E " << eff[ep2Plane]
596  << std::setprecision(2) << " P " << pur[ep2Plane] << " P:W:T " << ep2Plane << ":"
597  << (int)clBeg[0] << ":" << (int)clBeg[1] << "-" << ep2Plane << ":" << (int)clEnd[0]
598  << ":" << (int)clEnd[1] << " PDG " << PDG << " KE " << (int)KE << " MeV ";
599  if (fPrintLevel > 2) {
600  // print out the begin/end true hits
601  mf::LogVerbatim mfp("PFPAna");
602  mfp << " Truth P:W:T ";
603  for (unsigned int plane = 0; plane < wireReadoutGeom.Nplanes(); ++plane) {
604  unsigned short loW = 9999;
605  unsigned short loT = 0;
606  unsigned short hiW = 0;
607  unsigned short hiT = 0;
608  for (unsigned short ii = 0; ii < hlist2[ipl].size(); ++ii) {
609  if (ViewToPlane[hlist2[ipl][ii]->View()] == plane) {
610  art::Ptr<recob::Hit> theHit = hlist2[ipl][ii];
611  if (theHit->WireID().Wire < loW) {
612  loW = theHit->WireID().Wire;
613  loT = theHit->PeakTime();
614  }
615  if (theHit->WireID().Wire > hiW) {
616  hiW = theHit->WireID().Wire;
617  hiT = theHit->PeakTime();
618  }
619  } // correct view
620  } // ii
621  mfp << plane << ":" << loW << ":" << loT << "-" << plane << ":" << hiW << ":" << hiT
622  << " ";
623  } // plane
624  } // fPrintLevel > 2
625  } // !isCosmic
626  } // ipl
627 
628  float ave1 = -1.;
629  if (numNuEP2mu > 0.) ave1 = aveNuEP2mu / numNuEP2mu;
630 
631  float ave2 = -1.;
632  if (numNuEP2nm > 0.) ave2 = aveNuEP2nm / numNuEP2nm;
633 
634  float ave3 = -1.;
635  if (numCREP2 > 0.) ave3 = aveCREP2 / numCREP2;
636 
637  if (fPrintLevel > 0) {
638  std::string nuType = "Other";
639  if (neutCCNC == simb::kCC) {
640  if (neutIntType == 1001) nuType = "CCQE";
641  if (neutIntType == 1091) nuType = "DIS";
642  if (neutIntType == 1097) nuType = "COH";
643  if (neutIntType > 1002 && neutIntType < 1091) nuType = "RES";
644  }
645  else if (neutCCNC == simb::kNC) {
646  nuType = "NC";
647  }
648  else {
649  nuType = "Unknown";
650  }
651  mf::LogVerbatim("PFPAna") << "EvtEP2 " << evt.id().event() << " NuType " << nuType << " Enu "
652  << std::fixed << std::setprecision(0) << neutEnergy << std::right
653  << std::fixed << std::setprecision(2) << " NuMuons " << ave1
654  << " NuPiKp " << ave2 << " CosmicRays " << ave3 << " CCNC "
655  << neutCCNC << " IntType " << neutIntType;
656  }
657  } // analyze
658 
660 
661 } // end namespace
double E(const int i=0) const
Definition: MCParticle.h:234
void XYZ(double *xyz) const
Legacy method to access vertex position, preserved to avoid breaking code. Please try to use Vertex::...
Definition: Vertex.cxx:34
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::vector< float > fElecKERange
int PdgCode() const
Definition: MCParticle.h:213
int CCNC() const
Definition: MCNeutrino.h:148
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
Float_t mfp
Definition: plot.C:20
TH1F * fNuEP2_prot
Int_t ipart
Definition: Style.C:10
std::string fVertexModuleLabel
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
int Mother() const
Definition: MCParticle.h:214
Declaration of signal hit object.
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
list_type::const_iterator const_iterator
Definition: ParticleList.h:132
simb::Origin_t Origin() const
Definition: MCTruth.h:74
double Mass() const
Definition: MCParticle.h:240
TH1F * fNuP_kaon
std::string fHitsModuleLabel
TH1F * fNuKE_prot
std::vector< float > fProtKERange
constexpr auto abs(T v)
Returns the absolute value of the argument.
STL namespace.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
int PdgCode() const
Return the type of particle as a PDG ID.
Definition: PFParticle.h:79
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:430
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:290
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.cc:6
std::string Process() const
Definition: MCParticle.h:216
std::string fClusterModuleLabel
int TrackId() const
Definition: MCParticle.h:211
std::vector< float > fPionKERange
TProfile * fNuEP2_KE_kaon
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
TH1F * fNuVtx_dy
TH1F * fNuKE_elec
int InteractionType() const
Definition: MCNeutrino.h:150
TString part[npart]
Definition: Style.C:32
TProfile * fNuEP2_KE_elec
TH1F * fNuVtx_dz
TH1F * fNuEP2_muon
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
TH1F * fNuE_kaon
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
void beginJob() override
TH1F * fNuP_muon
short fTrackWeightOption
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int id) const
TProfile * fNuEP2_KE_pion
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:67
TH1F * fNuEP2_kaon
iterator begin()
Definition: ParticleList.h:305
TH1F * fNuP_elec
std::string fPFParticleModuleLabel
TH1F * fNHitInCluster
std::vector< std::vector< art::Ptr< recob::Hit > > > TrackIdsToHits_Ps(detinfo::DetectorClocksData const &clockData, std::vector< int > const &tkIds, std::vector< art::Ptr< recob::Hit >> const &hitsIn) const
TH1F * fNuKE_muon
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
Declaration of cluster object.
size_type size() const
Definition: PtrVector.h:302
TH1F * fNuP_prot
const sim::ParticleList & ParticleList() const
TH1F * fNuVtx_dx
TH1F * fNuE_prot
TProfile * fNuEP2_KE_prot
double Vx(const int i=0) const
Definition: MCParticle.h:222
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:226
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
TH1F * fNuE_muon
TH1F * fNClusters
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Encapsulate the construction of a single detector plane .
TProfile * fNuEP2_KE_muon
TH1F * fNuE_pion
TH1F * fNuE_elec
PFPAna(fhicl::ParameterSet const &pset)
TH1F * fNuEP2_elec
TH1F * fNuKE_pion
double Vz(const int i=0) const
Definition: MCParticle.h:224
TH1F * fNuP_pion
EventNumber_t event() const
Definition: EventID.h:116
std::vector< float > fMuonKERange
std::string fTrackModuleLabel
TH1F * fNuKE_kaon
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
Particle list in DetSim contains Monte Carlo particle information.
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
Event generator information.
Definition: MCNeutrino.h:18
Definition: fwd.h:26
void analyze(const art::Event &evt) override
EventID id() const
Definition: Event.cc:23
double Vy(const int i=0) const
Definition: MCParticle.h:223
art framework interface to geometry description
std::vector< float > fKaonKERange
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
Cosmic rays.
Definition: MCTruth.h:24
TH1F * fNuEP2_pion
Beam neutrinos.
Definition: MCTruth.h:23
vertex reconstruction