Classes | Functions
quad Namespace Reference

Classes

class  EvalVtx
 
class  HeatMap
 
struct  Line2D
 
struct  Pt2D
 
class  QuadVtx
 

Functions

recob::Vertex GetFirstVertex (const std::string &label, const art::Event &evt)
 
recob::Vertex GetVtxByAssns (const std::string &label, const art::Event &evt)
 
bool IntersectsCircle (float m, float c, float z0, float x0, float R, float &z1, float &z2)
 
void LinesFromPoints (const std::vector< Pt2D > &pts, std::vector< Line2D > &lines, float z0=0, float x0=0, float R=-1)
 
bool CloseAngles (float ma, float mb)
 
void MapFromLines (const std::vector< Line2D > &lines, HeatMap &hm)
 
TVector3 FindPeak3D (const std::vector< HeatMap > &hs, const std::vector< TVector3 > &dirs) noexcept
 
void GetPts2D (const detinfo::DetectorPropertiesData &detProp, const std::vector< recob::Hit > &hits, std::vector< std::vector< Pt2D >> &pts, std::vector< TVector3 > &dirs, const geo::GeometryCore *geom)
 

Function Documentation

bool quad::CloseAngles ( float  ma,
float  mb 
)
inline

Definition at line 188 of file QuadVtx_module.cc.

189  {
190  const float cosCrit = cos(10 * M_PI / 180.);
191  const float dot = 1 + ma * mb; // (1, ma)*(1, mb)
192  return cet::square(dot) > (1 + cet::square(ma)) * (1 + cet::square(mb)) * cet::square(cosCrit);
193  }
constexpr T square(T x)
Definition: pow.h:21
static constexpr double mb
Definition: Units.h:79
#define M_PI
Definition: includeROOT.h:54
TVector3 quad::FindPeak3D ( const std::vector< HeatMap > &  hs,
const std::vector< TVector3 > &  dirs 
)
noexcept

Definition at line 259 of file QuadVtx_module.cc.

260  {
261  assert(hs.size() == 3);
262  assert(dirs.size() == 3);
263 
264  const int Nx = hs[0].Nx;
265 
266  TMatrixD M(2, 2);
267  M(0, 0) = dirs[0].Y();
268  M(0, 1) = dirs[0].Z();
269  M(1, 0) = dirs[1].Y();
270  M(1, 1) = dirs[1].Z();
271 
272  // Singular, and stupid setup of exceptions means we can't test any other way
273  if (M(0, 0) * M(1, 1) - M(1, 0) * M(0, 1) == 0) return TVector3(0, 0, 0);
274 
275  M.Invert();
276 
277  float bestscore = -1;
278  TVector3 bestr;
279 
280  // Accumulate some statistics up front that will enable us to optimize
281  std::vector<float> colMax[3];
282  for (int view = 0; view < 3; ++view) {
283  colMax[view].resize(hs[view].Nz);
284  for (int iz = 0; iz < hs[view].Nz; ++iz) {
285  colMax[view][iz] = *std::max_element(&hs[view].map[Nx * iz], &hs[view].map[Nx * (iz + 1)]);
286  }
287  }
288 
289  for (int iz = 0; iz < hs[0].Nz; ++iz) {
290  const float z = hs[0].ZBinCenter(iz);
291  const float bonus = 1; // works badly... exp((hs[0].maxz-z)/1000.);
292 
293  for (int iu = 0; iu < hs[1].Nz; ++iu) {
294  const float u = hs[1].ZBinCenter(iu);
295  // r.Dot(d0) = z && r.Dot(d1) = u
296  TVectorD p(2);
297  p(0) = z;
298  p(1) = u;
299  const TVectorD r = M * p;
300  const float v = r[0] * dirs[2].Y() + r[1] * dirs[2].Z();
301  const int iv = hs[2].ZToBin(v);
302  if (iv < 0 || iv >= hs[2].Nz) continue;
303  const double y = r(0);
304 
305  // Even if the maxes were all at the same x we couldn't beat the record
306  if (colMax[0][iz] + colMax[1][iu] + colMax[2][iv] < bestscore) continue;
307 
308  // Attempt to micro-optimize the dx loop below
309  const float* __restrict__ h0 = &hs[0].map[Nx * iz];
310  const float* __restrict__ h1 = &hs[1].map[Nx * iu];
311  const float* __restrict__ h2 = &hs[2].map[Nx * iv];
312 
313  int bestix = -1;
314  for (int ix = 1; ix < Nx - 1; ++ix) {
315  const float score = bonus * (h0[ix] + h1[ix] + h2[ix]);
316 
317  if (score > bestscore) {
318  bestscore = score;
319  bestix = ix;
320  }
321  } // end for dx
322 
323  if (bestix != -1) { bestr = TVector3(hs[0].XBinCenter(bestix), y, z); }
324  } // end for u
325  } // end for z
326 
327  return bestr;
328  }
p
Definition: test.py:223
recob::Vertex quad::GetFirstVertex ( const std::string label,
const art::Event evt 
)

