All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
cluster::ClusterTrackAna Class Reference
Inheritance diagram for cluster::ClusterTrackAna:

Public Member Functions

 ClusterTrackAna (fhicl::ParameterSet const &p)
 
 ClusterTrackAna (ClusterTrackAna const &)=delete
 
 ClusterTrackAna (ClusterTrackAna &&)=delete
 
ClusterTrackAna & operator= (ClusterTrackAna const &)=delete
 
ClusterTrackAna & operator= (ClusterTrackAna &&)=delete
 
void analyze (art::Event const &e) override
 

Private Member Functions

void endJob () override
 
void FirstLastHitInPlane (unsigned int tpc, unsigned int plane, unsigned int mcpi, unsigned int &firstHitIndex, unsigned int &lastHitIndex)
 
std::string HitLocation (unsigned int iht)
 packs hit WireID and PeakTime into a compact format More...
 

Private Attributes

art::InputTag fHitModuleLabel
 
art::InputTag fClusterModuleLabel
 
art::InputTag fTrackModuleLabel
 
simb::Origin_t fTruthOrigin
 
std::vector< int > fSkipPDGCodes
 
short fPrintLevel
 
unsigned int fInTPC
 
float fBadEP
 
unsigned int fEventCnt {0}
 
unsigned int fNBadEP {0}
 
std::array< float, 5 > Cnts {{0}}
 
std::array< float, 5 > EPSums {{0}}
 
std::array< float, 5 > ESums {{0}}
 
std::array< float, 5 > PSums {{0}}
 
art::Handle< std::vector
< recob::Hit > > 		fHitHandle
 
std::vector< unsigned int > fHitMCPIndex
 
bool fCompareProductIDs
 compare Hit and Cluster-> Hit art product IDs on the first event More...
 
bool fFirstPrint {true}
 
unsigned int fCurrentRun {0}
 

Detailed Description

Definition at line 45 of file ClusterTrackAna_module.cc.

Constructor & Destructor Documentation

cluster::ClusterTrackAna::ClusterTrackAna ( fhicl::ParameterSet const &  p)
explicit

Definition at line 98 of file ClusterTrackAna_module.cc.

98  : EDAnalyzer{pset}
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);
113  fTruthOrigin = (simb::Origin_t)tmp;
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
cluster::ClusterTrackAna::PSums
std::array< float, 5 > PSums
Definition: ClusterTrackAna_module.cc:79
cluster::ClusterTrackAna::Cnts
std::array< float, 5 > Cnts
Definition: ClusterTrackAna_module.cc:74
cluster::ClusterTrackAna::EPSums
std::array< float, 5 > EPSums
Definition: ClusterTrackAna_module.cc:75
cluster::ClusterTrackAna::ESums
std::array< float, 5 > ESums
Definition: ClusterTrackAna_module.cc:77
cluster::ClusterTrackAna::ClusterTrackAna ( ClusterTrackAna const &  )
delete
cluster::ClusterTrackAna::ClusterTrackAna ( ClusterTrackAna &&  )
delete

Member Function Documentation

void cluster::ClusterTrackAna::analyze ( art::Event const &  e)
override

Definition at line 129 of file ClusterTrackAna_module.cc.

