LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
QuadVtx_module.cc
Go to the documentation of this file.
1 // Chris Backhouse - c.backhouse@ucl.ac.uk - Oct 2019
2 
4 
5 // C/C++ standard libraries
6 #include <iostream>
7 #include <random>
8 #include <string>
9 
10 // framework libraries
16 #include "art_root_io/TFileService.h"
18 #include "cetlib/pow.h"
19 #include "fhiclcpp/ParameterSet.h"
21 
22 // LArSoft libraries
27 
28 #include "TGraph.h"
29 #include "TH2F.h"
30 #include "TMatrixD.h"
31 #include "TVectorD.h"
32 
33 namespace quad {
34 
35  // ---------------------------------------------------------------------------
36  struct Pt2D {
37  Pt2D(double _x, double _z, int _view, double _energy)
38  : x(_x), z(_z), view(_view), energy(_energy)
39  {}
40 
41  bool operator<(const Pt2D& p) const { return z < p.z; }
42 
43  double x, z;
44  int view;
45  double energy;
46  };
47 
48  // ---------------------------------------------------------------------------
49  struct Line2D {
50  Line2D(const Pt2D& a, const Pt2D& b)
51  : m((b.x - a.x) / (b.z - a.z))
52  , c(b.x - m * b.z)
53  , // w(a.energy * b.energy),
54  minz(std::min(a.z, b.z))
55  , maxz(std::max(a.z, b.z))
56  {
57  assert(a.z != b.z); // no vertical lines
58  }
59 
60  bool operator<(const Line2D& l) const { return m < l.m; }
61 
62  float m, c, /*w,*/ minz, maxz;
63  };
64 
65  // ---------------------------------------------------------------------------
66  class QuadVtx : public art::EDProducer {
67  public:
68  explicit QuadVtx(const fhicl::ParameterSet& pset);
69 
70  void beginJob() override;
71  void produce(art::Event& evt) override;
72 
73  private:
74  bool FindVtx(const detinfo::DetectorPropertiesData& detProp,
75  const std::vector<recob::Hit>& hits,
77  int evt) const;
78 
79  std::string fHitLabel;
80 
81  bool fSavePlots;
82 
84  };
85 
87 
88  // ---------------------------------------------------------------------------
89  QuadVtx::QuadVtx(const fhicl::ParameterSet& pset)
90  : EDProducer(pset)
91  , fHitLabel(pset.get<std::string>("HitLabel"))
92  , fSavePlots(pset.get<bool>("SavePlots"))
93  {
94  produces<std::vector<recob::Vertex>>();
95  }
96 
97  // ---------------------------------------------------------------------------
99  {
100  geom = art::ServiceHandle<const geo::Geometry>()->provider();
101  }
102 
103  // ---------------------------------------------------------------------------
104  // x = m*z+c. z1 and z2 are the two intercepts in case of returning true
105  bool IntersectsCircle(float m, float c, float z0, float x0, float R, float& z1, float& z2)
106  {
107  // Change to the frame where (z0, x0) = (0, 0)
108  c += m * z0 - x0;
109 
110  // z^2 + (m*z+c)^2 = R^2
111  const float A = 1 + cet::square(m);
112  const float B = 2 * m * c;
113  const float C = cet::square(c) - cet::square(R);
114 
115  double desc = cet::square(B) - 4 * A * C;
116 
117  if (desc < 0) return false;
118 
119  desc = sqrt(desc);
120 
121  z1 = (-B - desc) / (2 * A);
122  z2 = (-B + desc) / (2 * A);
123 
124  // Back to the original frame
125  z1 += z0;
126  z2 += z0;
127 
128  return true;
129  }
130 
131  // ---------------------------------------------------------------------------
132  void LinesFromPoints(const std::vector<Pt2D>& pts,
133  std::vector<Line2D>& lines,
134  float z0 = 0,
135  float x0 = 0,
136  float R = -1)
137  {
138  constexpr size_t kMaxLines = 10 * 1000 * 1000; // This is 150MB of lines...
139 
140  const size_t product = (pts.size() * (pts.size() - 1)) / 2;
141  const int stride = product / kMaxLines + 1;
142 
143  lines.reserve(std::min(product, kMaxLines));
144 
145  for (int offset = 0; offset < stride; ++offset) {
146  for (unsigned int i = 0; i < pts.size(); ++i) {
147  for (unsigned int j = i + offset + 1; j < pts.size(); j += stride) {
148  const Line2D l(pts[i], pts[j]);
149 
150  if (isinf(l.m) || isnan(l.m) || isinf(l.c) || isnan(l.c)) continue;
151 
152  if (R > 0) {
153  float z1, z2;
154  if (!IntersectsCircle(l.m, l.c, z0, x0, 2.5, z1, z2)) continue;
155  if (l.minz < z1 && l.minz < z2 && l.maxz > z1 && l.maxz > z2) continue;
156  }
157 
158  lines.push_back(std::move(l));
159  if (lines.size() == kMaxLines) goto end; // break out of 3 loops
160  }
161  }
162  }
163 
164  end:
165 
166  lines.shrink_to_fit();
167 
168  mf::LogInfo() << "Made " << lines.size() << " lines using stride " << stride
169  << " to fit under cap of " << kMaxLines << std::endl;
170 
171  // Lines are required to be sorted by gradient for a later optimization
172  std::sort(lines.begin(), lines.end());
173  }
174 
175  // ---------------------------------------------------------------------------
176  inline bool CloseAngles(float ma, float mb)
177  {
178  const float cosCrit = cos(10 * M_PI / 180.);
179  const float dot = 1 + ma * mb; // (1, ma)*(1, mb)
180  return cet::square(dot) > (1 + cet::square(ma)) * (1 + cet::square(mb)) * cet::square(cosCrit);
181  }
182 
183  // ---------------------------------------------------------------------------
184  void MapFromLines(const std::vector<Line2D>& lines, HeatMap& hm)
185  {
186  // This maximum is driven by runtime
187  constexpr size_t kMaxPts = 10 * 1000 * 1000;
188 
189  unsigned int j0 = 0;
190  unsigned int jmax = 0;
191 
192  long npts = 0;
193  for (unsigned int i = 0; i + 1 < lines.size(); ++i) {
194  const Line2D a = lines[i];
195 
196  j0 = std::max(j0, i + 1);
197  while (j0 < lines.size() && CloseAngles(a.m, lines[j0].m))
198  ++j0;
199  jmax = std::max(jmax, j0);
200  while (jmax < lines.size() && !CloseAngles(a.m, lines[jmax].m))
201  ++jmax;
202 
203  npts += jmax - j0;
204  }
205 
206  const size_t product = (lines.size() * (lines.size() - 1)) / 2;
207  const int stride = npts / kMaxPts + 1;
208 
209  mf::LogInfo() << "Combining lines to points with stride " << stride << std::endl;
210 
211  mf::LogInfo() << npts << " cf " << product << " ie " << double(npts) / product << std::endl;
212 
213  j0 = 0;
214  jmax = 0;
215 
216  for (unsigned int i = 0; i + 1 < lines.size(); ++i) {
217  const Line2D a = lines[i];
218 
219  j0 = std::max(j0, i + 1);
220  while (j0 < lines.size() && CloseAngles(a.m, lines[j0].m))
221  ++j0;
222  jmax = std::max(jmax, j0);
223  while (jmax < lines.size() && !CloseAngles(a.m, lines[jmax].m))
224  ++jmax;
225 
226  for (unsigned int j = j0; j < jmax; j += stride) {
227  const Line2D& b = lines[j];
228 
229  // x = mA * z + cA = mB * z + cB
230  const float z = (b.c - a.c) / (a.m - b.m);
231  const float x = a.m * z + a.c;
232 
233  // No solutions within a line
234  if ((z < a.minz || z > a.maxz) && (z < b.minz || z > b.maxz)) {
235  const int iz = hm.ZToBin(z);
236  const int ix = hm.XToBin(x);
237  if (iz >= 0 && iz < hm.Nz && ix >= 0 && ix < hm.Nx) { hm.map[iz * hm.Nx + ix] += stride; }
238  }
239  } // end for i
240  }
241  }
242 
243  // ---------------------------------------------------------------------------
244  // Assumes that all three maps have the same vertical stride
245  recob::tracking::Point_t FindPeak3D(const std::vector<HeatMap>& hs,
246  const std::vector<recob::tracking::Vector_t>& dirs) noexcept
247  {
248  assert(hs.size() == 3);
249  assert(dirs.size() == 3);
250 
251  const int Nx = hs[0].Nx;
252 
253  TMatrixD M(2, 2);
254  M(0, 0) = dirs[0].Y();
255  M(0, 1) = dirs[0].Z();
256  M(1, 0) = dirs[1].Y();
257  M(1, 1) = dirs[1].Z();
258 
259  // Singular, and stupid setup of exceptions means we can't test any other way
260  if (M(0, 0) * M(1, 1) - M(1, 0) * M(0, 1) == 0) return {};
261 
262  M.Invert();
263 
264  float bestscore = -1;
266 
267  // Accumulate some statistics up front that will enable us to optimize
268  std::vector<float> colMax[3];
269  for (int view = 0; view < 3; ++view) {
270  colMax[view].resize(hs[view].Nz);
271  for (int iz = 0; iz < hs[view].Nz; ++iz) {
272  colMax[view][iz] = *std::max_element(&hs[view].map[Nx * iz], &hs[view].map[Nx * (iz + 1)]);
273  }
274  }
275 
276  for (int iz = 0; iz < hs[0].Nz; ++iz) {
277  const float z = hs[0].ZBinCenter(iz);
278  const float bonus = 1; // works badly... exp((hs[0].maxz-z)/1000.);
279 
280  for (int iu = 0; iu < hs[1].Nz; ++iu) {
281  const float u = hs[1].ZBinCenter(iu);
282  // r.Dot(d0) = z && r.Dot(d1) = u
283  TVectorD p(2);
284  p(0) = z;
285  p(1) = u;
286  const TVectorD r = M * p;
287  const float v = r[0] * dirs[2].Y() + r[1] * dirs[2].Z();
288  const int iv = hs[2].ZToBin(v);
289  if (iv < 0 || iv >= hs[2].Nz) continue;
290  const double y = r(0);
291 
292  // Even if the maxes were all at the same x we couldn't beat the record
293  if (colMax[0][iz] + colMax[1][iu] + colMax[2][iv] < bestscore) continue;
294 
295  const float* h0 = &hs[0].map[Nx * iz];
296  const float* h1 = &hs[1].map[Nx * iu];
297  const float* h2 = &hs[2].map[Nx * iv];
298 
299  int bestix = -1;
300  for (int ix = 1; ix < Nx - 1; ++ix) {
301  const float score = bonus * (h0[ix] + h1[ix] + h2[ix]);
302 
303  if (score > bestscore) {
304  bestscore = score;
305  bestix = ix;
306  }
307  } // end for dx
308 
309  if (bestix != -1) { bestr.SetXYZ(hs[0].XBinCenter(bestix), y, z); }
310  } // end for u
311  } // end for z
312 
313  return bestr;
314  }
315 
316  // ---------------------------------------------------------------------------
318  const std::vector<recob::Hit>& hits,
319  std::vector<std::vector<Pt2D>>& pts,
320  std::vector<recob::tracking::Vector_t>& dirs,
321  const geo::GeometryCore* geom)
322  {
323  pts.resize(3); // 3 views
324 
325  recob::tracking::Vector_t dirZ(0, 0, 1);
326  recob::tracking::Vector_t dirU, dirV;
327 
328  for (const recob::Hit& hit : hits) {
329  const geo::WireID wire = hit.WireID();
330 
331  const double xpos = detProp.ConvertTicksToX(hit.PeakTime(), wire);
332 
333  auto const r0 = geom->WireEndPoints(wire).start();
334  auto const r1 = geom->WireEndPoints(wire).end();
335 
336  const double energy = hit.Integral();
337 
338  if (geom->View(hit.Channel()) == geo::kZ) {
339  pts[0].emplace_back(xpos, r0.z(), 0, energy);
340  continue;
341  }
342 
343  // Compute the direction perpendicular to the wires
344  auto const unit = (r1 - r0).Unit();
345  recob::tracking::Vector_t perp{0, -unit.z(), unit.y()};
346  // We want to ultimately have a positive z component ion "perp"
347  if (perp.z() < 0) perp *= -1;
348 
349  // TODO check we never get a 4th view a-la the bug we had in the 3D version
350 
351  // The "U" direction is the first one we see
352  if (dirU.Mag2() == 0) { dirU = perp; }
353  else if (dirV.Mag2() == 0 && fabs(dirU.Dot(perp)) < 0.99) {
354  // If we still need a "V" and this direction differs from "U"
355  dirV = perp;
356  }
357 
358  // Hits belong to whichever view their perpendicular vector aligns with
359  if (fabs(dirU.Dot(perp)) > 0.99) { pts[1].emplace_back(xpos, r0.Dot(dirU), 1, energy); }
360  else {
361  pts[2].emplace_back(xpos, r0.Dot(dirV), 2, energy);
362  }
363  } // end for hits
364 
365  dirs = {dirZ, dirU, dirV};
366 
367  std::default_random_engine gen;
368 
369  // In case we need to sub-sample they should be shuffled
370  for (int view = 0; view < 3; ++view) {
371  std::shuffle(pts[view].begin(), pts[view].end(), gen);
372  }
373  }
374 
375  // ---------------------------------------------------------------------------
377  const std::vector<recob::Hit>& hits,
379  int evt) const
380  {
381  if (hits.empty()) return false;
382 
383  std::vector<std::vector<Pt2D>> pts;
384  std::vector<recob::tracking::Vector_t> dirs;
385 
386  GetPts2D(detProp, hits, pts, dirs, geom);
387 
388  double minx = +1e9;
389  double maxx = -1e9;
390  double minz[3] = {+1e9, +1e9, +1e9};
391  double maxz[3] = {-1e9, -1e9, -1e9};
392  for (int view = 0; view < 3; ++view) {
393  for (const Pt2D& p : pts[view]) {
394  minx = std::min(minx, p.x);
395  maxx = std::max(maxx, p.x);
396  minz[view] = std::min(minz[view], p.z);
397  maxz[view] = std::max(maxz[view], p.z);
398  }
399  }
400 
401  // Add some padding
402  for (int view = 0; view < 3; ++view) {
403  minz[view] -= 100;
404  maxz[view] += 100;
405  }
406  minx -= 20;
407  maxx += 20;
408 
409  // Don't allow the vertex further downstream in z (view 0) than 25% of the
410  // hits.
411  std::vector<float> zs;
412  zs.reserve(pts[0].size());
413  for (const Pt2D& p : pts[0])
414  zs.push_back(p.z);
415  auto mid = zs.begin() + zs.size() / 4;
416  if (mid != zs.end()) {
417  std::nth_element(zs.begin(), mid, zs.end());
418  maxz[0] = *mid;
419  }
420 
421  std::vector<HeatMap> hms;
422  hms.reserve(3);
423  for (int view = 0; view < 3; ++view) {
424  if (pts[view].empty()) return false;
425 
426  std::vector<Line2D> lines;
427  LinesFromPoints(pts[view], lines);
428 
429  if (lines.empty()) return false;
430 
431  // Approximately cm bins
432  hms.emplace_back(maxz[view] - minz[view], minz[view], maxz[view], maxx - minx, minx, maxx);
433  MapFromLines(lines, hms.back());
434  } // end for view
435 
436  vtx = FindPeak3D(hms, dirs);
437 
438  std::vector<HeatMap> hms_zoom;
439  hms_zoom.reserve(3);
440  for (int view = 0; view < 3; ++view) {
441  const double x0 = vtx.X();
442  const double z0 = vtx.Dot(dirs[view]);
443 
444  std::vector<Line2D> lines;
445  LinesFromPoints(pts[view], lines, z0, x0, 2.5);
446 
447  if (lines.empty()) return false; // How does this happen??
448 
449  // mm granularity
450  hms_zoom.emplace_back(50, z0 - 2.5, z0 + 2.5, 50, x0 - 2.5, x0 + 2.5);
451 
452  MapFromLines(lines, hms_zoom.back());
453  }
454 
455  vtx = FindPeak3D(hms_zoom, dirs);
456 
457  if (fSavePlots) {
458  art::TFileDirectory evt_dir =
459  art::ServiceHandle<art::TFileService>()->mkdir(TString::Format("evt%d", evt).Data());
460 
461  for (int view = 0; view < 3; ++view) {
462  art::TFileDirectory view_dir = evt_dir.mkdir(TString::Format("view%d", view).Data());
463 
464  TGraph* gpts = view_dir.makeAndRegister<TGraph>("hits", "");
465  for (const Pt2D& p : pts[view])
466  gpts->SetPoint(gpts->GetN(), p.z, p.x);
467 
468  view_dir.makeAndRegister<TH2F>("hmap", "", *hms[view].AsTH2());
469 
470  view_dir.makeAndRegister<TH2F>("hmap_zoom", "", *hms_zoom[view].AsTH2());
471 
472  const double x = vtx.X();
473  const double z = vtx.Dot(dirs[view]);
474  view_dir.makeAndRegister<TGraph>("vtx3d", "", 1, &z, &x);
475  } // end for view
476  } // end if saving plots
477 
478  return true;
479  }
480 
481  // ---------------------------------------------------------------------------
483  {
484  auto vtxcol = std::make_unique<std::vector<recob::Vertex>>();
485 
487  evt.getByLabel(fHitLabel, hits);
488 
489  auto const detProp =
492  if (FindVtx(detProp, *hits, vtx, evt.event())) {
493  vtxcol->emplace_back(vtx, recob::Vertex::SMatrixSym33(), 0, 0);
494  }
495 
496  evt.put(std::move(vtxcol));
497  }
498 
499 } // end namespace quad
TRandom r
Definition: spectrum.C:23
int XToBin(double x) const
Definition: HeatMap.h:19
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Declaration of signal hit object.
Float_t y
Definition: compare.C:6
bool FindVtx(const detinfo::DetectorPropertiesData &detProp, const std::vector< recob::Hit > &hits, recob::tracking::Point_t &vtx, int evt) const
tracking::SMatrixSym33 SMatrixSym33
Definition: Vertex.h:39
int ZToBin(double z) const
Definition: HeatMap.h:18
STL namespace.
Planes which measure Z direction.
Definition: geo_types.h:138
void GetPts2D(const detinfo::DetectorPropertiesData &detProp, const std::vector< recob::Hit > &hits, std::vector< std::vector< Pt2D >> &pts, std::vector< recob::tracking::Vector_t > &dirs, const geo::GeometryCore *geom)
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
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
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
tracking::Point_t Point_t
Definition: Vertex.h:38
void hits()
Definition: readHits.C:15
const int Nx
Definition: HeatMap.h:22
Pt2D(double _x, double _z, int _view, double _energy)
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void beginJob()
Definition: Breakpoints.cc:14
void MapFromLines(const std::vector< Line2D > &lines, HeatMap &hm)
parameter set interface
std::string fHitLabel
bool IntersectsCircle(float m, float c, float z0, float x0, float R, float &z1, float &z2)
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< Coord_t >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space. See recob::tracking::Coord_t for more details on the ...
Definition: TrackingTypes.h:31
bool operator<(const Pt2D &p) const
recob::tracking::Point_t FindPeak3D(const std::vector< HeatMap > &hs, const std::vector< recob::tracking::Vector_t > &dirs) noexcept
EventNumber_t event() const
Definition: Event.cc:41
constexpr auto dot(Vector const &a, OtherVector const &b)
Return cross product of two vectors.
Description of geometry of one entire detector.
Definition: GeometryCore.h:119
Detector simulation of raw signals on wires.
TH1F * h2
Definition: plot.C:44
double ConvertTicksToX(double ticks, int p, int t, int c) const
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
bool CloseAngles(float ma, float mb)
void beginJob() override
bool operator<(const Line2D &l) const
View_t View(PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
void LinesFromPoints(const std::vector< Pt2D > &pts, std::vector< Line2D > &lines, float z0=0, float x0=0, float R=-1)
TH1F * h1
Definition: plot.C:41
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
TCEvent evt
Definition: DataStructs.cxx:8
const geo::GeometryCore * geom
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< Coord_t >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space. See recob::tracking::Coord_t for more detai...
Definition: TrackingTypes.h:27
art framework interface to geometry description
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
void WireEndPoints(WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
void produce(art::Event &evt) override
Line2D(const Pt2D &a, const Pt2D &b)
std::vector< float > map
Definition: HeatMap.h:24