LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
TripletFinder.cxx
Go to the documentation of this file.
2 
4 
5 #include "TVector3.h"
6 
10 
11 namespace reco3d {
12  // -------------------------------------------------------------------------
14  const std::vector<art::Ptr<recob::Hit>>& xhits,
15  const std::vector<art::Ptr<recob::Hit>>& uhits,
16  const std::vector<art::Ptr<recob::Hit>>& vhits,
17  const std::vector<raw::ChannelID_t>& xbad,
18  const std::vector<raw::ChannelID_t>& ubad,
19  const std::vector<raw::ChannelID_t>& vbad,
20  double distThresh,
21  double distThreshDrift,
22  double xhitOffset)
23  : geom(art::ServiceHandle<geo::Geometry const>()->provider())
24  , fDistThresh(distThresh)
25  , fDistThreshDrift(distThreshDrift)
26  , fXHitOffset(xhitOffset)
27  {
28  FillHitMap(detProp, xhits, fX_by_tpc);
29  FillHitMap(detProp, uhits, fU_by_tpc);
30  FillHitMap(detProp, vhits, fV_by_tpc);
31 
32  FillBadMap(xbad, fXbad_by_tpc);
33  FillBadMap(ubad, fUbad_by_tpc);
34  FillBadMap(vbad, fVbad_by_tpc);
35  }
36 
37  // -------------------------------------------------------------------------
40  std::map<geo::TPCID, std::vector<HitOrChan>>& out)
41  {
42  for (const art::Ptr<recob::Hit>& hit : hits) {
43  for (geo::TPCID tpc : geom->ROPtoTPCs(geom->ChannelToROP(hit->Channel()))) {
44  double xpos = 0;
45  for (geo::WireID wire : geom->ChannelToWire(hit->Channel())) {
46  if (geo::TPCID(wire) == tpc) {
47  xpos = detProp.ConvertTicksToX(hit->PeakTime(), wire);
48  if (geom->SignalType(wire) == geo::kCollection) xpos += fXHitOffset;
49  }
50  }
51 
52  out[tpc].emplace_back(hit.get(), xpos);
53  }
54  }
55  for (auto& it : out)
56  std::sort(it.second.begin(), it.second.end(), [](auto a, auto b) { return a.xpos < b.xpos; });
57  }
58 
59  // -------------------------------------------------------------------------
60  void TripletFinder::FillBadMap(const std::vector<raw::ChannelID_t>& bads,
61  std::map<geo::TPCID, std::vector<raw::ChannelID_t>>& out)
62  {
63  for (raw::ChannelID_t chan : bads) {
64  for (geo::TPCID tpc : geom->ROPtoTPCs(geom->ChannelToROP(chan))) {
65  out[tpc].push_back(chan);
66  }
67  }
68  }
69 
70  // -------------------------------------------------------------------------
72  public:
74  : geom(art::ServiceHandle<geo::Geometry const>()->provider()), fTPC(tpc)
75  {}
76 
78  {
79  const auto key = std::make_pair(a, b);
80 
81  auto it = fMap.find(key);
82  if (it != fMap.end()) {
83  pt = fPtMap[key];
84  return it->second;
85  }
86 
87  const bool res = ISect(a, b, pt);
88  fMap.insert({key, res});
89  fPtMap.insert({key, pt});
90  return res;
91  }
92 
93  protected:
95  {
96  for (geo::WireID awire : geom->ChannelToWire(chanA)) {
97  if (geo::TPCID(awire) != fTPC) continue;
98  for (geo::WireID bwire : geom->ChannelToWire(chanB)) {
99  if (geo::TPCID(bwire) != fTPC) continue;
100 
101  if (geom->WireIDsIntersect(awire, bwire, pt)) return true;
102  }
103  }
104 
105  return false;
106  }
107 
109 
110  std::map<std::pair<raw::ChannelID_t, raw::ChannelID_t>, bool> fMap;
111  std::map<std::pair<raw::ChannelID_t, raw::ChannelID_t>, geo::WireIDIntersection> fPtMap;
112 
114  };
115 
116  // -------------------------------------------------------------------------
117  bool TripletFinder::CloseDrift(double xa, double xb) const
118  {
119  return fabs(xa - xb) < fDistThreshDrift;
120  }
121 
122  // -------------------------------------------------------------------------
124  {
125  const TVector3 pa(ra.y, ra.z, 0);
126  const TVector3 pb(rb.y, rb.z, 0);
127 
128  return (pa - pb).Mag() < fDistThresh;
129  }
130 
131  bool LessThanXHit(const ChannelDoublet& a, const ChannelDoublet& b)
132  {
133  // Make sure the bad hits get sorted too
134  if (a.a.hit == 0 && b.a.hit == 0) return a.a.chan < b.a.chan;
135  // But mostly just order the real hits (in some arbitrary order)
136  return a.a.hit < b.a.hit;
137  }
138 
139  bool SameXHit(const ChannelDoublet& a, const ChannelDoublet& b)
140  {
141  if (a.a.hit == 0 && b.a.hit == 0) return a.a.chan == b.a.chan;
142  return a.a.hit == b.a.hit;
143  }
144 
145  // -------------------------------------------------------------------------
146  std::vector<HitTriplet> TripletFinder::Triplets()
147  {
148  std::vector<HitTriplet> ret;
149 
150  for (const auto& it : fX_by_tpc) {
151  const geo::TPCID& tpc = it.first;
152 
153  std::vector<ChannelDoublet> xus = DoubletsXU(tpc);
154  std::vector<ChannelDoublet> xvs = DoubletsXV(tpc);
155 
156  // Cache to prevent repeating the same questions
157  IntersectionCache isectUV(tpc);
158 
159  // For the efficient looping below to work we need to sort the doublet
160  // lists so the X hits occur in the same order.
161  std::sort(xus.begin(), xus.end(), LessThanXHit);
162  std::sort(xvs.begin(), xvs.end(), LessThanXHit);
163 
164  auto xvit_begin = xvs.begin();
165 
166  int nxuv = 0;
167  for (const ChannelDoublet& xu : xus) {
168  const HitOrChan& x = xu.a;
169  const HitOrChan& u = xu.b;
170 
171  // Catch up until we're looking at the same X hit in XV
172  while (xvit_begin != xvs.end() && LessThanXHit(*xvit_begin, xu))
173  ++xvit_begin;
174 
175  // Loop through all those matching hits
176  for (auto xvit = xvit_begin; xvit != xvs.end() && SameXHit(*xvit, xu); ++xvit) {
177  const HitOrChan& v = xvit->b;
178 
179  // Only allow one bad channel per triplet
180  if (!x.hit && !u.hit) continue;
181  if (!x.hit && !v.hit) continue;
182  if (!u.hit && !v.hit) continue;
183 
184  if (u.hit && v.hit && !CloseDrift(u.xpos, v.xpos)) continue;
185 
187  if (!isectUV(u.chan, v.chan, ptUV)) continue;
188 
189  if (!CloseSpace(xu.pt, xvit->pt) || !CloseSpace(xu.pt, ptUV) ||
190  !CloseSpace(xvit->pt, ptUV))
191  continue;
192 
193  double xavg = 0;
194  int nx = 0;
195  if (x.hit) {
196  xavg += x.xpos;
197  ++nx;
198  }
199  if (u.hit) {
200  xavg += u.xpos;
201  ++nx;
202  }
203  if (v.hit) {
204  xavg += v.xpos;
205  ++nx;
206  }
207  xavg /= nx;
208 
209  const XYZ pt{
210  xavg, (xu.pt.y + xvit->pt.y + ptUV.y) / 3, (xu.pt.z + xvit->pt.z + ptUV.z) / 3};
211 
212  ret.emplace_back(HitTriplet{x.hit, u.hit, v.hit, pt});
213  ++nxuv;
214  } // end for xv
215  } // end for xu
216 
217  std::cout << tpc << " " << xus.size() << " XUs and " << xvs.size() << " XVs -> " << nxuv
218  << " XUVs" << std::endl;
219 
220  } // end for tpc
221 
222  std::cout << ret.size() << " XUVs total" << std::endl;
223 
224  return ret;
225  }
226 
227  // -------------------------------------------------------------------------
228  std::vector<HitTriplet> TripletFinder::TripletsTwoView()
229  {
230  std::vector<HitTriplet> ret;
231 
232  for (const auto& it : fX_by_tpc) {
233  const geo::TPCID& tpc = it.first;
234 
235  std::vector<ChannelDoublet> xus = DoubletsXU(tpc);
236 
237  for (const ChannelDoublet& xu : xus) {
238  const HitOrChan& x = xu.a;
239  const HitOrChan& u = xu.b;
240 
241  double xavg = x.xpos;
242  int nx = 1;
243  if (u.hit) {
244  xavg += u.xpos;
245  ++nx;
246  }
247  xavg /= nx;
248 
249  const XYZ pt{xavg, xu.pt.y, xu.pt.z};
250 
251  ret.emplace_back(HitTriplet{x.hit, u.hit, 0, pt});
252  } // end for xu
253  } // end for tpc
254 
255  std::cout << ret.size() << " XUs total" << std::endl;
256 
257  return ret;
258  }
259 
260  // -------------------------------------------------------------------------
261  std::vector<ChannelDoublet> TripletFinder::DoubletsXU(geo::TPCID tpc)
262  {
263  std::vector<ChannelDoublet> ret =
264  DoubletHelper(tpc, fX_by_tpc[tpc], fU_by_tpc[tpc], fUbad_by_tpc[tpc]);
265 
266  // Find X(bad)+U(good) doublets, have to flip them for the final result
267  for (auto it : DoubletHelper(tpc, fU_by_tpc[tpc], {}, fXbad_by_tpc[tpc])) {
268  ret.push_back({it.b, it.a, it.pt});
269  }
270 
271  return ret;
272  }
273 
274  // -------------------------------------------------------------------------
275  std::vector<ChannelDoublet> TripletFinder::DoubletsXV(geo::TPCID tpc)
276  {
277  std::vector<ChannelDoublet> ret =
278  DoubletHelper(tpc, fX_by_tpc[tpc], fV_by_tpc[tpc], fVbad_by_tpc[tpc]);
279 
280  // Find X(bad)+V(good) doublets, have to flip them for the final result
281  for (auto it : DoubletHelper(tpc, fV_by_tpc[tpc], {}, fXbad_by_tpc[tpc])) {
282  ret.push_back({it.b, it.a, it.pt});
283  }
284 
285  return ret;
286  }
287 
288  // -------------------------------------------------------------------------
289  std::vector<ChannelDoublet> TripletFinder::DoubletHelper(
290  geo::TPCID tpc,
291  const std::vector<HitOrChan>& ahits,
292  const std::vector<HitOrChan>& bhits,
293  const std::vector<raw::ChannelID_t>& bbads) const
294  {
295  std::vector<ChannelDoublet> ret;
296 
297  IntersectionCache isect(tpc);
298 
299  auto b_begin = bhits.begin();
300 
301  for (const HitOrChan& a : ahits) {
302  // Bad channels are easy because there's no timing constraint
303  for (raw::ChannelID_t b : bbads) {
305  if (isect(a.chan, b, pt)) { ret.emplace_back(a, b, pt); }
306  }
307 
308  while (b_begin != bhits.end() && b_begin->xpos < a.xpos && !CloseDrift(b_begin->xpos, a.xpos))
309  ++b_begin;
310 
311  for (auto bit = b_begin; bit != bhits.end(); ++bit) {
312  const HitOrChan& b = *bit;
313 
314  if (b.xpos > a.xpos && !CloseDrift(b.xpos, a.xpos)) break;
315 
317  if (!isect(a.chan, b.chan, pt)) continue;
318 
319  ret.emplace_back(a, b, pt);
320  } // end for b
321  } // end for a
322 
323  return ret;
324  }
325 }
Float_t x
Definition: compare.C:6
bool CloseSpace(geo::WireIDIntersection ra, geo::WireIDIntersection rb) const
raw::ChannelID_t chan
Definition: TripletFinder.h:30
double z
z position of intersection
Definition: geo_types.h:792
std::vector< WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
std::vector< ChannelDoublet > DoubletHelper(geo::TPCID tpc, const std::vector< HitOrChan > &ahits, const std::vector< HitOrChan > &bhits, const std::vector< raw::ChannelID_t > &bbads) const
IntersectionCache(geo::TPCID tpc)
std::map< std::pair< raw::ChannelID_t, raw::ChannelID_t >, geo::WireIDIntersection > fPtMap
std::map< geo::TPCID, std::vector< HitOrChan > > fX_by_tpc
Definition: TripletFinder.h:94
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, Point_t &intersection) const
Computes the intersection between two wires.
const geo::GeometryCore * geom
std::map< std::pair< raw::ChannelID_t, raw::ChannelID_t >, bool > fMap
std::map< geo::TPCID, std::vector< HitOrChan > > fV_by_tpc
Definition: TripletFinder.h:96
std::map< geo::TPCID, std::vector< raw::ChannelID_t > > fXbad_by_tpc
Definition: TripletFinder.h:99
std::vector< TPCID > ROPtoTPCs(readout::ROPID const &ropid) const
Returns a list of ID of TPCs the specified ROP spans.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
Access the description of detector geometry.
void hits()
Definition: readHits.C:15
std::map< geo::TPCID, std::vector< HitOrChan > > fU_by_tpc
Definition: TripletFinder.h:95
const geo::GeometryCore * geom
Definition: TripletFinder.h:69
TMarker * pt
Definition: egs.C:25
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
Description of geometry of one entire detector.
Definition: GeometryCore.h:119
bool SameXHit(const ChannelDoublet &a, const ChannelDoublet &b)
Detector simulation of raw signals on wires.
double ConvertTicksToX(double ticks, int p, int t, int c) const
readout::ROPID ChannelToROP(raw::ChannelID_t channel) const
TripletFinder(const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Hit >> &xhits, const std::vector< art::Ptr< recob::Hit >> &uhits, const std::vector< art::Ptr< recob::Hit >> &vhits, const std::vector< raw::ChannelID_t > &xbad, const std::vector< raw::ChannelID_t > &ubad, const std::vector< raw::ChannelID_t > &vbad, double distThresh, double distThreshDrift, double xhitOffset)
std::vector< HitTriplet > TripletsTwoView()
Only search for XU intersections.
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
const recob::Hit * hit
Definition: TripletFinder.h:31
std::vector< ChannelDoublet > DoubletsXV(geo::TPCID tpc)
std::vector< ChannelDoublet > DoubletsXU(geo::TPCID tpc)
double y
y position of intersection
Definition: geo_types.h:791
bool CloseDrift(double xa, double xb) const
std::vector< HitTriplet > Triplets()
bool LessThanXHit(const ChannelDoublet &a, const ChannelDoublet &b)
Definition: MVAAlg.h:12
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
std::map< geo::TPCID, std::vector< raw::ChannelID_t > > fUbad_by_tpc
void FillBadMap(const std::vector< raw::ChannelID_t > &bads, std::map< geo::TPCID, std::vector< raw::ChannelID_t >> &out)
Helper for constructor.
bool operator()(raw::ChannelID_t a, raw::ChannelID_t b, geo::WireIDIntersection &pt)
Namespace collecting geometry-related classes utilities.
bool ISect(raw::ChannelID_t chanA, raw::ChannelID_t chanB, geo::WireIDIntersection &pt) const
art framework interface to geometry description
std::map< geo::TPCID, std::vector< raw::ChannelID_t > > fVbad_by_tpc
void FillHitMap(const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits, std::map< geo::TPCID, std::vector< HitOrChan >> &out)
Helper for constructor.
Signal from collection planes.
Definition: geo_types.h:152