130 {
131  // Match hits to MCParticles, then consider reconstructed hits in each TPC and plane
132  // to calculate Efficiency, Purity and Efficiency * Purity (aka EP).
133 
134  ++fEventCnt;
135  auto const* geom = lar::providerFrom<geo::Geometry>();
136  // auto hitsHandle = art::Handle<std::vector<recob::Hit>>();
137  if (!evt.getByLabel(fHitModuleLabel, fHitHandle))
138  throw cet::exception("ClusterTrackAna")
139  << "Failed to get a handle to hit collection '" << fHitModuleLabel.label() << "'\n";
140  // get a reference to the MCParticles
141  auto mcpHandle = art::Handle<std::vector<simb::MCParticle>>();
142  if (!evt.getByLabel("largeant", mcpHandle))
143  throw cet::exception("ClusterTrackAna")
144  << "Failed to get a handle to MCParticles using largeant\n";
145 
146  // decide whether to consider cluster -> hit -> MCParticle for any MCParticle origin or for
147  // a specific user-specified origin
148  bool anySource = (fTruthOrigin == simb::kUnknown);
149 
150  if (fPrintLevel > 0 && evt.run() != fCurrentRun) {
151  mf::LogVerbatim("ClusterTrackAna") << "Run: " << evt.run();
152  fCurrentRun = evt.run();
153  }
154 
155  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
156  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
157  sim::ParticleList const& plist = pi_serv->ParticleList();
158  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
159 
160  // make a list of Hit -> MCParticle assns in all TPCs. The first step is
161  // to make a list of Geant TrackIDs whose origin was requested by the user.
162  std::vector<int> trkIDs;
163  // and a vector of MCParticle indices into the mcpHandle vector indexed by trkIDs
164  std::vector<unsigned int> MCPIs;
165  for (sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
166  const simb::MCParticle* part = (*ipart).second;
167  int trackID = part->TrackId();
168  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(trackID);
169  if (!anySource && theTruth->Origin() != fTruthOrigin) continue;
170  int pdg = abs(part->PdgCode());
171  bool isCharged = (pdg == 11) || (pdg == 13) || (pdg == 211) || (pdg == 321) || (pdg == 2212);
172  if (!isCharged) continue;
173  if (std::find(fSkipPDGCodes.begin(), fSkipPDGCodes.end(), pdg) != fSkipPDGCodes.end()) continue;
174  float TMeV = 1000 * (part->E() - part->Mass());
175  // ignore low energy particles
176  if (TMeV < 1) continue;
177  // find the MCParticle index and stash it in trkMCP
178  unsigned int mcpi = UINT_MAX;
179  for (mcpi = 0; mcpi < (*mcpHandle).size(); ++mcpi)
180  if ((*mcpHandle)[mcpi].TrackId() == trackID) break;
181  if (mcpi == UINT_MAX) {
182  std::cout << "Failed to find a MCParticle from ParticleList\n";
183  return;
184  }
185  trkIDs.push_back(trackID);
186  MCPIs.push_back(mcpi);
187  } // ipart
188  if (trkIDs.empty()) return;
189 
190  // next construct a companion vector of MCParticle indices indexed to the full hit collection
191  fHitMCPIndex.resize((*fHitHandle).size(), UINT_MAX);
192  // Make a list of the <first, last> MC-matched hit in each TPC. This will be used
193  // to only iterate through the range of hits that are interesting
194  std::vector<std::pair<unsigned int, unsigned int>> hitRange(geom->NTPC() + 1);
195  size_t noMatch = 0;
196  size_t nMatch = 0;
197  size_t nInTPC = 0;
198  for (auto& hr : hitRange)
199  hr = std::make_pair(UINT_MAX, UINT_MAX);
200  for (size_t iht = 0; iht < (*fHitHandle).size(); ++iht) {
201  auto& hit = (*fHitHandle)[iht];
202  // only consider hits in a single TPC?
203  unsigned int tpc = hit.WireID().TPC;
204  if (fInTPC != UINT_MAX && tpc != fInTPC) continue;
205  ++nInTPC;
206  auto tides = bt_serv->HitToTrackIDEs(clockData, hit);
207  if (tides.empty()) {
208  ++noMatch;
209  continue;
210  }
211  for (auto& tide : tides) {
212  // declare a match to a MCParticle if > 50% of energy is from it
213  if (tide.energyFrac < 0.5) continue;
214  int trackID = tide.trackID;
215  // find the MCParticle index and define the Hit -> MCParticle assn
216  for (size_t indx = 0; indx < trkIDs.size(); ++indx) {
217  if (trkIDs[indx] == trackID) {
218  fHitMCPIndex[iht] = MCPIs[indx];
219  break;
220  }
221  } // indx
222  break;
223  } // tide
224  if (fHitMCPIndex[iht] == UINT_MAX) continue;
225  ++nMatch;
226  // populate hitRange
227  if (hitRange[tpc].first == UINT_MAX) hitRange[tpc].first = iht;
228  hitRange[tpc].second = iht;
229  } // iht
230 
231  if (nMatch == 0) {
232  fHitMCPIndex.resize(0);
233  return;
234  }
235 
236  if (fPrintLevel > 3) {
237  // Print the gory details
238  mf::LogVerbatim myprt("ClusterTrackAna");
239  myprt << "Checking " << trkIDs.size() << " selected MCParticles with the requested TruthOrigin";
240  if (fInTPC == UINT_MAX) { myprt << " in all TPCs\n"; }
241  else {
242  myprt << " in TPC " << fInTPC;
243  }
244  myprt << " in Run " << evt.run() << " Event " << evt.event();
245  myprt << "\n";
246  myprt << "Found " << nMatch << " MC-matched hits with the requested TruthOrigin";
247  myprt << " out of " << nInTPC << " total hits.\n";
248  myprt << "Found " << noMatch << " hits not matched to ANY MCParticle.\n";
249  // count the number of hits matched to each MCParticle in the list
250  // mcpi count
251  std::vector<std::pair<unsigned int, unsigned int>> mcpCnt;
252  for (auto mcpi : fHitMCPIndex) {
253  if (mcpi == UINT_MAX) continue;
254  unsigned short mIndx = 0;
255  for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
256  if (mcpCnt[mIndx].first == mcpi) break;
257  if (mIndx == mcpCnt.size()) mcpCnt.push_back(std::make_pair(mcpi, 0));
258  // increment the count of MC-matched hits
259  ++mcpCnt[mIndx].second;
260  } // mcpi
261  myprt << " MCPI TrackID PDGCode T(MeV) nHits Process";
262  for (auto mcpcnt : mcpCnt) {
263  myprt << "\n";
264  unsigned int mcpi = mcpcnt.first;
265  auto& mcp = (*mcpHandle)[mcpi];
266  myprt << std::setw(5) << mcpi;
267  myprt << std::setw(8) << mcp.TrackId();
268  myprt << std::setw(8) << abs(mcp.PdgCode());
269  int TMeV = 1000 * (mcp.E() - mcp.Mass());
270  myprt << std::setw(9) << TMeV;
271  myprt << std::setw(8) << mcpcnt.second;
272  myprt << " " << mcp.Process();
273  } // mcpcnt
274  } // fPrintLevel > 3
275 
276  // fill a vector of hits indices from all clusters or all tracks
277  std::vector<std::vector<unsigned int>> recoHits;
278  // and a companion vector of indices into the cluster or track collection
279  std::vector<unsigned int> recoIndex;
280  // handles to these collections to enable printing the cluster or track ID. The
281  // handle will be valid for either clusters or tracks
282  art::Handle<std::vector<recob::Cluster>> inputClusters;
283  art::Handle<std::vector<recob::Track>> inputTracks;
284  if (fClusterModuleLabel.label() != "NA") {
285  // get MC-matched hits from clusters
286  if (!evt.getByLabel(fClusterModuleLabel, inputClusters))
287  throw cet::exception("ClusterTrackAna")
288  << "Failed to get a handle to the cluster collection '" << fClusterModuleLabel.label()
289  << "'\n";
290  art::FindManyP<recob::Hit> hitsFromCls(inputClusters, evt, fClusterModuleLabel);
291  if (!hitsFromCls.isValid())
292  throw cet::exception("ClusterTrackAna") << "Failed to get a handle to Cluster -> Hit assns\n";
293  for (unsigned int icl = 0; icl < (*inputClusters).size(); ++icl) {
294  std::vector<art::Ptr<recob::Hit>> cluhits = hitsFromCls.at(icl);
295  if (cluhits.empty()) continue;
296  if (fCompareProductIDs) {
297  if (cluhits[0].id() != fHitHandle.id())
298  throw cet::exception("ClusterTrackAna")
299  << "Hits associated with ClusterModuleLabel are in a different collection than "
300  "HitModuleLabel.\n";
301  fCompareProductIDs = false;
302  } // fCompareProductIDs
303  // only load MC-matched hits. Hits that are not MC-matched were either matched to a
304  // MCParticle that was not selected or were mis-reconstructed by the hit finder. Neither
305  // of these conditions warrant penalizing the reconstruction module performance
306  std::vector<unsigned int> hitsIndex;
307  for (auto& cluhit : cluhits) {
308  if (fHitMCPIndex[cluhit.key()] != UINT_MAX) hitsIndex.push_back(cluhit.key());
309  }
310  if (hitsIndex.empty()) continue;
311  recoIndex.push_back(icl);
312  recoHits.push_back(hitsIndex);
313  } // icl
314  }
315  else {
316  // get MC-matched hits from tracks
317  if (!evt.getByLabel(fTrackModuleLabel, inputTracks))
318  throw cet::exception("ClusterTrackAna")
319  << "Failed to get a handle to the track collection '" << fTrackModuleLabel.label() << "'\n";
320  art::FindManyP<recob::Hit> hitsFromTrk(inputTracks, evt, fTrackModuleLabel);
321  if (!hitsFromTrk.isValid())
322  throw cet::exception("ClusterTrackAna") << "Failed to get a handle to Track -> Hit assns\n";
323  for (unsigned int itk = 0; itk < (*inputTracks).size(); ++itk) {
324  std::vector<art::Ptr<recob::Hit>> trkhits = hitsFromTrk.at(itk);
325  if (trkhits.empty()) continue;
326  if (fCompareProductIDs) {
327  if (trkhits[0].id() != fHitHandle.id())
328  throw cet::exception("ClusterTrackAna")
329  << "Hits associated with TrackModuleLabel are in a different collection than "
330  "HitModuleLabel.\n";
331  fCompareProductIDs = false;
332  } // fCompareProductIDs
333  std::vector<unsigned int> hitsIndex;
334  for (auto& trkhit : trkhits) {
335  if (fHitMCPIndex[trkhit.key()] != UINT_MAX) hitsIndex.push_back(trkhit.key());
336  }
337  if (hitsIndex.empty()) continue;
338  recoIndex.push_back(itk);
339  recoHits.push_back(hitsIndex);
340  } // itk
341  } // get hits from tracks
342 
343  for (const auto& tpcid : geom->IterateTPCIDs()) {
344  unsigned int tpc = tpcid.TPC;
345  if (hitRange[tpc].first == UINT_MAX) continue;
346  // iterate over planes
347  for (unsigned short plane = 0; plane < geom->Nplanes(); ++plane) {
348  unsigned int tpcMatHitCnt = 0;
349  unsigned int tpcTotHitCnt = 0;
350  // create a list of (MCParticle index, matched hit count> pairs
351  // mcParticle plane
352  std::vector<std::pair<unsigned int, float>> mcpCnt;
353  // count MCParticle, Cluster/Track hit counts - size matched to mcpCnt
354  std::vector<std::vector<std::pair<unsigned int, float>>> mcpClsCnt;
355  for (unsigned int iht = hitRange[tpc].first; iht <= hitRange[tpc].second; ++iht) {
356  auto& hit = (*fHitHandle)[iht];
357  if (hit.WireID().TPC != tpc) continue;
358  if (hit.WireID().Plane != plane) continue;
359  ++tpcTotHitCnt;
360  // ignore hits that were either unmatched to the selected set of PDGCodes
361  // or have a not-requested origin
362  if (fHitMCPIndex[iht] == UINT_MAX) continue;
363  ++tpcMatHitCnt;
364  unsigned int mcpi = fHitMCPIndex[iht];
365  unsigned short mIndx = 0;
366  for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
367  if (mcpCnt[mIndx].first == mcpi) break;
368  if (mIndx == mcpCnt.size()) {
369  // found a MCParticle not in the list yet, so add it
370  mcpCnt.push_back(std::make_pair(mcpi, 0));
371  mcpClsCnt.resize(mcpCnt.size());
372  }
373  // increment the number of MC-matched hits
374  ++mcpCnt[mIndx].second;
375  } // iht
376  // ignore TPCs/planes with few MC-matched hits
377  if (fPrintLevel > 2) {
378  mf::LogVerbatim myprt("ClusterTrackAna");
379  myprt << "TPC:" << tpc << " Plane:" << plane << " has " << tpcMatHitCnt << "/"
380  << tpcTotHitCnt;
381  myprt << " (MC-matched hits) / (all hits)";
382  }
383  if (tpcMatHitCnt < 3) continue;
384  // next iterate over all clusters/tracks and count mc-matched hits that are in this TPC/plane
385  for (unsigned int ii = 0; ii < recoHits.size(); ++ii) {
386  float nRecoHitsInPlane = 0;
387  float nRecoMatHitsInPlane = 0;
388  for (auto iht : recoHits[ii]) {
389  auto& hit = (*fHitHandle)[iht];
390  if (hit.WireID().TPC != tpc) continue;
391  if (hit.WireID().Plane != plane) continue;
392  ++nRecoHitsInPlane;
393  if (fHitMCPIndex[iht] == UINT_MAX) continue;
394  ++nRecoMatHitsInPlane;
395  unsigned int mcpi = fHitMCPIndex[iht];
396  // find the MCParticle index in mcpCnt and use it to count the match entry
397  // in mcpClsCnt
398  unsigned short mIndx = 0;
399  for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
400  if (mcpCnt[mIndx].first == mcpi) break;
401  if (mIndx == mcpCnt.size()) {
402  std::cout << "Logic error: fHitMCPIndex = " << fHitMCPIndex[iht]
403  << " is valid but isn't in the list of MCParticles to use. Please send email "
404  "to baller@fnal.gov.\n";
405  continue;
406  }
407  unsigned short cIndx = 0;
408  for (cIndx = 0; cIndx < mcpClsCnt[mIndx].size(); ++cIndx)
409  if (mcpClsCnt[mIndx][cIndx].first == ii) break;
410  if (cIndx == mcpClsCnt[mIndx].size()) mcpClsCnt[mIndx].push_back(std::make_pair(ii, 0));
411  ++mcpClsCnt[mIndx][cIndx].second;
412  } // cluhit
413  if (nRecoMatHitsInPlane == 0) continue;
414  } // ii
415  // find the cluster that has the most hits matched to each MC Particle
416  for (unsigned short mIndx = 0; mIndx < mcpCnt.size(); ++mIndx) {
417  // require at least 3 MC-matched hits
418  if (mcpCnt[mIndx].second < 3) continue;
419  unsigned int mcpi = mcpCnt[mIndx].first;
420  auto& mcp = (*mcpHandle)[mcpi];
421  unsigned int pdgCode = abs(mcp.PdgCode());
422  unsigned short pIndx = USHRT_MAX;
423  if (pdgCode == 11) pIndx = 0;
424  if (pdgCode == 13) pIndx = 1;
425  if (pdgCode == 211) pIndx = 2;
426  if (pdgCode == 321) pIndx = 3;
427  if (pdgCode == 2212) pIndx = 4;
428  if (mcpClsCnt[mIndx].empty()) {
429  // Un-reconstructed MCParticle hits
430  ++Cnts[pIndx];
431  ++fNBadEP;
432  if (fPrintLevel > 0) {
433  mf::LogVerbatim myprt("ClusterTrackAna");
434  myprt << " MCPI " << mcpi;
435  int TMeV = 1000 * (mcp.E() - mcp.Mass());
436  myprt << " " << TMeV << " MeV";
437  std::string partName = std::to_string(pdgCode);
438  if (pdgCode == 11) partName = "El";
439  if (pdgCode == 13) partName = "Mu";
440  if (pdgCode == 211) partName = "Pi";
441  if (pdgCode == 311) partName = "Ka";
442  if (pdgCode == 2212) partName = "Pr";
443  myprt << std::setw(5) << partName;
444  // print out the range of truth-matched hits
445  unsigned int firstHitIndex = UINT_MAX;
446  unsigned int lastHitIndex = UINT_MAX;
447  FirstLastHitInPlane(tpc, plane, mcpi, firstHitIndex, lastHitIndex);
448  myprt << " Failed to reconstruct. Truth-matched hit range from ";
449  myprt << HitLocation(firstHitIndex);
450  myprt << " to ";
451  myprt << HitLocation(lastHitIndex);
452  myprt << " <- EP = 0!";
453  } // fPrintLevel > 0
454  continue;
455  } // (mcpClsCnt[mIndx].empty()
456  std::pair<unsigned int, float> big = std::make_pair(UINT_MAX, 0);
457  for (unsigned short cIndx = 0; cIndx < mcpClsCnt[mIndx].size(); ++cIndx) {
458  auto& mcc = mcpClsCnt[mIndx][cIndx];
459  if (mcc.second > big.second) big = mcc;
460  } // cIndx
461  if (big.first == UINT_MAX) {
462  if (fPrintLevel > 2) {
463  mf::LogVerbatim myprt("ClusterTrackAna");
464  unsigned int mcpi = mcpCnt[mIndx].first;
465  auto& mcp = (*mcpHandle)[mcpi];
466  myprt << "Match failed: mcpi " << mcpi << " pdg " << mcp.PdgCode();
467  }
468  std::cout << "match failed mIndx " << mIndx << "\n";
469  continue;
470  } // big.first == UINT_MAX
471  unsigned int ii = big.first;
472  float eff = big.second / mcpCnt[mIndx].second;
473  float nRecoHitsInPlane = 0;
474  // define some variables to print the range of the cluster/track
475  unsigned int firstRecoHitIndex = UINT_MAX;
476  unsigned int lastRecoHitIndex = UINT_MAX;
477  for (auto iht : recoHits[ii]) {
478  auto& hit = (*fHitHandle)[iht];
479  if (hit.WireID().TPC != tpc) continue;
480  if (hit.WireID().Plane != plane) continue;
481  ++nRecoHitsInPlane;
482  if (firstRecoHitIndex == UINT_MAX) {
483  firstRecoHitIndex = iht;
484  lastRecoHitIndex = iht;
485  }
486  unsigned int wire = (*fHitHandle)[iht].WireID().Wire;
487  if (wire < (*fHitHandle)[firstRecoHitIndex].WireID().Wire) firstRecoHitIndex = iht;
488  if (wire > (*fHitHandle)[lastRecoHitIndex].WireID().Wire) lastRecoHitIndex = iht;
489  } // iht
490  float pur = big.second / nRecoHitsInPlane;
491  ++Cnts[pIndx];
492  float ep = eff * pur;
493  EPSums[pIndx] += ep;
494  ESums[pIndx] += eff;
495  PSums[pIndx] += pur;
496  bool hasBadEP = (ep < fBadEP);
497  if (hasBadEP) ++fNBadEP;
498  bool prt = fPrintLevel > 1 || (fPrintLevel == 1 && hasBadEP);
499  if (prt) {
500  mf::LogVerbatim myprt("ClusterTrackAna");
501  if (fFirstPrint) {
502  myprt << "Hit location format is TPC:Plane:Wire:Tick\n";
503  myprt << " MCPI Ptcl T(MeV) nMCPHits ____From____ _____To_____";
504  if (inputClusters.isValid()) { myprt << " clsID"; }
505  else {
506  myprt << " trkID";
507  }
508  myprt
509  << " ____From____ _____To_____ nRecoHits nMCPRecoHits Eff Pur EP\n";
510  fFirstPrint = false;
511  } // fFirstPrint
512  myprt << std::setw(5) << mcpi;
513  // convert the PDG code into nicer format
514  std::string partName = std::to_string(pdgCode);
515  if (pdgCode == 11) partName = " Electron";
516  if (pdgCode == 13) partName = " Muon";
517  if (pdgCode == 211) partName = " Pion";
518  if (pdgCode == 311) partName = " Kaon";
519  if (pdgCode == 2212) partName = " Proton";
520  myprt << partName;
521  int TMeV = 1000 * (mcp.E() - mcp.Mass());
522  myprt << std::setw(7) << TMeV;
523  // nMCPHits
524  myprt << std::setw(10) << mcpCnt[mIndx].second;
525  unsigned int firstTruHitIndex = UINT_MAX;
526  unsigned int lastTruHitIndex = UINT_MAX;
527  FirstLastHitInPlane(tpc, plane, mcpi, firstTruHitIndex, lastTruHitIndex);
528  myprt << std::setw(14) << HitLocation(firstTruHitIndex);
529  myprt << std::setw(14) << HitLocation(lastTruHitIndex);
530  int id = -1;
531  if (inputClusters.isValid()) {
532  // print cluster info
533  auto& cls = (*inputClusters)[recoIndex[ii]];
534  id = cls.ID();
535  }
536  else if (inputTracks.isValid()) {
537  // print track info
538  auto& trk = (*inputTracks)[recoIndex[ii]];
539  id = trk.ID();
540  }
541  myprt << std::setw(6) << id;
542  myprt << std::setw(14) << HitLocation(firstRecoHitIndex);
543  myprt << std::setw(14) << HitLocation(lastRecoHitIndex);
544  myprt << std::setw(12) << (int)nRecoHitsInPlane;
545  myprt << std::setw(13) << (int)big.second;
546  myprt << std::fixed << std::setprecision(2);
547  myprt << std::setw(8) << eff;
548  myprt << std::setw(8) << pur;
549  myprt << std::setw(8) << ep;
550  if (hasBadEP) myprt << " <- BadEP";
551  myprt << " Evt: " << evt.event();
552  myprt << " Evt Cnt " << fEventCnt;
553  } // prt
554  } // mIndx
555  } // plane
556  } // tpcid
557 
558  fHitMCPIndex.resize(0);
559 
560 } // analyze
second
second
Definition: detectorproperties_icarus.fcl:8
tca::TCEvent::event
unsigned int event
Definition: DataStructs.h:634
pdg
var pdg
Definition: selectors.fcl:14
tca::TCEvent::run
unsigned int run
Definition: DataStructs.h:635
cluster::ClusterTrackAna::PSums
std::array< float, 5 > PSums
Definition: ClusterTrackAna_module.cc:79
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
hit
process_name hit
Definition: cheaterreco.fcl:51
cluster::ClusterTrackAna::HitLocation
std::string HitLocation(unsigned int iht)
packs hit WireID and PeakTime into a compact format
Definition: ClusterTrackAna_module.cc:564
cluster::ClusterTrackAna::Cnts
std::array< float, 5 > Cnts
Definition: ClusterTrackAna_module.cc:74
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
cluster::ClusterTrackAna::EPSums
std::array< float, 5 > EPSums
Definition: ClusterTrackAna_module.cc:75
cluster::ClusterTrackAna::FirstLastHitInPlane
void FirstLastHitInPlane(unsigned int tpc, unsigned int plane, unsigned int mcpi, unsigned int &firstHitIndex, unsigned int &lastHitIndex)
Definition: ClusterTrackAna_module.cc:575
icarus::trigger::to_string
std::string to_string(WindowPattern const &pattern)
Definition: WindowPattern.h:115
cluster::ClusterTrackAna::fCompareProductIDs
bool fCompareProductIDs
compare Hit and Cluster-&gt; Hit art product IDs on the first event
Definition: ClusterTrackAna_module.cc:84
cluster::ClusterTrackAna::fHitHandle
art::Handle< std::vector< recob::Hit > > fHitHandle
Definition: ClusterTrackAna_module.cc:81
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
cluster::ClusterTrackAna::ESums
std::array< float, 5 > ESums
Definition: ClusterTrackAna_module.cc:77
icarus::ns::util::empty
bool empty(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:555
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
cluster::ClusterTrackAna::fHitMCPIndex
std::vector< unsigned int > fHitMCPIndex
Definition: ClusterTrackAna_module.cc:82
void cluster::ClusterTrackAna::endJob ( )
overrideprivate

Definition at line 602 of file ClusterTrackAna_module.cc.

603 {
604  // output results
605  mf::LogVerbatim myprt("ClusterTrackAna");
606  myprt << "ClusterTrackAna summary results for " << fEventCnt;
607  if (fClusterModuleLabel.label() != "NA") {
608  myprt << " events using ClusterModuleLabel: " << fClusterModuleLabel.label();
609  }
610  else {
611  myprt << " events using TrackModuleLabel: " << fTrackModuleLabel.label();
612  }
613  myprt << " Origin: " << fTruthOrigin;
614  if (fInTPC != UINT_MAX) { myprt << " in TPC " << fInTPC; }
615  else {
616  myprt << " in all TPCs";
617  }
618  myprt << "\n";
619  float cnts = 0;
620  for (unsigned short pIndx = 0; pIndx < 5; ++pIndx)
621  cnts += Cnts[pIndx];
622  if (cnts == 0) {
623  myprt << "No ClusterTrackAna results";
624  return;
625  }
626  float sumEP = 0;
627  float sumE = 0;
628  float sumP = 0;
629  myprt << "Efficiency (Eff), Purity (Pur) and Eff * Pur (EP) by selected truth particle types\n";
630  std::array<std::string, 5> pName = {{"El", "Mu", "Pi", "K ", "P "}};
631  myprt << "particle Eff Pur EP\n";
632  myprt << "-------------------------------";
633  for (unsigned short pIndx = 0; pIndx < 5; ++pIndx) {
634  if (Cnts[pIndx] == 0) continue;
635  float ave;
636  myprt << "\n " << pName[pIndx] << " ";
637  myprt << std::fixed << std::setprecision(3);
638  ave = ESums[pIndx] / Cnts[pIndx];
639  myprt << std::setw(8) << ave;
640  ave = PSums[pIndx] / Cnts[pIndx];
641  myprt << std::setw(8) << ave;
642  ave = EPSums[pIndx] / Cnts[pIndx];
643  myprt << std::setw(8) << ave;
644  if (pIndx == 0) continue;
645  sumEP += EPSums[pIndx];
646  sumE += ESums[pIndx];
647  sumP += PSums[pIndx];
648  } // pIndx
649  if (cnts == 0) return;
650  myprt << "\n";
651  myprt << "Averages for all selected truth particles\n";
652  myprt << "Ave Eff " << sumE / cnts;
653  myprt << " Ave Pur " << sumP / cnts;
654  myprt << " Ave EP " << sumEP / cnts;
655  myprt << " nBadEP " << fNBadEP;
656  myprt << " (EP < " << std::fixed << std::setprecision(2) << fBadEP << ")";
657  myprt << "\n";
658  myprt << "MCParticle counts in all TPCs and Planes:";
659  for (unsigned short pIndx = 0; pIndx < 5; ++pIndx) {
660  if (Cnts[pIndx] == 0) continue;
661  myprt << " " << pName[pIndx] << " " << (int)Cnts[pIndx];
662  } // pIndx
663 } // endJob
cluster::ClusterTrackAna::PSums
std::array< float, 5 > PSums
Definition: ClusterTrackAna_module.cc:79
cluster::ClusterTrackAna::Cnts
std::array< float, 5 > Cnts
Definition: ClusterTrackAna_module.cc:74
cluster::ClusterTrackAna::EPSums
std::array< float, 5 > EPSums
Definition: ClusterTrackAna_module.cc:75
cluster::ClusterTrackAna::ESums
std::array< float, 5 > ESums
Definition: ClusterTrackAna_module.cc:77
void cluster::ClusterTrackAna::FirstLastHitInPlane ( unsigned int  tpc,
unsigned int  plane,
unsigned int  mcpi,
unsigned int &  firstHitIndex,
unsigned int &  lastHitIndex 
)
private

Definition at line 575 of file ClusterTrackAna_module.cc.

580 {
581  // Returns the index of the first hit (lowest wire number) and last hit (highest wire number)
582  // matched to the MCParticle indexed by mcpi in the requested tpc, plane
583  firstHitIndex = UINT_MAX;
584  lastHitIndex = UINT_MAX;
585  for (unsigned int iht = 0; iht < (*fHitHandle).size(); ++iht) {
586  if (fHitMCPIndex[iht] != mcpi) continue;
587  auto& hit = (*fHitHandle)[iht];
588  if (hit.WireID().TPC != tpc) continue;
589  if (hit.WireID().Plane != plane) continue;
590  if (firstHitIndex == UINT_MAX) {
591  firstHitIndex = iht;
592  lastHitIndex = iht;
593  }
594  unsigned int wire = (*fHitHandle)[iht].WireID().Wire;
595  if (wire < (*fHitHandle)[firstHitIndex].WireID().Wire) firstHitIndex = iht;
596  if (wire > (*fHitHandle)[lastHitIndex].WireID().Wire) lastHitIndex = iht;
597  } // iht
598 } // FirstLastHitInPlane
hit
process_name hit
Definition: cheaterreco.fcl:51
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
cluster::ClusterTrackAna::fHitHandle
art::Handle< std::vector< recob::Hit > > fHitHandle
Definition: ClusterTrackAna_module.cc:81
cluster::ClusterTrackAna::fHitMCPIndex
std::vector< unsigned int > fHitMCPIndex
Definition: ClusterTrackAna_module.cc:82
std::string cluster::ClusterTrackAna::HitLocation ( unsigned int  iht)
private

packs hit WireID and PeakTime into a compact format

Definition at line 564 of file ClusterTrackAna_module.cc.

565 {
566  // Put the hit location into a compact human-readable format
567  if (iht >= (*fHitHandle).size()) return "NA";
568  auto& hit = (*fHitHandle)[iht];
569  return std::to_string(hit.WireID().TPC) + ":" + std::to_string(hit.WireID().Plane) + ":" +
570  std::to_string(hit.WireID().Wire) + ":" + std::to_string((int)hit.PeakTime());
571 } // HitLocation
hit
process_name hit
Definition: cheaterreco.fcl:51
icarus::trigger::to_string
std::string to_string(WindowPattern const &pattern)
Definition: WindowPattern.h:115
ClusterTrackAna& cluster::ClusterTrackAna::operator= ( ClusterTrackAna const &  )
delete
ClusterTrackAna& cluster::ClusterTrackAna::operator= ( ClusterTrackAna &&  )
delete

Member Data Documentation

std::array<float, 5> cluster::ClusterTrackAna::Cnts {{0}}
private

Definition at line 74 of file ClusterTrackAna_module.cc.

std::array<float, 5> cluster::ClusterTrackAna::EPSums {{0}}
private

Definition at line 75 of file ClusterTrackAna_module.cc.

std::array<float, 5> cluster::ClusterTrackAna::ESums {{0}}
private

Definition at line 77 of file ClusterTrackAna_module.cc.

float cluster::ClusterTrackAna::fBadEP
private

Definition at line 67 of file ClusterTrackAna_module.cc.

art::InputTag cluster::ClusterTrackAna::fClusterModuleLabel
private

Definition at line 61 of file ClusterTrackAna_module.cc.

bool cluster::ClusterTrackAna::fCompareProductIDs
private
Initial value:
{
true}

compare Hit and Cluster-> Hit art product IDs on the first event

Definition at line 84 of file ClusterTrackAna_module.cc.

unsigned int cluster::ClusterTrackAna::fCurrentRun {0}
private

Definition at line 87 of file ClusterTrackAna_module.cc.

unsigned int cluster::ClusterTrackAna::fEventCnt {0}
private

Definition at line 70 of file ClusterTrackAna_module.cc.

bool cluster::ClusterTrackAna::fFirstPrint {true}
private

Definition at line 86 of file ClusterTrackAna_module.cc.

art::Handle<std::vector<recob::Hit> > cluster::ClusterTrackAna::fHitHandle
private

Definition at line 81 of file ClusterTrackAna_module.cc.

std::vector<unsigned int> cluster::ClusterTrackAna::fHitMCPIndex
private

Definition at line 82 of file ClusterTrackAna_module.cc.

art::InputTag cluster::ClusterTrackAna::fHitModuleLabel
private

Definition at line 60 of file ClusterTrackAna_module.cc.

unsigned int cluster::ClusterTrackAna::fInTPC
private

Definition at line 66 of file ClusterTrackAna_module.cc.

unsigned int cluster::ClusterTrackAna::fNBadEP {0}
private

Definition at line 72 of file ClusterTrackAna_module.cc.

short cluster::ClusterTrackAna::fPrintLevel
private

Definition at line 65 of file ClusterTrackAna_module.cc.

std::vector<int> cluster::ClusterTrackAna::fSkipPDGCodes
private

Definition at line 64 of file ClusterTrackAna_module.cc.

art::InputTag cluster::ClusterTrackAna::fTrackModuleLabel
private

Definition at line 62 of file ClusterTrackAna_module.cc.

simb::Origin_t cluster::ClusterTrackAna::fTruthOrigin
private

Definition at line 63 of file ClusterTrackAna_module.cc.

std::array<float, 5> cluster::ClusterTrackAna::PSums {{0}}
private

Definition at line 79 of file ClusterTrackAna_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.5