Definition at line 70 of file EvalVtx_module.cc.

71 {
72  const auto& vtxs = *evt.getValidHandle<std::vector<recob::Vertex>>(label);
73 
74  if(vtxs.empty()) return recob::Vertex();
75 
76  return vtxs[0];
77 }
Definition of vertex object for LArSoft.
Definition: Vertex.h:35
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
void quad::GetPts2D ( const detinfo::DetectorPropertiesData detProp,
const std::vector< recob::Hit > &  hits,
std::vector< std::vector< Pt2D >> &  pts,
std::vector< TVector3 > &  dirs,
const geo::GeometryCore geom 
)

Definition at line 332 of file QuadVtx_module.cc.

337  {
338  pts.resize(3); // 3 views
339 
340  TVector3 dirZ(0, 0, 1);
341  TVector3 dirU, dirV;
342 
343  for (const recob::Hit& hit : hits) {
344  const geo::WireID wire = hit.WireID();
345 
346  const double xpos = detProp.ConvertTicksToX(hit.PeakTime(), wire);
347 
348  const TVector3 r0 = geom->WireEndPoints(wire).start();
349  const TVector3 r1 = geom->WireEndPoints(wire).end();
350 
351  const double energy = hit.Integral();
352 
353  if (geom->View(hit.Channel()) == geo::kZ) {
354  pts[0].emplace_back(xpos, r0.z(), 0, energy);
355  continue;
356  }
357 
358  // Compute the direction perpendicular to the wires
359  TVector3 perp = (r1 - r0).Unit();
360  perp = TVector3(0, -perp.z(), perp.y());
361  // We want to ultimately have a positive z component in "perp"
362  if (perp.z() < 0) perp *= -1;
363 
364  // TODO check we never get a 4th view a-la the bug we had in the 3D version
365 
366  // The "U" direction is the first one we see
367  if (dirU.Mag2() == 0) { dirU = perp; }
368  else if (dirV.Mag2() == 0 && fabs(dirU.Dot(perp)) < 0.99) {
369  // If we still need a "V" and this direction differs from "U"
370  dirV = perp;
371  }
372 
373  // Hits belong to whichever view their perpendicular vector aligns with
374  if (fabs(dirU.Dot(perp)) > 0.99) { pts[1].emplace_back(xpos, r0.Dot(dirU), 1, energy); }
375  else {
376  pts[2].emplace_back(xpos, r0.Dot(dirV), 2, energy);
377  }
378  } // end for hits
379 
380  dirs = {dirZ, dirU, dirV};
381 
382  std::default_random_engine gen;
383 
384  // In case we need to sub-sample they should be shuffled
385  for (int view = 0; view < 3; ++view) {
386  std::shuffle(pts[view].begin(), pts[view].end(), gen);
387  }
388  }
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
TH3F * xpos
Definition: doAna.cpp:29
Planes which measure Z direction.
Definition: geo_types.h:132
gen
Definition: demo.py:24
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
Detector simulation of raw signals on wires.
double ConvertTicksToX(double ticks, int p, int t, int c) const
void WireEndPoints(geo::WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
recob::Vertex quad::GetVtxByAssns ( const std::string label,
const art::Event evt 
)

Definition at line 80 of file EvalVtx_module.cc.

81 {
83  evt.getByLabel(label, vtxs);
84 
85  if(vtxs->empty()) return recob::Vertex();
86 
88  evt.getByLabel(label, parts);
89 
90  art::FindManyP<recob::Vertex> fm(parts, evt, label);
91 
92  for(unsigned int i = 0; i < parts->size(); ++i){
93  const int pdg = abs((*parts)[i].PdgCode());
94  if(pdg == 12 || pdg == 14 || pdg == 16){
95  const std::vector<art::Ptr<recob::Vertex>>& vtxs = fm.at(i);
96  if(vtxs.size() == 1) return *vtxs[0];
97 
98  if(vtxs.empty()){
99  std::cout << "Warning: vertex list empty!" << std::endl;
100  }
101  else{
102  std::cout << "Warning: " << vtxs.size() << " vertices for daughter?" << std::endl;
103  }
104  }
105  }
106 
107  return recob::Vertex();
108 }
Definition of vertex object for LArSoft.
Definition: Vertex.h:35
T abs(T value)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
static constexpr double fm
Definition: Units.h:75
QTextStream & endl(QTextStream &s)
bool quad::IntersectsCircle ( float  m,
float  c,
float  z0,
float  x0,
float  R,
float &  z1,
float &  z2 
)

Definition at line 115 of file QuadVtx_module.cc.

116  {
117  // Change to the frame where (z0, x0) = (0, 0)
118  c += m * z0 - x0;
119 
120  // z^2 + (m*z+c)^2 = R^2
121  const float A = 1 + cet::square(m);
122  const float B = 2 * m * c;
123  const float C = cet::square(c) - cet::square(R);
124 
125  double desc = cet::square(B) - 4 * A * C;
126 
127  if (desc < 0) return false;
128 
129  desc = sqrt(desc);
130 
131  z1 = (-B - desc) / (2 * A);
132  z2 = (-B + desc) / (2 * A);
133 
134  // Back to the original frame
135  z1 += z0;
136  z2 += z0;
137 
138  return true;
139  }
constexpr T square(T x)
Definition: pow.h:21
void quad::LinesFromPoints ( const std::vector< Pt2D > &  pts,
std::vector< Line2D > &  lines,
float  z0 = 0,
float  x0 = 0,
float  R = -1 
)

Definition at line 143 of file QuadVtx_module.cc.

148  {
149  constexpr size_t kMaxLines = 10 * 1000 * 1000; // This is 150MB of lines...
150 
151  const size_t product = (pts.size() * (pts.size() - 1)) / 2;
152  const int stride = product / kMaxLines + 1;
153 
154  lines.reserve(std::min(product, kMaxLines));
155 
156  for (int offset = 0; offset < stride; ++offset) {
157  for (unsigned int i = 0; i < pts.size(); ++i) {
158  for (unsigned int j = i + offset + 1; j < pts.size(); j += stride) {
159  const Line2D l(pts[i], pts[j]);
160 
161  if (isinf(l.m) || isnan(l.m) || isinf(l.c) || isnan(l.c)) continue;
162 
163  if (R > 0) {
164  float z1, z2;
165  if (!IntersectsCircle(l.m, l.c, z0, x0, 2.5, z1, z2)) continue;
166  if (l.minz < z1 && l.minz < z2 && l.maxz > z1 && l.maxz > z2) continue;
167  }
168 
169  lines.push_back(std::move(l));
170  if (lines.size() == kMaxLines) goto end; // break out of 3 loops
171  }
172  }
173  }
174 
175  end:
176 
177  lines.shrink_to_fit();
178 
179  mf::LogInfo() << "Made " << lines.size() << " lines using stride " << stride
180  << " to fit under cap of " << kMaxLines << std::endl;
181 
182  // Lines are required to be sorted by gradient for a later optimization
183  std::sort(lines.begin(), lines.end());
184  }
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
static QStrList * l
Definition: config.cpp:1044
bool IntersectsCircle(float m, float c, float z0, float x0, float R, float &z1, float &z2)
def move(depos, offset)
Definition: depos.py:107
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
QTextStream & endl(QTextStream &s)
void quad::MapFromLines ( const std::vector< Line2D > &  lines,
HeatMap hm 
)

Definition at line 197 of file QuadVtx_module.cc.

198  {
199  // This maximum is driven by runtime
200  constexpr size_t kMaxPts = 10 * 1000 * 1000;
201 
202  unsigned int j0 = 0;
203  unsigned int jmax = 0;
204 
205  long npts = 0;
206  for (unsigned int i = 0; i + 1 < lines.size(); ++i) {
207  const Line2D a = lines[i];
208 
209  j0 = std::max(j0, i + 1);
210  while (j0 < lines.size() && CloseAngles(a.m, lines[j0].m))
211  ++j0;
212  jmax = std::max(jmax, j0);
213  while (jmax < lines.size() && !CloseAngles(a.m, lines[jmax].m))
214  ++jmax;
215 
216  npts += jmax - j0;
217  }
218 
219  const size_t product = (lines.size() * (lines.size() - 1)) / 2;
220  const int stride = npts / kMaxPts + 1;
221 
222  mf::LogInfo() << "Combining lines to points with stride " << stride << std::endl;
223 
224  mf::LogInfo() << npts << " cf " << product << " ie " << double(npts) / product << std::endl;
225 
226  j0 = 0;
227  jmax = 0;
228 
229  for (unsigned int i = 0; i + 1 < lines.size(); ++i) {
230  const Line2D a = lines[i];
231 
232  j0 = std::max(j0, i + 1);
233  while (j0 < lines.size() && CloseAngles(a.m, lines[j0].m))
234  ++j0;
235  jmax = std::max(jmax, j0);
236  while (jmax < lines.size() && !CloseAngles(a.m, lines[jmax].m))
237  ++jmax;
238 
239  for (unsigned int j = j0; j < jmax; j += stride) {
240  const Line2D& b = lines[j];
241 
242  // x = mA * z + cA = mB * z + cB
243  const float z = (b.c - a.c) / (a.m - b.m);
244  const float x = a.m * z + a.c;
245 
246  // No solutions within a line
247  if ((z < a.minz || z > a.maxz) && (z < b.minz || z > b.maxz)) {
248  const int iz = hm.ZToBin(z);
249  const int ix = hm.XToBin(x);
250  if (iz >= 0 && iz < hm.Nz && ix >= 0 && ix < hm.Nx) { hm.map[iz * hm.Nx + ix] += stride; }
251  }
252  } // end for i
253  }
254  }
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
const double a
static int max(int a, int b)
bool CloseAngles(float ma, float mb)
static bool * b
Definition: config.cpp:1043
list x
Definition: train.py:276
QTextStream & endl(QTextStream &s)