LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ClusterTrackAna_module.cc
Go to the documentation of this file.
1 // Class: ClusterTrackAna
3 // Plugin Type: analyzer (art v3_03_01)
4 // File: ClusterTrackAna_module.cc
5 //
6 // Calculates the performance of 2D cluster and 3D track reconstruction modules.
7 // The metrics Efficiency, Purity, and Efficiency * Purity (EP) are calculated using
8 // MC truth matched hits in each TPC and Plane. Note that the metrics are calculated in
9 // each TPC and plane for clusters or tracks that span multiple TPCs. The metrics are summed and
10 // reported at the end of the job for each selected true particle type. The metrics are also
11 // summed over ALL selected true particle types. Note that the term "cluster" refers to
12 // a subset of hits in a recob::Track or recob::Cluster that reside in one TPC and one plane
13 // and may be matched to a simb::MCParticle
14 //
15 // Generated at Wed Jan 8 10:33:20 2020 by Bruce Baller using cetskelgen
16 // from cetlib version v3_08_00.
18 
27 #include "fhiclcpp/ParameterSet.h"
29 
31 
40 
41 namespace cluster {
42  class ClusterTrackAna;
43 }
44 
46 public:
47  explicit ClusterTrackAna(fhicl::ParameterSet const& p);
48 
49  ClusterTrackAna(ClusterTrackAna const&) = delete;
50  ClusterTrackAna(ClusterTrackAna&&) = delete;
51  ClusterTrackAna& operator=(ClusterTrackAna const&) = delete;
53 
54  // Required functions.
55  void analyze(art::Event const& e) override;
56 
57 private:
58  void endJob() override;
59 
64  std::vector<int> fSkipPDGCodes;
65  short fPrintLevel;
66  unsigned int fInTPC;
67  float fBadEP;
68 
69  // Count of the number of events considered
70  unsigned int fEventCnt{0};
71  // Count of the number of poorly reconstructed clusters
72  unsigned int fNBadEP{0};
73  // count of EP entries for electrons(0), muons(1), pions(2), kaons(3), protons(4)
74  std::array<float, 5> Cnts{{0}};
75  std::array<float, 5> EPSums{{0}};
76  // same for Efficiency
77  std::array<float, 5> ESums{{0}};
78  // and for Purity
79  std::array<float, 5> PSums{{0}};
80 
82  std::vector<unsigned int> fHitMCPIndex;
83 
85  true};
86  bool fFirstPrint{true};
87  unsigned int fCurrentRun{0};
88 
89  void FirstLastHitInPlane(unsigned int tpc,
90  unsigned int plane,
91  unsigned int mcpi,
92  unsigned int& firstHitIndex,
93  unsigned int& lastHitIndex);
94  std::string HitLocation(
95  unsigned int iht);
96 };
97 
99 {
100  fHitModuleLabel = pset.get<art::InputTag>("HitModuleLabel");
101  bool validModuleLabel = false;
102  fClusterModuleLabel = "NA";
103  fTrackModuleLabel = "NA";
104  fClusterModuleLabel = pset.get<art::InputTag>("ClusterModuleLabel", "NA");
105  if (fClusterModuleLabel != "NA") validModuleLabel = true;
106  fTrackModuleLabel = pset.get<art::InputTag>("TrackModuleLabel", "NA");
107  if (validModuleLabel && fTrackModuleLabel != "NA")
108  throw cet::exception("ClusterTrackAna")
109  << "You must specify either a ClusterModuleLabel OR a TrackModuleLabel\n";
110  if (!validModuleLabel && fTrackModuleLabel != "NA") validModuleLabel = true;
111  // origin = 0 (anything), 1(nu), 2(cosmics), 3(SN nu), 4(SingleParticle)
112  int tmp = pset.get<int>("TruthOrigin", 0);
114  fPrintLevel = pset.get<short>("PrintLevel", 0);
115  if (pset.has_key("SkipPDGCodes")) fSkipPDGCodes = pset.get<std::vector<int>>("SkipPDGCodes");
116  fBadEP = pset.get<float>("BadEP", 0.);
117  fInTPC = UINT_MAX;
118  int intpc = pset.get<int>("InTPC", -1);
119  if (intpc >= 0) fInTPC = intpc;
120  // do some initialization
121  Cnts.fill(0.);
122  EPSums.fill(0.);
123  ESums.fill(0.);
124  PSums.fill(0.);
125 } // ClusterTrackAna constructor
126 
129 {
130  // Match hits to MCParticles, then consider reconstructed hits in each TPC and plane
131  // to calculate Efficiency, Purity and Efficiency * Purity (aka EP).
132 
133  ++fEventCnt;
134  auto const* geom = lar::providerFrom<geo::Geometry>();
135  // auto hitsHandle = art::Handle<std::vector<recob::Hit>>();
137  throw cet::exception("ClusterTrackAna")
138  << "Failed to get a handle to hit collection '" << fHitModuleLabel.label() << "'\n";
139  // get a reference to the MCParticles
140  auto mcpHandle = art::Handle<std::vector<simb::MCParticle>>();
141  if (!evt.getByLabel("largeant", mcpHandle))
142  throw cet::exception("ClusterTrackAna")
143  << "Failed to get a handle to MCParticles using largeant\n";
144 
145  // decide whether to consider cluster -> hit -> MCParticle for any MCParticle origin or for
146  // a specific user-specified origin
147  bool anySource = (fTruthOrigin == simb::kUnknown);
148 
149  if (fPrintLevel > 0 && evt.run() != fCurrentRun) {
150  mf::LogVerbatim("ClusterTrackAna") << "Run: " << evt.run();
151  fCurrentRun = evt.run();
152  }
153 
156  sim::ParticleList const& plist = pi_serv->ParticleList();
157  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
158 
159  // make a list of Hit -> MCParticle assns in all TPCs. The first step is
160  // to make a list of Geant TrackIDs whose origin was requested by the user.
161  std::vector<int> trkIDs;
162  // and a vector of MCParticle indices into the mcpHandle vector indexed by trkIDs
163  std::vector<unsigned int> MCPIs;
164  for (sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
165  const simb::MCParticle* part = (*ipart).second;
166  int trackID = part->TrackId();
167  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(trackID);
168  if (!anySource && theTruth->Origin() != fTruthOrigin) continue;
169  int pdg = abs(part->PdgCode());
170  bool isCharged = (pdg == 11) || (pdg == 13) || (pdg == 211) || (pdg == 321) || (pdg == 2212);
171  if (!isCharged) continue;
172  if (std::find(fSkipPDGCodes.begin(), fSkipPDGCodes.end(), pdg) != fSkipPDGCodes.end()) continue;
173  float TMeV = 1000 * (part->E() - part->Mass());
174  // ignore low energy particles
175  if (TMeV < 1) continue;
176  // find the MCParticle index and stash it in trkMCP
177  unsigned int mcpi = UINT_MAX;
178  for (mcpi = 0; mcpi < (*mcpHandle).size(); ++mcpi)
179  if ((*mcpHandle)[mcpi].TrackId() == trackID) break;
180  if (mcpi == UINT_MAX) {
181  std::cout << "Failed to find a MCParticle from ParticleList\n";
182  return;
183  }
184  trkIDs.push_back(trackID);
185  MCPIs.push_back(mcpi);
186  } // ipart
187  if (trkIDs.empty()) return;
188 
189  // next construct a companion vector of MCParticle indices indexed to the full hit collection
190  fHitMCPIndex.resize((*fHitHandle).size(), UINT_MAX);
191  // Make a list of the <first, last> MC-matched hit in each TPC. This will be used
192  // to only iterate through the range of hits that are interesting
193  std::vector<std::pair<unsigned int, unsigned int>> hitRange(geom->NTPC() + 1);
194  size_t noMatch = 0;
195  size_t nMatch = 0;
196  size_t nInTPC = 0;
197  for (auto& hr : hitRange)
198  hr = std::make_pair(UINT_MAX, UINT_MAX);
199  for (size_t iht = 0; iht < (*fHitHandle).size(); ++iht) {
200  auto& hit = (*fHitHandle)[iht];
201  // only consider hits in a single TPC?
202  unsigned int tpc = hit.WireID().TPC;
203  if (fInTPC != UINT_MAX && tpc != fInTPC) continue;
204  ++nInTPC;
205  auto tides = bt_serv->HitToTrackIDEs(clockData, hit);
206  if (tides.empty()) {
207  ++noMatch;
208  continue;
209  }
210  for (auto& tide : tides) {
211  // declare a match to a MCParticle if > 50% of energy is from it
212  if (tide.energyFrac < 0.5) continue;
213  int trackID = tide.trackID;
214  // find the MCParticle index and define the Hit -> MCParticle assn
215  for (size_t indx = 0; indx < trkIDs.size(); ++indx) {
216  if (trkIDs[indx] == trackID) {
217  fHitMCPIndex[iht] = MCPIs[indx];
218  break;
219  }
220  } // indx
221  break;
222  } // tide
223  if (fHitMCPIndex[iht] == UINT_MAX) continue;
224  ++nMatch;
225  // populate hitRange
226  if (hitRange[tpc].first == UINT_MAX) hitRange[tpc].first = iht;
227  hitRange[tpc].second = iht;
228  } // iht
229 
230  if (nMatch == 0) {
231  fHitMCPIndex.resize(0);
232  return;
233  }
234 
235  if (fPrintLevel > 3) {
236  // Print the gory details
237  mf::LogVerbatim myprt("ClusterTrackAna");
238  myprt << "Checking " << trkIDs.size() << " selected MCParticles with the requested TruthOrigin";
239  if (fInTPC == UINT_MAX) { myprt << " in all TPCs\n"; }
240  else {
241  myprt << " in TPC " << fInTPC;
242  }
243  myprt << " in Run " << evt.run() << " Event " << evt.event();
244  myprt << "\n";
245  myprt << "Found " << nMatch << " MC-matched hits with the requested TruthOrigin";
246  myprt << " out of " << nInTPC << " total hits.\n";
247  myprt << "Found " << noMatch << " hits not matched to ANY MCParticle.\n";
248  // count the number of hits matched to each MCParticle in the list
249  // mcpi count
250  std::vector<std::pair<unsigned int, unsigned int>> mcpCnt;
251  for (auto mcpi : fHitMCPIndex) {
252  if (mcpi == UINT_MAX) continue;
253  unsigned short mIndx = 0;
254  for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
255  if (mcpCnt[mIndx].first == mcpi) break;
256  if (mIndx == mcpCnt.size()) mcpCnt.push_back(std::make_pair(mcpi, 0));
257  // increment the count of MC-matched hits
258  ++mcpCnt[mIndx].second;
259  } // mcpi
260  myprt << " MCPI TrackID PDGCode T(MeV) nHits Process";
261  for (auto mcpcnt : mcpCnt) {
262  myprt << "\n";
263  unsigned int mcpi = mcpcnt.first;
264  auto& mcp = (*mcpHandle)[mcpi];
265  myprt << std::setw(5) << mcpi;
266  myprt << std::setw(8) << mcp.TrackId();
267  myprt << std::setw(8) << abs(mcp.PdgCode());
268  int TMeV = 1000 * (mcp.E() - mcp.Mass());
269  myprt << std::setw(9) << TMeV;
270  myprt << std::setw(8) << mcpcnt.second;
271  myprt << " " << mcp.Process();
272  } // mcpcnt
273  } // fPrintLevel > 3
274 
275  // fill a vector of hits indices from all clusters or all tracks
276  std::vector<std::vector<unsigned int>> recoHits;
277  // and a companion vector of indices into the cluster or track collection
278  std::vector<unsigned int> recoIndex;
279  // handles to these collections to enable printing the cluster or track ID. The
280  // handle will be valid for either clusters or tracks
283  if (fClusterModuleLabel.label() != "NA") {
284  // get MC-matched hits from clusters
285  if (!evt.getByLabel(fClusterModuleLabel, inputClusters))
286  throw cet::exception("ClusterTrackAna")
287  << "Failed to get a handle to the cluster collection '" << fClusterModuleLabel.label()
288  << "'\n";
289  art::FindManyP<recob::Hit> hitsFromCls(inputClusters, evt, fClusterModuleLabel);
290  if (!hitsFromCls.isValid())
291  throw cet::exception("ClusterTrackAna") << "Failed to get a handle to Cluster -> Hit assns\n";
292  for (unsigned int icl = 0; icl < (*inputClusters).size(); ++icl) {
293  std::vector<art::Ptr<recob::Hit>> cluhits = hitsFromCls.at(icl);
294  if (cluhits.empty()) continue;
295  if (fCompareProductIDs) {
296  if (cluhits[0].id() != fHitHandle.id())
297  throw cet::exception("ClusterTrackAna")
298  << "Hits associated with ClusterModuleLabel are in a different collection than "
299  "HitModuleLabel.\n";
300  fCompareProductIDs = false;
301  } // fCompareProductIDs
302  // only load MC-matched hits. Hits that are not MC-matched were either matched to a
303  // MCParticle that was not selected or were mis-reconstructed by the hit finder. Neither
304  // of these conditions warrant penalizing the reconstruction module performance
305  std::vector<unsigned int> hitsIndex;
306  for (auto& cluhit : cluhits) {
307  if (fHitMCPIndex[cluhit.key()] != UINT_MAX) hitsIndex.push_back(cluhit.key());
308  }
309  if (hitsIndex.empty()) continue;
310  recoIndex.push_back(icl);
311  recoHits.push_back(hitsIndex);
312  } // icl
313  }
314  else {
315  // get MC-matched hits from tracks
316  if (!evt.getByLabel(fTrackModuleLabel, inputTracks))
317  throw cet::exception("ClusterTrackAna")
318  << "Failed to get a handle to the track collection '" << fTrackModuleLabel.label() << "'\n";
319  art::FindManyP<recob::Hit> hitsFromTrk(inputTracks, evt, fTrackModuleLabel);
320  if (!hitsFromTrk.isValid())
321  throw cet::exception("ClusterTrackAna") << "Failed to get a handle to Track -> Hit assns\n";
322  for (unsigned int itk = 0; itk < (*inputTracks).size(); ++itk) {
323  std::vector<art::Ptr<recob::Hit>> trkhits = hitsFromTrk.at(itk);
324  if (trkhits.empty()) continue;
325  if (fCompareProductIDs) {
326  if (trkhits[0].id() != fHitHandle.id())
327  throw cet::exception("ClusterTrackAna")
328  << "Hits associated with TrackModuleLabel are in a different collection than "
329  "HitModuleLabel.\n";
330  fCompareProductIDs = false;
331  } // fCompareProductIDs
332  std::vector<unsigned int> hitsIndex;
333  for (auto& trkhit : trkhits) {
334  if (fHitMCPIndex[trkhit.key()] != UINT_MAX) hitsIndex.push_back(trkhit.key());
335  }
336  if (hitsIndex.empty()) continue;
337  recoIndex.push_back(itk);
338  recoHits.push_back(hitsIndex);
339  } // itk
340  } // get hits from tracks
341 
342  for (const auto& tpcid : geom->Iterate<geo::TPCID>()) {
343  unsigned int tpc = tpcid.TPC;
344  if (hitRange[tpc].first == UINT_MAX) continue;
345  // iterate over planes
346  for (unsigned short plane = 0; plane < geom->Nplanes(); ++plane) {
347  unsigned int tpcMatHitCnt = 0;
348  unsigned int tpcTotHitCnt = 0;
349  // create a list of (MCParticle index, matched hit count> pairs
350  // mcParticle plane
351  std::vector<std::pair<unsigned int, float>> mcpCnt;
352  // count MCParticle, Cluster/Track hit counts - size matched to mcpCnt
353  std::vector<std::vector<std::pair<unsigned int, float>>> mcpClsCnt;
354  for (unsigned int iht = hitRange[tpc].first; iht <= hitRange[tpc].second; ++iht) {
355  auto& hit = (*fHitHandle)[iht];
356  if (hit.WireID().TPC != tpc) continue;
357  if (hit.WireID().Plane != plane) continue;
358  ++tpcTotHitCnt;
359  // ignore hits that were either unmatched to the selected set of PDGCodes
360  // or have a not-requested origin
361  if (fHitMCPIndex[iht] == UINT_MAX) continue;
362  ++tpcMatHitCnt;
363  unsigned int mcpi = fHitMCPIndex[iht];
364  unsigned short mIndx = 0;
365  for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
366  if (mcpCnt[mIndx].first == mcpi) break;
367  if (mIndx == mcpCnt.size()) {
368  // found a MCParticle not in the list yet, so add it
369  mcpCnt.push_back(std::make_pair(mcpi, 0));
370  mcpClsCnt.resize(mcpCnt.size());
371  }
372  // increment the number of MC-matched hits
373  ++mcpCnt[mIndx].second;
374  } // iht
375  // ignore TPCs/planes with few MC-matched hits
376  if (fPrintLevel > 2) {
377  mf::LogVerbatim myprt("ClusterTrackAna");
378  myprt << "TPC:" << tpc << " Plane:" << plane << " has " << tpcMatHitCnt << "/"
379  << tpcTotHitCnt;
380  myprt << " (MC-matched hits) / (all hits)";
381  }
382  if (tpcMatHitCnt < 3) continue;
383  // next iterate over all clusters/tracks and count mc-matched hits that are in this TPC/plane
384  for (unsigned int ii = 0; ii < recoHits.size(); ++ii) {
385  float nRecoHitsInPlane = 0;
386  float nRecoMatHitsInPlane = 0;
387  for (auto iht : recoHits[ii]) {
388  auto& hit = (*fHitHandle)[iht];
389  if (hit.WireID().TPC != tpc) continue;
390  if (hit.WireID().Plane != plane) continue;
391  ++nRecoHitsInPlane;
392  if (fHitMCPIndex[iht] == UINT_MAX) continue;
393  ++nRecoMatHitsInPlane;
394  unsigned int mcpi = fHitMCPIndex[iht];
395  // find the MCParticle index in mcpCnt and use it to count the match entry
396  // in mcpClsCnt
397  unsigned short mIndx = 0;
398  for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
399  if (mcpCnt[mIndx].first == mcpi) break;
400  if (mIndx == mcpCnt.size()) {
401  std::cout << "Logic error: fHitMCPIndex = " << fHitMCPIndex[iht]
402  << " is valid but isn't in the list of MCParticles to use. Please send email "
403  "to baller@fnal.gov.\n";
404  continue;
405  }
406  unsigned short cIndx = 0;
407  for (cIndx = 0; cIndx < mcpClsCnt[mIndx].size(); ++cIndx)
408  if (mcpClsCnt[mIndx][cIndx].first == ii) break;
409  if (cIndx == mcpClsCnt[mIndx].size()) mcpClsCnt[mIndx].push_back(std::make_pair(ii, 0));
410  ++mcpClsCnt[mIndx][cIndx].second;
411  } // cluhit
412  if (nRecoMatHitsInPlane == 0) continue;
413  } // ii
414  // find the cluster that has the most hits matched to each MC Particle
415  for (unsigned short mIndx = 0; mIndx < mcpCnt.size(); ++mIndx) {
416  // require at least 3 MC-matched hits
417  if (mcpCnt[mIndx].second < 3) continue;
418  unsigned int mcpi = mcpCnt[mIndx].first;
419  auto& mcp = (*mcpHandle)[mcpi];
420  unsigned int pdgCode = abs(mcp.PdgCode());
421  unsigned short pIndx = USHRT_MAX;
422  if (pdgCode == 11) pIndx = 0;
423  if (pdgCode == 13) pIndx = 1;
424  if (pdgCode == 211) pIndx = 2;
425  if (pdgCode == 321) pIndx = 3;
426  if (pdgCode == 2212) pIndx = 4;
427  if (mcpClsCnt[mIndx].empty()) {
428  // Un-reconstructed MCParticle hits
429  ++Cnts[pIndx];
430  ++fNBadEP;
431  if (fPrintLevel > 0) {
432  mf::LogVerbatim myprt("ClusterTrackAna");
433  myprt << " MCPI " << mcpi;
434  int TMeV = 1000 * (mcp.E() - mcp.Mass());
435  myprt << " " << TMeV << " MeV";
436  std::string partName = std::to_string(pdgCode);
437  if (pdgCode == 11) partName = "El";
438  if (pdgCode == 13) partName = "Mu";
439  if (pdgCode == 211) partName = "Pi";
440  if (pdgCode == 311) partName = "Ka";
441  if (pdgCode == 2212) partName = "Pr";
442  myprt << std::setw(5) << partName;
443  // print out the range of truth-matched hits
444  unsigned int firstHitIndex = UINT_MAX;
445  unsigned int lastHitIndex = UINT_MAX;
446  FirstLastHitInPlane(tpc, plane, mcpi, firstHitIndex, lastHitIndex);
447  myprt << " Failed to reconstruct. Truth-matched hit range from ";
448  myprt << HitLocation(firstHitIndex);
449  myprt << " to ";
450  myprt << HitLocation(lastHitIndex);
451  myprt << " <- EP = 0!";
452  } // fPrintLevel > 0
453  continue;
454  } // (mcpClsCnt[mIndx].empty()
455  std::pair<unsigned int, float> big = std::make_pair(UINT_MAX, 0);
456  for (unsigned short cIndx = 0; cIndx < mcpClsCnt[mIndx].size(); ++cIndx) {
457  auto& mcc = mcpClsCnt[mIndx][cIndx];
458  if (mcc.second > big.second) big = mcc;
459  } // cIndx
460  if (big.first == UINT_MAX) {
461  if (fPrintLevel > 2) {
462  mf::LogVerbatim myprt("ClusterTrackAna");
463  unsigned int mcpi = mcpCnt[mIndx].first;
464  auto& mcp = (*mcpHandle)[mcpi];
465  myprt << "Match failed: mcpi " << mcpi << " pdg " << mcp.PdgCode();
466  }
467  std::cout << "match failed mIndx " << mIndx << "\n";
468  continue;
469  } // big.first == UINT_MAX
470  unsigned int ii = big.first;
471  float eff = big.second / mcpCnt[mIndx].second;
472  float nRecoHitsInPlane = 0;
473  // define some variables to print the range of the cluster/track
474  unsigned int firstRecoHitIndex = UINT_MAX;
475  unsigned int lastRecoHitIndex = UINT_MAX;
476  for (auto iht : recoHits[ii]) {
477  auto& hit = (*fHitHandle)[iht];
478  if (hit.WireID().TPC != tpc) continue;
479  if (hit.WireID().Plane != plane) continue;
480  ++nRecoHitsInPlane;
481  if (firstRecoHitIndex == UINT_MAX) {
482  firstRecoHitIndex = iht;
483  lastRecoHitIndex = iht;
484  }
485  unsigned int wire = (*fHitHandle)[iht].WireID().Wire;
486  if (wire < (*fHitHandle)[firstRecoHitIndex].WireID().Wire) firstRecoHitIndex = iht;
487  if (wire > (*fHitHandle)[lastRecoHitIndex].WireID().Wire) lastRecoHitIndex = iht;
488  } // iht
489  float pur = big.second / nRecoHitsInPlane;
490  ++Cnts[pIndx];
491  float ep = eff * pur;
492  EPSums[pIndx] += ep;
493  ESums[pIndx] += eff;
494  PSums[pIndx] += pur;
495  bool hasBadEP = (ep < fBadEP);
496  if (hasBadEP) ++fNBadEP;
497  bool prt = fPrintLevel > 1 || (fPrintLevel == 1 && hasBadEP);
498  if (prt) {
499  mf::LogVerbatim myprt("ClusterTrackAna");
500  if (fFirstPrint) {
501  myprt << "Hit location format is TPC:Plane:Wire:Tick\n";
502  myprt << " MCPI Ptcl T(MeV) nMCPHits ____From____ _____To_____";
503  if (inputClusters.isValid()) { myprt << " clsID"; }
504  else {
505  myprt << " trkID";
506  }
507  myprt
508  << " ____From____ _____To_____ nRecoHits nMCPRecoHits Eff Pur EP\n";
509  fFirstPrint = false;
510  } // fFirstPrint
511  myprt << std::setw(5) << mcpi;
512  // convert the PDG code into nicer format
513  std::string partName = std::to_string(pdgCode);
514  if (pdgCode == 11) partName = " Electron";
515  if (pdgCode == 13) partName = " Muon";
516  if (pdgCode == 211) partName = " Pion";
517  if (pdgCode == 311) partName = " Kaon";
518  if (pdgCode == 2212) partName = " Proton";
519  myprt << partName;
520  int TMeV = 1000 * (mcp.E() - mcp.Mass());
521  myprt << std::setw(7) << TMeV;
522  // nMCPHits
523  myprt << std::setw(10) << mcpCnt[mIndx].second;
524  unsigned int firstTruHitIndex = UINT_MAX;
525  unsigned int lastTruHitIndex = UINT_MAX;
526  FirstLastHitInPlane(tpc, plane, mcpi, firstTruHitIndex, lastTruHitIndex);
527  myprt << std::setw(14) << HitLocation(firstTruHitIndex);
528  myprt << std::setw(14) << HitLocation(lastTruHitIndex);
529  int id = -1;
530  if (inputClusters.isValid()) {
531  // print cluster info
532  auto& cls = (*inputClusters)[recoIndex[ii]];
533  id = cls.ID();
534  }
535  else if (inputTracks.isValid()) {
536  // print track info
537  auto& trk = (*inputTracks)[recoIndex[ii]];
538  id = trk.ID();
539  }
540  myprt << std::setw(6) << id;
541  myprt << std::setw(14) << HitLocation(firstRecoHitIndex);
542  myprt << std::setw(14) << HitLocation(lastRecoHitIndex);
543  myprt << std::setw(12) << (int)nRecoHitsInPlane;
544  myprt << std::setw(13) << (int)big.second;
545  myprt << std::fixed << std::setprecision(2);
546  myprt << std::setw(8) << eff;
547  myprt << std::setw(8) << pur;
548  myprt << std::setw(8) << ep;
549  if (hasBadEP) myprt << " <- BadEP";
550  myprt << " Evt: " << evt.event();
551  myprt << " Evt Cnt " << fEventCnt;
552  } // prt
553  } // mIndx
554  } // plane
555  } // tpcid
556 
557  fHitMCPIndex.resize(0);
558 
559 } // analyze
560 
562 std::string cluster::ClusterTrackAna::HitLocation(unsigned int iht)
563 {
564  // Put the hit location into a compact human-readable format
565  if (iht >= (*fHitHandle).size()) return "NA";
566  auto& hit = (*fHitHandle)[iht];
567  return std::to_string(hit.WireID().TPC) + ":" + std::to_string(hit.WireID().Plane) + ":" +
568  std::to_string(hit.WireID().Wire) + ":" + std::to_string((int)hit.PeakTime());
569 } // HitLocation
570 
573  unsigned int plane,
574  unsigned int mcpi,
575  unsigned int& firstHitIndex,
576  unsigned int& lastHitIndex)
577 {
578  // Returns the index of the first hit (lowest wire number) and last hit (highest wire number)
579  // matched to the MCParticle indexed by mcpi in the requested tpc, plane
580  firstHitIndex = UINT_MAX;
581  lastHitIndex = UINT_MAX;
582  for (unsigned int iht = 0; iht < (*fHitHandle).size(); ++iht) {
583  if (fHitMCPIndex[iht] != mcpi) continue;
584  auto& hit = (*fHitHandle)[iht];
585  if (hit.WireID().TPC != tpc) continue;
586  if (hit.WireID().Plane != plane) continue;
587  if (firstHitIndex == UINT_MAX) {
588  firstHitIndex = iht;
589  lastHitIndex = iht;
590  }
591  unsigned int wire = (*fHitHandle)[iht].WireID().Wire;
592  if (wire < (*fHitHandle)[firstHitIndex].WireID().Wire) firstHitIndex = iht;
593  if (wire > (*fHitHandle)[lastHitIndex].WireID().Wire) lastHitIndex = iht;
594  } // iht
595 } // FirstLastHitInPlane
596 
599 {
600  // output results
601  mf::LogVerbatim myprt("ClusterTrackAna");
602  myprt << "ClusterTrackAna summary results for " << fEventCnt;
603  if (fClusterModuleLabel.label() != "NA") {
604  myprt << " events using ClusterModuleLabel: " << fClusterModuleLabel.label();
605  }
606  else {
607  myprt << " events using TrackModuleLabel: " << fTrackModuleLabel.label();
608  }
609  myprt << " Origin: " << fTruthOrigin;
610  if (fInTPC != UINT_MAX) { myprt << " in TPC " << fInTPC; }
611  else {
612  myprt << " in all TPCs";
613  }
614  myprt << "\n";
615  float cnts = 0;
616  for (unsigned short pIndx = 0; pIndx < 5; ++pIndx)
617  cnts += Cnts[pIndx];
618  if (cnts == 0) {
619  myprt << "No ClusterTrackAna results";
620  return;
621  }
622  float sumEP = 0;
623  float sumE = 0;
624  float sumP = 0;
625  myprt << "Efficiency (Eff), Purity (Pur) and Eff * Pur (EP) by selected truth particle types\n";
626  std::array<std::string, 5> pName = {{"El", "Mu", "Pi", "K ", "P "}};
627  myprt << "particle Eff Pur EP\n";
628  myprt << "-------------------------------";
629  for (unsigned short pIndx = 0; pIndx < 5; ++pIndx) {
630  if (Cnts[pIndx] == 0) continue;
631  float ave;
632  myprt << "\n " << pName[pIndx] << " ";
633  myprt << std::fixed << std::setprecision(3);
634  ave = ESums[pIndx] / Cnts[pIndx];
635  myprt << std::setw(8) << ave;
636  ave = PSums[pIndx] / Cnts[pIndx];
637  myprt << std::setw(8) << ave;
638  ave = EPSums[pIndx] / Cnts[pIndx];
639  myprt << std::setw(8) << ave;
640  if (pIndx == 0) continue;
641  sumEP += EPSums[pIndx];
642  sumE += ESums[pIndx];
643  sumP += PSums[pIndx];
644  } // pIndx
645  if (cnts == 0) return;
646  myprt << "\n";
647  myprt << "Averages for all selected truth particles\n";
648  myprt << "Ave Eff " << sumE / cnts;
649  myprt << " Ave Pur " << sumP / cnts;
650  myprt << " Ave EP " << sumEP / cnts;
651  myprt << " nBadEP " << fNBadEP;
652  myprt << " (EP < " << std::fixed << std::setprecision(2) << fBadEP << ")";
653  myprt << "\n";
654  myprt << "MCParticle counts in all TPCs and Planes:";
655  for (unsigned short pIndx = 0; pIndx < 5; ++pIndx) {
656  if (Cnts[pIndx] == 0) continue;
657  myprt << " " << pName[pIndx] << " " << (int)Cnts[pIndx];
658  } // pIndx
659 } // endJob
double E(const int i=0) const
Definition: MCParticle.h:234
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
ClusterTrackAna(fhicl::ParameterSet const &p)
int PdgCode() const
Definition: MCParticle.h:213
Utilities related to art service access.
Int_t ipart
Definition: Style.C:10
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
ClusterTrackAna & operator=(ClusterTrackAna const &)=delete
Declaration of signal hit object.
std::array< float, 5 > PSums
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
enum simb::_ev_origin Origin_t
event origin types
constexpr auto abs(T v)
Returns the absolute value of the argument.
Float_t tmp
Definition: plot.C:35
Cluster finding and building.
std::string HitLocation(unsigned int iht)
packs hit WireID and PeakTime into a compact format
std::array< float, 5 > Cnts
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.cc:6
int TrackId() const
Definition: MCParticle.h:211
bool isValid() const noexcept
Definition: Handle.h:203
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
std::string const & label() const noexcept
Definition: InputTag.cc:79
TString part[npart]
Definition: Style.C:32
void analyze(art::Event const &e) override
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
decltype(auto) constexpr to_string(T &&obj)
ADL-aware version of std::to_string.
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int id) const
std::array< float, 5 > EPSums
iterator begin()
Definition: ParticleList.h:305
Provides recob::Track data product.
EventNumber_t event() const
Definition: Event.cc:41
void FirstLastHitInPlane(unsigned int tpc, unsigned int plane, unsigned int mcpi, unsigned int &firstHitIndex, unsigned int &lastHitIndex)
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
Declaration of cluster object.
Detector simulation of raw signals on wires.
const sim::ParticleList & ParticleList() const
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
bool fCompareProductIDs
compare Hit and Cluster-> Hit art product IDs on the first event
art::Handle< std::vector< recob::Hit > > fHitHandle
TCEvent evt
Definition: DataStructs.cxx:8
std::array< float, 5 > ESums
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.cc:29
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
ProductID id() const
Definition: Handle.h:224
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
std::vector< unsigned int > fHitMCPIndex