Public Member Functions | Private Member Functions | Private Attributes | List of all members
cluster::ClusterTrackAna Class Reference
Inheritance diagram for cluster::ClusterTrackAna:
art::EDAnalyzer art::detail::Analyzer art::detail::LegacyModule art::Observer art::ModuleBase

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
 
- Public Member Functions inherited from art::EDAnalyzer
 EDAnalyzer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDAnalyzer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Analyzer
virtual ~Analyzer () noexcept
 
 Analyzer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 Analyzer (Table< Config > const &config)
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Observer
 ~Observer () noexcept
 
 Observer (Observer const &)=delete
 
 Observer (Observer &&)=delete
 
Observer & operator= (Observer const &)=delete
 
Observer & operator= (Observer &&)=delete
 
void registerProducts (ProductDescriptions &, ModuleDescription const &)
 
void fillDescriptions (ModuleDescription const &)
 
fhicl::ParameterSetID selectorConfig () const
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

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}
 

Additional Inherited Members

- Public Types inherited from art::EDAnalyzer
using WorkerType = WorkerT< EDAnalyzer >
 
using ModuleType = EDAnalyzer
 
- Protected Member Functions inherited from art::Observer
std::string const & processName () const
 
bool wantAllEvents () const noexcept
 
bool wantEvent (ScheduleID id, Event const &e) const
 
Handle< TriggerResults > getTriggerResults (Event const &e) const
 
 Observer (fhicl::ParameterSet const &config)
 
 Observer (std::vector< std::string > const &select_paths, std::vector< std::string > const &reject_paths, fhicl::ParameterSet const &config)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 43 of file ClusterTrackAna_module.cc.

Constructor & Destructor Documentation

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

Definition at line 96 of file ClusterTrackAna_module.cc.

96  : EDAnalyzer{pset}
97 {
98  fHitModuleLabel = pset.get<art::InputTag>("HitModuleLabel");
99  bool validModuleLabel = false;
100  fClusterModuleLabel = "NA";
101  fTrackModuleLabel = "NA";
102  fClusterModuleLabel = pset.get<art::InputTag>("ClusterModuleLabel", "NA");
103  if (fClusterModuleLabel != "NA") validModuleLabel = true;
104  fTrackModuleLabel = pset.get<art::InputTag>("TrackModuleLabel", "NA");
105  if (validModuleLabel && fTrackModuleLabel != "NA")
106  throw cet::exception("ClusterTrackAna")
107  << "You must specify either a ClusterModuleLabel OR a TrackModuleLabel\n";
108  if (!validModuleLabel && fTrackModuleLabel != "NA") validModuleLabel = true;
109  // origin = 0 (anything), 1(nu), 2(cosmics), 3(SN nu), 4(SingleParticle)
110  int tmp = pset.get<int>("TruthOrigin", 0);
111  fTruthOrigin = (simb::Origin_t)tmp;
112  fPrintLevel = pset.get<short>("PrintLevel", 0);
113  if (pset.has_key("SkipPDGCodes")) fSkipPDGCodes = pset.get<std::vector<int>>("SkipPDGCodes");
114  fBadEP = pset.get<float>("BadEP", 0.);
115  fInTPC = UINT_MAX;
116  int intpc = pset.get<int>("InTPC", -1);
117  if (intpc >= 0) fInTPC = intpc;
118  // do some initialization
119  Cnts.fill(0.);
120  EPSums.fill(0.);
121  ESums.fill(0.);
122  PSums.fill(0.);
123 } // ClusterTrackAna constructor
cluster::ClusterTrackAna::PSums
std::array< float, 5 > PSums
Definition: ClusterTrackAna_module.cc:77
simb::Origin_t
enum simb::_ev_origin Origin_t
event origin types
cluster::ClusterTrackAna::Cnts
std::array< float, 5 > Cnts
Definition: ClusterTrackAna_module.cc:72
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
cluster::ClusterTrackAna::EPSums
std::array< float, 5 > EPSums
Definition: ClusterTrackAna_module.cc:73
languages.tmp
string tmp
Definition: languages.py:63
cluster::ClusterTrackAna::ESums
std::array< float, 5 > ESums
Definition: ClusterTrackAna_module.cc:75
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
cluster::ClusterTrackAna::ClusterTrackAna ( ClusterTrackAna const &  )
delete
cluster::ClusterTrackAna::ClusterTrackAna ( ClusterTrackAna &&  )
delete

Member Function Documentation

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

Implements art::EDAnalyzer.

Definition at line 127 of file ClusterTrackAna_module.cc.

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

Reimplemented from art::EDAnalyzer.

Definition at line 600 of file ClusterTrackAna_module.cc.

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

Definition at line 573 of file ClusterTrackAna_module.cc.

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

packs hit WireID and PeakTime into a compact format

Definition at line 562 of file ClusterTrackAna_module.cc.

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
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
art::to_string
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
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 72 of file ClusterTrackAna_module.cc.

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

Definition at line 73 of file ClusterTrackAna_module.cc.

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

Definition at line 75 of file ClusterTrackAna_module.cc.

float cluster::ClusterTrackAna::fBadEP
private

Definition at line 65 of file ClusterTrackAna_module.cc.

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

Definition at line 59 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 82 of file ClusterTrackAna_module.cc.

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

Definition at line 85 of file ClusterTrackAna_module.cc.

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

Definition at line 68 of file ClusterTrackAna_module.cc.

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

Definition at line 84 of file ClusterTrackAna_module.cc.

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

Definition at line 79 of file ClusterTrackAna_module.cc.

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

Definition at line 80 of file ClusterTrackAna_module.cc.

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

Definition at line 58 of file ClusterTrackAna_module.cc.

unsigned int cluster::ClusterTrackAna::fInTPC
private

Definition at line 64 of file ClusterTrackAna_module.cc.

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

Definition at line 70 of file ClusterTrackAna_module.cc.

short cluster::ClusterTrackAna::fPrintLevel
private

Definition at line 63 of file ClusterTrackAna_module.cc.

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

Definition at line 62 of file ClusterTrackAna_module.cc.

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

Definition at line 60 of file ClusterTrackAna_module.cc.

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

Definition at line 61 of file ClusterTrackAna_module.cc.

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

Definition at line 77 of file ClusterTrackAna_module.cc.

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