LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
DisambigAlg.cxx
Go to the documentation of this file.
1 //
3 // DisambigAlg.cxx
4 //
5 // tylerdalion@gmail.com
6 //
7 // description
8 //
9 //
10 // Since we are letting the module make the final column of hits, a disambiguated
11 // hit in this algorithm is a pair< art::Ptr<recob::Hit>, geo::WireID >
12 //
13 //
15 
16 //Framework includes:
18 #include "fhiclcpp/ParameterSet.h"
20 
28 
29 #include <cmath>
30 #include <cstdlib>
31 #include <map>
32 
33 #include "range/v3/view.hpp"
34 
35 using lar::to_element;
36 using ranges::views::filter;
37 using ranges::views::transform;
38 
39 namespace apa {
40 
42  : fWireReadoutGeom{&art::ServiceHandle<geo::WireReadout const>()->Get()}
43  , fCrawl{p.get<bool>("Crawl")}
44  , fUseEndP{p.get<bool>("UseEndP")}
45  , fCompareViews{p.get<bool>("CompareViews")}
46  , fNChanJumps{p.get<unsigned int>("NChanJumps")}
47  , fCloseHitsRadius{p.get<double>("CloseHitsRadius")}
48  , fMaxEndPDegRange{p.get<double>("MaxEndPDegRange")}
49  {}
50 
51  //----------------------------------------------------------
52  //----------------------------------------------------------
54  detinfo::DetectorPropertiesData const& detProp,
55  art::Handle<std::vector<recob::Hit>> ChannelHits)
56  {
57  fUeffSoFar.clear();
58  fVeffSoFar.clear();
59  fnUSoFar.clear();
60  fnVSoFar.clear();
61  fnDUSoFar.clear();
62  fnDVSoFar.clear();
63  fChannelToHits.clear();
64  fAPAToUVHits.clear();
65  fAPAToZHits.clear();
66  fAPAToHits.clear();
67  fAPAToEndPHits.clear();
68  fAPAToDHits.clear();
69  fDisambigHits.clear();
70  fChanTimeToWid.clear();
71  fHasBeenDisambiged.clear();
72 
73  std::vector<art::Ptr<recob::Hit>> ChHits;
74  art::fill_ptr_vector(ChHits, ChannelHits);
75 
76  fHasBeenDisambiged.clear();
77  unsigned int skipNoise(0);
78  // Map hits by channel/APA, initialize the disambiguation status map
79  for (size_t h = 0; h < ChHits.size(); h++) {
80  art::Ptr<recob::Hit> const& hit = ChHits[h];
81 
82  // **temporary** option to skip noise hits
83  try {
84  bt_serv->HitToXYZ(clockData, hit);
85  }
86  catch (...) {
87  skipNoise++;
88  continue;
89  }
90 
91  geo::View_t view = hit->View();
92  unsigned int apa(0), cryo(0);
93  fAPAGeo.ChannelToAPA(hit->Channel(), apa, cryo);
94  fAPAToHits[apa].push_back(hit);
95  if (view == geo::kZ) {
96  fAPAToZHits[apa].push_back(hit);
97  continue;
98  }
99  else if (view == geo::kU || view == geo::kV) {
100  std::pair<double, double> ChanTime(hit->Channel() * 1., hit->PeakTime() * 1.);
101  fHasBeenDisambiged[apa][ChanTime] = false;
102  fChannelToHits[hit->Channel()].push_back(hit);
103  fAPAToUVHits[apa].push_back(hit);
104  }
105  }
106 
107  if (skipNoise > 0)
108  mf::LogWarning("DisambigAlg")
109  << "\nSkipped " << skipNoise << " induction noise hits using the BackTrackerService.\n"
110  << "This is only to temporarily deal with the excessive amount of noise due to the bad "
111  "deconvolution.\n";
112 
113  mf::LogVerbatim("RunDisambig") << "\n~~~~~~~~~~~ Running Disambiguation ~~~~~~~~~~~\n";
114 
115  std::map<unsigned int, std::vector<art::Ptr<recob::Hit>>>::iterator APA_it;
116  for (APA_it = fAPAToUVHits.begin(); APA_it != fAPAToUVHits.end(); APA_it++) {
117  unsigned int apa = APA_it->first;
118 
119  mf::LogVerbatim("RunDisambig") << "APA " << apa << ":";
120 
121  fUeffSoFar[apa] = 0.;
122  fVeffSoFar[apa] = 0.;
123  fnUSoFar[apa] = 0;
124  fnVSoFar[apa] = 0;
125  fnDUSoFar[apa] = 0;
126  fnDVSoFar[apa] = 0;
127 
128  // Always run this...
129  TrivialDisambig(clockData, detProp, apa);
130  AssessDisambigSoFar(apa);
131  mf::LogVerbatim("RunDisambig")
132  << " Trivial Disambig --> " << fnDUSoFar[apa] << " / " << fnUSoFar[apa] << " U, "
133  << fnDVSoFar[apa] << " / " << fnVSoFar[apa] << " V";
134 
135  // ... and pick the rest with the configurations.
136  if (fCrawl) {
137  Crawl(apa);
138  AssessDisambigSoFar(apa);
139  mf::LogVerbatim("RunDisambig")
140  << " Crawl --> " << fnDUSoFar[apa] << " / " << fnUSoFar[apa] << " U, "
141  << fnDVSoFar[apa] << " / " << fnVSoFar[apa] << " V";
142  }
143 
144  if (fUseEndP) {
145  FindChanTimeEndPts(detProp, apa);
146  UseEndPts(detProp, apa); // does the crawl from inside
147  AssessDisambigSoFar(apa);
148  mf::LogVerbatim("RunDisambig")
149  << " Endpoint Crawl --> " << fnDUSoFar[apa] << " / " << fnUSoFar[apa] << " U, "
150  << fnDVSoFar[apa] << " / " << fnVSoFar[apa] << " V";
151  }
152 
153  if (fCompareViews) {
154  unsigned int nDisambig(1);
155  while (nDisambig > 0) {
156  nDisambig = CompareViews(detProp, apa);
157  Crawl(apa);
158  }
159  AssessDisambigSoFar(apa);
160  mf::LogVerbatim("RunDisambig")
161  << " Compare Views --> " << fnDUSoFar[apa] << " / " << fnUSoFar[apa] << " U, "
162  << fnDVSoFar[apa] << " / " << fnVSoFar[apa] << " V";
163  }
164 
165  // For now just buld a simple list to get from the module
166  for (size_t i = 0; i < fAPAToDHits[apa].size(); i++)
167  fDisambigHits.push_back(fAPAToDHits[apa][i]);
168 
169  } // end loop through APA
170  }
171 
172  //-------------------------------------------------
173  //-------------------------------------------------
175  geo::WireID wid,
176  unsigned int apa)
177  {
178  std::pair<double, double> ChanTime(hit->Channel() * 1., hit->PeakTime() * 1.);
179  if (fHasBeenDisambiged[apa][ChanTime]) return;
180 
181  if (!wid.isValid) {
182  mf::LogWarning("InvalidWireID") << "wid is invalid, hit not being made\n";
183  return;
184  }
185 
186  fAPAToDHits[apa].emplace_back(hit, wid);
187  fHasBeenDisambiged[apa][ChanTime] = true;
188  fChanTimeToWid[ChanTime] = wid;
189  }
190 
191  //----------------------------------------------------------
192  //----------------------------------------------------------
194  recob::Hit const& hitA,
195  recob::Hit const& hitB)
196  {
197  double AsT = hitA.PeakTimeMinusRMS();
198  double AeT = hitA.PeakTimePlusRMS();
199  double BsT = hitB.PeakTimeMinusRMS();
200  double BeT = hitB.PeakTimePlusRMS();
201 
202  if (hitA.View() == geo::kU) {
203  AsT -= detProp.TimeOffsetU();
204  AeT -= detProp.TimeOffsetU();
205  }
206  else if (hitA.View() == geo::kV) {
207  AsT -= detProp.TimeOffsetV();
208  AeT -= detProp.TimeOffsetV();
209  }
210  else if (hitA.View() == geo::kZ) {
211  AsT -= detProp.TimeOffsetZ();
212  AeT -= detProp.TimeOffsetZ();
213  }
214 
215  if (hitB.View() == geo::kU) {
216  BsT += detProp.TimeOffsetU();
217  BeT -= detProp.TimeOffsetU();
218  }
219  else if (hitB.View() == geo::kV) {
220  BsT -= detProp.TimeOffsetV();
221  BeT -= detProp.TimeOffsetV();
222  }
223  else if (hitA.View() == geo::kZ) { // FIXME: Shouldn't this be hitB, BsT, and BeT?
224  AsT -= detProp.TimeOffsetZ();
225  AeT -= detProp.TimeOffsetZ();
226  }
227 
228  return (AsT <= BsT && BsT <= AeT) || (AsT <= BeT && BeT <= AeT) || (BsT <= AsT && AsT <= BeT) ||
229  (BsT <= AeT && AeT <= BeT);
230  }
231 
232  //----------------------------------------------------------
233  //----------------------------------------------------------
235  detinfo::DetectorPropertiesData const& detProp,
236  unsigned int apa)
237  {
238  // Loop through ambiguous hits (U/V) in this APA
239  for (auto const& hitPtr : fAPAToUVHits[apa]) {
240  auto const& hit = *hitPtr;
241  raw::ChannelID_t chan = hit.Channel();
242  unsigned int peakT = hit.PeakTime();
243 
244  std::vector<geo::WireID> hitwids = fWireReadoutGeom->ChannelToWire(chan);
245  std::vector<bool> IsReasonableWid(hitwids.size(), false);
246  unsigned short nPossibleWids(0);
247  for (size_t w = 0; w < hitwids.size(); w++) {
248  geo::WireID wid = hitwids[w];
249 
250  auto const [start, end] = fWireReadoutGeom->WireEndPoints(wid);
251  unsigned int side(wid.TPC % 2), cryo(wid.Cryostat);
252  double zminPos(start.Z()), zmaxPos(end.Z());
253 
254  // get appropriate x and y with tpc center
255  unsigned int tpc =
256  2 * apa + side - cryo * geom->NTPC(); // apa number does not reset per cryo
257  auto const tpcCenter = geom->TPC(geo::TPCID(cryo, tpc)).GetCenter();
258 
259  // get channel range
260  auto Min = tpcCenter;
261  Min.SetZ(zminPos);
262  auto Max = tpcCenter;
263  Max.SetZ(zmaxPos);
264  raw::ChannelID_t ZminChan =
265  fWireReadoutGeom->NearestChannel(Min, geo::PlaneID{cryo, tpc, 2});
266  raw::ChannelID_t ZmaxChan =
267  fWireReadoutGeom->NearestChannel(Max, geo::PlaneID{cryo, tpc, 2});
268 
269  for (auto const& zhit : fAPAToZHits[apa] | transform(to_element)) {
270  raw::ChannelID_t chan = zhit.Channel();
271  if (chan <= ZminChan || ZmaxChan <= chan) continue;
272 
273  if (HitsOverlapInTime(detProp, hit, zhit)) {
274  IsReasonableWid[w] = true;
275  nPossibleWids++;
276  break;
277  }
278  }
279 
280  } // end hit chan-wid loop
281 
282  if (nPossibleWids == 0) {
283  std::vector<double> xyz;
284  try {
285  xyz = bt_serv->HitToXYZ(clockData, hit);
286  } // TEMPORARY
287  catch (...) {
288  continue;
289  }
291  mf::LogWarning("UniqueTimeSeg")
292  << "U/V hit inconsistent with Z info; peak time is " << peakT << " in APA " << apa
293  << " on channel " << hit.Channel();
294  }
295  else if (nPossibleWids == 1) {
296  for (size_t d = 0; d < hitwids.size(); d++)
297  if (IsReasonableWid[d]) MakeDisambigHit(hitPtr, hitwids[d], apa);
298  }
299  else if (nPossibleWids == 2) {
301  }
302 
303  } // end ambig hits loop
304  }
305 
306  //----------------------------------------------------------
307  //----------------------------------------------------------
308  unsigned int DisambigAlg::MakeCloseHits(int ext, geo::WireID Dwid, double Dmin, double Dmax)
309  {
310  // Function to look, on a channel *ext* channels away from a
311  // disambiguated hit channel, for hits with time windows touching
312  // range *Dmin to Dmax*. If found, make such a hit to have a
313  // wireID adjacent to supplied *wid*. Returns number of NEW hits
314  // made.
315 
317  geo::View_t view = fWireReadoutGeom->View(Dchan);
318  if (view == geo::kZ)
319  throw cet::exception("MakeCloseHits") << "Function not meant for non-wrapped channels.\n";
320 
321  // Account for wrapping
322  raw::ChannelID_t firstChan = fAPAGeo.FirstChannelInView(view, Dchan);
323  unsigned int ChanPerView = fAPAGeo.ChannelsInView(view);
324  int tempchan = Dchan + ext; // need sign for the set of channels starting with channel 0
325  if (tempchan < (int)firstChan) tempchan += ChanPerView;
326  if (tempchan > (int)(firstChan + ChanPerView - 1)) tempchan -= ChanPerView;
327  raw::ChannelID_t chan = (raw::ChannelID_t)(tempchan);
328 
329  // There may just be no hits
330  if (fChannelToHits.count(chan) == 0) return 0;
331 
332  // There are close channel hits, so for each
333  unsigned int apa(0), cryo(0);
334  fAPAGeo.ChannelToAPA(chan, apa, cryo);
335  unsigned int MakeCount(0);
336  for (size_t i = 0; i < fChannelToHits[chan].size(); i++) {
337  art::Ptr<recob::Hit> closeHit = fChannelToHits[chan][i];
338  double st = closeHit->PeakTimeMinusRMS();
339  double et = closeHit->PeakTimePlusRMS();
340  std::vector<geo::WireID> wids = fWireReadoutGeom->ChannelToWire(chan);
341 
342  if (!(Dmin <= st && st <= Dmax) && !(Dmin <= et && et <= Dmax)) continue;
343 
344  // Found hit with window overlapping given range,
345  // now find the only reasonable wireID.
346  for (size_t w = 0; w < wids.size(); w++) {
347  if (wids[w].TPC != Dwid.TPC) continue;
348  if ((int)(wids[w].Wire) - (int)(Dwid.Wire) != ext) continue;
349 
350  // In this case, we have a unique wireID.
351  // Check to see if it has already been made - if so, do not incriment count
352  std::pair<double, double> ChanTime(closeHit->Channel() * 1., closeHit->PeakTime() * 1.);
353  if (!fHasBeenDisambiged[apa][ChanTime]) {
354  MakeDisambigHit(closeHit, wids[w], apa);
355  MakeCount++;
356  }
357  break;
358  } // end find right wireID
359 
360  } // end loop through all hits on chan
361 
362  return MakeCount;
363  }
364 
365  //----------------------------------------------------------
366  //----------------------------------------------------------
367  void DisambigAlg::Crawl(unsigned int apa)
368  {
369 
370  std::vector<art::Ptr<recob::Hit>> hits = fAPAToUVHits[apa];
371 
372  // repeat this method until stable
373  unsigned int nExtended(1);
374  while (nExtended > 0) {
375  nExtended = 0;
376 
377  // Look for any disambiguated hit ...
378  for (size_t h = 0; h < hits.size(); h++) {
379  std::pair<double, double> ChanTime(hits[h]->Channel() * 1., hits[h]->PeakTime() * 1.);
380  if (!fHasBeenDisambiged[apa][ChanTime]) continue;
381  double stD = hits[h]->PeakTimePlusRMS(-1.);
382  double etD = hits[h]->PeakTimePlusRMS(+1.);
383  double hitWindow = etD - stD;
384  geo::WireID Dwid = fChanTimeToWid[ChanTime];
385 
386  // ... and if any neighboring-channel hits are close enough in time,
387  // extend the disambiguation to the neighboring wire.
388  unsigned int extensions = 0;
389  for (unsigned int ext = 1; ext < fNChanJumps + 1; ext++) {
391  unsigned int N(0);
392  double timeExt = hitWindow * ext;
393  N += MakeCloseHits((int)(-ext), Dwid, stD - 5 - timeExt, etD + 5 + timeExt);
394  N += MakeCloseHits((int)(ext), Dwid, stD - 5 - timeExt, etD + 5 + timeExt);
395  extensions += N;
396  }
397  nExtended += extensions;
398 
399  } // end UV hits loop
400 
401  } // end while still disambiguating***
402 
403  // *** nested while loops allow disambigauation to hop the channel-wrap-boundary
404  }
405 
406  //----------------------------------------------------------
407  //----------------------------------------------------------
409  unsigned int apa)
410  {
413 
414  double pi = 3.14159265;
415  double fMaxEndPRadRange = fMaxEndPDegRange / 180. * (2 * pi);
416 
417  for (size_t h = 0; h < fAPAToHits[apa].size(); h++) {
418  art::Ptr<recob::Hit> centhit = fAPAToHits[apa][h];
419  geo::View_t const view = centhit->View();
420  unsigned int plane = 0;
421  if (view == geo::kV) { plane = 1; }
422  else if (view == geo::kZ)
423  plane = 2;
424  std::vector<double> ChanTimeCenter(2, 0.);
425  unsigned int relchan = centhit->Channel() - fAPAGeo.FirstChannelInView(centhit->Channel());
426  auto const wire_pitch = fWireReadoutGeom->Plane({0, 0, view}).WirePitch();
427  ChanTimeCenter[0] = relchan * wire_pitch;
428  ChanTimeCenter[1] = detProp.ConvertTicksToX(centhit->PeakTime(),
429  plane,
430  apa * 2, // tpc doesnt matter
431  centhit->WireID().Cryostat);
432  //std::vector< art::Ptr<recob::Hit> > CloseHits;
433  std::vector<std::vector<double>> CloseHitsChanTime;
434  std::vector<double> FurthestCloseChanTime(2, 0.); //double maxDist = 0.;
435  std::vector<double> ClosestChanTime(2, 0.);
436  double minDist = fCloseHitsRadius + 1.;
437  double ChanDistRange = fAPAGeo.ChannelsInView(view) * wire_pitch;
438 
439  for (size_t c = 0; c < fAPAToHits[apa].size(); c++) {
440  art::Ptr<recob::Hit> closehit = fAPAToHits[apa][c];
441  if (view != closehit->View()) continue;
442  if (view == geo::kZ && centhit->WireID().TPC != closehit->WireID().TPC) continue;
443  unsigned int plane = 0;
444  if (view == geo::kV) { plane = 1; }
445  else if (view == geo::kZ)
446  plane = 2;
447  std::vector<double> ChanTimeClose(2, 0.);
448  unsigned int relchanclose =
449  closehit->Channel() - fAPAGeo.FirstChannelInView(closehit->Channel());
450  ChanTimeClose[0] = relchanclose * wire_pitch;
451  ChanTimeClose[1] = detProp.ConvertTicksToX(closehit->PeakTime(),
452  plane,
453  apa * 2, // tpc doesnt matter
454  closehit->WireID().Cryostat);
455  if (ChanTimeClose == ChanTimeCenter) continue; // move on if the same one
456 
457  double ChanDist = ChanTimeClose[0] - ChanTimeCenter[0];
458  if (ChanDist > ChanDistRange / 2) ChanDist = ChanDistRange - ChanDist;
459 
460  double distance = std::hypot(ChanDist, ChanTimeClose[1] - ChanTimeCenter[1]);
461 
462  if (distance <= fCloseHitsRadius) CloseHitsChanTime.push_back(ChanTimeClose);
463 
464  if (distance < minDist) {
465  ClosestChanTime = ChanTimeClose;
466  minDist = distance;
467  }
468 
469  } // end close-by hit loop
470 
471  if (CloseHitsChanTime.size() < 5) continue; // quick fix, to-be improved
472 
473  double minRad(2 * pi + 1.), maxRad(0.);
474  bool CloseToNegPi(false), CloseToPosPi(false);
475  for (size_t i = 0; i < CloseHitsChanTime.size(); i++) {
476  std::vector<double> ThisChanTime(CloseHitsChanTime[i]);
477  double ChanDist = ThisChanTime[0] - ChanTimeCenter[0];
478  if (ChanDist > ChanDistRange / 2) ChanDist = ChanDistRange - ChanDist;
479  double hitrad = std::atan2(ThisChanTime[1] - ChanTimeCenter[1], ChanDist);
480  if (hitrad > maxRad) maxRad = hitrad;
481  if (hitrad < minRad) minRad = hitrad;
482  if (hitrad + fMaxEndPRadRange > pi)
483  CloseToPosPi = true;
484  else if (hitrad - fMaxEndPRadRange < -pi)
485  CloseToNegPi = true;
486  }
487 
488  // activity at this boundary automatically kills the test, move boundary and redo
489  if (CloseToPosPi && CloseToNegPi) {
490  for (size_t i = 0; i < CloseHitsChanTime.size(); i++) {
491  std::vector<double> ThisChanTime(CloseHitsChanTime[i]);
492  double ChanDist = ThisChanTime[0] - ChanTimeCenter[0];
493  if (ChanDist > ChanDistRange / 2) ChanDist = ChanDistRange - ChanDist;
494  double hitrad = std::atan2(ThisChanTime[1] - ChanTimeCenter[1], ChanDist);
495  if (hitrad > 0) hitrad = pi - hitrad; // reflec across pi/2 line
496  if (hitrad < 0) hitrad = -pi - hitrad;
497  if (hitrad > maxRad) maxRad = hitrad;
498  if (hitrad < minRad) minRad = hitrad;
499  }
500  }
501 
502  if (maxRad - minRad < fMaxEndPRadRange) fAPAToEndPHits[apa].push_back(centhit);
503 
504  } // end UV hit loop
505 
506  if (fAPAToEndPHits[apa].size() == 0) return 0;
507  mf::LogVerbatim("FindChanTimeEndPts") << " Found " << fAPAToEndPHits[apa].size()
508  << " endpoint hits in apa " << apa << std::endl;
509  for (size_t ep = 0; ep < fAPAToEndPHits[apa].size(); ep++) {
510  art::Ptr<recob::Hit> epHit = fAPAToEndPHits[apa][ep];
511  mf::LogVerbatim("FindChanTimeEndPts") << " endP on channel " << epHit->Channel()
512  << " at time " << epHit->PeakTime() << std::endl;
513  }
514 
515  return fAPAToEndPHits[apa].size();
516  }
517 
518  //----------------------------------------------------------
519  //----------------------------------------------------------
520  void DisambigAlg::UseEndPts(detinfo::DetectorPropertiesData const& detProp, unsigned int apa)
521  {
523 
524  if (fAPAToEndPHits[apa].size() == 0) {
525  mf::LogVerbatim("UseEndPts") << " APA " << apa << " has no endpoints.";
526  return;
527  }
528  std::vector<art::Ptr<recob::Hit>> const& endPts = fAPAToEndPHits[apa];
529 
530  std::vector<std::vector<art::Ptr<recob::Hit>>> EndPMatch;
531  unsigned short nZendPts(0);
532 
533  auto on_z_plane = [](art::Ptr<recob::Hit> const& hit) { return hit->View() == geo::kZ; };
534  auto not_on_z_plane = [](art::Ptr<recob::Hit> const& hit) { return hit->View() != geo::kZ; };
535  for (auto const& ZHitPtr : endPts | filter(on_z_plane)) {
536  auto const& ZHit = *ZHitPtr;
537  art::Ptr<recob::Hit> Uhit = ZHitPtr;
538  art::Ptr<recob::Hit> Vhit = ZHitPtr;
539  unsigned short Umatch(0), Vmatch(0);
540  ++nZendPts;
541 
542  // look for U and V hits overlapping in time
543  for (auto const& hitPtr : endPts | filter(not_on_z_plane)) {
544  auto const& hit = *hitPtr;
545  if (not HitsOverlapInTime(detProp, ZHit, hit)) continue;
546 
547  if (hit.View() == geo::kU) {
548  Uhit = hitPtr;
549  Umatch++;
550  }
551  else if (hit.View() == geo::kV) {
552  Vhit = hitPtr;
553  Vmatch++;
554  }
555  }
556 
557  unsigned int tpc(ZHit.WireID().TPC), cryo(ZHit.WireID().Cryostat);
558  auto const tpcCenter = geom->TPC(geo::TPCID{cryo, tpc}).GetCenter();
559 
560  if (Umatch == 1 && Vmatch == 1) {
561 
562  std::vector<double> yzEndPt =
563  fAPAGeo.ThreeChanPos(Uhit->Channel(), Vhit->Channel(), ZHit.Channel());
564 
565  geo::Point_t const intersect{tpcCenter.X(), yzEndPt[0], yzEndPt[1]};
566  geo::WireID Uwid = fAPAGeo.NearestWireIDOnChan(intersect, Uhit->Channel(), {cryo, tpc, 0});
567  geo::WireID Vwid = fAPAGeo.NearestWireIDOnChan(intersect, Vhit->Channel(), {cryo, tpc, 1});
568  MakeDisambigHit(Uhit, Uwid, apa);
569  MakeDisambigHit(Vhit, Vwid, apa);
570  }
571  else if (Umatch == 1 && Vmatch != 1) {
572 
573  std::vector<geo::WireIDIntersection> widIntersects;
574  fAPAGeo.APAChannelsIntersect(Uhit->Channel(), ZHit.Channel(), widIntersects);
575  if (widIntersects.size() == 0)
576  continue;
577  else if (widIntersects.size() == 1) {
578  geo::Point_t const intersect{tpcCenter.X(), widIntersects[0].y, widIntersects[0].z};
579  geo::WireID Uwid =
580  fAPAGeo.NearestWireIDOnChan(intersect, Uhit->Channel(), {cryo, tpc, 0});
581  MakeDisambigHit(Uhit, Uwid, apa);
582  }
583  else {
584  for (size_t i = 0; i < widIntersects.size(); i++) {
585  // compare to V hit times, see if only one makes sense
586  }
587  }
588  }
589  else if (Umatch == 1 && Vmatch != 1) {
590 
591  std::vector<geo::WireIDIntersection> widIntersects;
592  fAPAGeo.APAChannelsIntersect(Vhit->Channel(), ZHit.Channel(), widIntersects);
593  if (widIntersects.size() == 0)
594  continue;
595  else if (widIntersects.size() == 1) {
596  geo::Point_t const intersect{tpcCenter.X(), widIntersects[0].y, widIntersects[0].z};
597  geo::WireID Vwid =
598  fAPAGeo.NearestWireIDOnChan(intersect, Vhit->Channel(), {cryo, tpc, 0});
599  MakeDisambigHit(Vhit, Vwid, apa);
600  }
601  }
602  }
603 
604  if (nZendPts == 0 && endPts.size() == 2 && HitsOverlapInTime(detProp, *endPts[0], *endPts[1])) {
605  std::vector<geo::WireIDIntersection> widIntersects;
606  fAPAGeo.APAChannelsIntersect(endPts[0]->Channel(), endPts[1]->Channel(), widIntersects);
607  if (widIntersects.size() == 1) {
608  unsigned int cryo = endPts[0]->WireID().Cryostat;
609  unsigned int tpc = widIntersects[0].TPC;
610  auto const tpcCenter = geom->TPC(geo::TPCID(cryo, tpc)).GetCenter();
611  geo::Point_t const intersect{tpcCenter.X(), widIntersects[0].y, widIntersects[0].z};
612  unsigned int plane0(0), plane1(0);
613  if (endPts[0]->View() == geo::kV) plane0 = 1;
614  if (endPts[1]->View() == geo::kV) plane1 = 1;
615  geo::WireID wid0 =
616  fAPAGeo.NearestWireIDOnChan(intersect, endPts[0]->Channel(), {cryo, tpc, plane0});
617  MakeDisambigHit(endPts[0], wid0, apa);
618  geo::WireID wid1 =
619  fAPAGeo.NearestWireIDOnChan(intersect, endPts[1]->Channel(), {cryo, tpc, plane1});
620  MakeDisambigHit(endPts[1], wid1, apa);
621  }
622  }
623 
624  Crawl(apa);
625  }
626 
627  //----------------------------------------------------------
628  //----------------------------------------------------------
630  {
631  unsigned int nU(0), nV(0);
632  for (size_t h = 0; h < fAPAToUVHits[apa].size(); h++) {
634  if (hit->View() == geo::kU)
635  nU++;
636  else if (hit->View() == geo::kV)
637  nV++;
638  }
639 
640  unsigned int nDU(0), nDV(0);
641  for (size_t h = 0; h < fAPAToDHits[apa].size(); h++) {
642  art::Ptr<recob::Hit> hit = fAPAToDHits[apa][h].first;
643  if (hit->View() == geo::kU)
644  nDU++;
645  else if (hit->View() == geo::kV)
646  nDV++;
647  }
648 
649  fUeffSoFar[apa] = (nDU * 1.) / (nU * 1.);
650  fVeffSoFar[apa] = (nDV * 1.) / (nV * 1.);
651  fnUSoFar[apa] = nU;
652  fnVSoFar[apa] = nV;
653  fnDUSoFar[apa] = nDU;
654  fnDVSoFar[apa] = nDV;
655  }
656 
657  //----------------------------------------------------------
658  //----------------------------------------------------------
660  unsigned int apa)
661  {
662  unsigned int nDisambiguations(0);
663 
664  // loop through all hits that are still ambiguous
665  for (auto const& ambighitPtr : fAPAToUVHits[apa]) {
666  auto const& ambighit = *ambighitPtr;
667  raw::ChannelID_t ambigchan = ambighit.Channel();
668  std::pair<double, double> ambigChanTime(ambigchan * 1., ambighit.PeakTime());
669  if (fHasBeenDisambiged[apa][ambigChanTime]) continue;
670  geo::View_t view = ambighit.View();
671  std::vector<geo::WireID> ambigwids = fWireReadoutGeom->ChannelToWire(ambigchan);
672  std::vector<unsigned int> widDcounts(ambigwids.size(), 0);
673  std::vector<unsigned int> widAcounts(ambigwids.size(), 0);
674 
675  // loop through hits in the other view which are close in time
676  for (auto const& hit : fAPAToUVHits[apa] | transform(to_element)) {
677  if (hit.View() == view || !HitsOverlapInTime(detProp, ambighit, hit)) continue;
678 
679  // An other-view-hit overlaps in time, see what
680  // wids of the ambiguous hit's channels it overlaps
681  raw::ChannelID_t chan = hit.Channel();
682  std::vector<geo::WireID> wids = fWireReadoutGeom->ChannelToWire(chan);
683  std::pair<double, double> ChanTime(chan * 1., hit.PeakTime());
684  if (fHasBeenDisambiged[apa][ChanTime]) {
685  for (size_t a = 0; a < ambigwids.size(); a++)
686  if (ambigwids[a].TPC == fChanTimeToWid[ChanTime].TPC &&
687  fWireReadoutGeom->WireIDsIntersect(ambigwids[a], fChanTimeToWid[ChanTime]))
688  widDcounts[a]++;
689  }
690  else {
691  // still might be able to glean disambiguation
692  // from the ambiguous hits at this time
693  for (size_t a = 0; a < ambigwids.size(); a++)
694  for (size_t w = 0; w < wids.size(); w++)
695  if (ambigwids[a].TPC == wids[w].TPC &&
696  fWireReadoutGeom->WireIDsIntersect(ambigwids[a], wids[w]))
697  widAcounts[a]++;
698  }
699  } // end loop through close-time hits
700 
701  // For now, just make a hit if either ambig or disambig hits
702  // unanimously intersect a single wireID
703  unsigned int Dcount(0), Acount(0);
704  for (size_t d = 0; d < widDcounts.size(); d++)
705  Dcount += widDcounts[d];
706  for (size_t a = 0; a < widAcounts.size(); a++)
707  Acount += widAcounts[a];
708  for (size_t d = 0; d < widDcounts.size(); d++) {
709  if (Dcount == widDcounts[d] && Dcount > 0 && Acount == 0) {
710  MakeDisambigHit(ambighitPtr, ambigwids[d], apa);
711  nDisambiguations++;
712  }
713  }
714  } // end loop through still ambiguous hits
715 
716  return nDisambiguations;
717  }
718 
719 } //end namespace apa
std::map< unsigned int, std::map< std::pair< double, double >, bool > > fHasBeenDisambiged
Convenient way to keep track of disambiguation so far.
Definition: DisambigAlg.h:84
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::map< raw::ChannelID_t, std::vector< art::Ptr< recob::Hit > > > fChannelToHits
Definition: DisambigAlg.h:73
constexpr to_element_t to_element
Definition: ToElement.h:25
Encapsulate the construction of a single cyostat .
double fCloseHitsRadius
Definition: DisambigAlg.h:103
bool HitsOverlapInTime(detinfo::DetectorPropertiesData const &detProp, recob::Hit const &hitA, recob::Hit const &hitB)
std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > fAPAToZHits
Definition: DisambigAlg.h:74
std::map< unsigned int, unsigned int > fnDVSoFar
Definition: DisambigAlg.h:59
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Planes which measure V.
Definition: geo_types.h:132
unsigned int NTPC(CryostatID const &cryoid=details::cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:416
apa::APAGeometryAlg fAPAGeo
Definition: DisambigAlg.h:66
Declaration of signal hit object.
std::map< unsigned int, unsigned int > fnVSoFar
Definition: DisambigAlg.h:57
The data type to uniquely identify a Plane.
Definition: geo_types.h:364
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:194
uint32_t FirstChannelInView(geo::View_t geoview, unsigned int apa, unsigned int cryo) const
void Crawl(unsigned int apa)
Extend what we disambiguation we do have in apa.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:195
Planes which measure Z direction.
Definition: geo_types.h:134
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
geo::View_t View() const
View for the plane of the hit.
Definition: Hit.h:286
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:430
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:290
void RunDisambig(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, art::Handle< std::vector< recob::Hit >> GausHits)
Definition: DisambigAlg.cxx:53
std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > fAPAToEndPHits
\ todo: Channel/APA to hits can be done in a unified way
Definition: DisambigAlg.h:77
std::vector< double > HitToXYZ(detinfo::DetectorClocksData const &clockData, const recob::Hit &hit) const
bool APAChannelsIntersect(uint32_t chan1, uint32_t chan2, std::vector< geo::WireIDIntersection > &IntersectVector) const
If the channels intersect, get all intersections.
art::ServiceHandle< geo::Geometry const > geom
Definition: DisambigAlg.h:67
Framework includes.
void TrivialDisambig(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, unsigned int apa)
Make the easiest and safest disambiguations in apa.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
Planes which measure U.
Definition: geo_types.h:131
Point_t GetCenter() const
Returns the center of the TPC volume in world coordinates [cm].
Definition: TPCGeo.h:132
unsigned int ChannelToAPA(uint32_t chan) const
Get number of the APA containing the given channel.
void hits()
Definition: readHits.C:15
void AssessDisambigSoFar(unsigned int apa)
See how much disambiguation has been done in this apa so far.
std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > fAPAToUVHits
Definition: DisambigAlg.h:74
std::map< unsigned int, std::vector< art::Ptr< recob::Hit > > > fAPAToHits
Definition: DisambigAlg.h:75
DisambigAlg(fhicl::ParameterSet const &pset)
Definition: DisambigAlg.cxx:41
Float_t d
Definition: plot.C:235
unsigned int ChannelsInView(geo::View_t geoview) const
virtual raw::ChannelID_t PlaneWireToChannel(WireID const &wireID) const =0
Returns the channel ID a wire is connected to.
std::map< unsigned int, double > fVeffSoFar
Definition: DisambigAlg.h:55
void MakeDisambigHit(art::Ptr< recob::Hit > const &hit, geo::WireID, unsigned int apa)
Makes a disambiguated hit while keeping track of what has already been disambiguated.
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
unsigned int fNChanJumps
Number of channels the crawl can jump over.
Definition: DisambigAlg.h:102
void UseEndPts(detinfo::DetectorPropertiesData const &detProp, unsigned int apa)
float PeakTimeMinusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Definition: Hit.h:300
Detector simulation of raw signals on wires.
raw::ChannelID_t NearestChannel(Point_t const &worldLoc, PlaneID const &planeid) const
Returns the ID of the channel nearest to the specified position.
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, Point_t &intersection) const
Computes the intersection between two wires.
double ConvertTicksToX(double ticks, int p, int t, int c) const
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
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:226
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
Contains all timing reference information for the detector.
View_t View(raw::ChannelID_t const channel) const
Returns the view (wire orientation) on the specified TPC channel.
geo::WireID NearestWireIDOnChan(geo::Point_t const &WorldLoc, uint32_t chan, geo::PlaneID const &planeID) const
bool fCrawl
\ todo: Write function that compares hits more detailedly
Definition: DisambigAlg.h:99
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
art::ServiceHandle< cheat::BackTrackerService const > bt_serv
Definition: DisambigAlg.h:69
std::vector< std::pair< art::Ptr< recob::Hit >, geo::WireID > > fDisambigHits
The final list of hits to pass back to be made.
Definition: DisambigAlg.h:61
double fMaxEndPDegRange
Definition: DisambigAlg.h:105
unsigned int CompareViews(detinfo::DetectorPropertiesData const &detProp, unsigned int apa)
Compare U and V to see if one says something about the other.
std::map< unsigned int, double > fUeffSoFar
Definition: DisambigAlg.h:54
geo::WireReadoutGeom const * fWireReadoutGeom
Definition: DisambigAlg.h:70
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
float PeakTimePlusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Definition: Hit.h:295
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:448
std::vector< double > ThreeChanPos(uint32_t u, uint32_t v, uint32_t z) const
Find the center of the 3 intersections, choose best if multiple.
std::map< unsigned int, unsigned int > fnUSoFar
Definition: DisambigAlg.h:56
Float_t w
Definition: plot.C:20
unsigned int FindChanTimeEndPts(detinfo::DetectorPropertiesData const &detProp, unsigned int apa)
Basic endpoint-hit finder per apa.
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
Definition: Hit.h:278
void WireEndPoints(WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
std::map< unsigned int, std::vector< std::pair< art::Ptr< recob::Hit >, geo::WireID > > > fAPAToDHits
Hold the disambiguations per APA.
Definition: DisambigAlg.h:78
unsigned int MakeCloseHits(int ext, geo::WireID wid, double Dmin, double Dmax)
Having disambiguated a time range on a wireID, extend to neighboring channels.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Encapsulate the construction of a single detector plane .
std::map< std::pair< double, double >, geo::WireID > fChanTimeToWid
If a hit is disambiguated, map its chan and peak time to the chosen wireID.
Definition: DisambigAlg.h:82
std::map< unsigned int, unsigned int > fnDUSoFar
Definition: DisambigAlg.h:58