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