LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
PhotonBackTracker.cc
Go to the documentation of this file.
1 //
3 // \file PhotonBackTracker.cc
4 // \brief The functions needed for the PhotonBackTracker class needed by the
5 // PhotonBackTrackerService in order to connect truth information with
6 // reconstruction.
7 // \author jason.stock@mines.sdsmt.edu
8 //
9 // Based on the original BackTracker by brebel@fnal.gov
10 //
12 //
13 //TODO: Impliment alternate backtracking scheme developed by T. Usher
14 //TODO: OpChanToOpDetSDPs (Expanded Clone of OpDetNumToOpDetSDPs
15 //
17 
18 // LArSoft
23 
24 // Framework
26 
27 // STL
28 #include <map>
29 
30 namespace cheat {
31 
32  //----------------------------------------------------------------
33  PhotonBackTracker::PhotonBackTracker(fhiclConfig const& config,
34  cheat::ParticleInventory const* partInv,
35  geo::WireReadoutGeom const* wireReadoutGeom)
36  : fPartInv(partInv)
37  , fWireReadoutGeom(wireReadoutGeom)
38  , fDelay(config.Delay())
39  , fG4ModuleLabel(config.G4ModuleLabel())
40  , fG4ModuleLabels(config.G4ModuleLabels())
41  , fOpHitLabel(config.OpHitLabel())
42  , fOpFlashLabel(config.OpFlashLabel())
43  , fMinOpHitEnergyFraction(config.MinOpHitEnergyFraction())
44  {}
45 
46  //----------------------------------------------------------------
48  cheat::ParticleInventory const* partInv,
49  geo::WireReadoutGeom const* wireReadoutGeom)
50  : fPartInv(partInv)
51  , fWireReadoutGeom(wireReadoutGeom)
52  , fDelay(pSet.get<double>("Delay"))
53  , fG4ModuleLabel(pSet.get<art::InputTag>("G4ModuleLabel", "largeant"))
54  , fG4ModuleLabels(pSet.get<std::vector<art::InputTag>>("G4ModuleLabels", {}))
55  , fOpHitLabel(pSet.get<art::InputTag>("OpHitLabel", "ophit"))
56  , fOpFlashLabel(pSet.get<art::InputTag>("OpFlashLabel", "opflash"))
57  , fMinOpHitEnergyFraction(pSet.get<double>("MinimumOpHitEnergyFraction", 0.1))
58  {}
59 
60  //----------------------------------------------------------------
62  {
63  return fDelay;
64  }
65 
66  //----------------------------------------------------------------
68  {
69  priv_OpDetBTRs.clear();
70  priv_OpFlashToOpHits.clear();
71  }
72 
73  //----------------------------------------------------------------
75  {
76  return !priv_OpDetBTRs.empty();
77  }
78 
79  //----------------------------------------------------------------
81  {
82  return !priv_OpFlashToOpHits.empty();
83  }
84 
85  //----------------------------------------------------------------
86  const std::vector<art::Ptr<sim::OpDetBacktrackerRecord>>& PhotonBackTracker::OpDetBTRs() const
87  {
88  return priv_OpDetBTRs;
89  }
90 
91  //----------------------------------------------------------------
92  std::vector<const sim::SDP*> PhotonBackTracker::TrackIdToSimSDPs_Ps(int const id) const
93  {
94  std::vector<const sim::SDP*> sdp_Ps;
95  for (size_t odet = 0; odet < priv_OpDetBTRs.size(); ++odet) {
96  const auto& pdTimeSDPmap = priv_OpDetBTRs[odet]->timePDclockSDPsMap();
97  for (auto mapitr = pdTimeSDPmap.begin(); mapitr != pdTimeSDPmap.end(); mapitr++) {
98  std::vector<sim::SDP> const& sdpvec = (*mapitr).second;
99  for (size_t iv = 0; iv < sdpvec.size(); ++iv) {
100  if (abs(sdpvec[iv].trackID) == id) sdp_Ps.push_back(&(sdpvec[iv]));
101  }
102  } // end loop over map from sim::OpDetBacktrackerRecord
103  } // end loop over sim::OpDetBacktrackerRecords
104  return sdp_Ps;
105  }
106 
107  //----------------------------------------------------------------
109  {
111  for (size_t sc = 0; sc < priv_OpDetBTRs.size(); ++sc) {
112  // This could become a bug. What if it occurs twice (shouldn't happen in correct
113  // records, but still, no error handeling included for the situation
114  if (priv_OpDetBTRs.at(sc)->OpDetNum() == opDetNum) opDet = priv_OpDetBTRs.at(sc);
115  }
116  if (!opDet) {
117  throw cet::exception("PhotonBackTracker2") << "No sim:: OpDetBacktrackerRecord corresponding "
118  << "to opDetNum: " << opDetNum << "\n";
119  }
120  return opDet;
121  }
122 
123  //----------------------------------------------------------------
124  std::vector<sim::TrackSDP> PhotonBackTracker::OpDetToTrackSDPs(int const OpDetNum,
125  double const opHit_start_time,
126  double const opHit_end_time) const
127  {
128  std::vector<sim::TrackSDP> tSDPs;
129  double totalE = 0;
130  try {
131  const art::Ptr<sim::OpDetBacktrackerRecord> opDetBTR = FindOpDetBTR(OpDetNum);
132  std::vector<sim::SDP> simSDPs =
133  opDetBTR->TrackIDsAndEnergies(opHit_start_time, opHit_end_time);
134  for (size_t e = 0; e < simSDPs.size(); ++e)
135  totalE += simSDPs[e].energy;
136  if (totalE < 1.e-5) totalE = 1.;
137  for (size_t e = 0; e < simSDPs.size(); ++e) {
138  if (simSDPs[e].trackID == sim::NoParticleId) continue;
139  sim::TrackSDP info;
140  info.trackID = std::abs(simSDPs[e].trackID);
141  info.energyFrac = simSDPs[e].energy / totalE;
142  info.energy = simSDPs[e].energy;
143  tSDPs.push_back(info);
144  }
145  }
146  catch (cet::exception const& e) {
147  // This needs to go. Make it specific if there is a really an exception we would
148  // like to catch.
149  mf::LogWarning("PhotonBackTracker") << "Exception caught\n" << e;
150  }
151  return tSDPs;
152  }
153 
154  //----------------------------------------------------------------
155  std::vector<sim::TrackSDP> PhotonBackTracker::OpHitToTrackSDPs(
156  art::Ptr<recob::OpHit> const& opHit_P) const
157  {
158  auto OpDetNum = fWireReadoutGeom->OpDetFromOpChannel(opHit_P->OpChannel());
159  const double pTime = opHit_P->PeakTime();
160  const double pWidth = opHit_P->Width();
161  const double start = (pTime - pWidth) * 1000 - fDelay;
162  const double end = (pTime + pWidth) * 1000 - fDelay;
163  return OpDetToTrackSDPs(OpDetNum, start, end);
164  }
165 
166  //----------------------------------------------------------------
167  std::vector<sim::TrackSDP> PhotonBackTracker::OpHitToTrackSDPs(recob::OpHit const& opHit) const
168  {
169  auto OpDetNum = fWireReadoutGeom->OpDetFromOpChannel(opHit.OpChannel());
170  const double pTime = opHit.PeakTime();
171  const double pWidth = opHit.Width();
172  const double start = (pTime - pWidth) * 1000 - fDelay;
173  const double end = (pTime + pWidth) * 1000 - fDelay;
174  return OpDetToTrackSDPs(OpDetNum, start, end);
175  }
176 
177  //----------------------------------------------------------------
178  std::vector<int> PhotonBackTracker::OpHitToTrackIds(recob::OpHit const& opHit) const
179  {
180  std::vector<int> retVec;
181  for (auto const trackSDP : OpHitToTrackSDPs(opHit)) {
182  retVec.push_back(trackSDP.trackID);
183  }
184  return retVec;
185  }
186 
187  //----------------------------------------------------------------
188  std::vector<int> PhotonBackTracker::OpHitToTrackIds(art::Ptr<recob::OpHit> const& opHit) const
189  {
190  return OpHitToTrackIds(*opHit);
191  }
192 
193  //----------------------------------------------------------------
195  {
196  std::vector<int> retVec;
197  for (auto const trackSDP : OpHitToEveTrackSDPs(opHit)) {
198  retVec.push_back(trackSDP.trackID);
199  }
200  return retVec;
201  }
202 
203  //----------------------------------------------------------------
205  {
206  return OpHitToEveTrackIds(*opHit_P);
207  }
208 
209  //----------------------------------------------------------------
210  std::vector<sim::TrackSDP> PhotonBackTracker::OpHitToEveTrackSDPs(
211  art::Ptr<recob::OpHit> const& opHit_P) const
212  {
213  return OpHitToEveTrackSDPs(*opHit_P);
214  }
215 
216  //----------------------------------------------------------------
217  std::vector<sim::TrackSDP> PhotonBackTracker::OpHitToEveTrackSDPs(recob::OpHit const& opHit) const
218  {
219  std::vector<sim::TrackSDP> trackSDPs = OpHitToTrackSDPs(opHit);
220 
221  // make a map of evd ID values and fraction of energy represented by that eve id in
222  // this opHit
223  std::map<int, float> eveIDtoEfrac;
224 
225  double totalE = 0.;
226  for (size_t t = 0; t < trackSDPs.size(); ++t) {
227  eveIDtoEfrac[(fPartInv->ParticleList()).EveId(trackSDPs[t].trackID)] += trackSDPs[t].energy;
228  totalE += trackSDPs[t].energy;
229  }
230 
231  // now fill the eveSDPs vector from the map
232  std::vector<sim::TrackSDP> eveSDPs;
233  eveSDPs.reserve(eveIDtoEfrac.size());
234  for (auto itr = eveIDtoEfrac.begin(); itr != eveIDtoEfrac.end(); itr++) {
235  sim::TrackSDP temp;
236  temp.trackID = (*itr).first;
237  temp.energyFrac = (*itr).second / totalE;
238  temp.energy = (*itr).second;
239  eveSDPs.push_back(std::move(temp));
240  }
241  return eveSDPs;
242  }
243 
244  //----------------------------------------------------------------
245  //TODO: Make a copy of this function that uses an allOpHits list.
246  std::vector<art::Ptr<recob::OpHit>> PhotonBackTracker::TrackIdToOpHits_Ps(
247  int const tkId,
248  std::vector<art::Ptr<recob::OpHit>> const& hitsIn) const
249  {
250  // One would think we would want to have this function defined, and call this function
251  // in the std::vector<tkids> to opHits, but that would require more loops (and a
252  // higher overhead.) Instead, to provide this, we will just call the existing
253  // std::vector<tkids>ToOpHits with an input of 1.
254  std::vector<int> tkidFake(1, tkId);
255  return TrackIdsToOpHits_Ps(tkidFake, hitsIn).at(0);
256  }
257 
258  //----------------------------------------------------------------
259  std::vector<std::vector<art::Ptr<recob::OpHit>>> PhotonBackTracker::TrackIdsToOpHits_Ps(
260  std::vector<int> const& tkIds,
261  std::vector<art::Ptr<recob::OpHit>> const& hitsIn) const
262  {
263  std::vector<std::pair<int, art::Ptr<recob::OpHit>>> opHitList;
264  for (auto itr = hitsIn.begin(); itr != hitsIn.end(); ++itr) {
265  art::Ptr<recob::OpHit> const& opHit = *itr;
266  auto OpDetNum = fWireReadoutGeom->OpDetFromOpChannel(opHit->OpChannel());
267  const double pTime = opHit->PeakTime(), pWidth = opHit->Width();
268  const double start = (pTime - pWidth) * 1000.0 - fDelay,
269  end = (pTime + pWidth) * 1000.0 - fDelay;
270  std::vector<sim::TrackSDP> tids = OpDetToTrackSDPs(OpDetNum, start, end);
271  for (auto itid = tids.begin(); itid != tids.end(); ++itid) {
272  for (auto itkid = tkIds.begin(); itkid != tkIds.end(); ++itkid) {
273  if (itid->trackID == *itkid and itid->energyFrac > fMinOpHitEnergyFraction)
274  opHitList.push_back(std::make_pair(*itkid, opHit));
275  } // itkid
276  } // itid
277  } // itr
278  // now build the truOpHits vector that will be returned to the caller
279  std::vector<std::vector<art::Ptr<recob::OpHit>>> truOpHits;
280  // temporary vector containing opHits assigned to one MC particle
281  std::vector<art::Ptr<recob::OpHit>> tmpOpHits;
282  for (auto itkid = tkIds.begin(); itkid != tkIds.end(); ++itkid) {
283  tmpOpHits.clear();
284  for (auto itr = opHitList.begin(); itr != opHitList.end(); ++itr) {
285  if (*itkid == (*itr).first) tmpOpHits.push_back((*itr).second);
286  }
287  truOpHits.push_back(tmpOpHits);
288  }
289  return truOpHits;
290  }
291 
292  //----------------------------------------------------------------
293  std::vector<const sim::SDP*> PhotonBackTracker::OpHitToSimSDPs_Ps(recob::OpHit const& opHit) const
294  {
295  std::vector<const sim::SDP*> retVec;
296  double fPeakTime = opHit.PeakTime();
297  double fWidth = opHit.Width();
299  ((fPeakTime - fWidth) * 1000.0) - fDelay;
300  sim::OpDetBacktrackerRecord::timePDclock_t end_time = ((fPeakTime + fWidth) * 1000.0) - fDelay;
301  if (start_time > end_time) { throw; }
302 
303  auto const& timeSDPMap = FindOpDetBTR(fWireReadoutGeom->OpDetFromOpChannel(opHit.OpChannel()))
304  ->timePDclockSDPsMap(); //Not guranteed to be sorted.
305 
306  std::vector<const std::pair<double, std::vector<sim::SDP>>*> timePDclockSDPMap_SortedPointers;
307  for (auto& pair : timeSDPMap) {
308  timePDclockSDPMap_SortedPointers.push_back(&pair);
309  }
310  auto pairSort = [](auto& a, auto& b) { return a->first < b->first; };
311  if (!std::is_sorted(timePDclockSDPMap_SortedPointers.begin(),
312  timePDclockSDPMap_SortedPointers.end(),
313  pairSort)) {
314  std::sort(
315  timePDclockSDPMap_SortedPointers.begin(), timePDclockSDPMap_SortedPointers.end(), pairSort);
316  }
317 
318  //This section is a hack to make comparisons work right.
319  std::vector<sim::SDP> dummyVec;
320  std::pair<double, std::vector<sim::SDP>> start_timePair = std::make_pair(start_time, dummyVec);
321  std::pair<double, std::vector<sim::SDP>> end_timePair = std::make_pair(end_time, dummyVec);
322  auto start_timePair_P = &start_timePair;
323  auto end_timePair_P = &end_timePair;
324  //First interesting iterator.
325  auto mapFirst = std::lower_bound(timePDclockSDPMap_SortedPointers.begin(),
326  timePDclockSDPMap_SortedPointers.end(),
327  start_timePair_P,
328  pairSort);
329  //Last interesting iterator.
330  auto mapLast =
331  std::upper_bound(mapFirst, timePDclockSDPMap_SortedPointers.end(), end_timePair_P, pairSort);
332 
333  for (auto& mapitr = mapFirst; mapitr != mapLast; ++mapitr)
334  for (auto& sdp : (*mapitr)->second)
335  retVec.push_back(&sdp);
336 
337  return retVec;
338  }
339 
340  //----------------------------------------------------------------
341  std::vector<const sim::SDP*> PhotonBackTracker::OpHitToSimSDPs_Ps(
342  art::Ptr<recob::OpHit> const& opHit_P) const
343  {
344  return OpHitToSimSDPs_Ps(*opHit_P);
345  }
346 
347  //----------------------------------------------------------------
348  std::vector<double> PhotonBackTracker::SimSDPsToXYZ(std::vector<sim::SDP> const& sdps) const&
349  {
350  std::vector<double> xyz(3, -999.);
351  double x = 0.;
352  double y = 0.;
353  double z = 0.;
354  double w = 0.;
355  // loop over photons.
356  for (auto const& sdp : sdps) {
357  double weight = sdp.numPhotons;
358  w += weight;
359  x += weight * sdp.x;
360  y += weight * sdp.y;
361  z += weight * sdp.z;
362  } // end loop over sim::SDPs
363  // If the sum of the weights is still zero, then fail to return a value. A hit with
364  // no contributing photons does't make sense.
365  if (w < 1.e-5)
366  throw cet::exception("PhotonBackTracker")
367  << "No sim::SDPs providing non-zero number of photons"
368  << " can't determine originating location from truth\n";
369  xyz[0] = x / w;
370  xyz[1] = y / w;
371  xyz[2] = z / w;
372  return xyz;
373  }
374 
375  //----------------------------------------------------------------
376  std::vector<double> PhotonBackTracker::SimSDPsToXYZ(
377  std::vector<const sim::SDP*> const& sdps_Ps) const&
378  {
379  std::vector<double> xyz(3, -999.);
380  double x = 0.;
381  double y = 0.;
382  double z = 0.;
383  double w = 0.;
384  // loop over photons.
385  for (const sim::SDP* sdp_P : sdps_Ps) {
386  auto& sdp = *sdp_P;
387  double weight = sdp.numPhotons;
388  w += weight;
389  x += weight * sdp.x;
390  y += weight * sdp.y;
391  z += weight * sdp.z;
392  } // end loop over sim::SDPs
393  // If the sum of the weights is still zero, then fail to return a value. A hit with
394  // no contributing photons does't make sense.
395  if (w < 1.e-5)
396  throw cet::exception("PhotonBackTracker")
397  << "No sim::SDPs providing non-zero number of photons"
398  << " can't determine originating location from truth\n";
399  xyz[0] = x / w;
400  xyz[1] = y / w;
401  xyz[2] = z / w;
402  return xyz;
403  }
404 
405  //----------------------------------------------------------------
406  std::vector<double> PhotonBackTracker::OpHitToXYZ(recob::OpHit const& opHit)
407  {
408  return SimSDPsToXYZ(OpHitToSimSDPs_Ps(opHit));
409  }
410 
411  //----------------------------------------------------------------
412  std::vector<double> PhotonBackTracker::OpHitToXYZ(art::Ptr<recob::OpHit> const& opHit)
413  {
414  return SimSDPsToXYZ(OpHitToSimSDPs_Ps(*opHit));
415  }
416 
417  //----------------------------------------------------------------
418  std::vector<const sim::SDP*> PhotonBackTracker::OpHitsToSimSDPs_Ps(
419  std::vector<art::Ptr<recob::OpHit>> const& opHits_Ps) const
420  {
421  std::vector<const sim::SDP*> sdps_Ps;
422  for (auto opHit_P : opHits_Ps) {
423  std::vector<const sim::SDP*> to_add_sdps_Ps = OpHitToSimSDPs_Ps(opHit_P);
424  sdps_Ps.insert(sdps_Ps.end(), to_add_sdps_Ps.begin(), to_add_sdps_Ps.end());
425  }
426  return sdps_Ps;
427  }
428 
429  //----------------------------------------------------------------
430  std::vector<double> PhotonBackTracker::OpHitsToXYZ(
431  std::vector<art::Ptr<recob::OpHit>> const& opHits_Ps) const
432  {
433  const std::vector<const sim::SDP*> SDPs_Ps = OpHitsToSimSDPs_Ps(opHits_Ps);
434  return SimSDPsToXYZ(SDPs_Ps);
435  }
436 
437  //----------------------------------------------------------------
438  std::unordered_set<const sim::SDP*> PhotonBackTracker::OpHitToEveSimSDPs_Ps(
439  recob::OpHit const& opHit_P)
440  {
441  const std::vector<int> ids = OpHitToEveTrackIds(opHit_P);
442  std::unordered_set<const sim::SDP*> sdps;
443  for (auto const& id : ids) {
444  std::vector<const sim::SDP*> tmp_sdps = TrackIdToSimSDPs_Ps(id);
445  for (const sim::SDP* tmp_sdp : tmp_sdps) {
446  sdps.insert(tmp_sdp); //emplace not needed here.
447  }
448  }
449  return sdps;
450  }
451 
452  //----------------------------------------------------------------
453  std::unordered_set<const sim::SDP*> PhotonBackTracker::OpHitToEveSimSDPs_Ps(
454  art::Ptr<recob::OpHit>& opHit)
455  {
456  const std::vector<int> ids = OpHitToEveTrackIds(opHit);
457  std::unordered_set<const sim::SDP*> sdps;
458  for (auto const& id : ids) {
459  std::vector<const sim::SDP*> tmp_sdps = TrackIdToSimSDPs_Ps(id);
460  for (const sim::SDP* tmp_sdp : tmp_sdps) {
461  sdps.insert(tmp_sdp); //emplace not needed here.
462  }
463  }
464  return sdps;
465  }
466 
467  //----------------------------------------------------------------
468  std::set<int> PhotonBackTracker::GetSetOfEveIds() const
469  {
470  return fPartInv->GetSetOfEveIds();
471  }
472 
473  //----------------------------------------------------------------
475  {
476  return fPartInv->GetSetOfTrackIds();
477  }
478 
479  //----------------------------------------------------------------
481  std::vector<art::Ptr<recob::OpHit>> const& opHits_Ps) const
482  {
483  std::set<int> eveIds;
484  for (auto const& opHit_P : opHits_Ps) {
485  const std::vector<sim::TrackSDP> sdps = OpHitToEveTrackSDPs(opHit_P);
486  for (auto const& sdp : sdps) {
487  eveIds.insert(sdp.trackID);
488  } //end sdps
489  } //End for hits
490  return eveIds;
491  }
492 
493  //----------------------------------------------------------------
494  std::set<int> PhotonBackTracker::GetSetOfEveIds(std::vector<recob::OpHit> const& opHits) const
495  {
496  std::set<int> eveIds;
497  for (auto const& opHit : opHits) {
498  const std::vector<sim::TrackSDP> sdps = OpHitToEveTrackSDPs(opHit);
499  for (auto const& sdp : sdps) {
500  eveIds.insert(sdp.trackID);
501  } //end sdps
502  } //End for hits
503  return eveIds;
504  }
505 
506  //----------------------------------------------------------------
508  std::vector<art::Ptr<recob::OpHit>> const& opHits) const
509  {
510  std::set<int> tids;
511  for (auto const& opHit : opHits) {
512  for (auto const& sdp : OpHitToTrackSDPs(opHit)) {
513  tids.insert(sdp.trackID);
514  } //End for TrackSDPs
515  } //End for hits
516  return tids;
517  }
518 
519  //----------------------------------------------------------------
520  std::set<int> PhotonBackTracker::GetSetOfTrackIds(std::vector<recob::OpHit> const& opHits) const
521  {
522  std::set<int> tids;
523  for (auto const& opHit : opHits) {
524  for (auto const& sdp : OpHitToTrackSDPs(opHit)) {
525  tids.insert(sdp.trackID);
526  } //End for TrackSDPs
527  } //End for hits
528  return tids;
529  }
530 
531  //----------------------------------------------------------------
532  double PhotonBackTracker::OpHitCollectionPurity(std::set<int> const& tkIds,
533  std::vector<art::Ptr<recob::OpHit>> const& opHits)
534  {
535  // get the list of EveIDs that correspond to the opHits in this collection if the
536  // EveID shows up in the input list of tkIds, then it counts
537  float total = 1. * opHits.size();
538  float desired = 0.;
539  for (size_t h = 0; h < opHits.size(); ++h) {
540  art::Ptr<recob::OpHit> opHit = opHits[h];
541  std::vector<sim::TrackSDP> opHitTrackSDPs = OpHitToTrackSDPs(opHit);
542  // don't double count if this opHit has more than one of the desired track IDs
543  // associated with it
544  for (size_t e = 0; e < opHitTrackSDPs.size(); ++e) {
545  if (tkIds.find(opHitTrackSDPs[e].trackID) != tkIds.end()) {
546  desired += 1.;
547  break;
548  }
549  }
550  } // end loop over opHits
551  double purity = 0.;
552  if (total > 0) purity = desired / total;
553  return purity;
554  }
555 
556  //----------------------------------------------------------------
558  std::set<int> const& tkIds,
559  std::vector<art::Ptr<recob::OpHit>> const& opHits)
560  {
561  // get the list of EveIDs that correspond to the opHits in this collection if the
562  // EveID shows up in the input list of tkIds, then it counts
563  float total = 0;
564  float desired = 0.;
565  // don't have to check the view in the opHits collection because those are assumed to
566  // be from the object we are testing and will the correct view by definition then.
567  for (size_t h = 0; h < opHits.size(); ++h) {
568  art::Ptr<recob::OpHit> opHit = opHits[h];
569  std::vector<sim::TrackSDP> opHitTrackIDs = OpHitToTrackSDPs(opHit);
570  total += opHit->Area(); // sum up the charge in the cluster
571  // don't double count if this opHit has more than one of the desired track IDs
572  // associated with it
573  for (size_t e = 0; e < opHitTrackIDs.size(); ++e) {
574  if (tkIds.find(opHitTrackIDs[e].trackID) != tkIds.end()) {
575  desired += opHit->Area();
576  break;
577  }
578  }
579  } // end loop over opHits
580  double purity = 0.;
581  if (total > 0) purity = desired / total;
582  return purity;
583  }
584 
585  //----------------------------------------------------------------
587  std::set<int> const& tkIds,
588  std::vector<art::Ptr<recob::OpHit>> const& opHits,
589  std::vector<art::Ptr<recob::OpHit>> const& opHitsIn)
590  {
591  float desired = 0.;
592  float total = 0.;
593  for (size_t h = 0; h < opHits.size(); ++h) {
594  art::Ptr<recob::OpHit> opHit = opHits[h];
595  std::vector<sim::TrackSDP> opHitTrackSDPs = OpHitToTrackSDPs(opHit);
596  // also don't double count if this opHit has more than one of the desired track IDs
597  // associated with it
598  for (size_t e = 0; e < opHitTrackSDPs.size(); ++e) {
599  if (tkIds.find(opHitTrackSDPs[e].trackID) != tkIds.end() &&
600  opHitTrackSDPs[e].energyFrac >= fMinOpHitEnergyFraction) {
601  desired += 1.;
602  break;
603  }
604  }
605  } // end loop over opHits
606  // now figure out how many opHits in the whole collection are associated with this id
607  for (size_t h = 0; h < opHitsIn.size(); ++h) {
608  art::Ptr<recob::OpHit> opHit = opHitsIn[h];
609  std::vector<sim::TrackSDP> opHitTrackSDPs = OpHitToTrackSDPs(opHit);
610  for (size_t e = 0; e < opHitTrackSDPs.size(); ++e) {
611  // don't worry about opHits where the energy fraction for the chosen trackID is <
612  // 0.1; also don't double count if this opHit has more than one of the desired
613  // track IDs associated with it
614  if (tkIds.find(opHitTrackSDPs[e].trackID) != tkIds.end() &&
615  opHitTrackSDPs[e].energyFrac >= fMinOpHitEnergyFraction) {
616  total += 1.;
617  break;
618  }
619  }
620  } // end loop over all opHits
621  double efficiency = 0.;
622  if (total > 0.) efficiency = desired / total;
623  return efficiency;
624  }
625 
626  //----------------------------------------------------------------
628  std::set<int> const& tkIds,
629  std::vector<art::Ptr<recob::OpHit>> const& opHits,
630  std::vector<art::Ptr<recob::OpHit>> const& opHitsIn)
631  {
632  float desired = 0.;
633  float total = 0.;
634 
635  // don't have to check the view in the opHits collection because those are assumed to
636  // be from the object we are testing and will the correct view by definition then.
637  for (size_t h = 0; h < opHits.size(); ++h) {
638 
639  art::Ptr<recob::OpHit> opHit = opHits[h];
640  std::vector<sim::TrackSDP> opHitTrackIDs = OpHitToTrackSDPs(opHit);
641 
642  // don't worry about opHits where the energy fraction for the chosen trackID is <
643  // 0.1; also don't double count if this opHit has more than one of the desired track
644  // IDs associated with it
645  for (size_t e = 0; e < opHitTrackIDs.size(); ++e) {
646  if (tkIds.find(opHitTrackIDs[e].trackID) != tkIds.end() &&
647  opHitTrackIDs[e].energyFrac >= fMinOpHitEnergyFraction) {
648  desired += opHit->Area();
649  break;
650  }
651  }
652  } // end loop over opHits
653  for (size_t h = 0; h < opHitsIn.size(); ++h) {
654  art::Ptr<recob::OpHit> opHit = opHitsIn[h];
655  // check that we are looking at the appropriate view here in the case of 3D objects
656  // we take all opHits
657  std::vector<sim::TrackSDP> opHitTrackIDs = OpHitToTrackSDPs(opHit);
658  for (size_t e = 0; e < opHitTrackIDs.size(); ++e) {
659  // don't worry about opHits where the energy fraction for the chosen trackID is <
660  // 0.1; also don't double count if this opHit has more than one of the desired
661  // track IDs associated with it
662  if (tkIds.find(opHitTrackIDs[e].trackID) != tkIds.end() &&
663  opHitTrackIDs[e].energyFrac >= fMinOpHitEnergyFraction) {
664  total += opHit->Area();
665  break;
666  }
667  }
668  } // end loop over all opHits
669  double efficiency = 0.;
670  if (total > 0.) efficiency = desired / total;
671  return efficiency;
672  }
673  //--------------------------------------------------
674  std::vector<art::Ptr<recob::OpHit>> PhotonBackTracker::OpFlashToOpHits_Ps(
675  art::Ptr<recob::OpFlash>& flash_P) const
676  {
677  // There is not "non-pointer" version of this because the art::Ptr is needed to look
678  // up the assn. One could loop the Ptrs and dereference them, but I will not encourage
679  // the behavior by building the tool to do it.
680  return priv_OpFlashToOpHits.at(flash_P);
681  }
682 
683  //--------------------------------------------------
684  std::vector<double> PhotonBackTracker::OpFlashToXYZ(art::Ptr<recob::OpFlash>& flash_P) const
685  {
686  return OpHitsToXYZ(OpFlashToOpHits_Ps(flash_P));
687  }
688 
689  //--------------------------------------------------
691  {
692  std::vector<art::Ptr<recob::OpHit>> opHits_Ps = OpFlashToOpHits_Ps(flash_P);
693  std::set<int> ids;
694  for (auto& opHit_P : opHits_Ps) {
695  for (const int& id : OpHitToTrackIds(opHit_P)) {
696  ids.insert(id);
697  } // end for ids
698  } // end for opHits
699  return ids;
700  }
701 } // namespace
Float_t x
Definition: compare.C:6
std::vector< double > OpHitsToXYZ(std::vector< art::Ptr< recob::OpHit >> const &opHits_Ps) const
double OpHitLightCollectionPurity(std::set< int > const &tkIds, std::vector< art::Ptr< recob::OpHit >> const &opHits)
std::vector< double > OpHitToXYZ(art::Ptr< recob::OpHit > const &hit)
std::vector< double > OpFlashToXYZ(art::Ptr< recob::OpFlash > &flash_P) const
std::vector< sim::TrackSDP > OpDetToTrackSDPs(int OpDetNum, double opHit_start_time, double opHit_end_time) const
std::vector< const sim::SDP * > TrackIdToSimSDPs_Ps(int id) const
art::Ptr< sim::OpDetBacktrackerRecord > FindOpDetBTR(int opDetNum) const
Float_t y
Definition: compare.C:6
Double_t z
Definition: plot.C:276
std::vector< double > SimSDPsToXYZ(std::vector< sim::SDP > const &sdps)
constexpr auto abs(T v)
Returns the absolute value of the argument.
double PeakTime() const
Definition: OpHit.h:94
STL namespace.
std::vector< int > OpHitToEveTrackIds(recob::OpHit const &opHit)
double OpHitLightCollectionEfficiency(std::set< int > const &tkIds, std::vector< art::Ptr< recob::OpHit >> const &opHits, std::vector< art::Ptr< recob::OpHit >> const &opHitsIn)
std::set< int > OpFlashToTrackIds(art::Ptr< recob::OpFlash > &flash_P) const
std::vector< art::Ptr< sim::OpDetBacktrackerRecord > > const & OpDetBTRs() const
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
std::vector< const sim::SDP * > OpHitToSimSDPs_Ps(recob::OpHit const &opHit) const
double Width() const
Definition: OpHit.h:110
int trackID
Geant4 supplied trackID.
std::vector< const sim::SDP * > OpHitsToSimSDPs_Ps(std::vector< art::Ptr< recob::OpHit >> const &opHits_Ps) const
static const int NoParticleId
Definition: sim.h:21
Interface for a class providing readout channel mapping to geometry.
float energyFrac
fraction of OpHit energy from the particle with this trackID
double timePDclock_t
Type for iTimePDclock tick used in the interface.
double energy
Definition: plottest35.C:25
std::set< int > GetSetOfTrackIds() const
back track the reconstruction to the simulation
Ionization photons from a Geant4 track.
double weight
Definition: plottest35.C:25
std::vector< sim::TrackSDP > OpHitToTrackSDPs(art::Ptr< recob::OpHit > const &hit)
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
PhotonBackTracker(fhicl::ParameterSet const &pset, art::ActivityRegistry &reg)
code to link reconstructed objects back to the MC truth information
Definition: BackTracker.cc:24
bool OpFlashToOpHitsReady() const
std::vector< int > OpHitToTrackIds(recob::OpHit const &opHit) const
Float_t sc
Definition: plot.C:23
double OpHitCollectionEfficiency(std::set< int > trackIDs, std::vector< art::Ptr< recob::OpHit >> const &hits, std::vector< art::Ptr< recob::OpHit >> const &allhits)
std::vector< std::vector< art::Ptr< recob::OpHit > > > TrackIdsToOpHits_Ps(std::vector< int > const &tkIds, std::vector< art::Ptr< recob::OpHit >> const &hitsIn) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Header for the ParticleInvenotry Service Provider.
double OpHitCollectionPurity(std::set< int > trackIDs, std::vector< art::Ptr< recob::OpHit >> const &hits)
Definition: MVAAlg.h:12
double Area() const
Definition: OpHit.h:114
std::unordered_set< const sim::SDP * > OpHitToEveSimSDPs_Ps(recob::OpHit const &opHit)
std::vector< sim::SDP > TrackIDsAndEnergies(timePDclock_t startTimePDclock, timePDclock_t endTimePDclock) const
Return all the recorded energy deposition within a time interval.
std::set< int > GetSetOfEveIds() const
float energy
energy from the particle with this trackID [MeV]
Float_t e
Definition: plot.C:35
std::vector< art::Ptr< recob::OpHit > > OpFlashToOpHits_Ps(art::Ptr< recob::OpFlash > &flash_P) const
std::vector< art::Ptr< recob::OpHit > > TrackIdToOpHits_Ps(int tkId, std::vector< art::Ptr< recob::OpHit >> const &hitsIn) const
std::vector< sim::TrackSDP > OpHitToEveTrackSDPs(art::Ptr< recob::OpHit > const &opHit_P) const
Float_t w
Definition: plot.C:20
int OpChannel() const
Definition: OpHit.h:86
Tools and modules for checking out the basics of the Monte Carlo.
Interface to geometry for wire readouts .
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33