LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
BackTracker.cc
Go to the documentation of this file.
1 //
3 //
4 // \file BackTracker.cc
5 // \brief The functions needed for the BackTracker class needed by the BackTracker service in order to connect truth information with reconstruction.
6 // \author jason.stock@mines.sdsmt.edu
7 //
8 // Based on the original BackTracker by brebel@fnal.gov
9 //
10 //
12 
14 
15 #include "BackTracker.h"
17 
21 
22 namespace cheat {
23 
24  //-----------------------------------------------------------------------
26  const cheat::ParticleInventory* partInv,
27  const geo::GeometryCore* geom)
28  : fPartInv(partInv)
29  , fGeom(geom)
30  , fG4ModuleLabel(config.G4ModuleLabel())
31  , fSimChannelModuleLabel(config.SimChannelModuleLabel())
32  , fHitLabel(config.DefaultHitModuleLabel())
33  , fMinHitEnergyFraction(config.MinHitEnergyFraction())
34  , fOverrideRealData(config.OverrideRealData())
35  , fHitTimeRMS(config.HitTimeRMS())
36  {}
37 
38  //-----------------------------------------------------------------------
40  const cheat::ParticleInventory* partInv,
41  const geo::GeometryCore* geom)
42  : fPartInv(partInv)
43  , fGeom(geom)
44  , fG4ModuleLabel(pSet.get<art::InputTag>("G4ModuleLabel", "largeant"))
45  , fSimChannelModuleLabel(pSet.get<art::InputTag>("SimChannelModuleLabel", fG4ModuleLabel))
46  , // -- D.R. if not provided, default behavior is to use the G4ModuleLabel
47  fHitLabel(pSet.get<art::InputTag>("DefaultHitModuleLabel", "hitfd"))
48  , fMinHitEnergyFraction(pSet.get<double>("MinHitEnergyFraction", 0.010))
49  , fOverrideRealData(pSet.get<bool>("OverrideRealData", false))
50  , fHitTimeRMS(pSet.get<double>("HitTimeRMS", 1.0))
51  {}
52 
53  //-----------------------------------------------------------------------
55  {
56  fSimChannels.clear();
57  // fAllHitList.clear();
58  }
59 
60  //-----------------------------------------------------------------------
61  std::vector<const sim::IDE*> BackTracker::TrackIdToSimIDEs_Ps(int const& id) const
62  {
63  std::vector<const sim::IDE*> ideps;
64  for (size_t sc = 0; sc < fSimChannels.size(); ++sc) {
65  const auto& tdcidemap = fSimChannels[sc]->TDCIDEMap(); // This returns a reference.
66  // loop over the IDEMAP
67  for (auto mapitr = tdcidemap.begin(); mapitr != tdcidemap.end(); mapitr++) {
68  // loop over the vector of IDE objects.
69  const std::vector<sim::IDE>& idevec =
70  (*mapitr).second; // note, mapitr.second returns the actual data from
71  // the map, not a copy
72  for (size_t iv = 0; iv < idevec.size(); ++iv) {
73  // const sim::IDE* const idep = &idevec[iv];
74  // if( abs(idevec[iv].trackID) == id) continue;
75  // ideps.push_back(idep);
76  if (abs(idevec[iv].trackID) == id) ideps.push_back(&(idevec[iv]));
77  } // end for index in idevec
78  } // end loop over map from sim::SimChannel
79  } // end loop over sim::SimChannels
80  return ideps;
81  }
82 
83  //-----------------------------------------------------------------------
84  std::vector<const sim::IDE*> BackTracker::TrackIdToSimIDEs_Ps(int const& id,
85  const geo::View_t view) const
86  {
87  std::vector<const sim::IDE*> ide_Ps;
89  if (fGeom->View(sc->Channel()) != view) continue;
90 
91  // loop over the IDEMAP
92  for (const auto& item : sc->TDCIDEMap()) {
93 
94  // loop over the vector of IDE objects.
95  for (const sim::IDE& ide : item.second) {
96  if (abs(ide.trackID) == id) ide_Ps.push_back(&ide);
97  }
98  } // end loop over map from sim::SimChannel
99  } // end loop over sim::SimChannels
100 
101  return ide_Ps;
102  }
103 
104  //-----------------------------------------------------------------------
106  {
107  auto ilb = std::lower_bound(fSimChannels.begin(),
108  fSimChannels.end(),
109  channel,
110  [](art::Ptr<sim::SimChannel> const& a, raw::ChannelID_t channel) {
111  return (a->Channel() < channel);
112  });
113  return ((ilb != fSimChannels.end()) && ((*ilb)->Channel() == channel)) ?
114  *ilb :
116  }
117 
118  //-----------------------------------------------------------------------
120  {
121  if (auto const chan = FindSimChannelPtr(channel)) return chan;
122  throw cet::exception("BackTracker") << "No sim::SimChannel corresponding "
123  << "to channel: " << channel << "\n";
124  }
125 
126  //-----------------------------------------------------------------------
127  std::vector<sim::TrackIDE> BackTracker::ChannelToTrackIDEs(
128  detinfo::DetectorClocksData const& clockData,
129  raw::ChannelID_t channel,
130  const double hit_start_time,
131  const double hit_end_time) const
132  {
133  art::Ptr<sim::SimChannel> schannel = this->FindSimChannelPtr(channel);
134  if (!schannel) return {};
135 
136  std::vector<sim::TrackIDE> trackIDEs;
137  double totalE = 0.;
138 
139  // loop over the electrons in the channel and grab those that are in time
140  // with the identified hit start and stop times
141  int start_tdc = clockData.TPCTick2TDC(hit_start_time);
142  int end_tdc = clockData.TPCTick2TDC(hit_end_time);
143  if (start_tdc < 0) start_tdc = 0;
144  if (end_tdc < 0) end_tdc = 0;
145  std::vector<sim::IDE> simides = schannel->TrackIDsAndEnergies(start_tdc, end_tdc);
146 
147  // first get the total energy represented by all track ids for
148  // this channel and range of tdc values
149  for (size_t e = 0; e < simides.size(); ++e)
150  totalE += simides[e].energy;
151 
152  // protect against a divide by zero below
153  if (totalE < 1.e-5) totalE = 1.;
154 
155  // loop over the entries in the map and fill the input vectors
156 
157  for (size_t e = 0; e < simides.size(); ++e) {
158 
159  if (simides[e].trackID == sim::NoParticleId) continue;
160 
161  sim::TrackIDE info;
162  info.trackID = simides[e].trackID;
163  info.energyFrac = simides[e].energy / totalE;
164  info.energy = simides[e].energy;
165  info.numElectrons = simides[e].numElectrons;
166 
167  trackIDEs.push_back(info);
168  }
169 
170  return trackIDEs;
171  }
172 
173  //-----------------------------------------------------------------------
174  std::vector<sim::TrackIDE> BackTracker::HitToTrackIDEs(
175  detinfo::DetectorClocksData const& clockData,
176  recob::Hit const& hit) const
177  {
178  const double start = hit.PeakTimeMinusRMS(fHitTimeRMS);
179  const double end = hit.PeakTimePlusRMS(fHitTimeRMS);
180  return this->ChannelToTrackIDEs(clockData, hit.Channel(), start, end);
181  }
182 
183  //-----------------------------------------------------------------------
184  std::vector<int> BackTracker::HitToTrackIds(detinfo::DetectorClocksData const& clockData,
185  recob::Hit const& hit) const
186  {
187  std::vector<int> retVec;
188  for (auto const trackIDE : this->HitToTrackIDEs(clockData, hit)) {
189  retVec.push_back(trackIDE.trackID);
190  }
191  return retVec;
192  }
193 
194  // These don't exist in the event. They are generated on the fly.
195  //----------------------------------------------------------------------------
196  std::vector<sim::TrackIDE> BackTracker::HitToEveTrackIDEs(
197  detinfo::DetectorClocksData const& clockData,
198  recob::Hit const& hit) const
199  {
200  std::vector<sim::TrackIDE> eveIDEs;
201  std::vector<sim::TrackIDE> trackIDEs = this->HitToTrackIDEs(clockData, hit);
202  std::map<int, std::pair<double, double>> eveToEMap;
203  double totalE = 0.0;
204  for (const auto& trackIDE : trackIDEs) {
205  auto check = eveToEMap.emplace(fPartInv->TrackIdToEveTrackId(trackIDE.trackID),
206  std::make_pair(trackIDE.energy, trackIDE.numElectrons));
207  if (check.second == false) {
208  check.first->second.first += trackIDE.energy;
209  check.first->second.second += trackIDE.numElectrons;
210  }
211  // eveToEMap[fPartInv->TrackIdToEveTrackId(trackIDE.trackID)].first
212  // += trackIDE.energy;
213  totalE += trackIDE.energy;
214  } // End for trackIDEs
215  eveIDEs.reserve(eveToEMap.size());
216  for (const auto& eveToE : eveToEMap) {
217  sim::TrackIDE eveTrackIDE_tmp;
218 
219  eveTrackIDE_tmp.trackID = eveToE.first;
220  eveTrackIDE_tmp.energy = eveToE.second.first;
221  eveTrackIDE_tmp.energyFrac = (eveTrackIDE_tmp.energy) / (totalE);
222  eveTrackIDE_tmp.numElectrons = eveToE.second.second;
223 
224  eveIDEs.push_back(eveTrackIDE_tmp);
225  } // END eveToEMap loop
226  return eveIDEs;
227  }
228 
229  //-----------------------------------------------------------------------
230  std::vector<art::Ptr<recob::Hit>> BackTracker::TrackIdToHits_Ps(
231  detinfo::DetectorClocksData const& clockData,
232  const int tkId,
233  std::vector<art::Ptr<recob::Hit>> const& hitsIn) const
234  {
235  // returns a subset of the hits in the hitsIn collection that are matched
236  // to the given track
237 
238  // temporary vector of TrackIds and Ptrs to hits so only one
239  // loop through the (possibly large) hitsIn collection is needed
240  std::vector<art::Ptr<recob::Hit>> hitList;
241  std::vector<sim::TrackIDE> trackIDE;
242  for (auto itr = hitsIn.begin(); itr != hitsIn.end(); ++itr) {
243  trackIDE.clear();
244  art::Ptr<recob::Hit> const& hit = *itr;
245  trackIDE = this->ChannelToTrackIDEs(clockData,
246  hit->Channel(),
249  for (auto itr_trakIDE = trackIDE.begin(); itr_trakIDE != trackIDE.end(); ++itr_trakIDE) {
250  if (itr_trakIDE->trackID == tkId && itr_trakIDE->energyFrac > fMinHitEnergyFraction)
251  hitList.push_back(hit);
252  } // itr_trakIDE
253  } // itr
254  return hitList;
255  }
256 
257  //-----------------------------------------------------------------------
258  //This function could clearly be made by calling TrackIdToHits for each trackId, but that would be significantly slower because we would loop through all hits many times.
259  std::vector<std::vector<art::Ptr<recob::Hit>>> BackTracker::TrackIdsToHits_Ps(
260  detinfo::DetectorClocksData const& clockData,
261  std::vector<int> const& tkIds,
262  std::vector<art::Ptr<recob::Hit>> const& hitsIn) const
263  {
264  // returns a subset of the hits in the hitsIn collection that are matched
265  // to MC particles listed in tkIds
266 
267  // temporary vector of TrackIds and Ptrs to hits so only one
268  // loop through the (possibly large) hitsIn collection is needed
269  std::vector<std::pair<int, art::Ptr<recob::Hit>>> hitList;
270  std::vector<sim::TrackIDE> tids;
271  for (auto itr = hitsIn.begin(); itr != hitsIn.end(); ++itr) {
272  tids.clear();
273  art::Ptr<recob::Hit> const& hit = *itr;
274  tids = this->ChannelToTrackIDEs(clockData,
275  hit->Channel(),
278  for (auto itid = tids.begin(); itid != tids.end(); ++itid) {
279  for (auto itkid = tkIds.begin(); itkid != tkIds.end(); ++itkid) {
280  if (itid->trackID == *itkid) {
281  if (itid->energyFrac > fMinHitEnergyFraction)
282  hitList.push_back(std::make_pair(*itkid, hit));
283  }
284  } // itkid
285  } // itid
286  } // itr
287 
288  // now build the truHits vector that will be returned to the caller
289  std::vector<std::vector<art::Ptr<recob::Hit>>> truHits;
290  // temporary vector containing hits assigned to one MC particle
291  std::vector<art::Ptr<recob::Hit>> tmpHits;
292  for (auto itkid = tkIds.begin(); itkid != tkIds.end(); ++itkid) {
293  tmpHits.clear();
294  for (auto itr = hitList.begin(); itr != hitList.end(); ++itr) {
295  if (*itkid == (*itr).first) tmpHits.push_back((*itr).second);
296  }
297  truHits.push_back(tmpHits);
298  }
299  return truHits;
300  }
301 
302  //-----------------------------------------------------------------------
303  //Cannot be returned as a reference, as these IDEs do not exist in the event. They are constructed on the fly.
304  std::vector<sim::IDE> BackTracker::HitToAvgSimIDEs(detinfo::DetectorClocksData const& clockData,
305  recob::Hit const& hit) const
306  {
307  // Get services.
308 
309  int start_tdc = clockData.TPCTick2TDC(hit.PeakTimeMinusRMS(fHitTimeRMS));
310  int end_tdc = clockData.TPCTick2TDC(hit.PeakTimePlusRMS(fHitTimeRMS));
311  if (start_tdc < 0) start_tdc = 0;
312  if (end_tdc < 0) end_tdc = 0;
313 
314  return (this->FindSimChannel(hit.Channel()))->TrackIDsAndEnergies(start_tdc, end_tdc);
315  }
316 
317  //-----------------------------------------------------------------------
318  std::vector<const sim::IDE*> BackTracker::HitToSimIDEs_Ps(
319  detinfo::DetectorClocksData const& clockData,
320  recob::Hit const& hit) const
321  {
322  std::vector<const sim::IDE*> retVec;
323  int start_tdc = clockData.TPCTick2TDC(hit.PeakTimeMinusRMS(fHitTimeRMS));
324  int end_tdc = clockData.TPCTick2TDC(hit.PeakTimePlusRMS(fHitTimeRMS));
325  if (start_tdc < 0) start_tdc = 0;
326  if (end_tdc < 0) end_tdc = 0;
327 
328  if (start_tdc > end_tdc) { throw; }
329 
330  const std::vector<std::pair<unsigned short, std::vector<sim::IDE>>>& tdcIDEMap =
331  (this->FindSimChannel(hit.Channel()))
332  ->TDCIDEMap(); // This in fact does not return a std::map. It returns a
333  // vector... with no guarantee that it is sorted...
334  std::vector<const std::pair<unsigned short, std::vector<sim::IDE>>*> tdcIDEMap_SortedPointers;
335  for (auto& pair : tdcIDEMap) {
336  tdcIDEMap_SortedPointers.push_back(&pair);
337  }
338 
339  // This is a bunch of extra steps, due to needing a vector we can sort, and
340  // needing those items in the sorted vector to be the items from the sim
341  // channels (so a pointer to the IDEs inside the sim channels can be made).
342  // The work around is to make a vector of pointers to IDEs inside the
343  // TDCIDEMap (which is a constant reference to the fTDCIDEs in the
344  // SimChannel.)
345  auto pairSort = [](auto& a, auto& b) { return a->first < b->first; };
346  if (!std::is_sorted(
347  tdcIDEMap_SortedPointers.begin(), tdcIDEMap_SortedPointers.end(), pairSort)) {
348  std::sort(tdcIDEMap_SortedPointers.begin(), tdcIDEMap_SortedPointers.end(), pairSort);
349  }
350 
351  std::vector<sim::IDE> dummyVec; // I need something to stick in a pair to compare pair<tdcVal,
352  // IDE>. This is an otherwise useless "hack".
353  std::pair<double, std::vector<sim::IDE>> start_tdcPair =
354  std::make_pair(start_tdc, dummyVec); // This pair is a "hack" to make my comparison work
355  // for lower and upper bound.
356  std::pair<double, std::vector<sim::IDE>> end_tdcPair = std::make_pair(end_tdc, dummyVec);
357  auto start_tdcPair_P = &start_tdcPair;
358  auto end_tdcPair_P = &end_tdcPair;
359  auto mapFirst = std::lower_bound(
360  tdcIDEMap_SortedPointers.begin(), tdcIDEMap_SortedPointers.end(), start_tdcPair_P, pairSort);
361 
362  // iterator to just after the last interesting IDE
363  auto mapLast = std::upper_bound(
364  tdcIDEMap_SortedPointers.begin(), tdcIDEMap_SortedPointers.end(), end_tdcPair_P, pairSort);
365  for (auto& mapitr = mapFirst; mapitr != mapLast; ++mapitr) {
366  for (auto& ide : (*mapitr)->second) {
367  retVec.push_back(&ide);
368  } // Add all interesting IDEs to the retVec
369  }
370  return retVec;
371  }
372 
373  //------------------------------------------------------------------------------
374  std::vector<double> BackTracker::SimIDEsToXYZ(std::vector<sim::IDE> const& ides) const
375  {
376  std::vector<double> xyz(3, 0.0);
377  double w = 0.0;
378  for (auto const& ide : ides) {
379  double weight = ide.numElectrons;
380  w += weight;
381  xyz[0] += (weight * ide.x);
382  xyz[1] += (weight * ide.y);
383  xyz[2] += (weight * ide.z);
384  }
385  if (w < 1.e-5)
386  throw cet::exception("BackTracker") << "No sim::IDEs providing non-zero number of electrons"
387  << " can't determine originating location from truth\n";
388  xyz[0] = xyz[0] / w;
389  xyz[1] = xyz[1] / w;
390  xyz[2] = xyz[2] / w;
391  return xyz;
392  }
393 
394  //-------------------------------------------------------------------------------
395  std::vector<double> BackTracker::SimIDEsToXYZ(std::vector<const sim::IDE*> const& ide_Ps) const
396  {
397  std::vector<sim::IDE> ides;
398  for (auto ide_P : ide_Ps) {
399  ides.push_back(*ide_P);
400  }
401  return this->SimIDEsToXYZ(ides);
402  }
403 
404  //--------------------------------------------------------------------------------
405  std::vector<double> BackTracker::HitToXYZ(detinfo::DetectorClocksData const& clockData,
406  const recob::Hit& hit) const
407  {
408  std::vector<const sim::IDE*> ide_Ps = this->HitToSimIDEs_Ps(clockData, hit);
409  return this->SimIDEsToXYZ(ide_Ps);
410  }
411 
412  //-----------------------------------------------------------------------------------
414  std::set<int> const& trackIds,
415  std::vector<art::Ptr<recob::Hit>> const& hits) const
416  {
417  int desired = 0;
418  for (const auto& hit : hits) {
419  std::vector<sim::TrackIDE> hitTrackIDEs = this->HitToTrackIDEs(clockData, hit);
420  for (const auto& tIDE : hitTrackIDEs) {
421  if (trackIds.find(tIDE.trackID) != trackIds.end()) {
422  ++desired;
423  break;
424  } // End if TID Found
425  } // END for trackIDE in TrackIDEs
426  } // End for hit in hits
427  if (hits.size() > 0) { return double(double(desired) / double(hits.size())); }
428  return 0;
429  }
430 
431  //-----------------------------------------------------------------------------------
433  std::set<int> const& trackIds,
434  std::vector<art::Ptr<recob::Hit>> const& hits) const
435  {
436  double totalCharge = 0., desired = 0.;
437  for (const auto& hit : hits) {
438  totalCharge += hit->Integral();
439  std::vector<sim::TrackIDE> trackIDEs = this->HitToTrackIDEs(clockData, hit);
440  for (const auto& trackIDE : trackIDEs) {
441  if (trackIds.find(trackIDE.trackID) != trackIds.end()) {
442  desired += hit->Integral();
443  break;
444  } // End if trackId in trackIds.
445  } // End for trackIDE in trackIDEs
446  } // End for Hit in Hits
447  if (totalCharge > 0.0) { return (desired / totalCharge); }
448  return 0.0;
449  }
450 
451  //-----------------------------------------------------------------------------------
453  std::set<int> const& trackIds,
455  std::vector<art::Ptr<recob::Hit>> const& allHits,
456  geo::View_t const& view) const
457  {
458 
459  int desired = 0, total = 0;
460 
461  for (const auto& hit : hits) {
462  std::vector<sim::TrackIDE> hitTrackIDEs = this->HitToTrackIDEs(clockData, hit);
463  for (const auto& trackIDE : hitTrackIDEs) {
464  if (trackIds.find(trackIDE.trackID) != trackIds.end() &&
465  trackIDE.energyFrac >= fMinHitEnergyFraction) {
466  ++desired;
467  break;
468  } // End if trackID in trackIds.
469  } // end for trackIDE in TrackIDEs
470  } // end for hit in hits
471 
472  for (const auto& hit : allHits) {
473  if (hit->View() != view && view != geo::k3D) {
474  continue;
475  } // End if hit.view = view or view = geo::k3D
476  std::vector<sim::TrackIDE> hitTrackIDEs = this->HitToTrackIDEs(clockData, hit);
477  for (const auto& hitIDE : hitTrackIDEs) {
478  if (trackIds.find(hitIDE.trackID) != trackIds.end() &&
479  hitIDE.energyFrac >= fMinHitEnergyFraction) {
480  ++total;
481  break;
482  }
483  } // END for all IDEs in HitTrackIDEs.
484  } // end for hit in allHits.
485  if (total >= 0) { return double(double(desired) / double(total)); }
486  return 0.;
487  }
488 
489  //-----------------------------------------------------------------------------------
491  detinfo::DetectorClocksData const& clockData,
492  std::set<int> const& trackIds,
494  std::vector<art::Ptr<recob::Hit>> const& allHits,
495  geo::View_t const& view) const
496  {
497  double desired = 0., total = 0.;
498  for (const auto& hit : hits) {
499  std::vector<sim::TrackIDE> hitTrackIDEs = this->HitToTrackIDEs(clockData, hit);
500  for (const auto& hitIDE : hitTrackIDEs) {
501  if (trackIds.find(hitIDE.trackID) != trackIds.end() &&
502  hitIDE.energyFrac >= fMinHitEnergyFraction) {
503  desired += hit->Integral();
504  break;
505  } // end if hit id matches and energy sufficient.
506  } // End for IDE in HitTrackIDEs.
507  } // End for hit in hits.
508 
509  for (const auto& hit : allHits) {
510  if (hit->View() != view && view != geo::k3D) { continue; }
511  std::vector<sim::TrackIDE> hitTrackIDEs = this->HitToTrackIDEs(clockData, hit);
512  for (const auto& hitIDE : hitTrackIDEs) {
513  if (trackIds.find(hitIDE.trackID) != trackIds.end() &&
514  hitIDE.energyFrac >= fMinHitEnergyFraction) {
515  total += hit->Integral();
516  break;
517  } // end if hit matches
518  } // end for ide in ides
519  } // End for hit in allHits
520 
521  if (total > 0.) { return desired / total; }
522  return 0.;
523  }
524 
525  //-----------------------------------------------------------------------------------
527  std::vector<art::Ptr<recob::Hit>> const& hits) const
528  {
529  std::set<int> tids;
530  for (const auto& hit : hits) {
531  const double start = hit->PeakTimeMinusRMS(fHitTimeRMS);
532  const double end = hit->PeakTimePlusRMS(fHitTimeRMS);
533  std::vector<sim::TrackIDE> trackIDEs =
534  this->ChannelToTrackIDEs(clockData, hit->Channel(), start, end);
535  for (const auto& ide : trackIDEs) {
536  tids.insert(ide.trackID);
537  } // End for TrackIDEs
538  } // End for hits
539  return tids;
540  } // End GetSetOfTrackIds
541 
542  //-----------------------------------------------------------------------------------
544  std::vector<art::Ptr<recob::Hit>> const& hits) const
545  {
546  std::set<int> eveIds;
547  for (const auto& hit : hits) {
548  const std::vector<sim::TrackIDE> ides = this->HitToEveTrackIDEs(clockData, hit);
549  for (const auto& ide : ides) {
550  eveIds.insert(ide.trackID);
551  } // end ides
552  } // End for hits
553  return eveIds;
554  }
555 
556  // This function definitely needs a new implementation. There must
557  // be a better way than so many loops.
559  detinfo::DetectorClocksData const& clockData,
560  std::vector<art::Ptr<recob::Hit>> const& hits) const
561  {
562  std::vector<std::vector<std::vector<int>>> numHits(fGeom->Ncryostats());
563  std::vector<std::vector<std::vector<double>>> hitWeight(fGeom->Ncryostats());
564  std::vector<std::vector<std::vector<std::vector<double>>>> hitPos(fGeom->Ncryostats());
565  // Do we need to resize everything...
566  for (size_t c = 0; c < numHits.size(); ++c) {
567  geo::CryostatID const cid(c);
568  numHits[c].resize(fGeom->NTPC(cid));
569  hitWeight[c].resize(fGeom->NTPC(cid));
570  hitPos[c].resize(fGeom->NTPC(cid));
571  for (size_t t = 0; t < numHits[c].size(); ++t) {
572  geo::TPCID const tpcid(cid, t);
573  numHits[c][t].resize(fGeom->Nplanes(tpcid));
574  hitWeight[c][t].resize(fGeom->Nplanes(tpcid));
575  hitPos[c][t].resize(fGeom->Nplanes(tpcid));
576  }
577  }
578 
579  for (auto const& hit_ptr : hits) {
580  const recob::Hit& hit = *hit_ptr;
581 
582  // use the HitToXYZ and Geometry::PositionToTPC
583  // to figure out which drift volume the hit originates from
584  std::vector<double> hitOrigin = this->HitToXYZ(clockData, hit_ptr);
585  geo::Point_t const worldLoc{hitOrigin[0], hitOrigin[1], hitOrigin[2]};
586  auto const tpcid = fGeom->PositionToTPCID(worldLoc);
587 
588  auto const [cstat, tpc] = std::make_tuple(tpcid.Cryostat, tpcid.TPC);
589 
590  if (hit.WireID().Cryostat == cstat && hit.WireID().TPC == tpc) {
591  ++numHits[cstat][tpc][hit.WireID().Plane];
592  hitWeight[cstat][tpc][hit.WireID().Plane] = hit.Integral();
593  hitPos[cstat][tpc][hit.WireID().Plane] = hitOrigin;
594  }
595  }
596 
597  // loop over the vectors we made and find the average position for the hits
598  // in the future we might take a weighted average
599  int nhits = 0;
600  std::vector<double> xyz(3);
601  for (size_t c = 0; c < numHits.size(); ++c) {
602  for (size_t t = 0; t < numHits[c].size(); ++t) {
603  for (size_t p = 0; p < numHits[c][t].size(); ++p) {
604 
605  if (numHits[c][t][p] == 1) {
606  ++nhits;
607  xyz[0] += hitPos[c][t][p][0];
608  xyz[1] += hitPos[c][t][p][1];
609  xyz[2] += hitPos[c][t][p][2];
610  }
611 
612  } // end loop over planes
613  } // end loop over tpcs
614  } // end loop over cryostats
615 
616  // get the average position
617  if (nhits < 1)
618  throw cet::exception("BackTracker")
619  << "No hits to determine originating location from truth\n";
620 
621  xyz[0] /= nhits;
622  xyz[1] /= nhits;
623  xyz[2] /= nhits;
624 
625  // Done.
626  return xyz;
627  }
628 
629 } // End namespace cheat
TrackID_t trackID
Geant4 supplied track ID.
Definition: SimChannel.h:107
const art::InputTag fHitLabel
Definition: BackTracker.h:244
unsigned int NTPC(CryostatID const &cryoid=cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:686
const geo::GeometryCore * fGeom
Definition: BackTracker.h:241
std::vector< art::Ptr< sim::SimChannel > > fSimChannels
Definition: BackTracker.h:249
std::vector< std::vector< art::Ptr< recob::Hit > > > TrackIdsToHits_Ps(detinfo::DetectorClocksData const &clockData, std::vector< int > const &tkIds, std::vector< art::Ptr< recob::Hit >> const &hitsIn) const
Definition: BackTracker.cc:259
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
art::Ptr< sim::SimChannel > FindSimChannelPtr(raw::ChannelID_t channel) const
Returns the cached sim::SimChannel on the specified channel.
Definition: BackTracker.cc:105
std::vector< sim::TrackIDE > HitToEveTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Definition: BackTracker.cc:196
BackTracker(const fhiclConfig &config, const cheat::ParticleInventory *partInv, const geo::GeometryCore *geom)
Definition: BackTracker.cc:25
const art::InputTag fSimChannelModuleLabel
Definition: BackTracker.h:243
float numElectrons
number of electrons from the particle detected on the wires
Definition: SimChannel.h:30
constexpr auto abs(T v)
Returns the absolute value of the argument.
const double fMinHitEnergyFraction
Definition: BackTracker.h:245
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:244
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
float energy
energy from the particle with this trackID [MeV]
Definition: SimChannel.h:29
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:430
pure virtual base interface for detector clocks
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Definition: BackTracker.cc:174
3-dimensional objects, potentially hits, clusters, prongs, etc.
Definition: geo_types.h:141
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
Access the description of detector geometry.
void hits()
Definition: readHits.C:15
std::vector< double > HitToXYZ(detinfo::DetectorClocksData const &clockData, const recob::Hit &hit) const
Definition: BackTracker.cc:405
Ionization at a point of the TPC sensitive volume.
Definition: SimChannel.h:85
static const int NoParticleId
Definition: sim.h:21
art::Ptr< sim::SimChannel > FindSimChannel(raw::ChannelID_t channel) const
Returns the cached sim::SimChannel on the specified channel.
Definition: BackTracker.cc:119
std::vector< art::Ptr< recob::Hit > > TrackIdToHits_Ps(detinfo::DetectorClocksData const &clockData, int tkId, std::vector< art::Ptr< recob::Hit >> const &hitsIn) const
Definition: BackTracker.cc:230
std::set< int > GetSetOfTrackIds() const
Definition: BackTracker.h:227
double HitCollectionPurity(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits) const
Definition: BackTracker.cc:413
double energy
Definition: plottest35.C:25
float energyFrac
fraction of hit energy from the particle with this trackID
Definition: SimChannel.h:28
const cheat::ParticleInventory * fPartInv
Definition: BackTracker.h:240
const double fHitTimeRMS
Definition: BackTracker.h:247
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
Description of geometry of one entire detector.
Definition: GeometryCore.h:119
Definition of data types for geometry description.
float PeakTimeMinusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Definition: Hit.h:290
Detector simulation of raw signals on wires.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
std::vector< sim::TrackIDE > ChannelToTrackIDEs(detinfo::DetectorClocksData const &clockData, raw::ChannelID_t channel, const double hit_start_time, const double hit_end_time) const
Definition: BackTracker.cc:127
double weight
Definition: plottest35.C:25
std::set< int > GetSetOfEveIds() const
Definition: BackTracker.h:228
double HitChargeCollectionPurity(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits) const
Definition: BackTracker.cc:432
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
code to link reconstructed objects back to the MC truth information
Definition: BackTracker.cc:22
int trackID
Geant4 supplied trackID.
Definition: SimChannel.h:27
const bool fOverrideRealData
Definition: BackTracker.h:246
Contains all timing reference information for the detector.
std::vector< sim::IDE > TrackIDsAndEnergies(TDC_t startTDC, TDC_t endTDC) const
Return all the recorded energy deposition within a time interval.
Definition: SimChannel.cxx:153
std::vector< int > HitToTrackIds(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Definition: BackTracker.cc:184
Float_t sc
Definition: plot.C:23
View_t View(PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
double HitCollectionEfficiency(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits, std::vector< art::Ptr< recob::Hit >> const &allhits, geo::View_t const &view) const
Definition: BackTracker.cc:452
std::vector< const sim::IDE * > TrackIdToSimIDEs_Ps(int const &id) const
Definition: BackTracker.cc:61
const art::InputTag fG4ModuleLabel
Definition: BackTracker.h:242
Definition: MVAAlg.h:12
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
float PeakTimePlusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Definition: Hit.h:285
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
std::vector< sim::IDE > HitToAvgSimIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Definition: BackTracker.cc:304
Float_t e
Definition: plot.C:35
int TrackIdToEveTrackId(const int &tid) const
Float_t w
Definition: plot.C:20
Ionization energy from a Geant4 track.
Definition: SimChannel.h:26
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
Definition: Hit.h:268
std::vector< double > SimIDEsToXYZ(std::vector< sim::IDE > const &ides) const
Definition: BackTracker.cc:374
double HitChargeCollectionEfficiency(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits, std::vector< art::Ptr< recob::Hit >> const &allhits, geo::View_t const &view) const
Definition: BackTracker.cc:490
double TPCTick2TDC(double const tick) const
Tools and modules for checking out the basics of the Monte Carlo.
TPCID PositionToTPCID(Point_t const &point) const
Returns the ID of the TPC at specified location.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
The data type to uniquely identify a cryostat.
Definition: geo_types.h:192
std::vector< const sim::IDE * > HitToSimIDEs_Ps(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Definition: BackTracker.cc:318
std::vector< double > SpacePointHitsToWeightedXYZ(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> const &hits) const
Definition: BackTracker.cc:558