BackTracker_service.cc
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////
2 // $Id: BackTracker_service.cc,v 1.7 2012-12-06 15:03:39 lein Exp $
3 //
4 //
5 // \file: BackTracker_service.cc
6 //
7 // brebel@fnal.gov
8 //
9 ////////////////////////////////////////////////////////////////////////
10 
11 #include "MCCheater/BackTracker.h"
12 
13 #include <map>
14 #include <algorithm>
15 
16 // Framework includes
22 
23 #include "Simulation/FLSHitList.h"
24 #include "Geometry/Geometry.h"
25 
26 // need a map of CellUniqueId to vector of FLSHits to make this work
27 // really should just have a CellHit/RawDigit to provide the Cell and Plane
28 // and then look through the event to gett the proper FLSHits
29 // have to know where the FLSHits are stored in the event if want to pass art::Event around
30 
31 namespace cheat{
32 
33  //----------------------------------------------------------------------
34  static bool sortFLSHit(sim::FLSHit const& a, sim::FLSHit const& b)
35  {
36  return a.GetEdep() > b.GetEdep();
37  }
38 
39  //----------------------------------------------------------------------
40  bool CompareByEnergy(const TrackIDE& a, const TrackIDE& b)
41  {
42  return a.energy < b.energy;
43  }
44 
45  //----------------------------------------------------------------------
47  {
48  return ep.efficiency;
49  }
50 
51  //----------------------------------------------------------------------
53  {
54  return ep.purity;
55  }
56 
57  //----------------------------------------------------------------------
59  {
60  return ep.efficiency*ep.purity;
61  }
62 
63  //----------------------------------------------------------------------
65  {
66  return ep.energySlice;
67  }
68 
69  //----------------------------------------------------------------------
72  : fParams(params()), fHaveTruthInfo(false)
73  {
74  reg.sPreProcessEvent.watch(this, &BackTracker::Rebuild);
76  }
77 
78  //----------------------------------------------------------------------
80  {
81  }
82 
83  //----------------------------------------------------------------------
85  {
86  // Clear out anything remaining from previous calls to Rebuild
87  fCellToFLSHit.clear();
88  fTrackIDToFLSHit.clear();
89  fCellToPhotonSignal.clear();
90  fTrackIDs.clear();
92  fMCTruthList.clear();
94  fTrueEnergyList.clear();
95  fTrackIDToTrueEIndex.clear();
96  fHaveTruthInfo = false; // Haven't got anything yet
97  fDetID = novadaq::cnv::kUNKNOWN_DET; // Haven't found detector type yet
98 
99  // do nothing if this is data
100  if(evt.isRealData() && !fParams.MRMode()
101  && !fParams.MCOverlayMode()) return;
102 
103  // get detector for CloseInTime function
105  fDetID = geom->DetId();
106 
107  // fill the map of OfflineChan to std::vector<FLSHit>
108  // first get the FLSHitList from the event
110  evt.getByLabel(fParams.GeantModuleLabel(), flslistHandle);
111 
112  // Probably this is a MC generation job and this is the preEvent callback
113  // which is too early. That's OK, we'll get called at the right time by
114  // ReadoutSim in a moment. Or we're mock data or something.
115  if(flslistHandle.failedToGet()) return;
116 
117  fHaveTruthInfo = true; // We got FLSHits, the rest should go well too
118 
119  // loop over the FLSHits and fill the std::vector for each
120  // channel
121  for( size_t l = 0; l < flslistHandle->size(); ++l){
122  const std::vector<sim::FLSHit> &hits = flslistHandle->at(l).fHits;
123  for( size_t f = 0; f < hits.size(); ++f){
124  const sim::FLSHit& hit = hits[f];
125  const geo::OfflineChan chan(hit.GetPlaneID(), hit.GetCellID());
126  fCellToFLSHit[chan].push_back(hit);
127  fTrackIDToFLSHit[hit.GetTrackID()].push_back(hit);
128  }// end loop over flshits
129  }// end loop over flshit lists
130 
131  // loop over the map and sort the vectors of flshits
132  for(auto itr = fCellToFLSHit.begin(); itr != fCellToFLSHit.end(); itr++){
133  itr->second.sort(sortFLSHit);
134  }
135 
136  // loop over the map and sort the vectors of flshits
137  for(auto pfitr = fTrackIDToFLSHit.begin(); pfitr != fTrackIDToFLSHit.end(); pfitr++){
138  pfitr->second.sort(sortFLSHit);
139  }
140 
141  // Similarly, PhotonSignals
143  evt.getByLabel(fParams.PhotonModuleLabel(), phots);
144 
145  if(!phots.failedToGet()){
146  for(size_t i = 0; i < phots->size(); ++i){
147  const sim::PhotonSignal& phot = (*phots)[i];
148  const geo::OfflineChan chan(phot.Plane(), phot.Cell());
149  fCellToPhotonSignal[chan].push_back(phot);
150  } // end for i
151 
152  for(auto it = fCellToPhotonSignal.begin(); it != fCellToPhotonSignal.end(); ++it){
153  // Sort in reverse order, ie from most photons to least
154  std::sort(it->second.rbegin(), it->second.rend(), sim::CompareByNPhoton);
155  }
156  }
157  else{
158  mf::LogWarning("BackTracker") << "unable to find PhotonTransport "
159  << "created in " << fParams.PhotonModuleLabel() << " "
160  << "any attempt to get the PhotonTransport objects from "
161  << "the backtracker will fail";
162  }
163 
165  evt.getByLabel(fParams.GeantModuleLabel(), pHandle);
166  if(pHandle.failedToGet()){
167  mf::LogWarning("BackTracker") << "unable to find sim::Particles "
168  << "created in " << fParams.GeantModuleLabel() << " "
169  << "any attempt to get the Particles from "
170  << "the backtracker will fail";
171  }
172  else{
173  // Keep track of all the MCTruths we discover associated to the particles
174  std::set<art::Ptr<simb::MCTruth>> truths;
175  // And which truth each particle points at
176  std::map<int, art::Ptr<simb::MCTruth>> trackIDToMCTruthPtr;
177 
179  for(size_t p = 0; p < pHandle->size(); ++p){
180 
181  sim::Particle *part = new sim::Particle(pHandle->at(p));
182  fParticleNav.Add(part);
183 
184  // get the simb::MCTruth associated to this vector<Particle>
185  art::Ptr<simb::MCTruth> mct = fo.at(p);
186  truths.insert(mct);
187  try{
188  trackIDToMCTruthPtr.emplace(pHandle->at(p).TrackId(), mct);
189  }
190  catch(cet::exception &ex){
191  mf::LogWarning("BackTracker") << "unable to find MCTruth from Particle list "
192  << "created in " << fParams.GeantModuleLabel() << " "
193  << "any attempt to get the MCTruth objects from "
194  << "the backtracker will fail\n"
195  << "message from caught exception:\n" << ex;
196  }
197  } // end loop over particles to get MCTruthList
198 
199  // Flatten the MCTruth set into a vector
200  fMCTruthList.insert(fMCTruthList.end(), truths.begin(), truths.end());
201 
202  // Translate all the truth pointers into indices
203  std::map<art::Ptr<simb::MCTruth>, int> truthIdxs;
204  for(unsigned int i = 0; i < fMCTruthList.size(); ++i){
205  truthIdxs.emplace(fMCTruthList[i], i);
206  }
207 
208  // Convert trackIDToMCTruthPtr from Ptrs to indices
209  for(auto it: trackIDToMCTruthPtr){
210  fTrackIDToMCTruthIndex[it.first] = truthIdxs[it.second];
211  }
212  }
213 
214 
215  // fill the set of track ids from the ParticleNavigator
216  fTrackIDs.clear();
217  for(auto& it: fParticleNav) fTrackIDs.insert(it.first);
218 
219  LOG_DEBUG("BackTracker") << fParticleNav;
220 
221  // Get the sim::TrueEnergy data products out of the event. I don't think that I'll be able to utilise the Assn
222  // between sim::Particle and sim::TrueEnergy here, as I would have to somehow have to pass the Assn around.
224  evt.getByLabel(fParams.TrueEnergyModuleLabel(), teHandle);
225  if(!teHandle.failedToGet()) {
226  // Get the TrueEnergy Assn, and loop through the particles.
228  for(size_t p = 0; p < pHandle->size(); ++p){
229  int parTrId = pHandle->at(p).TrackId();
230  // There should just be a 1 to 1 correspondance, however for some reason
231  // the first few particles may have multiple assns, and the last few none.
232  bool match_truep = false;
233  for(size_t te = 0; te < teAssn.at(p).size(); ++te){
234  if (parTrId == teAssn.at(p)[te]->TrackId()) {
235  match_truep = true;
236  fTrueEnergyList.push_back ( teAssn.at(p)[te] );
237  fTrackIDToTrueEIndex.emplace( parTrId, p );
238  }
239  }
240  if (!match_truep) {
241  // If didn't match, then it's one of the last particles, and so has an assn with a particle at the.
242  // start of the particle loop. Lets restart the loop over the particles, pulling those assns...
243  for(size_t par2 = 0; par2 < pHandle->size(); ++par2){
244  for(size_t te = 0; te < teAssn.at(par2).size(); ++te){
245  if (parTrId == teAssn.at(par2)[te]->TrackId()) {
246  fTrueEnergyList.push_back ( teAssn.at(par2)[te] );
247  fTrackIDToTrueEIndex.emplace( parTrId, p );
248  match_truep = true;
249  }
250  if (match_truep) break;
251  }
252  if (match_truep) break;
253  }
254  }
255  }
256 
257  /*
258  // I think that if using the true energy module you may have to use this technique.
259  // Totally my fault, I made the assns in different way for some reason...
260  art::FindOneP<sim::TrueEnergy> teAssn(pHandle, evt, fParams.TrueEnergyModuleLabel());
261  for(size_t p = 0; p < pHandle->size(); ++p){
262  art::Ptr<sim::TrueEnergy > MyEns = teAssn.at( p );
263  fTrueEnergyList.push_back( MyEns );
264  fTrackIDToTrueEIndex.emplace( pHandle->at(p).TrackId(), p );
265  }
266  */
267  } else if (fNumTrueEnergyWarnings<10) {
268  mf::LogWarning("BackTracker") << "unable to find sim::TrueEnergy "
269  << "created in " << fParams.TrueEnergyModuleLabel() << " "
270  << "any attempt to get the sim::TrueEnergy's from "
271  << "the backtracker will fail";
273  }
274  }
275 
276  //----------------------------------------------------------------------
277  const std::vector<sim::PhotonSignal> BackTracker::CellToPhotonSignal(unsigned int const& plane,
278  unsigned int const& cell) const
279  {
280  auto it = fCellToPhotonSignal.find(geo::OfflineChan(plane, cell));
281 
282  if(it == fCellToPhotonSignal.end()){
283  LOG_DEBUG("BackTracker") << "no collection of sim::PhotonSignal for "
284  << " plane = " << plane << " cell = " << cell
285  << " probably this was a noise hit, returning an"
286  << " empty vector";
287  const std::vector<sim::PhotonSignal> empty;
288  return empty;
289  }
290  return it->second;
291  }
292 
293  //----------------------------------------------------------------------
294  const std::vector<sim::PhotonSignal> BackTracker::HitToPhotonSignal(const art::Ptr<rb::CellHit>& hit) const
295  {
296  auto it = fCellToPhotonSignal.find(geo::OfflineChan(hit->Plane(), hit->Cell()));
297 
298  // Nothing to accumulate from this cell
299  if(it == fCellToPhotonSignal.end()){
300  LOG_DEBUG("BackTracker") << "no collection of sim::PhotonSignal for "
301  << " plane = " << hit->Plane() << " cell = " << hit->Cell()
302  << " probably this was a noise hit, returning an"
303  << " empty vector";
304  const std::vector<sim::PhotonSignal> empty;
305  return empty;
306  }
307 
308  // Candidate photons, in that they're in the right cell.
309  const std::vector<sim::PhotonSignal>& phots = it->second;
310 
311  std::vector<sim::PhotonSignal> hitPhots;
312 
313  // For every photon that's sufficiently in time with the hit,
314  // accumulate its energy into the relevant trackID.
315  for(unsigned int i = 0; i < phots.size(); ++i){
316  const sim::PhotonSignal& phot = phots[i];
317  if(CloseInTime(*hit.get(), phot)){
318  hitPhots.push_back(phot);
319  }
320  }
321  return hitPhots;
322  }
323 
324  //----------------------------------------------------------------------
325  const std::vector<sim::FLSHit> BackTracker::CellToFLSHit(unsigned int const& plane,
326  unsigned int const& cell) const
327  {
328  auto it = fCellToFLSHit.find(geo::OfflineChan(plane, cell));
329 
330  if(it == fCellToFLSHit.end()){
331  LOG_DEBUG("BackTracker") << "no collection of sim::FLSHit for "
332  << " plane = " << plane << " cell = " << cell
333  << " probably this was a noise hit, returning an"
334  << " empty vector";
335  const std::vector<sim::FLSHit> empty;
336  return empty;
337  }
338  return std::vector<sim::FLSHit>(it->second.begin(), it->second.end());
339  }
340 
341  //----------------------------------------------------------------------
342  std::vector<sim::FLSHit> BackTracker::HitToFLSHit(const rb::CellHit& hit) const
343  {
344  std::map<int, int> idToNPhot;
345  AccumulateHitContributions(hit, idToNPhot);
346 
347  // Next, sort by number of photons. This turns out to be a bit fiddly.
348 
349  std::vector<std::pair<int, int> > nphotToId;
350  nphotToId.reserve(idToNPhot.size());
351  for(std::map<int, int>::iterator it = idToNPhot.begin(); it != idToNPhot.end(); ++it){
352  nphotToId.push_back(std::make_pair(it->second, it->first));
353  }
354 
355  if(nphotToId.empty()){
356  LOG_DEBUG("BackTracker") << "No matching photons, this is probably a noise hit. Returning empty vector";
357  return std::vector<sim::FLSHit>();
358  }
359 
360  // Sort in reverse order (ie decreasing number of photons).
361  // Tuples compare by first element first, so this works
362  std::sort(nphotToId.rbegin(), nphotToId.rend());
363 
364  auto it = fCellToFLSHit.find(geo::OfflineChan(hit.Plane(), hit.Cell()));
365  if ( it == fCellToFLSHit.end() ) {
366  throw cet::exception("BackTracker")
367  << "PhotonSignal with no matching FLSHit. That shouldn't be possible\n"
368  << __FILE__ << ":" << __LINE__ << "\n";
369  }
370 
371  std::vector<sim::FLSHit> ret;
372  const std::list<sim::FLSHit>& flshits = it->second;
373 
374  ret.reserve(nphotToId.size()*flshits.size()); // Overestimate
375 
376  for(unsigned int idIdx = 0; idIdx < nphotToId.size(); ++idIdx){
377  const int trackid = nphotToId[idIdx].second;
378  for(const sim::FLSHit& fls: flshits){
379  if(fls.GetTrackID() == trackid){
380  ret.push_back(fls);
381  }
382  } // end for flsIdx
383  } // end for idIdx
384 
385  return ret;
386  }
387 
388  //----------------------------------------------------------------------
389  const double BackTracker::EnergyFromTrackId(const art::Ptr<rb::CellHit>& hit, int trackId, bool useBirksE) const
390  {
391  auto it = fCellToFLSHit.find(geo::OfflineChan(hit->Plane(), hit->Cell()));
392  if(it == fCellToFLSHit.end()){
393  mf::LogWarning("BackTracker") << "No FLSHits in this cell";
394  return 0;
395  }
396 
397  const std::list<sim::FLSHit>& flshits = it->second;
398 
399  double energy = 0;
400 
401  for(const sim::FLSHit& fls: flshits){
402  if(fls.GetTrackID() == trackId){
403  energy += ((useBirksE) ? fls.GetEdepBirks() : fls.GetEdep());
404  }
405  } // end for flsIdx
406 
407  return energy;
408  }
409 
410 
411  //----------------------------------------------------------------------
413  bool quiet) const
414  {
415  return HitToParticle( *hit.get(), quiet);
416  }
417 
418  //----------------------------------------------------------------------
420  bool quiet) const
421  {
422  std::map<int, int> idToNPhot;
423  AccumulateHitContributions(hit, idToNPhot);
424 
425  // Go through the map and find the particle ID that contributed the most light
426  int bestId = -1;
427  int mostPhot = -1;
428  for(std::map<int, int>::iterator it = idToNPhot.begin(); it != idToNPhot.end(); ++it){
429  if(it->second > mostPhot){
430  mostPhot = it->second;
431  bestId = it->first;
432  }
433  }
434 
435  if(bestId == -1){
436  if(!quiet){
437  mf::LogWarning("BackTracker") << "No photons in this cell are in time with the hit. "
438  << "Returning null pointer";
439  }
440  return 0;
441  }
442 
443  return TrackIDToParticle(bestId);
444  }
445 
446  //----------------------------------------------------------------------
448  {
449  return (HitToParticle(hit, true) == 0);
450  }
451 
452  //----------------------------------------------------------------------
453  bool BackTracker::IsNoise(const rb::CellHit& hit) const
454  {
455  return (HitToParticle(hit, true) == 0);
456  }
457 
458  //----------------------------------------------------------------------
459  bool BackTracker::IsHitsNoise(std::vector< art::Ptr<rb::CellHit> > const& hits) const
460  {
461  int totalNoise = 0;
462  int total = hits.size();
463 
464  for(size_t h = 0; h < hits.size(); ++h){
466  if (IsNoise(hit)) totalNoise += 1;
467  }
468 
469  int totalPhysics = total - totalNoise;
470  double physicsFrac = (double)totalPhysics / (double)total;
471 
472  // If at least MinPhysicsHits physics hits in the collection, not a noise group.
473  // If at least MinPhysicsFrac of the hits are physics, not a noise group.
474  // Otherwise, call it a noise group.
475 
476  return((totalPhysics < fParams.MinPhysicsHits()) && (physicsFrac < fParams.MinPhysicsFrac()));
477  }
478 
479 
480  //----------------------------------------------------------------------
481  const sim::Particle* BackTracker::TrackIDToParticle(int const& id) const
482  {
484 
485  if(part_it == fParticleNav.end()){
486  mf::LogWarning("BackTracker") << "can't find particle with track id "
487  << id << " in sim::ParticleNavigator"
488  << " returning null pointer";
489  return 0;
490  }
491 
492  return part_it->second;
493  }
494 
495  //----------------------------------------------------------------------
497  {
498  // get the mother id from the particle navigator
499  // the EveId was adopted in the Rebuild method
500 
501  return this->TrackIDToParticle(fParticleNav.EveId(abs(id)));
502  }
503 
504  //----------------------------------------------------------------------
506  {
507  // find the entry in the MCTruth collection for this track id
508  size_t mct = fTrackIDToMCTruthIndex.find(abs(id))->second;
509 
510  if (mct > fMCTruthList.size() ) {
511  throw cet::exception("BackTracker")
512  << "attempting to find MCTruth index for "
513  << "out of range value: " << mct
514  << "/" << fMCTruthList.size() << "\n"
515  << __FILE__ << ":" << __LINE__ << "\n";
516  }
517  return fMCTruthList[mct];
518  }
519 
520  //----------------------------------------------------------------------
522  {
523  return this->TrackIDToMCTruth(p->TrackId());
524  }
525 
526  //----------------------------------------------------------------------
527  std::vector<const sim::Particle*> BackTracker::MCTruthToParticles(art::Ptr<simb::MCTruth> const& mct) const
528  {
529  std::vector<const sim::Particle*> ret;
530 
531  for(std::set<int>::const_iterator itr = fTrackIDs.begin(); itr != fTrackIDs.end(); itr++){
532  if( this->TrackIDToMCTruth(*itr) == mct ) ret.push_back(fParticleNav.find(*itr)->second);
533  }
534 
535  return ret;
536  }
537 
538  //----------------------------------------------------------------------
539  // this method orders the vector of particles from contributing most
540  // to the collection to least
541  std::vector<const sim::Particle*> BackTracker::
542  HitsToParticle(const std::vector<const rb::CellHit*>& hits) const
543  {
544  std::map<int, int> idToNPhot;
545  for(size_t h = 0; h < hits.size(); ++h){
546  AccumulateHitContributions(*hits[h], idToNPhot);
547  }
548 
549  // Now "all" we need to do is sort by number of photons, and discard any IDs
550  // that we can't find particles for. This turns out to be a bit fiddly.
551 
552  std::vector<std::pair<int, int> > nphotToId;
553  for(std::map<int, int>::iterator it = idToNPhot.begin(); it != idToNPhot.end(); ++it){
554  nphotToId.push_back(std::make_pair(it->second, it->first));
555  }
556 
557  // Sort in reverse order (ie increasing number of photons).
558  // Tuples compare by first element first, so this works
559  std::sort(nphotToId.rbegin(), nphotToId.rend());
560 
561  std::vector<const sim::Particle*> ret;
562  for(unsigned int n = 0; n < nphotToId.size(); ++n){
563  const sim::Particle* part = TrackIDToParticle(nphotToId[n].second);
564  if(part) ret.push_back(part);
565  }
566 
567  return ret;
568  }
569 
570  //----------------------------------------------------------------------
571  // this method orders the vector of particles from contributing most
572  // to the collection to least
573  std::vector<TrackIDE> BackTracker::
574  HitsToTrackIDE(const std::vector<const rb::CellHit*>& hits, bool useBirksE) const
575  {
576 
577  std::map<int, double> sigtracks;
578 
579  double totalE = 0;
580  for(size_t h = 0; h < hits.size(); ++h){
581  const std::vector<sim::FLSHit> fls = HitToFLSHit(*hits[h]);
582  for(size_t f = 0; f < fls.size(); ++f){
583  if( useBirksE){
584  sigtracks[fls[f].GetTrackID()] += fls[f].GetEdepBirks();
585  totalE += fls[f].GetEdepBirks();
586  }
587  else{
588  sigtracks[fls[f].GetTrackID()] += fls[f].GetEdep();
589  totalE += fls[f].GetEdep();
590  }
591  }
592  } // end for h
593 
594  // Make trackIDE objects from our collected energies
595  std::vector<TrackIDE> trackIDE;
596  for(std::map<int, double>::iterator it = sigtracks.begin(); it != sigtracks.end(); ++it){
597  TrackIDE ide;
598  ide.trackID = it->first;
599  ide.energy = it->second;
600  ide.energyFrac = it->second/totalE;
601  trackIDE.push_back(ide);
602  }
603 
604  // Sort TrackIDE objects from most energy contributed to least
605  std::sort(trackIDE.rbegin(), trackIDE.rend(), CompareByEnergy);
606 
607  return trackIDE;
608  }
609 
610  //----------------------------------------------------------------------
611  std::vector<TrackIDE> BackTracker::
612  HitToTrackIDE(const art::Ptr<rb::CellHit>& hit, bool useBirksE) const
613  {
614  return HitsToTrackIDE(std::vector<const rb::CellHit*>(1, hit.get()), useBirksE);
615  }
616 
617  //----------------------------------------------------------------------
618  std::vector<TrackIDE> BackTracker::
619  HitToTrackIDE(const rb::CellHit& hit, bool useBirksE) const
620  {
621  return HitsToTrackIDE(std::vector<const rb::CellHit*>(1, &hit), useBirksE);
622  }
623 
624  //----------------------------------------------------------------------
625  TVector3 BackTracker::HitToXYZ(art::Ptr<rb::CellHit> const& hit, bool useBirksE) const
626  {
627  double x = 0.;
628  double y = 0.;
629  double z = 0.;
630  double w = 0.;
631 
632  const std::vector<sim::FLSHit> fhits = this->HitToFLSHit(hit);
633 
634  // loop over the FLSHits for this cell and get the weighted central
635  // position of each track id in it.
636  for(const sim::FLSHit& fhit: fhits){
637  double weight = ((useBirksE) ? fhit.GetEdepBirks() : fhit.GetEdep());
638 
639  w += weight;
640  x += weight * 0.5*(fhit.GetExitX() + fhit.GetEntryX());
641  y += weight * 0.5*(fhit.GetExitY() + fhit.GetEntryY());
642  z += weight * 0.5*(fhit.GetExitZ() + fhit.GetEntryZ());
643 
644  }// end loop over sim::FLSHits
645 
646  // if the sum of the weights is still 0, then return
647  // the obviously stupid default values
648  if(w < 1.e-5) return TVector3(-999, -999, -999);
649 
650  return TVector3(x/w, y/w, z/w);
651  }
652 
653  //----------------------------------------------------------------------
654  double BackTracker::
655  HitCollectionPurity(const std::set<int>& trackIDs,
656  const std::vector<rb::WeightedHit>& whits,
657  std::map<int, double>* purMap,
658  std::map<int, int>* parents,
659  bool energyPur) const
660  {
661  const std::vector<RawWeightedHit> rwh(whits.begin(), whits.end());
662  return HitCollectionPurity(trackIDs, rwh, purMap, parents, energyPur);
663  }
664 
665  //----------------------------------------------------------------------
666  double BackTracker::
667  HitCollectionPurity(const std::set<int>& trackIDs,
668  const std::vector<RawWeightedHit>& whits,
669  std::map<int, double>* purMap,
670  std::map<int, int>* parents,
671  bool energyPur) const
672 
673  {
674  // get the list of EveIDs that correspond to the hits in this collection
675  // if the EveID shows up in the input list of trackIDs, then it counts
676  double total = whits.size();
677  double desired = 0;
678 
679  double desiredE = 0;
680  double totalE = 0;
681 
682  if(purMap) purMap->clear();
683 
684  for(size_t h = 0; h < whits.size(); ++h){
685 
686  double totalEInCell = 0;
687  double desiredEInCell = 0;
688 
689  const rb::CellHit* hit = whits[h].hit;
690  double weight = whits[h].weight;
691 
692  assert( weight <= 1);
693 
694  if(IsNoise(*hit)) continue;
695 
696  std::vector<TrackIDE> tracks = HitToTrackIDE(*hit);
697 
698  // The set of trackIDs (or parents) that contributed to this hit
699  std::set<int> tally;
700 
701  for(size_t t = 0; t < tracks.size(); ++t){
702  const int trackID = tracks[t].trackID;
703  if(trackIDs.find(trackID) != trackIDs.end()){
704  tally.insert(trackID);
705  }
706  if(parents){
707  std::map<int, int>::iterator parentsIt = parents->find(trackID);
708  if(parentsIt != parents->end()){
709  tally.insert(parentsIt->second);
710  }
711  }
712  } // end for t
713 
714 
715  if (energyPur){
716  for(size_t t = 0; t < tracks.size(); ++t){
717  totalEInCell += tracks[t].energy * weight;
718  if(trackIDs.find(tracks[t].trackID) != trackIDs.end()){
719  desiredE += tracks[t].energy;
720  desiredEInCell += tracks[t].energy;
721  }
722  }
723  }//End of loop for counting energy when doing purity based on FLS hits
724 
725  // Purity in a cell can sometimes be greater than one if the weight of
726  // a hit is less than the fraction of energy from desired track IDs
727  // in that cell.
728  if( desiredEInCell >= totalEInCell )
729  totalE += desiredEInCell;
730  else
731  totalE += totalEInCell;
732 
733  if(purMap){
734  for(std::set<int>::iterator it = tally.begin(); it != tally.end(); ++it){
735  ++(*purMap)[*it];
736  }
737  }
738 
739  if(!tally.empty()) ++desired;
740 
741 
742  }// end loop over whits
743 
744  double purity = 0.;
745 
746  if((total > 0) && !energyPur) purity = desired/total;
747  if((totalE > 0) && energyPur) purity = desiredE/totalE;
748 
749  assert(purity >= 0);
750  assert(purity <= 1);
751 
752  if(total > 0 && purMap){
753  for(std::map<int, double>::iterator it = purMap->begin(); it != purMap->end(); ++it){
754  it->second /= total;
755  // This is a purity, it should be between 0-100%
756  assert(it->second >= 0);
757  assert(it->second <= 1);
758  }
759  }
760 
761  return purity;
762 
763  }
764 
765 
766  //----------------------------------------------------------------------
767 
768 
769  double BackTracker::
770  HitCollectionPurity(const std::set<int>& trackIDs,
771  const std::vector<const rb::CellHit*>& hits,
772  std::map<int, double>* purMap,
773  std::map<int, int>* parents,
774  bool energyPur) const
775  {
776  // convert all cell hits to weighted hits with weight 1.
777  std::vector<RawWeightedHit> wHits;
778  for(const rb::CellHit* h: hits) wHits.emplace_back(h, 1);
779 
780  return HitCollectionPurity(trackIDs, wHits, purMap, parents, energyPur);
781 
782  }
783 
784  //----------------------------------------------------------------------
785  double BackTracker::
786  HitCollectionEfficiency(const std::set<int>& trackIDs,
787  const std::vector<const rb::CellHit*>& hits,
788  const std::vector<const rb::CellHit*>& allhits,
789  const geo::View_t& view,
790  std::map<int, double>* effMap,
791  bool energyEff,
792  double* desiredEnergy,
793  double* totalEnergy,
794  int* desiredHits,
795  int* totHits) const
796  {
797 
798  // convert cell hits to weighted hits with weight 1
799  std::vector<RawWeightedHit> wHits;
800  for(const rb::CellHit* h: hits) wHits.emplace_back(h, 1);
801 
802  return HitCollectionEfficiency(trackIDs, wHits, allhits, view,
803  effMap, energyEff, desiredEnergy,
804  totalEnergy, desiredHits, totHits);
805  }
806  //----------------------------------------------------------------------
807  double BackTracker::
808  HitCollectionEfficiency(const std::set<int>& trackIDs,
809  const std::vector<rb::WeightedHit>& whits,
810  const std::vector<const rb::CellHit*>& allhits,
811  const geo::View_t& view,
812  std::map<int, double>* effMap,
813  bool energyEff,
814  double* desiredEnergy,
815  double* totalEnergy,
816  int* desiredHits,
817  int* totHits) const
818  {
819  const std::vector<RawWeightedHit> rwh(whits.begin(), whits.end());
820  return HitCollectionEfficiency(trackIDs, rwh, allhits, view,
821  effMap, energyEff, desiredEnergy, totalEnergy, desiredHits, totHits);
822  }
823 
824  //----------------------------------------------------------------------
825  double BackTracker::
826  HitCollectionEfficiency(const std::set<int>& trackIDs,
827  const std::vector<RawWeightedHit>& whits,
828  const std::vector<const rb::CellHit*>& allhits,
829  const geo::View_t& view,
830  std::map<int, double>* effMap,
831  bool energyEff,
832  double* desiredEnergy,
833  double* totalEnergy,
834  int* desiredHits,
835  int* totHits) const
836 
837  {
838 
839  // allhits should be a superset of hits
840  assert(allhits.size() >= whits.size());
841 
842  // to check it's a superset more thoroughly, optional
843  if(fParams.DoHitsCheck()){
844  for(unsigned int n = 0; n < whits.size(); ++n){
845  // TODO: should be checking actual equality not the pointers
846  // assert(std::find(allhits.begin(), allhits.end(), whits[n].hit) != allhits.end());
847  }
848  }
849  std::vector<const rb::CellHit*> hits_copy;
850  for(int i = 0; i < (int)whits.size(); i++)
851  hits_copy.push_back(whits[i].hit);
852 
853  // Check for duplicates
854  std::sort(hits_copy.begin(), hits_copy.end());
855  for(unsigned int n = 0; n+1 < hits_copy.size(); ++n)
856  assert(hits_copy[n] < hits_copy[n+1]);
857 
858  // get the list of EveIDs that correspond to the hits in this collection
859  // and the energy associated with the desired trackID
860 
861  double desired = 0.;
862  double total = 0.;
863 
864  double desiredE = 0.;
865  double totalE = 0.;
866 
867  if(effMap) effMap->clear();
868  std::map<int, double> totMap;
869 
870  for(size_t h = 0; h < whits.size(); ++h){
871 
872  double totalEInCell = 0.;
873  double desiredEInCell = 0.;
874 
875  const rb::CellHit* hit = whits[h].hit;
876  double weight = whits[h].weight;
877  assert( weight <= 1);
878 
879  if(IsNoise(*hit)) continue;
880 
881  bool hitCounted = false;
882 
883  std::vector<TrackIDE> tracks = HitToTrackIDE(*hit);
884 
885  assert(hit->View() == view || view == geo::kXorY);
886 
887  // don't double count if this hit has more than one of the
888  // desired track IDs associated with it
889  for(size_t t = 0; t < tracks.size(); ++t){
890  totalEInCell += tracks[t].energy;
891  if(trackIDs.find(tracks[t].trackID) != trackIDs.end() ){
892  desiredEInCell += tracks[t].energy;
893  if(tracks[t].energyFrac >= fParams.MinContribFrac()){
894  if(!hitCounted){
895  desired += 1.;
896  hitCounted = true;
897  }
898  if(effMap) ++(*effMap)[tracks[t].trackID];
899  }
900  }
901  }// end loop over track IDEs
902 
903  if( (weight*totalEInCell) >= desiredEInCell)
904  desiredE += desiredEInCell;
905  else
906  desiredE += (weight*totalEInCell);
907 
908  }// end loop over hits
909 
910  // now figure out how many hits in the whole collection are associated with this id
911  for(size_t h = 0; h < allhits.size(); ++h){
912  const rb::CellHit* hit = allhits[h];
913  if(IsNoise(*hit)) continue;
914 
915  std::vector<TrackIDE> tracks = HitToTrackIDE(*hit);
916 
917  bool hitCounted = false;
918 
919  // check that we are looking at the appropriate view here
920  // in the case of 3D objects we take all hits
921  if(hit->View() != view && view != geo::kXorY) continue;
922 
923  for(size_t t = 0; t < tracks.size(); ++t){
924  // don't worry about hits where the energy fraction for the chosen
925  // trackID is < MinContribFrac
926  // also don't double count if this hit has more than one of the
927  // desired track IDs associated with it
928  if(trackIDs.find(tracks[t].trackID) != trackIDs.end()){
929  totalE += tracks[t].energy;
930  if(tracks[t].energyFrac >= fParams.MinContribFrac()){
931  if(!hitCounted){
932  total += 1.;
933  hitCounted = true;
934  }
935  if(effMap) ++totMap[tracks[t].trackID];
936  }// end if energy fraction is above threshold
937  }// end if the current track ID is in the list of ones to check
938  }// end loop over tracks
939  }// end loop over all hits
940 
941  double efficiency = 0.;
942  if((total > 0.) && !energyEff) efficiency = desired/total;
943  if((totalE > 0.) && energyEff) efficiency = desiredE/totalE;
944 
945  if((total > 0.) && desiredHits ) *desiredHits = desired;
946  if((total > 0.) && totHits ) *totHits = total;
947  if((totalE > 0.) && desiredEnergy) *desiredEnergy = desiredE;
948  if((totalE > 0.) && totalEnergy ) *totalEnergy = totalE;
949 
950  if(effMap){
951  for(std::map<int, double>::iterator it = effMap->begin(); it != effMap->end(); ++it){
952  if(totMap[it->first] > 0)
953  it->second /= totMap[it->first];
954  // This is an efficiency. It should be between 0-100%
955  assert(it->second >= 0);
956  assert(it->second <= 1);
957  }
958  }
959 
960  return efficiency;
961  }
962 
963 
964  //----------------------------------------------------------------------
966  art::PtrVector<rb::CellHit> const& hits) const
967  {
968  int nxcells = 0;
969  int nycells = 0;
970 
971  for(size_t h = 0; h < hits.size(); ++h){
972 
973  art::Ptr<rb::CellHit> hit = hits[h];
974  std::vector<TrackIDE> tracks = HitToTrackIDE(hit);
975  for(size_t t = 0; t < tracks.size(); ++t){
976  if(tracks[t].trackID == trackID){
977 
978  if (hit->View() == geo::kX){
979  nxcells++;
980  }// end of loop over x hits
981 
982  if (hit->View() == geo::kY){
983  nycells++;
984  }// end of loop over y hits
985 
986  }// end of loop over hits from the muon
987  }// end of loop over tracks
988  }// end loop over hits
989 
990 
991  if ((nxcells >= 3) && (nycells >= 3)) return 1;
992  else return 0;
993  }
994 
995 
996  //----------------------------------------------------------------------
998  return a.efficiency>b.efficiency;
999  }
1000 
1001  //----------------------------------------------------------------------
1003  return a.purity>b.purity;
1004  }
1005  //----------------------------------------------------------------------
1006  bool BackTracker::sort_idxeff(const std::pair<int,NeutrinoEffPur>& a, const std::pair<int,NeutrinoEffPur>& b){
1007  return a.second.efficiency>b.second.efficiency;
1008  }
1009  //----------------------------------------------------------------------
1010  bool BackTracker::sort_idxpur(const std::pair<int,NeutrinoEffPur>& a, const std::pair<int,NeutrinoEffPur>& b){
1011  return a.second.purity>b.second.purity;
1012  }
1013 
1014  //----------------------------------------------------------------------
1015  std::vector<int> BackTracker::
1016  SliceToNeutrinoIndex(const std::vector<const rb::CellHit*>& sliceHits,
1017  const std::vector<const rb::CellHit*>& allHits,
1018  bool sortPur) const
1019  {
1020  std::vector < std::pair<unsigned int,NeutrinoEffPur> > sliceNeutrinoInteractions;
1021 
1022  for (unsigned int s = 0; s < fMCTruthList.size(); s++){
1023  const art::Ptr<simb::MCTruth> neutrinoInteraction = fMCTruthList[s];
1024  if(!neutrinoInteraction->NeutrinoSet())
1025  continue;
1026 
1027  std::vector<const sim::Particle*> particleList = MCTruthToParticles(neutrinoInteraction);
1028  std::set<int> neutrinoTrackIDs;
1029  for (unsigned int i = 0; i < particleList.size(); ++i){
1030  neutrinoTrackIDs.insert(particleList[i]->TrackId());
1031  }
1032 
1033  double energySlice = -1.0;
1034  double energyTotal = -1.0;
1035  int nSliceHits = 0;
1036  int nTotalHits = 0;
1037 
1038  double sliceEff = HitCollectionEfficiency(neutrinoTrackIDs, sliceHits, allHits, geo::kXorY, 0, true, &energySlice, &energyTotal, &nSliceHits, &nTotalHits);
1039  double slicePur = HitCollectionPurity(neutrinoTrackIDs, sliceHits, 0, 0, true);
1040 
1041  if (sliceEff > 0){
1042  NeutrinoEffPur sliceNeutrinoInteraction = {neutrinoInteraction, sliceEff, slicePur, s, energySlice, energyTotal, nSliceHits, nTotalHits};
1043  sliceNeutrinoInteractions.push_back( std::make_pair(s,sliceNeutrinoInteraction) );
1044  }
1045  }//End of loop over neutrino interactions
1046 
1047  if (!sortPur) std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_idxeff);
1048  else std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_idxpur);
1049 
1050  std::vector<int> nuIndex;
1051 
1052  for (unsigned int i = 0; i<sliceNeutrinoInteractions.size(); ++i ){
1053  nuIndex.push_back( int(sliceNeutrinoInteractions[i].first) );
1054  }
1055 
1056  return nuIndex;
1057  }
1058 
1059 
1060  //----------------------------------------------------------------------
1061  std::vector<NeutrinoEffPur> BackTracker::
1062  SliceToNeutrinoInteractions(const std::vector<const rb::CellHit*>& sliceHits,
1063  const std::vector<const rb::CellHit*>& allHits,
1064  bool sortPur) const
1065  {
1066  std::vector <NeutrinoEffPur> sliceNeutrinoInteractions;
1067 
1068  for (unsigned int s = 0; s < fMCTruthList.size(); s++){
1069  const art::Ptr<simb::MCTruth> neutrinoInteraction = fMCTruthList[s];
1070  if(!neutrinoInteraction->NeutrinoSet())
1071  continue;
1072 
1073  std::vector<const sim::Particle*> particleList = MCTruthToParticles(neutrinoInteraction);
1074  std::set<int> neutrinoTrackIDs;
1075  for (unsigned int i = 0; i < particleList.size(); ++i){
1076  neutrinoTrackIDs.insert(particleList[i]->TrackId());
1077  }
1078 
1079  double energySlice = -1.0;
1080  double energyTotal = -1.0;
1081  int nSliceHits = 0;
1082  int nTotalHits = 0;
1083 
1084  double sliceEff = HitCollectionEfficiency(neutrinoTrackIDs, sliceHits, allHits, geo::kXorY, 0, true, &energySlice, &energyTotal, &nSliceHits, &nTotalHits);
1085  double slicePur = HitCollectionPurity(neutrinoTrackIDs, sliceHits, 0, 0, true);
1086 
1087  if (sliceEff > 0){
1088  NeutrinoEffPur sliceNeutrinoInteraction = {neutrinoInteraction, sliceEff, slicePur, s, energySlice, energyTotal, nSliceHits, nTotalHits};
1089  sliceNeutrinoInteractions.push_back(sliceNeutrinoInteraction);
1090  }
1091  }//End of loop over neutrino interactions
1092 
1093  if (!sortPur) std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_eff);
1094  else std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_pur);
1095 
1096  return sliceNeutrinoInteractions;
1097  }
1098 
1099  //----------------------------------------------------------------------
1100  std::vector<NeutrinoEffPur> BackTracker::
1102  const std::vector<const rb::CellHit*>& allHits,
1103  bool sortPur) const
1104  {
1105  return SliceToNeutrinoInteractions(PtrVecToVecRawPtr(sliceCluster->AllCells()),
1106  allHits, sortPur);
1107  }
1108 
1109  //......................................................................
1110  std::vector<int> BackTracker::
1111  SliceToOrderedNuIds(const std::vector<cheat::NeutrinoWithIndex>& nusWithIdx,
1112  const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable,
1113  std::function<double(const cheat::NeutrinoEffPur&)> slMetric,
1114  std::function<double(const cheat::NeutrinoEffPur&)> nuMetric ) const
1115  {
1116 
1117  for( auto slc: slTruthTable )
1118  assert( slc.size() == nusWithIdx.size() && "Badly shaped Slice/truth table!" );
1119 
1120  // Fill with one Nu Id per slice
1121  std::vector<int> nuIds( slTruthTable.size(), -1 );
1122  if( nusWithIdx.empty() ) return nuIds;
1123 
1124  // Find the best slice for each nu
1125  std::vector<int> bestSlices( nusWithIdx.size() );
1126  for( size_t nuIdx = 0; nuIdx < nusWithIdx.size(); ++nuIdx ){
1127  float maxValue = 0;
1128  int bestSlice = -1;
1129 
1130  for( size_t sliceIdx = 0; sliceIdx < slTruthTable.size(); ++sliceIdx ){
1131  const cheat::NeutrinoEffPur& effpur = slTruthTable[sliceIdx][nuIdx];
1132 
1133  if( !effpur.neutrinoInt.isNull() ){
1134  if( slMetric(effpur) > maxValue ){
1135  maxValue = slMetric(effpur);
1136  bestSlice = sliceIdx;
1137  }
1138  }
1139 
1140  } //sliceIdx
1141 
1142  bestSlices[nuIdx] = bestSlice;
1143 
1144  } // nuIdx
1145 
1146  // Find the best nu for each slice.
1147  for( size_t sliceIdx = 0; sliceIdx < slTruthTable.size(); ++sliceIdx ){
1148  float maxValue = 0;
1149  int bestNu = -1;
1150 
1151  for( size_t nuIdx = 0; nuIdx < slTruthTable[sliceIdx].size(); ++nuIdx ){
1152  const cheat::NeutrinoEffPur& effpur = slTruthTable[sliceIdx][nuIdx];
1153 
1154  if( !effpur.neutrinoInt.isNull() ){
1155  if( nuMetric(effpur) > maxValue ){
1156  maxValue = nuMetric(effpur);
1157  bestNu = nuIdx;
1158  }
1159  }
1160  } //nuIdx
1161 
1162  if( bestNu >= 0 && bestSlices[bestNu] == (int)sliceIdx )
1163  nuIds[sliceIdx] = bestNu;
1164 
1165  }//sliceIdx
1166 
1167  return nuIds;
1168  }
1169 
1170  //----------------------------------------------------------------------
1171  std::vector<int> BackTracker::
1172  SliceToOrderedNuIdsByEff(const std::vector<cheat::NeutrinoWithIndex>& nusWithIdx,
1173  const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const
1174  {
1175 
1176  return SliceToOrderedNuIds(nusWithIdx, slTruthTable, EffMetric, EffMetric);
1177 
1178  }
1179 
1180  //----------------------------------------------------------------------
1181  std::vector<int> BackTracker::
1182  SliceToOrderedNuIdsByPur(const std::vector<cheat::NeutrinoWithIndex>& nusWithIdx,
1183  const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const
1184  {
1185 
1186  return SliceToOrderedNuIds(nusWithIdx, slTruthTable, PurMetric, PurMetric);
1187 
1188  }
1189 
1190  //----------------------------------------------------------------------
1191  std::vector<int> BackTracker::
1192  SliceToOrderedNuIdsByEffPur(const std::vector<cheat::NeutrinoWithIndex>& nusWithIdx,
1193  const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const
1194  {
1195 
1196  return SliceToOrderedNuIds(nusWithIdx, slTruthTable, EffPurMetric, EffPurMetric);
1197 
1198  }
1199 
1200  //----------------------------------------------------------------------
1201  std::vector<int> BackTracker::
1202  SliceToOrderedNuIdsByEnergy(const std::vector<cheat::NeutrinoWithIndex>& nusWithIdx,
1203  const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const
1204  {
1205 
1206  return SliceToOrderedNuIds(nusWithIdx, slTruthTable, EnergyMetric, EnergyMetric);
1207 
1208  }
1209 
1210  //----------------------------------------------------------------------
1211  std::vector<NeutrinoEffPur> BackTracker::
1212  SliceToMCTruth(const std::vector<const rb::CellHit*>& sliceHits,
1213  const std::vector<const rb::CellHit*>& allHits,
1214  bool sortPur) const
1215  {
1216  std::vector <NeutrinoEffPur> sliceNeutrinoInteractions;
1217 
1218  for (unsigned int s = 0; s < fMCTruthList.size(); s++){
1219  const art::Ptr<simb::MCTruth> neutrinoInteraction = fMCTruthList[s];
1220 
1221  std::vector<const sim::Particle*> particleList = MCTruthToParticles(neutrinoInteraction);
1222  std::set<int> neutrinoTrackIDs;
1223  for (unsigned int i = 0; i < particleList.size(); ++i){
1224  neutrinoTrackIDs.insert(particleList[i]->TrackId());
1225  }
1226 
1227  double energySlice = -1.0;
1228  double energyTotal = -1.0;
1229  int nSliceHits = 0;
1230  int nTotalHits = 0;
1231 
1232  double sliceEff = HitCollectionEfficiency(neutrinoTrackIDs, sliceHits, allHits, geo::kXorY, 0, true, &energySlice, &energyTotal, &nSliceHits, &nTotalHits);
1233  double slicePur = HitCollectionPurity(neutrinoTrackIDs, sliceHits, 0, 0, true);
1234 
1235  if (sliceEff > 0){
1236  NeutrinoEffPur sliceNeutrinoInteraction = {neutrinoInteraction, sliceEff, slicePur, s, energySlice, energyTotal, nSliceHits, nTotalHits};
1237  sliceNeutrinoInteractions.push_back(sliceNeutrinoInteraction);
1238  }
1239  }//End of loop over neutrino interactions
1240 
1241  if (!sortPur) std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_eff);
1242  else std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_pur);
1243 
1244  return sliceNeutrinoInteractions;
1245  }
1246 
1247  //----------------------------------------------------------------------
1248  std::vector<NeutrinoEffPur> BackTracker::
1250  const std::vector<const rb::CellHit*>& allHits,
1251  bool sortPur) const
1252  {
1253  return SliceToMCTruth(PtrVecToVecRawPtr(sliceCluster->AllCells()),
1254  allHits, sortPur);
1255  }
1256 
1257  //----------------------------------------------------------------------
1258  std::vector< std::vector<cheat::NeutrinoEffPur> > BackTracker::
1259  SlicesToMCTruthsTable(const std::vector<const rb::Cluster*>& sliceList) const
1260  {
1261 
1262  unsigned int numSlices = sliceList.size();
1263  std::vector<NeutrinoWithIndex> truthList = allMCTruth();
1264  unsigned int numTruths = truthList.size();
1265 
1266  //This is what will eventually be returned
1267  std::vector< std::vector<cheat::NeutrinoEffPur> > table(numSlices,std::vector<cheat::NeutrinoEffPur>(numTruths));
1268 
1269  //Loop to fill NeutrinoEffPur's with proper truth information
1270  for (unsigned int s = 0; s < numSlices; s++){
1271  for (unsigned int t = 0; t < numTruths; t++){
1272  cheat::NeutrinoEffPur truthFill = {truthList[t].neutrinoInt,0.0,0.0,truthList[t].truthColIndex,0.0,0.0,0,0};
1273  table[s][t] = truthFill;
1274  }
1275  }//End of defining basic container with truths filled in
1276 
1277  for (unsigned int s = 0; s < numSlices; s++){
1278  //Vector to hold the total energy contributions to slice from each MC Truth
1279  std::vector<double> sliceEnergies(numTruths,0.0);
1280  //Vector to hold the total number of hits in the slice from each MC Truth
1281  std::vector<int> sliceHits(numTruths,0);
1282 
1283  for (unsigned int h = 0; h < sliceList[s]->NCell(); h++){
1284  //Don't worry about noise hits - they can't count towards the energy an interaction deposited
1285  if (IsNoise(sliceList[s]->Cell(h))) continue;
1286  //Vector to hold the total energy contributions to hit from each MC Truth
1287  std::vector<double> hitEnergies(numTruths,0.0);
1288 
1289  const std::vector<TrackIDE> trackIDs = HitToTrackIDE(sliceList[s]->Cell(h));
1290  for (unsigned int id = 0; id < trackIDs.size(); id++){
1291  int trackID = trackIDs[id].trackID;
1292  //Gets the MCTruth index for the particle trackID
1293  int truthID = fTrackIDToMCTruthIndex.find(abs(trackID))->second;
1294  hitEnergies[truthID] += trackIDs[id].energy;
1295  }//End of loop over track ids for the hit
1296 
1297  //Add hit values into the slice vectors
1298  for (unsigned int q = 0; q < numTruths; q++){
1299  if (hitEnergies[q] > 0.0){
1300  sliceEnergies[q] += hitEnergies[q];
1301  sliceHits[q] += 1;
1302  }//Only do for MCTruths with nonzero contributions to hit
1303  }//End of loop adding hit values to slice vectors
1304 
1305  }//End of loop over hits in slice
1306 
1307  double totalSliceE = std::accumulate(sliceEnergies.begin(),sliceEnergies.end(),0.0);
1308  for (unsigned int t = 0; t < numTruths; t++){
1309 
1310  double purity = 0.0;
1311  if (totalSliceE > 0){
1312  purity = (sliceEnergies[t])/totalSliceE;
1313  }
1314 
1315  if ((purity < 0) || (purity > 1)) {
1316  throw cet::exception("BackTrackerIllegalPurity")
1317  << "The purity value "<<purity<<" is not allowed. \n";
1318  }
1319 
1320  table[s][t].purity = purity;
1321  table[s][t].energySlice = sliceEnergies[t];
1322  table[s][t].nSliceHits = sliceHits[t];
1323 
1324  }//End of loop over MCTruths for the slice
1325  }//End of loop over slices
1326 
1327  for (unsigned int t = 0; t < numTruths; t++){
1328 
1329  double totalMCTruthE = 0.0;
1330  int totalMCTruthHits = 0;
1331 
1332  for (unsigned int s = 0; s < numSlices; s++){
1333  totalMCTruthE += table[s][t].energySlice;
1334  totalMCTruthHits += table[s][t].nSliceHits;
1335  }//Loop over slices to find summed values
1336 
1337  for (unsigned int s = 0; s < numSlices; s++){
1338 
1339  double efficiency = 0.0;
1340  if (totalMCTruthE > 0){
1341  efficiency = (table[s][t].energySlice)/totalMCTruthE;
1342  }
1343 
1344  if ((efficiency < 0) || (efficiency > 1)) {
1345  throw cet::exception("BackTrackerIllegalEfficiency")
1346  << "The efficiency value "<<efficiency<<" is not allowed. \n";
1347  }
1348 
1349  table[s][t].efficiency = efficiency;
1350  table[s][t].energyTotal = totalMCTruthE;
1351  table[s][t].nTotalHits = totalMCTruthHits;
1352  }//Loop over slices to insert final values
1353 
1354  }//End of final loop over MCTruths for summed values
1355 
1356  return table;
1357  }
1358 
1359  //----------------------------------------------------------------------
1360  std::vector<ParticleEffPur> BackTracker::
1361  TracksToParticles(const std::vector<const rb::Track*>& tracks,
1362  const std::set<int>& trackIDs,
1363  const std::set<int>& allowedDaughterIDs,
1364  const std::vector<const rb::CellHit*>& allHits) const
1365  {
1366  std::vector<ParticleEffPur> matches;
1367 
1368  std::set<int> trkIDs = trackIDs;
1369  std::set<int> allowedDaughters = allowedDaughterIDs;
1370 
1371  std::map<int, int> parents;
1372  for(std::set<int>::iterator it = trkIDs.begin(); it != trkIDs.end(); ++it){
1373  const sim::Particle* mother = TrackIDToParticle(*it);
1374  if(!mother) continue;
1375  for(int i = 0; i < mother->NumberDaughters(); ++i){
1376  const int id = mother->Daughter(i);
1378  if(allowedDaughters.count(TrackIDToParticle(id)->PdgCode())){
1379  parents[id] = *it;
1380  }
1381  } // end for i
1382  } // end for it
1383 
1384  // loop over all tracks
1385  for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack){
1386  // get the purity of the track for all track IDs
1387  std::map<int, double> purs;
1388  HitCollectionPurity(trkIDs,tracks[itrack]->AllCells(),&purs,&parents);
1389  // find the most pure trackID for this track
1390  int bestid = -1;
1391  double bestpur = 0.0;
1392  for(std::set<int>::iterator it = trkIDs.begin(); it != trkIDs.end(); ++it){
1393  const double p = purs[*it];
1394  if(p > bestpur){
1395  bestpur = p;
1396  bestid = *it;
1397  }
1398  } // end of loop over track ids
1399 
1400  // update this tracks info for the matched true particle
1401  double recoPur = bestpur;
1402  int recotrkid = bestid;
1403  int recopdg = 0;
1404  double recoEff = 0.0;
1405  if(bestpur > 0.0){
1406  // found a real particle that matched the track, update pdg and efficiency to match this
1407  recopdg = TrackIDToParticle(recotrkid)->PdgCode();
1408  std::set<int> ideff;
1409  ideff.insert(recotrkid);
1410  recoEff = HitCollectionEfficiency(ideff,
1411  PtrVecToVecRawPtr(tracks[itrack]->AllCells()),
1412  allHits, geo::kXorY);
1413  // remove this track id from the set of track ids so no other tracks can try to match to it
1414  trkIDs.erase(trkIDs.find(recotrkid));
1415  // remove any instance where this track id is the parent particle in the parents map
1416  std::map<int, int>::iterator parentsIt = parents.begin();
1417  while(parentsIt != parents.end()){
1418  if(parentsIt->second == recotrkid){ parents.erase(parentsIt++); }
1419  else{ ++parentsIt; }
1420  }
1421  }
1422  // record this tracks match
1423  ParticleEffPur match = {recotrkid, recoEff, recoPur, recopdg};
1424  matches.push_back(match);
1425  } // end of loop over tracks for matching to truth
1426 
1427  return matches;
1428 
1429  }
1430 
1431  //----------------------------------------------------------------------
1432  std::vector<NeutrinoWithIndex> BackTracker::allNeutrinoInteractions() const
1433  {
1434  std::vector <NeutrinoWithIndex> allNeutrinoInteractions;
1435 
1436  for (unsigned int s = 0; s < fMCTruthList.size(); s++){
1437  const art::Ptr<simb::MCTruth> neutrinoInteraction = fMCTruthList[s];
1438  if(!neutrinoInteraction->NeutrinoSet())
1439  continue;
1440 
1441  NeutrinoWithIndex neutrino = {neutrinoInteraction, s};
1442  allNeutrinoInteractions.push_back(neutrino);
1443  }
1444 
1445  return allNeutrinoInteractions;
1446  }
1447 
1448  //----------------------------------------------------------------------
1449  std::vector<NeutrinoWithIndex> BackTracker::allMCTruth() const
1450  {
1451  std::vector <NeutrinoWithIndex> allTruth;
1452 
1453  for (unsigned int s = 0; s < fMCTruthList.size(); s++){
1454  const art::Ptr<simb::MCTruth> neutrinoInteraction = fMCTruthList[s];
1455  NeutrinoWithIndex neutrino = {neutrinoInteraction, s};
1456  allTruth.push_back(neutrino);
1457  }
1458 
1459  return allTruth;
1460  }
1461 
1462 
1463  //----------------------------------------------------------------------
1464  bool BackTracker::CloseInTime(const rb::CellHit& chit, const sim::PhotonSignal& phot) const
1465  {
1466  //The constants used to determine this value change depending on detector
1467  //and hit timing reconstruction. See docdb 9092 to see how the current
1468  //constants were determined (May 2013). To redo the study and see if new
1469  //constants are required, use MCCheckOut/CloseInTime_module.cc.
1470 
1471  // Empirically determined as an over-generous width of the timing reconstruction.
1472  double kWindowNs = 1000;
1473  // Empirically determined as the approximate mean bias of the timing reconstruction.
1474  // This is the amount the CellHit is expected to be earlier than the photons.
1475  double kOffsetNs = 350;
1476 
1477  if ((fDetID == novadaq::cnv::kFARDET)&&(chit.GoodTiming())){
1478  kWindowNs = 500;
1479  kOffsetNs = 10;
1480  }
1481 
1482  if ((fDetID == novadaq::cnv::kNEARDET)&&(chit.GoodTiming())){
1483  kWindowNs = 200;
1484  kOffsetNs = 10;
1485  }
1486 
1487  const double dt = chit.TNS()-phot.TimeMean();
1488 
1489  return fabs(dt+kOffsetNs) <= kWindowNs;
1490  }
1491 
1492  //----------------------------------------------------------------------
1494  std::map<int, int> & idToNPhot) const
1495  {
1496  auto it = fCellToPhotonSignal.find(geo::OfflineChan(hit.Plane(), hit.Cell()));
1497  // Nothing to accumulate from this cell
1498  if(it == fCellToPhotonSignal.end()) return;
1499 
1500  // Candidate photons, in that they're in the right cell.
1501  const std::vector<sim::PhotonSignal>& phots = it->second;
1502 
1503  // For every photon that's sufficiently in time with the hit,
1504  // accumulate its energy into the relevant trackID.
1505  for(unsigned int i = 0; i < phots.size(); ++i){
1506  const sim::PhotonSignal& phot = phots[i];
1507  if(CloseInTime(hit, phot)){
1508  idToNPhot[phot.TrackId()] += phot.NPhoton();
1509  }
1510  }
1511  }
1512  //----------------------------------------------------------------------
1513 
1514  std::vector<sim::FLSHit> BackTracker::ParticleToFLSHit(const int& trackID) const{
1515  std::unordered_map< int, std::list<sim::FLSHit> >::const_iterator it = fTrackIDToFLSHit.find(trackID);
1516  if( it == fTrackIDToFLSHit.end() ){
1517  std::vector<sim::FLSHit> empty;
1518  LOG_DEBUG("BackTracker") << "no collection of sim::FLSHit for "
1519  <<"Track ID "<<trackID
1520  << ". Returning empty vector";
1521 
1522  return empty;
1523  }
1524  else
1525  return std::vector<sim::FLSHit>(it->second.begin(), it->second.end());
1526  }
1527 
1528 
1529  //----------------------------------------------------------------------
1530  std::vector<sim::FLSHit> BackTracker::ParticleToFLSHit(const int trackID, const int pdg) const{
1531 
1532  const std::vector<sim::FLSHit>& flshits_with_same_track_id = ParticleToFLSHit(trackID);
1533 
1534  std::vector<sim::FLSHit> out_flshits;
1535 
1536  const int nflshits = flshits_with_same_track_id.size();
1537 
1538  /// Fill only FLSHits with a given PDG
1539  for(int i=0; i<nflshits; ++i){
1540  const sim::FLSHit& flshit= flshits_with_same_track_id[i];
1541  if(flshit.GetPDG() == pdg){
1542  out_flshits.push_back(flshit);
1543  }
1544  }
1545 
1546  return out_flshits;
1547  }
1548 
1549  //----------------------------------------------------------------------
1550  std::vector<rb::Cluster> BackTracker::ClusterByTruth(std::vector< art::Ptr< rb::CellHit > > const& hits)
1551  {
1552  std::vector<rb::Cluster> clusters(fMCTruthList.size()+1);
1553 
1554  for(unsigned int hiti = 0; hiti < hits.size(); ++hiti) {
1555  art::Ptr<rb::CellHit> hit = hits[hiti];
1556 
1557  // If the hit is noise, put it into the last cluster.
1558  if(this->IsNoise(hit)) {
1559  clusters[fMCTruthList.size()].Add(hit);
1560  continue;
1561  }
1562 
1563  const sim::Particle *part = this->HitToParticle(hit);
1564  unsigned int idx = fTrackIDToMCTruthIndex[part->TrackId()];
1565  clusters[idx].Add(hit);
1566 
1567  } // end loop over hits
1568 
1569  // Label the noise slice.
1570  clusters[fMCTruthList.size()].SetNoise(true);
1571 
1572  return clusters;
1573  }
1574 
1575  //----------------------------------------------------------------------
1577  std::vector< art::Ptr< rb::CellHit > > const& hits)
1578  {
1579  rb::Cluster cluster;
1580 
1581  // figure out where in the MCTruth list mct belongs
1582  unsigned int MCTidx;
1583  for(MCTidx = 0; MCTidx < fMCTruthList.size(); ++MCTidx) {
1584  if(fMCTruthList[MCTidx] == mct) break;
1585  }
1586 
1587  // loop over hits and add to the cluster
1588  for(unsigned int hiti = 0; hiti < hits.size(); ++hiti) {
1589  art::Ptr<rb::CellHit> hit = hits[hiti];
1590 
1591  // If the hit is noise, skip it.
1592  if(this->IsNoise(hit)) continue;
1593 
1594  const sim::Particle *part = this->HitToParticle(hit);
1595  unsigned int idx = fTrackIDToMCTruthIndex[part->TrackId()];
1596  if(idx == MCTidx) cluster.Add(hit);
1597 
1598  } // end loop over hits
1599 
1600  return cluster;
1601  }
1602 
1603  //----------------------------------------------------------------------
1604 
1605  void BackTracker::HitToParticlesMap(const std::vector<const rb::Cluster*>& sliceList,
1606  const int& sliceIdx)
1607  {
1608  fCellToIDs.clear();
1609  fIdToEnergy.clear();
1610 
1611  unsigned int nsli = sliceList.size();
1612 
1613  for (unsigned int isli = 0; isli < nsli; isli++){
1614  unsigned int nh = sliceList[isli]->NCell();
1615  for (unsigned int ih = 0; ih < nh; ih++){
1616 
1617  art::Ptr<rb::CellHit> hit = sliceList[isli]->Cell(ih);
1618 
1619  if (IsNoise(hit)) continue;
1620  double chan = (hit->TDC()*1000000) + (hit->Plane()*1000) + hit->Cell();
1621 
1622  std::vector<cheat::TrackIDE> ides = HitToTrackIDE(hit);
1623  fCellToIDs[chan] = ides;
1624 
1625  for(auto ide : ides)
1626  fIdToEnergy[ide.trackID] += ide.energy;
1627  }// end of loop over slice hits
1628  }// end of loop over slices
1629  fSliceIDForHitParticleMap = sliceIdx;
1630 
1631  }// end of hits to particles table
1632 
1633  //----------------------------------------------------------------------
1634 
1635  cheat::ParticleEffPur BackTracker::ClusterToParticle(const std::vector<const rb::Cluster*>& sliceList,
1636  const std::vector<const rb::CellHit*>& hits,
1637  const int& sliceIdx)
1638  {
1639 
1640  if(sliceIdx != fSliceIDForHitParticleMap){
1641  HitToParticlesMap(sliceList, sliceIdx);
1642  }
1643 
1644  double totalE = 0.0;
1645  unsigned int nh = hits.size();
1646 
1647  // make a map of how much energy each track ID contributes to
1648  // this cluster. Also keep track of how much total true energy
1649  // there is in this cluster.
1650 
1651  std::map< int, double> idToEnergyInProng;
1652  for(unsigned int ih = 0; ih < nh; ih++){
1653  const rb::CellHit* hit = hits[ih];
1654  double chan = (hit->TDC()*1000000) + (hit->Plane()*1000) + hit->Cell();
1655  std::vector<cheat::TrackIDE> ides = fCellToIDs[chan];
1656  for(auto ide: ides){
1657  idToEnergyInProng[ide.trackID] += ide.energy;
1658  totalE += ide.energy;
1659  }// end loop over ides
1660  }// end loop over cluster hits
1661 
1662  // now find which trackID contributed the most energy to
1663  // the cluster
1664  int maxID = -1;
1665  double maxE = -1.0;
1666 
1667  cheat::ParticleEffPur puff = {-1,-1,-1,-1};
1668  if(idToEnergyInProng.size() == 0)
1669  return puff;
1670 
1671  for(auto idtoe: idToEnergyInProng){
1672  if(idtoe.second > maxE ){
1673  maxE = idtoe.second;
1674  maxID = idtoe.first;
1675  }
1676  }// end loop over trackID to energy map
1677 
1678  puff.trackID = maxID;
1679  puff.purity = maxE/totalE;
1680  puff.efficiency = maxE/fIdToEnergy[maxID];
1681  puff.pdg = TrackIDToParticle(maxID)->PdgCode();
1682 
1683  return puff;
1684  }// end of ClusterToParticle
1685 
1686  //----------------------------------------------------------------------
1687 
1688  float BackTracker::CalcEscapingEnergy( const sim::Particle &Par, bool UseMCPart, std::string G4Label ) {
1689  float TempVar = -5;
1690  return CalcEscapingEnergy( Par, TempVar, UseMCPart, G4Label );
1691  }// end of CalcEscapingEnergy
1692 
1693  //----------------------------------------------------------------------
1694 
1695  float BackTracker::CalcEscapingEnergy( const sim::Particle &Par, float& EnDep, bool UseMCPart, std::string G4Label ) {
1696  /// A function to calculate the amount of energy from a given MCParticle which is not contained in the detector.
1697  /// Will return the energy the particle had when it exited the detector.
1698  /// Should a particle stop in the detector then the "Escaping energy" is 0.
1699  /// It has the option of returning the energy that a particle deposited energy in the detector.
1700  /// This is calclated from EnEnt - EnEsc, this has the deficiency that depositions from daughters may be doubly counted,
1701  /// however this is what the function CalcDaughterEscapingEnergy is for.
1702  /// There is a flag "UseMCPart", this controls whether to use the saved MCParticles, or to use the more accurate GEANT4
1703  /// information, which takes into account all daughter particles.
1704  /// Currently it is defaulted to true, ie use the MCParticles, as the GEANT4 info is not usually in the event record yet.
1705 
1706  // Set to some default values.
1707  float EnEsc, EnEnt;
1708  EnEsc = EnEnt = EnDep = -5;
1709 
1710  // --- Do I want to use the information stored in MCParticle?
1711  if (UseMCPart) {
1712  // Make variables to call InterceptsDetector.
1713  int TrajEnter, TrajExit;
1714  TrajEnter = TrajExit = -1;
1715  // Call InterceptsDetector
1716  bool IntDet = InterceptsDetector( Par, TrajEnter, TrajExit );
1717 
1718  // If IntDet == false then the particle never entered the detector, so return 0's.
1719  if (!IntDet) {
1720  return EnEsc;
1721  } else {
1722  // Now work out the energy which the particle deposited...
1723  // Set EnEnt
1724  EnEnt = (float)Par.E( TrajEnter ) - Par.Mass();
1725  // If TrajExit isn't NTraj-1, the particle left the volume then some energy was lost!
1726  // If TrajExit == NTraj-1 then we can just leave EnEsc as 0.
1727  if ( TrajExit != (int)Par.NumberTrajectoryPoints()-1 ) {
1728  EnEsc = (float)Par.E( TrajExit ) - Par.Mass();
1729  }
1730  // The energy deposited is just the entering energy - the exiting energy.
1731  EnDep = EnEnt - EnEsc;
1732  }
1733  // now return
1734  return EnEsc;
1735  } else {
1736  std::map<int, int>::iterator TrEIt = fTrackIDToTrueEIndex.find( Par.TrackId() );
1737  if ( TrEIt != fTrackIDToTrueEIndex.end() ) {
1738  art::Ptr< sim::TrueEnergy > ThisEn = fTrueEnergyList[ TrEIt->second ];
1739  /*
1740  std::cout << "Looking at the sim::TrueEnergy from the data product it is, TrackID " << ThisEn->TrackID() << ", EntEn " << ThisEn->EnteringEnergy()
1741  << ", EscEn " << ThisEn->EscapingEnergy() << ", TotEDep " << ThisEn->TotalDepEnergy() << ", TotEscEn " << ThisEn->TotalEscEnergy()
1742  << std::endl;
1743  */
1744  EnDep = ThisEn->TotalDepEnergy();
1745  return ThisEn->TotalEscEnergy();
1746  } else {
1747  mf::LogWarning("BackTracker") << "Didn't find TrackId " << Par.TrackId() << " in my map of sim::TrueEnergy's, returning -5.";
1748  return -5;
1749  }
1750  }
1751  }// end of CalcEscapingEnergy
1752 
1753  //----------------------------------------------------------------------
1754 
1756  float DauEsc = 0.;
1757  int NumDaught = Par.NumberDaughters();
1758  for (int dau=0; dau<NumDaught; ++dau) {
1759  const sim::Particle* ThisDau = TrackIDToParticle( Par.Daughter( dau ) );
1760  float ThisEsc = CalcEscapingEnergy( (*ThisDau), true );
1761  DauEsc += std::max((float)0., ThisEsc);
1762  /*
1763  std::cout << "\tDaughter " << dau << " of " << NumDaught << " was a " << ThisDau->PdgCode() << ", TrackId " << ThisDau->TrackId()
1764  << ", it had " << ThisDau->NumberDaughters() << " daughters itself. "
1765  << " EscEn for this particle = " << ThisEsc << " ==> DauEsc = " << DauEsc
1766  << std::endl;
1767  */
1768  if ( ThisDau->NumberDaughters() ) {
1769  //std::cout << "\t\t****** This particle has daughters, so recall the function... ******" << std::endl;
1770  float DaughtEsc = CalcDaughterEscapingEnergy( (*ThisDau) );
1771  DauEsc += DaughtEsc;
1772  //std::cout << "\t\t****** The sum of all EEscs was " << DaughtEsc << " ==> DauEsc is " << DauEsc << std::endl;
1773  } // If this particle has daughters itself
1774  } // Loop through the daughters of my particle.
1775  //std::cout << "\n\tAt the end of all of that, the sum of all escaping daughters was " << DauEsc << std::endl;
1776  return DauEsc;
1777  } // end of CalcDaughterEscapingEnergy
1778 
1779  //----------------------------------------------------------------------
1780 
1781  float BackTracker::CalcTotalEscapingEnergy( const sim::Particle &Par, bool UseMCPart, std::string G4Label ) {
1782  float TempVar = 0.;
1783  return CalcTotalEscapingEnergy( Par, TempVar, UseMCPart, G4Label );
1784  }
1785 
1786  //----------------------------------------------------------------------
1787 
1788  float BackTracker::CalcTotalEscapingEnergy( const sim::Particle &Par, float& TotEnDep, bool UseMCPart, std::string G4Label ) {
1789  if ( UseMCPart ) {
1790  float ParEDep = 0;
1791  float ParEEsc = CalcEscapingEnergy( Par, ParEDep, UseMCPart, G4Label );
1792  float DauEEsc = CalcDaughterEscapingEnergy( Par );
1793 
1794  float TotEnEsc = std::max((float)0., ParEEsc) + std::max((float)0., DauEEsc);
1795  TotEnDep = ParEDep - DauEEsc;
1796  // --- If negative, then had particle had very little energy --> may need to consider rest mass energy.
1797  if ( TotEnDep < 0 ) {
1798  //std::cout << "\t\t ** TotEnDep is negative " << TotEnDep << ", initially set to 0, but will check Part type. " << std::endl;
1799  // --- Default TotEnDep to 0.
1800  TotEnDep = 0;
1801  // --- If this particle decayed in the detector (ParEdep != 0) then it's likely that the particles it produced left the detector,
1802  // and the sum of their energies will be more than the initial KE of the particle.
1803  // For now I will modify this for some choice particles, not sure how best to do this...
1804  if ( ParEDep != 0 &&
1805  ( abs(Par.PdgCode()) == 11 || // Electron
1806  abs(Par.PdgCode()) == 13 || // Muon
1807  abs(Par.PdgCode()) == 111 || // Pi0
1808  abs(Par.PdgCode()) == 211 || // Pion
1809  abs(Par.PdgCode()) == 321 ) // Kaon
1810  ) {
1811  /*
1812  std::cout << "\t\t ** This particle likely decayed, so adding RestMass of " << Par.Mass()
1813  << " --- TotEnDep = " << ParEDep << " - " << DauEEsc << " + " << Par.Mass() << " = " << ParEDep - DauEEsc + Par.Mass()
1814  << std::endl;
1815  */
1816  TotEnDep = std::max( 0., ParEDep - DauEEsc + Par.Mass() );
1817  }
1818  }
1819  /*
1820  std::cout << "\t***The end of CalcTotalEscapingEnergy -- ParEEsc " << ParEEsc << ", ParEDep " << ParEDep << ", DauEEsc " << DauEEsc << ".\n"
1821  << "\t TotEnEsc = " << ParEEsc << " + " << DauEEsc << " = " << TotEnEsc << ".\n"
1822  << "\t TotEnDep = " << ParEDep << " - " << DauEEsc << " = " << TotEnDep << "."
1823  << std::endl;
1824  */
1825  return TotEnEsc;
1826  } else {
1827  std::map<int, int>::iterator TrEIt = fTrackIDToTrueEIndex.find( Par.TrackId() );
1828  if ( TrEIt != fTrackIDToTrueEIndex.end() ) {
1829  art::Ptr< sim::TrueEnergy > ThisEn = fTrueEnergyList[ TrEIt->second ];
1830  /*
1831  std::cout << "Looking at the sim::TrueEnergy from the data product it is, TrackID " << ThisEn->TrackID() << ", EntEn " << ThisEn->EnteringEnergy()
1832  << ", EscEn " << ThisEn->EscapingEnergy() << ", TotEDep " << ThisEn->TotalDepEnergy() << ", TotEscEn " << ThisEn->TotalEscEnergy()
1833  << std::endl;
1834  */
1835  TotEnDep = ThisEn->TotalDepEnergy();
1836  return ThisEn->TotalEscEnergy();
1837  } else {
1838  mf::LogWarning("BackTracker") << "Didn't find TrackId " << Par.TrackId() << " in my map of sim::TrueEnergy's, returning -5.";
1839  return -5;
1840  }
1841  }
1842  }
1843 
1844  //----------------------------------------------------------------------
1845 
1846  bool BackTracker::InterceptsDetector( const sim::Particle &Par, int &TrajEnt, int &TrajEx ) {
1847  /// A function which returns the trajectory point at which the particle first enters and leave the detector volume.
1848  /// The trajectory point is returned so that the user can determine the energies etc which the particle had at this
1849  /// point.
1850  /// Should the particle have been generated in the detector then TrajEnt == 0 will be returned.
1851  /// Should the particle stop in the detector then TrajEx == NumberTrajectoryPoints()-1 will be returned.
1852  /// Should the particle never enter the detector then holder values TrajEnter == TrajEx == -5 will be returned.
1853 
1854  // Initially set some holder values.
1855  TrajEnt = TrajEx = -5;
1856  // Make a handle to geometry service
1858  // Get the Number of trajectory points.
1859  unsigned int NumTraj = Par.NumberTrajectoryPoints();
1860  bool EnteredVol = false;
1861  bool PrevHit = false;
1862  for (unsigned int trj=0; trj<NumTraj; ++trj) {
1863  // Is the hit in the detector?
1864  bool InDet = geom->isInsideDetectorBigBox( Par.Vx(trj), Par.Vy(trj), Par.Vz(trj) );
1865  if ( InDet ) {
1866  // If this is the first hit in the volume.
1867  if ( !EnteredVol ) {
1868  TrajEnt = trj;
1869  EnteredVol = true;
1870  }
1871  // If this is the last trajectory point.
1872  if ( EnteredVol && InDet && (trj==NumTraj-1) ) {
1873  TrajEx = trj;
1874  }
1875  }
1876  // If last hit was in volume, and this hit isn't.
1877  if ( PrevHit && !InDet) {
1878  TrajEx = trj-1;
1879  }
1880  // Set the PrevHit variable.
1881  PrevHit = InDet;
1882  } // Loop over TrajPoints
1883  // Finally, return
1884  return EnteredVol;
1885  } //end InterceptsDetector
1886 
1887  //----------------------------------------------------------------------
1888 
1889 } // namespace
1890 
1891 
1892 namespace cheat{
1894 }
double E(const int i=0) const
Definition: MCParticle.h:232
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
float TNS() const
Definition: CellHit.h:46
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
std::vector< const T * > PtrVecToVecRawPtr(const art::PtrVector< T > &xs) const
Helper function for implementing overloads.
Definition: BackTracker.h:861
std::unordered_map< geo::OfflineChan, std::vector< sim::PhotonSignal > > fCellToPhotonSignal
map of PhotonSignals in each cell in event
Definition: BackTracker.h:941
set< int >::iterator it
int GetPlaneID() const
Plane ID.
Definition: FLSHit.h:37
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
BackTracker(const Parameters &params, art::ActivityRegistry &reg)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::set< int > fTrackIDs
G4 track ids for all particles in the events.
Definition: BackTracker.h:943
const sim::Particle * TrackIDToMotherParticle(int const &id) const
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
static bool sort_idxeff(const std::pair< int, NeutrinoEffPur > &a, const std::pair< int, NeutrinoEffPur > &b)
Tool for function SliceToNeutrinoIndex, tells how to sort by putting structures with highest efficien...
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
int TrackId() const
Definition: PhotonSignal.h:33
const Var weight
void HitToParticlesMap(const std::vector< const rb::Cluster * > &sliceList, const int &sliceIdx)
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Atom< std::string > GeantModuleLabel
Definition: BackTracker.h:102
unsigned short Plane() const
Definition: CellHit.h:39
#define DEFINE_ART_SERVICE(svc)
Definition: ServiceMacros.h:93
float energyFrac
fraction of hit energy from the particle with this trackID
Definition: BackTracker.h:41
bool InterceptsDetector(const sim::Particle &Par, int &TrajEnt, int &TrajEx)
A function which returns the trajectory point at which the particle first enters and leave the detect...
int GetCellID() const
Cell ID.
Definition: FLSHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
std::vector< int > SliceToOrderedNuIdsByEff(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
const char * p
Definition: xmltok.h:285
std::vector< int > SliceToOrderedNuIdsByEffPur(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
std::vector< NeutrinoWithIndex > allMCTruth() const
const double EnergyFromTrackId(const art::Ptr< rb::CellHit > &hit, int trackId, bool useBirksE=false) const
The total energy deposited in this cell by all FLSHits with this trackId.
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
double Mass() const
Definition: MCParticle.h:238
std::map< double, std::vector< cheat::TrackIDE > > fCellToIDs
Definition: BackTracker.h:951
cheat::ParticleEffPur ClusterToParticle(const std::vector< const rb::Cluster * > &sliceList, const std::vector< const rb::CellHit * > &hits, const int &sliceIdx)
list_type::const_iterator const_iterator
Vertical planes which measure X.
Definition: PlaneGeo.h:28
int trackID
Truth particle track ID.
Definition: BackTracker.h:66
float energy
energy from the particle with this trackID
Definition: BackTracker.h:42
Atom< std::string > PhotonModuleLabel
Definition: BackTracker.h:103
static bool sort_idxpur(const std::pair< int, NeutrinoEffPur > &a, const std::pair< int, NeutrinoEffPur > &b)
Tool for function SliceToNeutrinoIndex, tells how to sort by putting structures with highest purity f...
A collection of associated CellHits.
Definition: Cluster.h:47
std::vector< NeutrinoWithIndex > allNeutrinoInteractions() const
Function primarily for CAFMaker - returns a vector of all MCTruth interactions along with their truth...
std::vector< rb::Cluster > ClusterByTruth(std::vector< art::Ptr< rb::CellHit > > const &hits)
:
float CalcEscapingEnergy(const sim::Particle &Par, bool UseMCPart=true, std::string G4Label="trueens")
A function to calculate the amount of energy from a given MCParticle which is not contained in the de...
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void abs(TH1 *hist)
float CalcDaughterEscapingEnergy(const sim::Particle &Par)
A function to be used in conjunction with CalcTotalEscapingEnergy to calculate the energy which the d...
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
double purity
Putity of particle relative to track.
Definition: BackTracker.h:68
bool isRealData() const
Definition: Event.h:83
int NPhoton() const
Definition: PhotonSignal.h:30
const sim::Particle * HitToParticle(art::Ptr< rb::CellHit > const &hit, bool quiet=false) const
Returns the sim::Particle object contributing the most light to the specified rb::CellHit.
std::vector< int > SliceToOrderedNuIdsByPur(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
void AccumulateHitContributions(rb::CellHit const &hit, std::map< int, int > &idToNPhot) const
Internal helper, adds the contributions from each particle in this cell to the map.
double efficiency
Efficiency (based on FLS energy) of neutrino interaction relative to slice.
Definition: BackTracker.h:48
int NumberDaughters() const
Definition: MCParticle.h:216
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
int maxID(double arrayInput[])
float TotalEscEnergy() const
Definition: TrueEnergy.h:40
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
int TrackId() const
Definition: MCParticle.h:209
std::vector< int > SliceToOrderedNuIds(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable, std::function< double(const cheat::NeutrinoEffPur &)> slMetric, std::function< double(const cheat::NeutrinoEffPur &)> nuMetric) const
Given a vector of indexed neutrino interaction and a vector of slice truth tables, returns a vector of neutrino interaction indices ordered by best match to the corresponding slice. Here, best match is determined according to the given cheat::NeutrinoEffPur functions for the slice and the nu.
int Daughter(const int i) const
Definition: MCParticle.cxx:112
std::unordered_map< geo::OfflineChan, std::list< sim::FLSHit > > fCellToFLSHit
map of FLSHits in each cell in event
Definition: BackTracker.h:939
Atom< std::string > TrueEnergyModuleLabel
Definition: BackTracker.h:110
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
std::vector< art::Ptr< simb::MCTruth > > fMCTruthList
all the MCTruths for this event
Definition: BackTracker.h:944
int GetPDG() const
PDG.
Definition: FLSHit.h:43
std::vector< int > SliceToOrderedNuIdsByEnergy(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
unsigned short Cell() const
Definition: CellHit.h:40
const XML_Char * s
Definition: expat.h:262
Definition: Cand.cxx:23
TString part[npart]
Definition: Style.C:32
std::unordered_map< int, std::list< sim::FLSHit > > fTrackIDToFLSHit
map of G4 track ids to FLSHits
Definition: BackTracker.h:940
std::map< int, float > fIdToEnergy
Definition: BackTracker.h:952
Far Detector at Ash River, MN.
void hits()
Definition: readHits.C:15
int trackID
Geant4 supplied trackID.
Definition: BackTracker.h:40
bool IsHitsNoise(std::vector< art::Ptr< rb::CellHit > > const &hits) const
Of this collection of hits, is it mostly noise or physics hits? True indicates noise.
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
std::vector< art::Ptr< sim::TrueEnergy > > fTrueEnergyList
all the TrueEnergy&#39;s for this event
Definition: BackTracker.h:946
sim::ParticleNavigator fParticleNav
Particle navigator to map track ID to Particle.
Definition: BackTracker.h:942
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
const double a
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
double efficiency
Efficiency of particle relative to track.
Definition: BackTracker.h:67
double energy
Definition: plottest35.C:25
bool CompareByEnergy(const TrackIDE &a, const TrackIDE &b)
Does a have less energy than b?
Near Detector in the NuMI cavern.
const art::Ptr< simb::MCTruth > & ParticleToMCTruth(const sim::Particle *p) const
double EnergyMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur&#39;s nu interaction to slice energy.
bool isInsideDetectorBigBox(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Big Box?
std::map< int, int > fTrackIDToTrueEIndex
map of track id values to TrueEnergy entry
Definition: BackTracker.h:947
int Cell() const
Definition: PhotonSignal.h:31
int pdg
Pdg of particle.
Definition: BackTracker.h:69
z
Definition: test.py:28
size_type size() const
Definition: PtrVector.h:308
static bool sort_pur(const NeutrinoEffPur &a, const NeutrinoEffPur &b)
Tool for function SliceToNeutrinoInteractions, tells how to sort by putting structures with highest p...
double HitCollectionPurity(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, std::map< int, double > *purMap=0, std::map< int, int > *parents=0, bool energyPur=false) const
Returns the fraction of hits in a collection that come from the specified Geant4 track ids...
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
bool CompareByNPhoton(const sim::PhotonSignal &a, const sim::PhotonSignal &b)
Helper for SortByNPhoton.
void Rebuild(const art::Event &evt)
std::vector< TrackIDE > HitsToTrackIDE(const std::vector< const rb::CellHit * > &hits, bool useBirksE=false) const
Returns vector of TrackIDE structs contributing to the given collection of hits.
double Vx(const int i=0) const
Definition: MCParticle.h:220
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double EffPurMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur&#39;s nu interaction to slice efficiency * purity.
const std::vector< sim::PhotonSignal > HitToPhotonSignal(const art::Ptr< rb::CellHit > &hit) const
Returns the PhotonSignals contributing the signal in the specified hit.
const std::vector< sim::PhotonSignal > CellToPhotonSignal(unsigned int const &plane, unsigned int const &cell) const
Returns the PhotonSignals contributing the signal in the specified cell. WARNING: Use with extreme ca...
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
code to link reconstructed objects back to the MC truth information
Definition: FillTruth.h:16
std::vector< int > SliceToNeutrinoIndex(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of neutrino indices corresponding to the v...
int GetTrackID() const
Track ID.
Definition: FLSHit.h:45
const art::Ptr< simb::MCTruth > & TrackIDToMCTruth(int const &id) const
bool fHaveTruthInfo
Set by Rebuild.
Definition: BackTracker.h:935
T const * get() const
Definition: Ptr.h:321
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
GlobalSignal< detail::SignalResponseType::FIFO, void(Event const &)> sPreProcessEvent
std::map< int, int > fTrackIDToMCTruthIndex
map of track id values to MCTruthList entry
Definition: BackTracker.h:945
Definition: event.h:1
A (plane, cell) pair.
Definition: OfflineChan.h:17
static bool sortFLSHit(sim::FLSHit const &a, sim::FLSHit const &b)
double energySlice
Sum of FLS hits from the neutrino contributing to hits included in the slice.
Definition: BackTracker.h:51
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
const hit & b
Definition: hits.cxx:21
double Vz(const int i=0) const
Definition: MCParticle.h:222
assert(nhit_max >=nhit_nbins)
double PurMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur&#39;s nu interaction to slice purity.
static bool sort_eff(const NeutrinoEffPur &a, const NeutrinoEffPur &b)
Tool for function SliceToNeutrinoInteractions, tells how to sort by putting structures with highest e...
double EffMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur&#39;s nu interaction to slice efficiency.
particleList
Definition: demo.py:41
const std::vector< sim::FLSHit > CellToFLSHit(unsigned int const &plane, unsigned int const &cell) const
Returns the FLSHits contributing the signal in the specified cell. WARNING: Use with extreme caution!...
bool PassMuonCuts(int trackID, art::PtrVector< rb::CellHit > const &hits) const
Tool for NumuEAna that allows one to see if primary muon (or any track ID you feed the function) cont...
bool isNull() const
Definition: Ptr.h:328
bool NeutrinoSet() const
Definition: MCTruth.h:77
static Var totalEnergy(const std::shared_ptr< CAFAnaModel > &model)
Definition: LSTMEVar.h:77
float GetEdep() const
Get total Energy deposited into the cell for the whole FLSHit.
Definition: FLSHit.h:31
bool GoodTiming() const
Definition: CellHit.h:48
ParticleNavigator Add(const int &offset) const
TVector3 HitToXYZ(art::Ptr< rb::CellHit > const &hit, bool useBirksE=false) const
Returns the XYZ position of the energy deposition for a given hit.
Int_t trackID
Definition: plot.C:84
double maxE
Definition: plot_hist.C:8
Float_t e
Definition: plot.C:35
Float_t w
Definition: plot.C:20
void efficiency()
Definition: efficiency.C:58
float TotalDepEnergy() const
Definition: TrueEnergy.h:39
std::vector< ParticleEffPur > TracksToParticles(const std::vector< const rb::Track * > &tracks, const std::set< int > &trkIDs, const std::set< int > &allowedDaughters, const std::vector< const rb::CellHit * > &allHits) const
Given a collection of reconstructed tracks, this method finds the best match for each reco track to t...
iterator find(const key_type &key)
novadaq::cnv::DetId fDetID
detector id
Definition: BackTracker.h:937
std::vector< NeutrinoEffPur > SliceToMCTruth(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of MCTruth, efficiency, and purity of that neutrino interaction relative to the slice. Efficiency is defined as FLS energy from neutrino interaction in slice / total FLS energy from neutrino interaction in event. This vector is sorted from the highest efficiency interaction to lowest. This function returns all MCTruth, including those without neutrino interactions, i.e. cosmics.
Atom< double > MinContribFrac
Definition: BackTracker.h:105
float CalcTotalEscapingEnergy(const sim::Particle &Par, bool UseMCPart=true, std::string G4Label="trueens")
An extension of CalcEscapingEnergy which also takes into account the energy which the daughters of gi...
double TimeMean() const
Definition: PhotonSignal.h:29
rb::Cluster MCTruthToCluster(art::Ptr< simb::MCTruth > const &truth, std::vector< art::Ptr< rb::CellHit > > const &hits)
:
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
Atom< double > MinPhysicsFrac
Definition: BackTracker.h:107
art::Ptr< simb::MCTruth > neutrinoInt
Truth information about neutrino interaction.
Definition: BackTracker.h:47
double ideff[10]
Definition: Xdiff_gwt.C:104
int Plane() const
Definition: PhotonSignal.h:32
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
double purity
Purity (based on FLS energy) of neutrino interaction relative to slice.
Definition: BackTracker.h:49
int EveId(const int trackID) const
bool CloseInTime(const rb::CellHit &chit, const sim::PhotonSignal &phot) const
Internal helper, decides if a PhotonSignal and CellHit are close enough to match. ...