LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
SpacePointAlg.cxx
Go to the documentation of this file.
1 
11 // LArSoft includes
12 //\todo Remove include of BackTrackerService.h once this algorithm is stripped of test for MC
27 
28 // Framework includes
30 #include "cetlib_except/exception.h"
32 
33 // C++ STL includes
34 #include <algorithm>
35 #include <cassert>
36 #include <cmath>
37 #include <iostream>
38 #include <map>
39 
40 using lar::to_element;
41 
42 //----------------------------------------------------------------------
43 // Constructor.
44 //
45 namespace trkf {
46 
48  : fMaxDT{pset.get<double>("MaxDT")}
49  , fMaxS{pset.get<double>("MaxS")}
50  , fMinViews{pset.get<int>("MinViews")}
51  , fEnableU{pset.get<bool>("EnableU")}
52  , fEnableV{pset.get<bool>("EnableV")}
53  , fEnableW{pset.get<bool>("EnableW")}
54  , fFilter{pset.get<bool>("Filter")}
55  , fMerge{pset.get<bool>("Merge")}
56  , fPreferColl{pset.get<bool>("PreferColl")}
57  , fTickOffsetU{pset.get<double>("TickOffsetU", 0.)}
58  , fTickOffsetV{pset.get<double>("TickOffsetV", 0.)}
59  , fTickOffsetW{pset.get<double>("TickOffsetW", 0.)}
60  {
61  // Only allow one of fFilter and fMerge to be true.
62 
63  if (fFilter && fMerge)
64  throw cet::exception("SpacePointAlg") << "Filter and Merge flags are both true.\n";
65 
66  // Report.
67 
68  std::cout << "SpacePointAlg configured with the following parameters:\n"
69  << " MaxDT = " << fMaxDT << "\n"
70  << " MaxS = " << fMaxS << "\n"
71  << " MinViews = " << fMinViews << "\n"
72  << " EnableU = " << fEnableU << "\n"
73  << " EnableV = " << fEnableV << "\n"
74  << " EnableW = " << fEnableW << "\n"
75  << " Filter = " << fFilter << "\n"
76  << " Merge = " << fMerge << "\n"
77  << " PreferColl = " << fPreferColl << "\n"
78  << " TickOffsetU = " << fTickOffsetU << "\n"
79  << " TickOffsetV = " << fTickOffsetV << "\n"
80  << " TickOffsetW = " << fTickOffsetW << std::endl;
81  }
82 
83  //----------------------------------------------------------------------
84  // Print geometry and properties constants.
85  //
87  {
88  // Generate info report on first call only.
89 
90  static bool first = true;
91  bool report = first;
92  first = false;
93 
94  // Get services.
95 
97  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
98 
99  // Calculate and print geometry information.
100 
101  if (report) mf::LogInfo("SpacePointAlg") << "Updating geometry constants.\n";
102 
103  for (auto const& plane : wireReadoutGeom.Iterate<geo::PlaneGeo>()) {
104 
105  // Fill view-dependent quantities.
106 
107  geo::View_t view = plane.View();
108  std::string viewname = "?";
109  if (view == geo::kU) { viewname = "U"; }
110  else if (view == geo::kV) {
111  viewname = "V";
112  }
113  else if (view == geo::kZ) {
114  viewname = "Z";
115  }
116  else
117  throw cet::exception("SpacePointAlg") << "Bad view = " << view << "\n";
118 
119  std::string sigtypename = "?";
120  geo::SigType_t sigtype = wireReadoutGeom.SignalType(plane.ID());
121  if (sigtype == geo::kInduction)
122  sigtypename = "Induction";
123  else if (sigtype == geo::kCollection)
124  sigtypename = "Collection";
125  else
126  throw cet::exception("SpacePointAlg") << "Bad signal type = " << sigtype << "\n";
127 
128  std::string orientname = "?";
129  geo::Orient_t orient = plane.Orientation();
130  if (orient == geo::kVertical)
131  orientname = "Vertical";
132  else if (orient == geo::kHorizontal)
133  orientname = "Horizontal";
134  else
135  throw cet::exception("SpacePointAlg") << "Bad orientation = " << orient << "\n";
136 
137  if (report) {
138  auto const xyz = plane.GetCenter();
139  auto const& planeid = plane.ID();
140  mf::LogInfo("SpacePointAlg")
141  << '\n'
142  << planeid << '\n'
143  << " View: " << viewname << "\n"
144  << " SignalType: " << sigtypename << "\n"
145  << " Orientation: " << orientname << "\n"
146  << " Plane location: " << xyz.X() << "\n"
147  << " Plane pitch: " << wireReadoutGeom.Plane0Pitch(planeid.parentID(), planeid.Plane)
148  << "\n"
149  << " Wire angle: " << plane.Wire(0).ThetaZ() << "\n"
150  << " Wire pitch: " << wireReadoutGeom.Plane(planeid).WirePitch() << "\n"
151  << " Time offset: " << detProp.GetXTicksOffset(planeid) << "\n";
152  }
153 
154  if (orient != geo::kVertical)
155  throw cet::exception("SpacePointAlg") << "Horizontal wire geometry not implemented.\n";
156  } // end loop over planes
157  }
158 
159  //----------------------------------------------------------------------
160  // Get corrected time for the specified hit.
162  const recob::Hit& hit) const
163  {
164  // Correct time for trigger offset and plane-dependent time offsets.
165 
166  double t = hit.PeakTime() - detProp.GetXTicksOffset(hit.WireID());
167  if (hit.View() == geo::kU)
168  t -= fTickOffsetU;
169  else if (hit.View() == geo::kV)
170  t -= fTickOffsetV;
171  else if (hit.View() == geo::kW)
172  t -= fTickOffsetW;
173 
174  return t;
175  }
176 
177  //----------------------------------------------------------------------
178  // Spatial separation of hits (zero if two or fewer).
180  {
181  // Trivial case - fewer than three hits.
182 
183  if (hits.size() < 3) return 0.;
184 
185  // Error case - more than three hits.
186 
187  if (hits.size() > 3) {
188  mf::LogError("SpacePointAlg") << "Method separation called with more than three htis.";
189  return 0.;
190  }
191 
192  // Got exactly three hits.
193 
194  // Calculate angles and distance of each hit from origin.
195 
196  double dist[3] = {0., 0., 0.};
197  double sinth[3] = {0., 0., 0.};
198  double costh[3] = {0., 0., 0.};
199  unsigned int cstats[3];
200  unsigned int tpcs[3];
201  unsigned int planes[3];
202 
203  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
204 
205  for (int i = 0; i < 3; ++i) {
206 
207  // Get tpc, plane, wire.
208 
209  const recob::Hit& hit = *(hits[i]);
210  const geo::WireGeo& wgeom = wireReadoutGeom.Wire(hit.WireID());
211  cstats[i] = hit.WireID().Cryostat;
212  tpcs[i] = hit.WireID().TPC;
213  planes[i] = hit.WireID().Plane;
214 
215  // Check tpc and plane errors.
216 
217  for (int j = 0; j < i; ++j) {
218 
219  if (cstats[j] != hit.WireID().Cryostat) {
220  mf::LogError("SpacePointAlg")
221  << "Method separation called with hits from multiple cryostats..";
222  return 0.;
223  }
224 
225  if (tpcs[j] != hit.WireID().TPC) {
226  mf::LogError("SpacePointAlg")
227  << "Method separation called with hits from multiple tpcs..";
228  return 0.;
229  }
230 
231  if (planes[j] == hit.WireID().Plane) {
232  mf::LogError("SpacePointAlg")
233  << "Method separation called with hits from the same plane..";
234  return 0.;
235  }
236  }
237 
238  // Get angles and distance of wire.
239 
240  double const hl = wgeom.HalfL();
241  auto const xyz = wgeom.GetCenter();
242  auto const xyz1 = wgeom.GetEnd();
243  double s = (xyz1.Y() - xyz.Y()) / hl;
244  double c = (xyz1.Z() - xyz.Z()) / hl;
245  sinth[hit.WireID().Plane] = s;
246  costh[hit.WireID().Plane] = c;
247  dist[hit.WireID().Plane] = xyz.Z() * s - xyz.Y() * c;
248  }
249 
250  double S = ((sinth[1] * costh[2] - costh[1] * sinth[2]) * dist[0] +
251  (sinth[2] * costh[0] - costh[2] * sinth[0]) * dist[1] +
252  (sinth[0] * costh[1] - costh[0] * sinth[1]) * dist[2]);
253  return S;
254  }
255 
256  //----------------------------------------------------------------------
257  // Check hits for compatibility.
258  // Check hits pairwise for different views and maximum time difference.
259  // Check three hits for spatial compatibility.
262  bool useMC) const
263  {
265 
266  int nhits = hits.size();
267 
268  // Fewer than two or more than three hits can never be compatible.
269 
270  bool result = nhits >= 2 && nhits <= 3;
271  bool mc_ok = true;
272  unsigned int tpc = 0;
273  unsigned int cstat = 0;
274 
275  if (result) {
276 
277  // First do pairwise tests.
278  // Do double loop over hits.
279 
280  for (int ihit1 = 0; result && ihit1 < nhits - 1; ++ihit1) {
281  const recob::Hit& hit1 = *(hits[ihit1]);
282  geo::WireID hit1WireID = hit1.WireID();
283  geo::View_t view1 = hit1.View();
284 
285  double t1 = hit1.PeakTime() -
286  detProp.GetXTicksOffset(hit1WireID.Plane, hit1WireID.TPC, hit1WireID.Cryostat);
287 
288  // If using mc information, get a collection of track ids for hit 1.
289  // If not using mc information, this section of code will trigger the
290  // insertion of a single invalid HitMCInfo object into fHitMCMap.
291 
292  const HitMCInfo& mcinfo1 = fHitMCMap[(useMC ? &hit1 : 0)];
293  const std::vector<int>& tid1 = mcinfo1.trackIDs;
294  bool only_neg1 = tid1.size() > 0 && tid1.back() < 0;
295 
296  // Loop over second hit.
297 
298  for (int ihit2 = ihit1 + 1; result && ihit2 < nhits; ++ihit2) {
299  const recob::Hit& hit2 = *(hits[ihit2]);
300  geo::WireID hit2WireID = hit2.WireID();
301  geo::View_t view2 = hit2.View();
302 
303  // Test for same tpc and different views.
304 
305  result = result && hit1WireID.TPC == hit2WireID.TPC && view1 != view2 &&
306  hit1WireID.Cryostat == hit2WireID.Cryostat;
307  if (result) {
308 
309  // Remember which tpc and cryostat we are in.
310 
311  tpc = hit1WireID.TPC;
312  cstat = hit1WireID.Cryostat;
313 
314  double t2 = hit2.PeakTime() - detProp.GetXTicksOffset(
315  hit2WireID.Plane, hit2WireID.TPC, hit2WireID.Cryostat);
316 
317  // Test maximum time difference.
318 
319  result = result && std::abs(t1 - t2) <= fMaxDT;
320 
321  // Test mc truth.
322 
323  if (result && useMC) {
324 
325  // Test whether hits have a common parent track id.
326 
327  const HitMCInfo& mcinfo2 = fHitMCMap[&hit2];
328  std::vector<int> tid2 = mcinfo2.trackIDs;
329  bool only_neg2 = tid2.size() > 0 && tid2.back() < 0;
330  std::vector<int>::iterator it = std::set_intersection(
331  tid1.begin(), tid1.end(), tid2.begin(), tid2.end(), tid2.begin());
332  tid2.resize(it - tid2.begin());
333 
334  // Hits are compatible if they have parents in common.
335  // If the only parent id in common is negative (-999),
336  // then hits are compatible only if both hits have only
337  // negative parent tracks.
338 
339  bool only_neg3 = tid2.size() > 0 && tid2.back() < 0;
340  mc_ok = tid2.size() > 0 && (!only_neg3 || (only_neg1 && only_neg2));
341  result = result && mc_ok;
342 
343  // If we are still OK, check that either hit is
344  // the nearest neighbor of the other.
345 
346  if (result) {
347  result = mcinfo1.pchit[hit2WireID.Plane] == &hit2 ||
348  mcinfo2.pchit[hit1WireID.Plane] == &hit1;
349  }
350  }
351  }
352  }
353  }
354 
355  // If there are exactly three hits, and they pass pairwise tests, check
356  // for spatial compatibility.
357 
358  if (result && nhits == 3) {
359 
360  // Loop over hits.
361 
362  double dist[3] = {0., 0., 0.};
363  double sinth[3] = {0., 0., 0.};
364  double costh[3] = {0., 0., 0.};
365 
366  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
367 
368  for (int i = 0; i < 3; ++i) {
369 
370  // Get tpc, plane, wire.
371 
372  const recob::Hit& hit = *(hits[i]);
373  geo::WireID hitWireID = hit.WireID();
374 
375  const geo::WireGeo& wgeom = wireReadoutGeom.Wire(hit.WireID());
376  if ((hitWireID.TPC != tpc) || (hitWireID.Cryostat != cstat))
377  throw cet::exception("SpacePointAlg") << "compatible(): geometry mismatch\n";
378 
379  // Get angles and distance of wire.
380 
381  double const hl = wgeom.HalfL();
382  auto const xyz = wgeom.GetCenter();
383  auto const xyz1 = wgeom.GetEnd();
384  double s = (xyz1.Y() - xyz.Y()) / hl;
385  double c = (xyz1.Z() - xyz.Z()) / hl;
386  sinth[hit.WireID().Plane] = s;
387  costh[hit.WireID().Plane] = c;
388  dist[hit.WireID().Plane] = xyz.Z() * s - xyz.Y() * c;
389  }
390 
391  // Do space cut.
392 
393  double S = ((sinth[1] * costh[2] - costh[1] * sinth[2]) * dist[0] +
394  (sinth[2] * costh[0] - costh[2] * sinth[0]) * dist[1] +
395  (sinth[0] * costh[1] - costh[0] * sinth[1]) * dist[2]);
396 
397  result = result && std::abs(S) < fMaxS;
398  }
399  }
400 
401  // Done.
402 
403  return result;
404  }
405 
406  //----------------------------------------------------------------------
407  // Fill one space point using a colleciton of hits.
408  // Assume points have already been tested for compatibility.
409  //
412  std::vector<recob::SpacePoint>& sptv,
413  int sptid) const
414  {
416 
417  double timePitch = detProp.GetXTicksCoefficient();
418 
419  int nhits = hits.size();
420 
421  // Remember associated hits internally.
422 
423  if (fSptHitMap.find(sptid) != fSptHitMap.end())
424  throw cet::exception("SpacePointAlg") << "fillSpacePoint(): hit already present!\n";
425  fSptHitMap[sptid] = hits;
426 
427  // Calculate position and error matrix.
428 
429  double xyz[3] = {0., 0., 0.};
430  double errxyz[6] = {0., 0., 0., 0., 0., 0.};
431 
432  // Calculate x using drift times.
433  // Loop over all hits and calculate the weighted average drift time.
434  // Also calculate time variance and chisquare.
435 
436  double sumt2w = 0.;
437  double sumtw = 0.;
438  double sumw = 0.;
439 
440  for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
441  ++ihit) {
442 
443  const recob::Hit& hit = **ihit;
444  geo::WireID hitWireID = hit.WireID();
445 
446  // Correct time for trigger offset and view-dependent time offsets.
447 
448  double t0 = detProp.GetXTicksOffset(hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat);
449  double t = hit.PeakTime() - t0;
450  double et = hit.SigmaPeakTime();
451  double w = 1. / (et * et);
452 
453  sumt2w += w * t * t;
454  sumtw += w * t;
455  sumw += w;
456  }
457 
458  double drift_time = 0.;
459  double var_time = 0.;
460  double chisq = 0.;
461  if (sumw != 0.) {
462  drift_time = sumtw / sumw;
463  //var_time = sumt2w / sumw - drift_time * drift_time;
464  var_time = 1. / sumw;
465  if (var_time < 0.) var_time = 0.;
466  chisq = sumt2w - sumtw * drift_time;
467  if (chisq < 0.) chisq = 0.;
468  }
469  xyz[0] = drift_time * timePitch;
470  errxyz[0] = var_time * timePitch * timePitch;
471 
472  // Calculate y, z using wires (need at least two hits).
473 
474  if (nhits >= 2) {
475 
476  // Calculate y and z by chisquare minimization of wire coordinates.
477 
478  double sus = 0.; // sum w_i u_i sin_th_i
479  double suc = 0.; // sum w_i u_i cos_th_i
480  double sc2 = 0.; // sum w_i cos2_th_i
481  double ss2 = 0.; // sum w_i sin2_th_i
482  double ssc = 0.; // sum w_i sin_th_i cos_th_i
483 
484  // Loop over points.
485 
486  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
487  for (auto const& hit : hits | ranges::views::transform(to_element)) {
488  geo::WireID hitWireID = hit.WireID();
489  const geo::WireGeo& wgeom = wireReadoutGeom.Wire(hit.WireID());
490 
491  // Calculate angle and wire coordinate in this view.
492 
493  double const hl = wgeom.HalfL();
494  auto const cen = wgeom.GetCenter();
495  auto const cen1 = wgeom.GetEnd();
496  double s = (cen1.Y() - cen.Y()) / hl;
497  double c = (cen1.Z() - cen.Z()) / hl;
498  double u = cen.Z() * s - cen.Y() * c;
499  double eu = wireReadoutGeom.Plane(hitWireID).WirePitch() / std::sqrt(12.);
500  double w = 1. / (eu * eu);
501 
502  // Summations
503 
504  sus += w * u * s;
505  suc += w * u * c;
506  sc2 += w * c * c;
507  ss2 += w * s * s;
508  ssc += w * s * c;
509  }
510 
511  // Calculate y,z
512 
513  double denom = sc2 * ss2 - ssc * ssc;
514  if (denom != 0.) {
515  xyz[1] = (-suc * ss2 + sus * ssc) / denom;
516  xyz[2] = (sus * sc2 - suc * ssc) / denom;
517  errxyz[2] = ss2 / denom;
518  errxyz[4] = ssc / denom;
519  errxyz[5] = sc2 / denom;
520  }
521 
522  // Make space point.
523 
524  recob::SpacePoint spt(xyz, errxyz, chisq, sptid);
525  sptv.push_back(spt);
526  }
527  return;
528  }
529 
549  std::vector<recob::SpacePoint>& spts,
550  std::multimap<double, KHitTrack> const& trackMap) const
551  {
552  // Loop over KHitTracks.
553 
555  art::PtrVector<recob::Hit> compatible_hits;
556  for (std::multimap<double, KHitTrack>::const_iterator it = trackMap.begin();
557  it != trackMap.end();
558  ++it) {
559  const KHitTrack& track = (*it).second;
560 
561  // Extrack Hit from track.
562 
563  const std::shared_ptr<const KHitBase>& hit = track.getHit();
564  const KHitWireX* phit = dynamic_cast<const KHitWireX*>(&*hit);
565  if (phit != 0) {
566  const art::Ptr<recob::Hit> prhit = phit->getHit();
567 
568  // Test this hit for compatibility.
569 
570  hits.push_back(prhit);
571  bool ok = this->compatible(detProp, hits);
572  if (!ok) {
573 
574  // The new hit is not compatible. Make a space point out of
575  // the last known compatible hits, provided there are at least
576  // two.
577 
578  if (compatible_hits.size() >= 2) {
579  this->fillSpacePoint(detProp, compatible_hits, spts, this->numHitMap());
580  compatible_hits.clear();
581  }
582 
583  // Forget about any previous hits.
584 
585  hits.clear();
586  hits.push_back(prhit);
587  }
588 
589  // Update the list of known compatible hits.
590 
591  compatible_hits = hits;
592  }
593  }
594 
595  // Maybe make one final space point.
596 
597  if (compatible_hits.size() >= 2) {
598  this->fillSpacePoint(detProp, compatible_hits, spts, this->numHitMap());
599  }
600  }
601 
602  //----------------------------------------------------------------------
603  // Fill one space point using a colleciton of hits.
604  // Assume points have already been tested for compatibility.
605  // This version assumes there can be multiple hits per view,
606  // and gives unequal weight to different hits.
607  //
610  std::vector<recob::SpacePoint>& sptv,
611  int sptid) const
612  {
614  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
615 
616  // Calculate time pitch.
617 
618  double timePitch = detProp.GetXTicksCoefficient(); // cm / tick
619 
620  // Figure out which tpc we are in.
621 
622  unsigned int tpc0 = 0;
623  unsigned int cstat0 = 0;
624  int nhits = hits.size();
625  if (nhits > 0) {
626  tpc0 = hits.front()->WireID().TPC;
627  cstat0 = hits.front()->WireID().Cryostat;
628  }
629 
630  // Remember associated hits internally.
631 
632  if (fSptHitMap.count(sptid) != 0)
633  throw cet::exception("SpacePointAlg") << "fillComplexSpacePoint(): hit already present!\n";
634  fSptHitMap[sptid] = hits;
635 
636  // Do a preliminary scan of hits.
637  // Determine weight given to hits in each view.
638 
639  unsigned int nplanes = wireReadoutGeom.Nplanes(geo::TPCID(cstat0, tpc0));
640  std::vector<int> numhits(nplanes, 0);
641  std::vector<double> weight(nplanes, 0.);
642 
643  for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
644  ++ihit) {
645 
646  const recob::Hit& hit = **ihit;
647  geo::WireID hitWireID = hit.WireID();
648  // kept as assertions for performance reasons
649  assert(hitWireID.Cryostat == cstat0);
650  assert(hitWireID.TPC == tpc0);
651  assert(hitWireID.Plane < nplanes);
652  ++numhits[hitWireID.Plane];
653  }
654 
655  for (unsigned int plane = 0; plane < nplanes; ++plane) {
656  double np = numhits[plane];
657  if (np > 0.) weight[plane] = 1. / (np * np * np);
658  }
659 
660  // Calculate position and error matrix.
661 
662  double xyz[3] = {0., 0., 0.};
663  double errxyz[6] = {0., 0., 0., 0., 0., 0.};
664 
665  // Calculate x using drift times.
666  // Loop over all hits and calculate the weighted average drift time.
667  // Also calculate time variance and chisquare.
668 
669  double sumt2w = 0.;
670  double sumtw = 0.;
671  double sumw = 0.;
672 
673  for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
674  ++ihit) {
675 
676  const recob::Hit& hit = **ihit;
677  geo::WireID hitWireID = hit.WireID();
678 
679  // Correct time for trigger offset and view-dependent time offsets.
680 
681  double t0 = detProp.GetXTicksOffset(hitWireID.Plane, hitWireID.TPC, hitWireID.Cryostat);
682  double t = hit.PeakTime() - t0;
683  double et = hit.SigmaPeakTime();
684  double w = weight[hitWireID.Plane] / (et * et);
685 
686  sumt2w += w * t * t;
687  sumtw += w * t;
688  sumw += w;
689  }
690 
691  double drift_time = 0.;
692  double var_time = 0.;
693  double chisq = 0.;
694  if (sumw != 0.) {
695  drift_time = sumtw / sumw;
696  var_time = sumt2w / sumw - drift_time * drift_time;
697  if (var_time < 0.) var_time = 0.;
698  chisq = sumt2w - sumtw * drift_time;
699  if (chisq < 0.) chisq = 0.;
700  }
701  xyz[0] = drift_time * timePitch;
702  errxyz[0] = var_time * timePitch * timePitch;
703 
704  // Calculate y, z using wires (need at least two hits).
705 
706  if (nhits >= 2) {
707 
708  // Calculate y and z by chisquare minimization of wire coordinates.
709 
710  double sus = 0.; // sum w_i u_i sin_th_i
711  double suc = 0.; // sum w_i u_i cos_th_i
712  double sc2 = 0.; // sum w_i cos2_th_i
713  double ss2 = 0.; // sum w_i sin2_th_i
714  double ssc = 0.; // sum w_i sin_th_i cos_th_i
715 
716  // Loop over points.
717 
718  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
719  for (auto const& hit : hits | ranges::views::transform(to_element)) {
720  geo::WireID hitWireID = hit.WireID();
721  const geo::WireGeo& wgeom = wireReadoutGeom.Wire(hit.WireID());
722 
723  // Calculate angle and wire coordinate in this view.
724 
725  double const hl = wgeom.HalfL();
726  auto const cen = wgeom.GetCenter();
727  auto const cen1 = wgeom.GetEnd();
728  double s = (cen1.Y() - cen.Y()) / hl;
729  double c = (cen1.Z() - cen.Z()) / hl;
730  double u = cen.Z() * s - cen.Y() * c;
731  double eu = wireReadoutGeom.Plane(hitWireID).WirePitch() / std::sqrt(12.);
732  double w = weight[hitWireID.Plane] / (eu * eu);
733 
734  // Summations
735 
736  sus += w * u * s;
737  suc += w * u * c;
738  sc2 += w * c * c;
739  ss2 += w * s * s;
740  ssc += w * s * c;
741  }
742 
743  // Calculate y,z
744 
745  double denom = sc2 * ss2 - ssc * ssc;
746  if (denom != 0.) {
747  xyz[1] = (-suc * ss2 + sus * ssc) / denom;
748  xyz[2] = (sus * sc2 - suc * ssc) / denom;
749  errxyz[2] = ss2 / denom;
750  errxyz[4] = ssc / denom;
751  errxyz[5] = sc2 / denom;
752  }
753 
754  // Make space point.
755 
756  recob::SpacePoint spt(xyz, errxyz, chisq, sptid);
757  sptv.push_back(spt);
758  }
759  }
760 
761  //----------------------------------------------------------------------
762  // Fill a vector of space points for all compatible combinations of hits
763  // from an input vector of hits (non-mc-truth version).
764  //
766  detinfo::DetectorPropertiesData const& detProp,
768  std::vector<recob::SpacePoint>& spts) const
769  {
770  makeSpacePoints(clockData, detProp, hits, spts, false);
771  }
772 
773  //----------------------------------------------------------------------
774  // Fill a vector of space points for all compatible combinations of hits
775  // from an input vector of hits (mc-truth version).
776  //
778  detinfo::DetectorPropertiesData const& detProp,
780  std::vector<recob::SpacePoint>& spts) const
781  {
782  makeSpacePoints(clockData, detProp, hits, spts, true);
783  }
784 
785  //----------------------------------------------------------------------
786  // Fill a vector of space points for all compatible combinations of hits
787  // from an input vector of hits (general version).
788  //
790  detinfo::DetectorPropertiesData const& detProp,
792  std::vector<recob::SpacePoint>& spts,
793  bool useMC) const
794  {
796 
797  fSptHitMap.clear();
798 
799  // Print diagnostic information.
800 
801  update(detProp);
802 
803  // First make sure result vector is empty.
804 
805  spts.erase(spts.begin(), spts.end());
806 
807  // Statistics.
808 
809  int n2 = 0; // Number of two-hit space points.
810  int n3 = 0; // Number of three-hit space points.
811  int n2filt = 0; // Number of two-hit space points after filtering/merging.
812  int n3filt = 0; // Number of three-hit space pointe after filtering/merging.
813 
814  // Sort hits into maps indexed by [cryostat][tpc][plane][wire].
815  // If using mc information, also generate maps of sim::IDEs and mc
816  // position indexed by hit.
817 
818  std::vector<std::vector<std::vector<std::multimap<unsigned int, art::Ptr<recob::Hit>>>>> hitmap;
819  fHitMCMap.clear();
820 
821  unsigned int ncstat = geom->Ncryostats();
822  hitmap.resize(ncstat);
823 
824  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
825 
826  for (auto const& cryoid : geom->Iterate<geo::CryostatID>()) {
827  unsigned int cstat = cryoid.Cryostat;
828  unsigned int ntpc = geom->Cryostat(cryoid).NTPC();
829  hitmap[cstat].resize(ntpc);
830  for (unsigned int tpc = 0; tpc < ntpc; ++tpc) {
831  unsigned int nplane = wireReadoutGeom.Nplanes({cstat, tpc});
832  hitmap[cstat][tpc].resize(nplane);
833  }
834  }
835 
836  for (art::PtrVector<recob::Hit>::const_iterator ihit = hits.begin(); ihit != hits.end();
837  ++ihit) {
838  const art::Ptr<recob::Hit>& phit = *ihit;
839  geo::View_t view = phit->View();
840  if ((view == geo::kU && fEnableU) || (view == geo::kV && fEnableV) ||
841  (view == geo::kZ && fEnableW)) {
842  geo::WireID phitWireID = phit->WireID();
843  hitmap[phitWireID.Cryostat][phitWireID.TPC][phitWireID.Plane].insert(
844  std::make_pair(phitWireID.Wire, phit));
845  }
846  }
847 
848  // Fill mc information, including IDEs and closest neighbors
849  // of each hit.
852  if (useMC) {
854 
855  // First loop over hits and fill track ids and mc position.
856  for (auto const& id : wireReadoutGeom.Iterate<geo::PlaneID>()) {
857  auto const [cstat, tpc, plane] = std::make_tuple(id.Cryostat, id.TPC, id.Plane);
858  int nplane = wireReadoutGeom.Nplanes(id);
859  for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator ihit =
860  hitmap[cstat][tpc][plane].begin();
861  ihit != hitmap[cstat][tpc][plane].end();
862  ++ihit) {
863  const art::Ptr<recob::Hit>& phit = ihit->second;
864  const recob::Hit& hit = *phit;
865  HitMCInfo& mcinfo = fHitMCMap[&hit]; // Default HitMCInfo.
866 
867  // Fill default nearest neighbor information (i.e. none).
868 
869  mcinfo.pchit.resize(nplane, 0);
870  mcinfo.dist2.resize(nplane, 1.e20);
871 
872  // Get sim::IDEs for this hit.
873 
874  std::vector<sim::IDE> ides = bt_serv->HitToAvgSimIDEs(clockData, phit);
875 
876  // Get sorted track ids. for this hit.
877 
878  mcinfo.trackIDs.reserve(ides.size());
879  for (std::vector<sim::IDE>::const_iterator i = ides.begin(); i != ides.end(); ++i)
880  mcinfo.trackIDs.push_back(i->trackID);
881  sort(mcinfo.trackIDs.begin(), mcinfo.trackIDs.end());
882 
883  // Get position of ionization for this hit.
884 
885  try {
886  mcinfo.xyz = bt_serv->SimIDEsToXYZ(ides);
887  }
888  catch (cet::exception& x) {
889  mcinfo.xyz.clear();
890  }
891  } // end loop over ihit
892  } // end loop oer planes
893 
894  // Loop over hits again and fill nearest neighbor information for real.
895  for (auto const& id : wireReadoutGeom.Iterate<geo::PlaneID>()) {
896  auto const [cstat, tpc, plane] = std::make_tuple(id.Cryostat, id.TPC, id.Plane);
897  int nplane = wireReadoutGeom.Nplanes(id);
898  for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator ihit =
899  hitmap[cstat][tpc][plane].begin();
900  ihit != hitmap[cstat][tpc][plane].end();
901  ++ihit) {
902  const art::Ptr<recob::Hit>& phit = ihit->second;
903  const recob::Hit& hit = *phit;
904  HitMCInfo& mcinfo = fHitMCMap[&hit];
905  if (mcinfo.xyz.size() != 0) {
906  assert(mcinfo.xyz.size() == 3);
907 
908  // Fill nearest neighbor information for this hit.
909 
910  for (int plane2 = 0; plane2 < nplane; ++plane2) {
911  for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator jhit =
912  hitmap[cstat][tpc][plane2].begin();
913  jhit != hitmap[cstat][tpc][plane2].end();
914  ++jhit) {
915  const art::Ptr<recob::Hit>& phit2 = jhit->second;
916  const recob::Hit& hit2 = *phit2;
917  const HitMCInfo& mcinfo2 = fHitMCMap[&hit2];
918 
919  if (mcinfo2.xyz.size() != 0) {
920  assert(mcinfo2.xyz.size() == 3);
921  double dx = mcinfo.xyz[0] - mcinfo2.xyz[0];
922  double dy = mcinfo.xyz[1] - mcinfo2.xyz[1];
923  double dz = mcinfo.xyz[2] - mcinfo2.xyz[2];
924  double dist2 = dx * dx + dy * dy + dz * dz;
925  if (dist2 < mcinfo.dist2[plane2]) {
926  mcinfo.dist2[plane2] = dist2;
927  mcinfo.pchit[plane2] = &hit2;
928  }
929  } // end if mcinfo2.xyz valid
930  } // end loop over jhit
931  } // end loop over plane2
932  } // end if mcinfo.xyz valid.
933  } // end loop over ihit
934  } // end loop over planes
935  } // end if MC
936 
937  // use mf::LogDebug instead of MF_LOG_DEBUG because we reuse it in many lines
938  // insertions are protected by mf::isDebugEnabled()
939  mf::LogDebug debug("SpacePointAlg");
940  if (mf::isDebugEnabled()) {
941  debug << "Total hits = " << hits.size() << "\n\n";
942  for (auto const& id : geom->Iterate<geo::TPCID>()) {
943  auto const [cstat, tpc] = std::make_tuple(id.Cryostat, id.TPC);
944  int nplane = hitmap[cstat][tpc].size();
945  for (int plane = 0; plane < nplane; ++plane) {
946  debug << "TPC, Plane: " << tpc << ", " << plane
947  << ", hits = " << hitmap[cstat][tpc][plane].size() << "\n";
948  }
949  } // end loop over TPCs
950  } // if debug
951 
952  // Make empty multimap from hit pointer on preferred
953  // (most-populated or collection) plane to space points that
954  // include that hit (used for sorting, filtering, and
955  // merging).
956 
957  typedef const recob::Hit* sptkey_type;
958  std::multimap<sptkey_type, recob::SpacePoint> sptmap;
959  std::set<sptkey_type> sptkeys; // Keys of multimap.
960 
961  // Loop over TPCs.
962  for (auto const& tpcid : geom->Iterate<geo::TPCID>()) {
963 
964  auto const [cstat, tpc] = std::make_tuple(tpcid.Cryostat, tpcid.TPC);
965  // Sort maps in increasing order of number of hits.
966  // This is so that we can do the outer loops over hits
967  // over the views with fewer hits.
968  //
969  // If config parameter PreferColl is true, treat the colleciton
970  // plane as if it had the most hits, regardless of how many
971  // hits it actually has. This will force space points to be
972  // filtered and merged with respect to the collection plane
973  // wires. It will also force space points to be sorted by
974  // collection plane wire.
975 
976  int nplane = hitmap[cstat][tpc].size();
977  std::vector<int> index(nplane);
978 
979  for (int i = 0; i < nplane; ++i)
980  index[i] = i;
981 
982  for (int i = 0; i < nplane - 1; ++i) {
983  for (int j = i + 1; j < nplane; ++j) {
984  bool icoll = fPreferColl && wireReadoutGeom.SignalType(geo::PlaneID(tpcid, index[i])) ==
986  bool jcoll = fPreferColl && wireReadoutGeom.SignalType(geo::PlaneID(tpcid, index[j])) ==
988  if ((hitmap[cstat][tpc][index[i]].size() > hitmap[cstat][tpc][index[j]].size() &&
989  !jcoll) ||
990  icoll) {
991  int temp = index[i];
992  index[i] = index[j];
993  index[j] = temp;
994  }
995  }
996  } // end loop over i
997 
998  // how many views with hits?
999  // This will allow for the special case where we might have only 2 planes of information and
1000  // still want space points even if a three plane TPC
1001  std::vector<std::multimap<unsigned int, art::Ptr<recob::Hit>>>& hitsByPlaneVec =
1002  hitmap[cstat][tpc];
1003  int nViewsWithHits(0);
1004 
1005  for (int i = 0; i < nplane; i++) {
1006  if (hitsByPlaneVec[index[i]].size() > 0) nViewsWithHits++;
1007  }
1008 
1009  // If two-view space points are allowed, make a double loop
1010  // over hits and produce space points for compatible hit-pairs.
1011 
1012  if ((nViewsWithHits == 2 || nplane == 2) && fMinViews <= 2) {
1013 
1014  // Loop over pairs of views.
1015  for (int i = 0; i < nplane - 1; ++i) {
1016  unsigned int plane1 = index[i];
1017 
1018  if (hitmap[cstat][tpc][plane1].empty()) continue;
1019 
1020  for (int j = i + 1; j < nplane; ++j) {
1021  unsigned int plane2 = index[j];
1022 
1023  if (hitmap[cstat][tpc][plane2].empty()) continue;
1024 
1025  // Get angle, pitch, and offset of plane2 wires.
1026  geo::PlaneID const plane2_id{tpcid, plane2};
1027  const geo::WireGeo& wgeo2 = wireReadoutGeom.Plane(plane2_id).Wire(0);
1028  double const hl2 = wgeo2.HalfL();
1029  auto const xyz21 = wgeo2.GetStart();
1030  auto const xyz22 = wgeo2.GetEnd();
1031  double s2 = (xyz22.Y() - xyz21.Y()) / (2. * hl2);
1032  double c2 = (xyz22.Z() - xyz21.Z()) / (2. * hl2);
1033  double dist2 = -xyz21.Y() * c2 + xyz21.Z() * s2;
1034  double pitch2 = wireReadoutGeom.Plane(plane2_id).WirePitch();
1035 
1036  if (!fPreferColl &&
1037  hitmap[cstat][tpc][plane1].size() > hitmap[cstat][tpc][plane2].size())
1038  throw cet::exception("SpacePointAlg")
1039  << "makeSpacePoints(): hitmaps with incompatible size\n";
1040 
1041  // Loop over pairs of hits.
1042 
1044  hitvec.reserve(2);
1045 
1046  for (std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator ihit1 =
1047  hitmap[cstat][tpc][plane1].begin();
1048  ihit1 != hitmap[cstat][tpc][plane1].end();
1049  ++ihit1) {
1050 
1051  const art::Ptr<recob::Hit>& phit1 = ihit1->second;
1052  geo::WireID phit1WireID = phit1->WireID();
1053  const geo::WireGeo& wgeo = wireReadoutGeom.Wire(phit1WireID);
1054 
1055  // Get endpoint coordinates of this wire.
1056  // (kept as assertions for performance reasons)
1057  assert(phit1WireID.Cryostat == cstat);
1058  assert(phit1WireID.TPC == tpc);
1059  assert(phit1WireID.Plane == plane1);
1060  auto const xyz1 = wgeo.GetStart();
1061  auto const xyz2 = wgeo.GetEnd();
1062 
1063  // Find the plane2 wire numbers corresponding to the endpoints.
1064 
1065  double wire21 = (-xyz1.Y() * c2 + xyz1.Z() * s2 - dist2) / pitch2;
1066  double wire22 = (-xyz2.Y() * c2 + xyz2.Z() * s2 - dist2) / pitch2;
1067 
1068  int wmin = std::max(0., std::min(wire21, wire22));
1069  int wmax = std::max(0., std::max(wire21, wire22) + 1.);
1070 
1071  std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
1072  ihit2 = hitmap[cstat][tpc][plane2].lower_bound(wmin),
1073  ihit2end = hitmap[cstat][tpc][plane2].upper_bound(wmax);
1074 
1075  for (; ihit2 != ihit2end; ++ihit2) {
1076 
1077  const art::Ptr<recob::Hit>& phit2 = ihit2->second;
1078 
1079  // Check current pair of hits for compatibility.
1080  // By construction, hits should always have compatible views
1081  // and times, but may not have compatible mc information.
1082 
1083  hitvec.clear();
1084  hitvec.push_back(phit1);
1085  hitvec.push_back(phit2);
1086  bool ok = compatible(detProp, hitvec, useMC);
1087  if (ok) {
1088 
1089  // Add a space point.
1090 
1091  ++n2;
1092 
1093  // make a dummy vector of recob::SpacePoints
1094  // as we are filtering or merging and don't want to
1095  // add the created SpacePoint to the final collection just yet
1096  // This dummy vector will hold just one recob::SpacePoint,
1097  // which will go into the multimap and then the vector
1098  // will go out of scope.
1099 
1100  std::vector<recob::SpacePoint> sptv;
1101  fillSpacePoint(detProp, hitvec, sptv, sptmap.size());
1102  sptkey_type key = &*phit2;
1103  sptmap.insert(std::pair<sptkey_type, recob::SpacePoint>(key, sptv.back()));
1104  sptkeys.insert(key);
1105  }
1106  }
1107  }
1108  }
1109  }
1110  } // end if fMinViews <= 2
1111 
1112  // If three-view space points are allowed, make a triple loop
1113  // over hits and produce space points for compatible triplets.
1114 
1115  if (nplane >= 3 && fMinViews <= 3) {
1116 
1117  // Loop over triplets of hits.
1118 
1120  hitvec.reserve(3);
1121 
1122  unsigned int plane1 = index[0];
1123  unsigned int plane2 = index[1];
1124  unsigned int plane3 = index[2];
1125 
1126  // Get angle, pitch, and offset of plane1 wires.
1127 
1128  geo::PlaneID const plane1_id{tpcid, plane1};
1129  const geo::WireGeo& wgeo1 = wireReadoutGeom.Plane(plane1_id).Wire(0);
1130  double const hl1 = wgeo1.HalfL();
1131  auto const xyz11 = wgeo1.GetStart();
1132  auto const xyz12 = wgeo1.GetEnd();
1133  double s1 = (xyz12.Y() - xyz11.Y()) / (2. * hl1);
1134  double c1 = (xyz12.Z() - xyz11.Z()) / (2. * hl1);
1135  double dist1 = -xyz11.Y() * c1 + xyz11.Z() * s1;
1136  double pitch1 = wireReadoutGeom.Plane(plane1_id).WirePitch();
1137  const double TicksOffset1 = detProp.GetXTicksOffset(plane1_id);
1138 
1139  // Get angle, pitch, and offset of plane2 wires.
1140 
1141  geo::PlaneID const plane2_id{tpcid, plane2};
1142  const geo::WireGeo& wgeo2 = wireReadoutGeom.Plane(plane2_id).Wire(0);
1143  double const hl2 = wgeo2.HalfL();
1144  auto const xyz21 = wgeo2.GetStart();
1145  auto const xyz22 = wgeo2.GetEnd();
1146  double s2 = (xyz22.Y() - xyz21.Y()) / (2. * hl2);
1147  double c2 = (xyz22.Z() - xyz21.Z()) / (2. * hl2);
1148  double dist2 = -xyz21.Y() * c2 + xyz21.Z() * s2;
1149  double pitch2 = wireReadoutGeom.Plane(plane2_id).WirePitch();
1150  const double TicksOffset2 = detProp.GetXTicksOffset(plane2_id);
1151 
1152  // Get angle, pitch, and offset of plane3 wires.
1153 
1154  geo::PlaneID const plane3_id{tpcid, plane3};
1155  const geo::WireGeo& wgeo3 = wireReadoutGeom.Plane(plane3_id).Wire(0);
1156  double const hl3 = wgeo3.HalfL();
1157  auto const xyz31 = wgeo3.GetStart();
1158  auto const xyz32 = wgeo3.GetEnd();
1159  double s3 = (xyz32.Y() - xyz31.Y()) / (2. * hl3);
1160  double c3 = (xyz32.Z() - xyz31.Z()) / (2. * hl3);
1161  double dist3 = -xyz31.Y() * c3 + xyz31.Z() * s3;
1162  double pitch3 = wireReadoutGeom.Plane(plane3_id).WirePitch();
1163  const double TicksOffset3 = detProp.GetXTicksOffset(plane3_id);
1164 
1165  // Get sine of angle differences.
1166 
1167  double s12 = s1 * c2 - s2 * c1; // sin(theta1 - theta2).
1168  double s23 = s2 * c3 - s3 * c2; // sin(theta2 - theta3).
1169  double s31 = s3 * c1 - s1 * c3; // sin(theta3 - theta1).
1170 
1171  // Loop over hits in plane1.
1172 
1173  std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
1174  ihit1 = hitmap[cstat][tpc][plane1].begin(),
1175  ihit1end = hitmap[cstat][tpc][plane1].end();
1176  for (; ihit1 != ihit1end; ++ihit1) {
1177 
1178  unsigned int wire1 = ihit1->first;
1179  const art::Ptr<recob::Hit>& phit1 = ihit1->second;
1180  geo::WireID phit1WireID = phit1->WireID();
1181  const geo::WireGeo& wgeo = wireReadoutGeom.Wire(phit1WireID);
1182 
1183  // Get endpoint coordinates of this wire from plane1.
1184  // (kept as assertions for performance reasons)
1185  assert(phit1WireID.Cryostat == cstat);
1186  assert(phit1WireID.TPC == tpc);
1187  assert(phit1WireID.Plane == plane1);
1188  assert(phit1WireID.Wire == wire1);
1189  auto const xyz1 = wgeo.GetStart();
1190  auto const xyz2 = wgeo.GetEnd();
1191 
1192  // Get corrected time and oblique coordinate of first hit.
1193 
1194  double t1 = phit1->PeakTime() - TicksOffset1;
1195  double u1 = wire1 * pitch1 + dist1;
1196 
1197  // Find the plane2 wire numbers corresponding to the endpoints.
1198 
1199  double wire21 = (-xyz1.Y() * c2 + xyz1.Z() * s2 - dist2) / pitch2;
1200  double wire22 = (-xyz2.Y() * c2 + xyz2.Z() * s2 - dist2) / pitch2;
1201 
1202  int wmin = std::max(0., std::min(wire21, wire22));
1203  int wmax = std::max(0., std::max(wire21, wire22) + 1.);
1204 
1205  std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
1206  ihit2 = hitmap[cstat][tpc][plane2].lower_bound(wmin),
1207  ihit2end = hitmap[cstat][tpc][plane2].upper_bound(wmax);
1208 
1209  for (; ihit2 != ihit2end; ++ihit2) {
1210 
1211  int wire2 = ihit2->first;
1212  const art::Ptr<recob::Hit>& phit2 = ihit2->second;
1213 
1214  // Get corrected time of second hit.
1215 
1216  double t2 = phit2->PeakTime() - TicksOffset2;
1217 
1218  // Check maximum time difference with first hit.
1219 
1220  bool dt12ok = std::abs(t1 - t2) <= fMaxDT;
1221  if (dt12ok) {
1222 
1223  // Test first two hits for compatibility before looping
1224  // over third hit.
1225 
1226  hitvec.clear();
1227  hitvec.push_back(phit1);
1228  hitvec.push_back(phit2);
1229  bool h12ok = compatible(detProp, hitvec, useMC);
1230  if (h12ok) {
1231 
1232  // Get oblique coordinate of second hit.
1233 
1234  double u2 = wire2 * pitch2 + dist2;
1235 
1236  // Predict plane3 oblique coordinate and wire number.
1237 
1238  double u3pred = (-u1 * s23 - u2 * s31) / s12;
1239  double w3pred = (u3pred - dist3) / pitch3;
1240  double w3delta = std::abs(fMaxS / (s12 * pitch3));
1241  int w3min = std::max(0., std::ceil(w3pred - w3delta));
1242  int w3max = std::max(0., std::floor(w3pred + w3delta));
1243 
1244  std::map<unsigned int, art::Ptr<recob::Hit>>::const_iterator
1245  ihit3 = hitmap[cstat][tpc][plane3].lower_bound(w3min),
1246  ihit3end = hitmap[cstat][tpc][plane3].upper_bound(w3max);
1247 
1248  for (; ihit3 != ihit3end; ++ihit3) {
1249 
1250  int wire3 = ihit3->first;
1251  const art::Ptr<recob::Hit>& phit3 = ihit3->second;
1252 
1253  // Get corrected time of third hit.
1254 
1255  double t3 = phit3->PeakTime() - TicksOffset3;
1256 
1257  // Check time difference of third hit compared to first two hits.
1258 
1259  bool dt123ok = std::abs(t1 - t3) <= fMaxDT && std::abs(t2 - t3) <= fMaxDT;
1260  if (dt123ok) {
1261 
1262  // Get oblique coordinate of third hit and check spatial separation.
1263 
1264  double u3 = wire3 * pitch3 + dist3;
1265  double S = s23 * u1 + s31 * u2 + s12 * u3;
1266  bool sok = std::abs(S) <= fMaxS;
1267  if (sok) {
1268 
1269  // Test triplet for compatibility.
1270 
1271  hitvec.clear();
1272  hitvec.push_back(phit1);
1273  hitvec.push_back(phit2);
1274  hitvec.push_back(phit3);
1275  bool h123ok = compatible(detProp, hitvec, useMC);
1276  if (h123ok) {
1277 
1278  // Add a space point.
1279 
1280  ++n3;
1281 
1282  // make a dummy vector of recob::SpacePoints
1283  // as we are filtering or merging and don't want to
1284  // add the created SpacePoint to the final collection just yet
1285  // This dummy vector will hold just one recob::SpacePoint,
1286  // which will go into the multimap and then the vector
1287  // will go out of scope.
1288 
1289  std::vector<recob::SpacePoint> sptv;
1290  fillSpacePoint(detProp, hitvec, sptv, sptmap.size() - 1);
1291  sptkey_type key = &*phit3;
1292  sptmap.insert(std::pair<sptkey_type, recob::SpacePoint>(key, sptv.back()));
1293  sptkeys.insert(key);
1294  }
1295  }
1296  }
1297  }
1298  }
1299  }
1300  }
1301  }
1302  } // end if fMinViews <= 3
1303 
1304  // Do Filtering.
1305 
1306  if (fFilter) {
1307 
1308  // Transfer (some) space points from sptmap to spts.
1309 
1310  spts.reserve(spts.size() + sptkeys.size());
1311 
1312  // Loop over keys of space point map.
1313  // Space points that have the same key are candidates for filtering.
1314 
1315  for (std::set<sptkey_type>::const_iterator i = sptkeys.begin(); i != sptkeys.end(); ++i) {
1316  sptkey_type key = *i;
1317 
1318  // Loop over space points corresponding to the current key.
1319  // Choose the single best space point from among this group.
1320 
1321  double best_chisq = 0.;
1322  const recob::SpacePoint* best_spt = 0;
1323 
1325  sptmap.lower_bound(key);
1326  j != sptmap.upper_bound(key);
1327  ++j) {
1328  const recob::SpacePoint& spt = j->second;
1329  if (best_spt == 0 || spt.Chisq() < best_chisq) {
1330  best_spt = &spt;
1331  best_chisq = spt.Chisq();
1332  }
1333  }
1334 
1335  // Transfer best filtered space point to result vector.
1336 
1337  if (!best_spt)
1338  throw cet::exception("SpacePointAlg") << "makeSpacePoints(): no best point\n";
1339  spts.push_back(*best_spt);
1340  if (fMinViews <= 2)
1341  ++n2filt;
1342  else
1343  ++n3filt;
1344  }
1345  } // end if filtering
1346 
1347  // Do merging.
1348 
1349  else if (fMerge) {
1350 
1351  // Transfer merged space points from sptmap to spts.
1352 
1353  spts.reserve(spts.size() + sptkeys.size());
1354 
1355  // Loop over keys of space point map.
1356  // Space points that have the same key are candidates for merging.
1357 
1358  for (std::set<sptkey_type>::const_iterator i = sptkeys.begin(); i != sptkeys.end(); ++i) {
1359  sptkey_type key = *i;
1360 
1361  // Loop over space points corresponding to the current key.
1362  // Make a collection of hits that is the union of the hits
1363  // from each candidate space point.
1364 
1366  sptmap.lower_bound(key),
1367  jSPTend =
1368  sptmap.upper_bound(key);
1369 
1370  art::PtrVector<recob::Hit> merged_hits;
1371  for (; jSPT != jSPTend; ++jSPT) {
1372  const recob::SpacePoint& spt = jSPT->second;
1373 
1374  // Loop over hits from this space points.
1375  // Add each hit to the collection of all hits.
1376 
1377  const art::PtrVector<recob::Hit>& spt_hits = getAssociatedHits(spt);
1378  merged_hits.reserve(merged_hits.size() +
1379  spt_hits.size()); // better than nothing, but not ideal
1381  k != spt_hits.end();
1382  ++k) {
1383  const art::Ptr<recob::Hit>& hit = *k;
1384  merged_hits.push_back(hit);
1385  }
1386  }
1387 
1388  // Remove duplicates.
1389 
1390  std::sort(merged_hits.begin(), merged_hits.end());
1392  std::unique(merged_hits.begin(), merged_hits.end());
1393  merged_hits.erase(it, merged_hits.end());
1394 
1395  // Construct a complex space points using merged hits.
1396 
1397  fillComplexSpacePoint(detProp, merged_hits, spts, sptmap.size() + spts.size() - 1);
1398 
1399  if (fMinViews <= 2)
1400  ++n2filt;
1401  else
1402  ++n3filt;
1403  }
1404  } // end if merging
1405 
1406  // No filter, no merge.
1407 
1408  else {
1409 
1410  // Transfer all space points from sptmap to spts.
1411 
1412  spts.reserve(spts.size() + sptkeys.size());
1413 
1414  // Loop over space points.
1415 
1417  j != sptmap.end();
1418  ++j) {
1419  const recob::SpacePoint& spt = j->second;
1420  spts.push_back(spt);
1421  }
1422 
1423  // Update statistics.
1424 
1425  n2filt = n2;
1426  n3filt = n3;
1427  }
1428  } // end loop over tpcs
1429 
1430  if (mf::isDebugEnabled()) {
1431  debug << "\n2-hit space points = " << n2 << "\n"
1432  << "3-hit space points = " << n3 << "\n"
1433  << "2-hit filtered/merged space points = " << n2filt << "\n"
1434  << "3-hit filtered/merged space points = " << n3filt;
1435  } // if debug
1436  }
1437 
1438  //----------------------------------------------------------------------
1439  // Get hits associated with a particular space point, based on most recent
1440  // call of any make*SpacePoints method.
1442  const recob::SpacePoint& spt) const
1443  {
1444  // It is an error if no hits are associated with this space point (throw exception).
1445 
1446  std::map<int, art::PtrVector<recob::Hit>>::const_iterator it = fSptHitMap.find(spt.ID());
1447  if (it == fSptHitMap.end()) {
1448  mf::LogWarning("SpacePointAlg")
1449  << "Looking for ID " << spt.ID() << " from " << fSptHitMap.size() << std::endl;
1450  throw cet::exception("SpacePointAlg") << "No Hits associated with space point.\n";
1451  }
1452  return (*it).second;
1453  }
1454 
1455 } // end namespace trkf
Float_t x
Definition: compare.C:6
code to link reconstructed objects back to the MC truth information
intermediate_table::iterator iterator
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:112
void reserve(size_type n)
Definition: PtrVector.h:337
Basic Kalman filter track class, plus one measurement on same surface.
const std::shared_ptr< const KHitBase > & getHit() const
Measurement.
Definition: KHitTrack.h:51
Point_t const & GetCenter() const
Returns the world coordinate of the center of the wire [cm].
Definition: WireGeo.h:219
std::vector< const recob::Hit * > pchit
Pointer to nearest neighbor hit (indexed by plane).
typename data_t::iterator iterator
Definition: PtrVector.h:54
std::map< const recob::Hit *, HitMCInfo > fHitMCMap
double fTickOffsetU
Tick offset for plane U.
constexpr to_element_t to_element
Definition: ToElement.h:25
TTree * t1
Definition: plottest35.C:26
Encapsulate the construction of a single cyostat .
double fMaxDT
Maximum time difference between planes.
double GetXTicksCoefficient(int t, int c) const
bool fFilter
Filter flag.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Planes which measure V.
Definition: geo_types.h:132
double GetXTicksOffset(int p, int t, int c) const
Float_t ssc
Definition: plot.C:23
enum geo::_plane_orient Orient_t
Enumerate the possible plane projections.
iterator begin()
Definition: PtrVector.h:217
Declaration of signal hit object.
Kalman filter wire-time measurement on a SurfWireX surface.
void update(detinfo::DetectorPropertiesData const &detProp) const
The data type to uniquely identify a Plane.
Definition: geo_types.h:364
std::vector< double > dist2
Distance to nearest neighbor hit (indexed by plane).
bool fEnableW
Enable flag (W).
iterator erase(iterator position)
Definition: PtrVector.h:504
std::vector< int > trackIDs
Parent trackIDs.
const art::Ptr< recob::Hit > & getHit() const
Get original hit.
Definition: KHitWireX.h:47
constexpr auto abs(T v)
Returns the absolute value of the argument.
Double32_t Chisq() const
Definition: SpacePoint.h:86
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
intermediate_table::const_iterator const_iterator
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:290
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:303
Point_t GetStart() const
Returns the world coordinate of one end of the wire [cm].
Definition: WireGeo.h:224
Point_t GetEnd() const
Returns the world coordinate of one end of the wire [cm].
Definition: WireGeo.h:229
const art::PtrVector< recob::Hit > & getAssociatedHits(const recob::SpacePoint &spt) const
void makeMCTruthSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
double separation(const art::PtrVector< recob::Hit > &hits) const
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
Planes that are in the horizontal plane.
Definition: geo_types.h:142
void hits()
Definition: readHits.C:15
SpacePointAlg(const fhicl::ParameterSet &pset)
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
TCanvas * c1
Definition: plotHisto.C:7
bool fEnableU
Enable flag (U).
TCanvas * c2
Definition: plot_hist.C:75
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
Signal from induction planes.
Definition: geo_types.h:147
DebugStuff debug
Definition: DebugStruct.cxx:4
std::map< int, art::PtrVector< recob::Hit > > fSptHitMap
Planes that are in the vertical plane (e.g. ArgoNeuT).
Definition: geo_types.h:143
std::vector< double > xyz
Location of ionization (all tracks).
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
void makeSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
T get(std::string const &key) const
Definition: ParameterSet.h:314
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:67
iterator end()
Definition: PtrVector.h:231
void fillComplexSpacePoint(detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:171
bool isDebugEnabled()
double fMaxS
Maximum space separation between wires.
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
double correctedTime(detinfo::DetectorPropertiesData const &detProp, const recob::Hit &hit) const
bool fPreferColl
Sort by collection wire.
CryostatGeo const & Cryostat(CryostatID const &cryoid=details::cryostat_zero) const
Returns the specified cryostat.
size_type size() const
Definition: PtrVector.h:302
void fillSpacePoints(detinfo::DetectorPropertiesData const &detProp, std::vector< recob::SpacePoint > &spts, std::multimap< double, KHitTrack > const &trackMap) const
Fill a collection of space points.
Detector simulation of raw signals on wires.
TTree * t2
Definition: plottest35.C:36
Class defining a plane for tracking. It provides various functionalities to convert track parameters ...
Definition: TrackingPlane.h:37
Encapsulate the geometry of a wire .
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:170
reference front()
Definition: PtrVector.h:373
double weight
Definition: plottest35.C:25
ID_t ID() const
Definition: SpacePoint.h:74
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:226
double fTickOffsetV
Tick offset for plane V.
Encapsulate the construction of a single detector plane .
double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Contains all timing reference information for the detector.
std::vector< sim::IDE > HitToAvgSimIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double fTickOffsetW
Tick offset for plane W.
std::vector< double > SimIDEsToXYZ(std::vector< sim::IDE > const &ides) const
bool compatible(detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, bool useMC=false) const
int fMinViews
Mininum number of views per space point.
int numHitMap() const
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
float SigmaPeakTime() const
Uncertainty for the signal peak, in tick units.
Definition: Hit.h:230
Planes which measure W (third view for Bo, MicroBooNE, etc).
Definition: geo_types.h:133
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
Algorithm for generating space points from hits.
bool fMerge
Merge flag.
void fillSpacePoint(detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &sptv, int sptid) const
Float_t track
Definition: plot.C:35
void clear()
Definition: PtrVector.h:533
Float_t w
Definition: plot.C:20
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
Encapsulate the construction of a single detector plane .
The data type to uniquely identify a cryostat.
Definition: geo_types.h:187
Signal from collection planes.
Definition: geo_types.h:148
bool fEnableV
Enable flag (V).