MultiHough2P.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Perform a "2 point" Hough transform on a collection of hits
3 /// \authors messier@indiana.edu, mibaird@indiana.edu
4 /// \date
5 ////////////////////////////////////////////////////////////////////////
7 
8 #include <algorithm>
9 #include <cmath>
10 #include <cstdlib>
11 #include <cassert>
12 #include <iostream>
13 
14 #include "TH2F.h"
15 #include "TMath.h"
16 
17 #include "GeometryObjects/Geo.h"
18 #include "RecoBase/HitList.h"
19 #include "Utilities/func/MathUtil.h"
20 
21 namespace hough
22 {
23  // TH2F is a bad temporary storage for a fast algorithm. There's a lot of
24  // overhead. Using this simple version is at least 2x faster.
25  class LiteTH2
26  {
27  public:
28  LiteTH2(int nx, float x0, float x1,
29  int ny, float y0, float y1)
30  : fNx(nx), fX0(x0), fX1(x1),
31  fNy(ny), fY0(y0), fY1(y1)
32  {
33  fBins = new float[(nx+2)*(ny+2)];
34  memset(fBins, 0, (nx+2)*(ny+2)*sizeof(float));
35  }
36  LiteTH2(const LiteTH2& h)
37  : fNx(h.fNx), fX0(h.fX0), fX1(h.fX1),
38  fNy(h.fNy), fY0(h.fY0), fY1(h.fY1)
39  {
40  fBins = new float[(fNx+2)*(fNy+2)];
41  memcpy(fBins, h.fBins, (fNx+2)*(fNy+2)*sizeof(float));
42  }
44  {
45  delete[] fBins;
46  }
47  int GetNbinsX() const {return fNx;}
48  int GetNbinsY() const {return fNy;}
49  float GetBinCenterX(int ix) const
50  {
51  assert(ix >= 0 && ix < fNx+2);
52  return fX0+(ix-.5)/fNx*(fX1-fX0);
53  }
54  float GetBinCenterY(int iy) const
55  {
56  assert(iy >= 0 && iy < fNy+2);
57  return fY0+(iy-.5)/fNy*(fY1-fY0);
58  }
59  float GetBinContent(int ix, int iy) const
60  {
61  const int bin = iy*(fNx+2)+ix;
62  if(bin < 0 || bin > (fNx+2)*(fNy+2)) return 0;
63  return fBins[bin];
64  }
65  int FindBinX(float x) const
66  {
67  int ret = fNx*(x-fX0)/(fX1-fX0)+1;
68  if(ret < 0) return 0;
69  if(ret > fNx+1) return fNx+1;
70  return ret;
71  }
72  int FindBinY(float y) const
73  {
74  int ret = fNy*(y-fY0)/(fY1-fY0)+1;
75  if(ret < 0) return 0;
76  if(ret > fNy+1) return fNy+1;
77  return ret;
78  }
79  void SetBinContent(int ix, int iy, float val)
80  {
81  const int bin = iy*(fNx+2)+ix;
82  assert(bin >= 0 && bin < (fNx+2)*(fNy+2));
83  fBins[iy*(fNx+2)+ix] = val;
84  }
85  void AddBinContent(int ix, int iy, float val)
86  {
87  const int bin = iy*(fNx+2)+ix;
88  assert(bin >= 0 && bin < (fNx+2)*(fNy+2));
89  fBins[iy*(fNx+2)+ix] += val;
90  }
91  TH2F* ToTH2F() const
92  {
93  TH2F* ret = new TH2F("", "", fNx, fX0, fX1, fNy, fY0, fY1);
94  for(int ix = 0; ix < fNx+2; ++ix)
95  for(int iy = 0; iy < fNy+2; ++iy)
96  ret->SetBinContent(ix, iy, GetBinContent(ix, iy));
97  return ret;
98  }
99  protected:
100  float* fBins;
101  int fNx;
102  float fX0, fX1;
103  int fNy;
104  float fY0, fY1;
105  };
106 
107  //......................................................................
108 
110  double rhosz,
111  double thetasz,
112  int rhosm,
113  int thetasm,
114  double rsqr,
115  double pco,
116  int loops,
117  double rmdst) :
118  fView(v),
119  fRhoSz(rhosz),
120  fThetaSz(thetasz),
121  fXwinSz(rhosm),
122  fYwinSz(thetasm),
123  fSigma(3.0/sqrt(12.)),
124  fRsqr(rsqr),
125  fPco(pco),
126  fLoops(loops),
127  fRmDst(rmdst),
128  fHspaceW(0),
129  fH(v, 0, 0, 0, 0)
130  { }
131 
132  //......................................................................
133 
134  void MultiHough2P::RhoTheta(double x1, double y1,
135  double x2, double y2,
136  double* rho, double* theta,
137  double* sigmarho, double* sigmatheta)
138  {
139  *theta = atan2(x1-x2,y2-y1);
140  *rho = 0.5*(cos(*theta)*(x1+x2)+sin(*theta)*(y1+y2));
141  //
142  // Keep rho and theta within range [0,pi) (-R,+R)
143  //
144  while (*theta<0) {
145  *theta += M_PI;
146  *rho = -(*rho);
147  }
148  while (*theta>=M_PI) {
149  *theta -= M_PI;
150  *rho = -(*rho);
151  }
152 
153  // Compute the error estimates
154  double d = util::pythag((x2-x1),(y2-y1));
155  *sigmarho = fSigma;
156  *sigmatheta = M_SQRT2*fSigma/d;
157 
158  // these were failed attempts at making sigmatheta level off at a fixed value
159  //*sigmatheta = M_SQRT2*fSigma/d + 0.1;
160  /*
161  if(d < 20.0) *sigmatheta = M_SQRT2*fSigma/d;
162  else *sigmatheta = M_SQRT2*fSigma/20.0;
163  */
164  }
165 
166  //......................................................................
167 
169  {
170  if (fHspaceW) { delete fHspaceW; fHspaceW = 0;}
171  }
172 
173  //......................................................................
174 
175  void MultiHough2P::BuildMap(double xlo, double xhi,
176  double ylo, double yhi)
177  {
178  double r; // Size of hough map in rho
179  int nr; // Number of bins to use in rho
180  int nth; // Number of bins to use in theta
181 
182  // 0.55 gives a little buffer beyond the minimum required 0.5
183  r = 0.55*util::pythag((xhi-xlo),(yhi-ylo));
184  nr = (int)rint(r/fRhoSz+0.5);
185 
186  nth = (int)rint(180.0/fThetaSz+0.5);
187 
188  // Allocate the Hough space histogram
189  if (fHspaceW!=0) { delete fHspaceW; fHspaceW=0; }
190  fHspaceW = new LiteTH2(nr, -r, r,
191  nth,-0.5*M_PI/180.0,179.5*M_PI/180.0);
192  }
193 
194  //......................................................................
195 
197  {
198  if(!fHspaceW) return 0;
199  TH2F* ret = fHspaceW->ToTH2F();
200  ret->SetTitle(";rho [cm];theta [radians]");
201 
202  const char* tiw = 0;
203  if(fView == geo::kX) tiw = "fHough2P_Hx";
204  else if (fView == geo::kY) tiw = "fHough2P_Hy";
205  else abort();
206  ret->SetTitle(tiw);
207 
208  return ret;
209  }
210 
211  //......................................................................
212 
213  // This is special-case version of TMath::Gaus that saves time by assuming
214  // the mean is zero and not including a switch for whether to normalize.
215  static Double_t Gaus(const Double_t x, const Double_t sigma)
216  {
217  const Double_t arg = x/sigma;
218  // for |arg| > 39 result is zero in double precision
219  if (arg < -39.0 || arg > 39.0) return 0.0;
220  const Double_t res = exp(-0.5*arg*arg);
221  return res/(2.50662827463100024*sigma); //sqrt(2*Pi)=2.50662827463100024
222  }
223 
224 
226  {
227  unsigned int i;
228  unsigned int j;
229 
230  //
231  // Find a box that completely contains all the hits. A note on
232  // coordinates: Since the literature on Hough transforms uses an
233  // "x-y" coodinate system I map NOvA z to Hough x and NOvA X/Y to
234  // Hough y so that the formalizm matches what one would find in the
235  // literature.
236  //
237  double tlo=9E9, thi=-9e9;
238  double xlo=9e9, xhi=-9e9;
239  double ylo=9e9, yhi=-9e9;
240  for (i=0; i<h.size(); ++i) {
241  if (h[i].fView == fView) {
242  if (h[i].fZ<xlo) xlo = h[i].fZ;
243  if (h[i].fZ>xhi) xhi = h[i].fZ;
244  if (h[i].fXY<ylo) ylo = h[i].fXY;
245  if (h[i].fXY>yhi) yhi = h[i].fXY;
246  if (h[i].fT<tlo) tlo = h[i].fT;
247  if (h[i].fT>tlo) thi = h[i].fT;
248  }
249  }
250  fH.fTNSlo = tlo;
251  fH.fTNShi = thi;
252  fH.fZ0 = 0.5*(xlo+xhi);
253  fH.fXY0 = 0.5*(ylo+yhi);
254 
255 
256 
257  //
258  // Fill the hough map
259  //
260  this->BuildMap(xlo,xhi,ylo,yhi);
261 
262  for (i=0; i<h.size(); ++i) {
263  if (h[i].fView != fView) continue;
264  for (j=i+1; j<h.size(); ++j) {
265  if (h[j].fView != fView) continue;
266  const bool hasweight = h[i].fW != 0.0 && h[j].fW != 0.0;
267  if(!hasweight) continue;
268 
269  // Check distance between hit pair
270  double rsqr = util::sqr(h[i].fZ-h[j].fZ) + util::sqr(h[i].fXY-h[j].fXY);
271 
272  if (rsqr>fRsqr) continue;
273 
274  // don't let points in the same column vote unless they are close together
275  // if ( abs(h[i].fZ - h[j].fZ) < 5.9 && rsqr > fRsqr/4.0 ) continue;
276 
277  // don't let points in the same row vote unless they are far apart
278  if ( abs(h[i].fXY - h[j].fXY) < 3.9 && rsqr < fRsqr/4.0 ) continue;
279 
280  double rho, theta;
281  double sigmarho, sigmatheta;
282  this->RhoTheta(h[i].fZ-fH.fZ0, h[i].fXY-fH.fXY0,
283  h[j].fZ-fH.fZ0, h[j].fXY-fH.fXY0,
284  &rho, &theta,
285  &sigmarho, &sigmatheta);
286  //
287  // Fill region weighted by errors
288  //
289  int ixlo = fHspaceW->FindBinX(rho-3.0*sigmarho);
290  int ixhi = fHspaceW->FindBinX(rho+3.0*sigmarho);
291  int iylo = fHspaceW->FindBinY(theta-3.0*sigmatheta);
292  int iyhi = fHspaceW->FindBinY(theta+3.0*sigmatheta);
293 
294  if (ixlo<1) ixlo = 1;
295  if (iylo<1) iylo = 1;
296  if (ixhi>fHspaceW->GetNbinsX()) ixhi = fHspaceW->GetNbinsX();
297  if (iyhi>fHspaceW->GetNbinsY()) iyhi = fHspaceW->GetNbinsY();
298 
299  int ix, iy;
300 
301  // Precalcuate gaussian weights since they will be the same every loop
302  double yprob[iyhi]; // won't use first elements: use memory to save time
303  for (iy=iylo; iy<=iyhi; ++iy) {
304  const double y = fHspaceW->GetBinCenterY(iy);
305  yprob[iy] = Gaus(y-theta, sigmatheta);
306  }
307 
308  for (ix=ixlo; ix<=ixhi; ++ix) {
309  double x = fHspaceW->GetBinCenterX(ix);
310  const double xprob = Gaus(x-rho, sigmarho);
311  for (iy=iylo; iy<=iyhi; ++iy) {
312  double w = 0.0;
313 
314  // if(h[i].fW > 0.1 && h[j].fW > 0.1) w = 2.0*h[i].fW*h[j].fW/(h[i].fW+h[j].fW);
315  w = 2.0*h[i].fW*h[j].fW/(h[i].fW+h[j].fW) *
316  xprob * yprob[iy];
317  fHspaceW->AddBinContent(ix,iy,w);
318 
319  }
320  }
321  }
322  }
323  if (fXwinSz > 0 || fYwinSz > 0) this->SmoothMap();
324 
325  }
326 
327  //......................................................................
328 
330  {
331 
332  int flag = 0; // flag to jump out of multi hough loop
333  int mh = 1; // count of number of multihough loops
334 
335  // make the initial Hough Map
336  this->Map(h);
337 
338  // calculate the threshold (needs to be done once at the beginning so that it is
339  // the same for all subsequent maps)
340  double Pave = 0.0; // average peak height
341  double Ps = 0.0; // standard deviation of average peak height
342  double th = 0.0; // threshold to use for peak cutoff
343  int nb = 0; // number of bins counted for avergae peak height
344 
345  // calculate average peak height to determine threshold
346  // empty bins are currently counted
347  for (int a = 1; a <= fHspaceW->GetNbinsX(); ++a) {
348  for (int b = 1; b <= fHspaceW->GetNbinsY(); ++b) {
349  double hab = fHspaceW->GetBinContent(a, b);
350  if (hab != 0) {
351  Pave += hab;
352  Ps += hab*hab;
353  } // end if
354  nb++;
355  } // end for b
356  } // end for a
357 
358  Pave = Pave/nb;
359  Ps = sqrt(Ps/nb - Pave*Pave);
360  th = Pave + fPco*Ps;
361 
362  // If the threshold is zero, then the initial Hough map is empty so there
363  // is no reason to proceed.
364  if(th == 0.0) return;
365 
366  // find the tallest peak
367  int ix, iy;
368  int xpeak = 0, ypeak = 0;
369  double peak = 0.0;
370  double r = 0.0;
371  double theta = 0.0;
372  double p;
373  for(ix = 1; ix <= fHspaceW->GetNbinsX(); ++ix) {
374  for(iy = 1; iy <= fHspaceW->GetNbinsY(); ++iy) {
375  p = fHspaceW->GetBinContent(ix,iy);
376  if (p>peak) {
377  peak = p;
378  xpeak = ix;
379  ypeak = iy;
380  } // end if p>peak
381  } // loop on y bins
382  } // loop on x bins
383 
384  // Average bins around the peak to get a better rho and theta value
385  this->RefinePeak(r, theta, xpeak, ypeak);
386 
387  // AS OF RIGHT NOW, SIGMA-RHO AND SIGMA-THETA ARE NOT CALCULATED!!!
388  double sr = 0.0;
389  double st = 0.0;
390 
391  // put tallest peak into Hough Results
392  if(peak > 0.0) fH.fPeak.push_back(rb::HoughPeak(r,theta,peak,th,sr,st));
393 
394  // remove hits on the line by reweighting them
395  this->ReweightHits(r, theta, h);
396 
397 
398 
399 
400  //
401  // start the MultiHough loop
402  //
403 
404  // if either the number of desired loops is only one or the tallest peak was already below threshold,
405  // then skip the multi-hough loop because we are already done
406  if(peak < th) flag = 1;
407  if(mh >= fLoops) flag = 1;
408 
409  while(flag == 0) {
410 
411  mh++;
412  if(mh >= fLoops) flag = 1;
413 
414  // repeat everything done above EXCEPT: calculating the threshold
415 
416  // build the new map
417  this->Map(h);
418 
419 
420 
421  // find the tallest peak
422  xpeak = 0;
423  ypeak = 0;
424  peak = 0.0;
425  r = 0.0;
426  theta = 0.0;
427  p = 0.0;
428 
429  for (ix=1; ix<=fHspaceW->GetNbinsX(); ++ix) {
430  for (iy=1; iy<=fHspaceW->GetNbinsY(); ++iy) {
431  p = fHspaceW->GetBinContent(ix,iy);
432  if (p>peak) {
433  peak = p;
434  xpeak = ix;
435  ypeak = iy;
436  } // end if p>peak
437  } // loop on y bins
438  } // loop on x bins
439 
440  this->RefinePeak(r, theta, xpeak, ypeak);
441 
442  sr = 0.0;
443  st = 0.0;
444 
445 
446 
447  // Run parallel line test: (If two lines have almost the same rho and theta values,
448  // they are assumed to be part of the same line and that line is NOT put into the
449  // Hough results.)
450  //
451  // NOTE: This is a little sensitive to the numbers 15.0 and 0.02.
452  int parflag = 0;
453  for(unsigned int i = 0; i < fH.fPeak.size(); ++i) {
454  if(fabs(r - fH.fPeak[i].fRho) < 15.0 && fabs(theta - fH.fPeak[i].fTheta) < 0.02) parflag = 1;
455  }
456 
457  // put only the tallest peak into Houghresults
458  if(parflag != 1 && peak >= th) fH.fPeak.push_back(rb::HoughPeak(r,theta,peak,th,sr,st));
459  if(peak < th) flag = 1;
460 
461  this->ReweightHits(r, theta, h);
462 
463  } // end MultiHough loop
464 
465  fH.SortPeaks();
466 
467  }
468 
469  //......................................................................
470 
472  {
473 
474  // NOTE: Since most of the Hough space map is empty, we only apply smoothing to the non-zero
475  // bins in order to save lots of calculation time.
476 
477  int i, j, m, n;
478  int nx, ny;
479  int ix, iy;
480  double sum;
481  int wx, wy;
482 
483  wx = 2*fXwinSz + 1;
484  wy = 2*fYwinSz + 1;
485 
486  double w[wx][wy];
487 
488  // fill the smoothing weight map
489  for(int a = 0; a < wx; a++) {
490  for(int b = 0; b < wy; b++) {
491 
492  w[a][b] = Gaus(a - fXwinSz, fXwinSz/3.0)*Gaus(b - fYwinSz, fYwinSz/3.0);
493 
494  } // end for b
495  } // end for a
496 
497  LiteTH2* htmp = new LiteTH2(*fHspaceW);
498 
499  nx = htmp->GetNbinsX();
500  ny = htmp->GetNbinsY();
501  for (i=1; i<=nx; ++i) {
502  for(j=1; j<=ny; ++j) {
503  sum = 0.0;
504  // skip smoothing if bin is empty
505  if(htmp->GetBinContent(i,j) == 0.0) continue;
506  for (m=-fXwinSz; m<=fXwinSz; m++) {
507  ix = i+m;
508  // Don't wrap around in rho
509  if (ix>nx) continue;
510  if (ix<1) continue;
511  for (n=-fYwinSz; n<=fYwinSz; ++n) {
512  iy = j+n;
513  // Wrap in theta
514  if (iy>ny) iy = iy%ny;
515  if (iy<1) iy = ny+iy;
516 
517  sum += w[m + fXwinSz][n + fYwinSz]*htmp->GetBinContent(ix,iy);
518  }
519  }
520  fHspaceW->SetBinContent(i,j,sum);
521  }
522  }
523  delete htmp;
524  htmp = 0;
525 
526  }
527 
528 
529 
530  //......................................................................
531 
532  void MultiHough2P::RefinePeak(double &rho, double &theta, int xp, int yp)
533  {
534  // make these fcl parameters???
535  int nx = 3, ny = 3; // +/- number of bins around peak to use for averaging
536  double peaktotal = 0.0;
537  double BinValue = 0.0;
538  int nXbins = fHspaceW->GetNbinsX();
539  int nYbins = fHspaceW->GetNbinsY();
540 
541 
542  // determine if bins to be averaged over will be off the map and constrain nx and ny accordingly
543  if((xp+nx) >= nXbins) nx = (nXbins-1) - xp;
544  if((xp-nx) < 0) nx = xp;
545 
546  if((yp+ny) >= nYbins) ny = (nYbins-1) - yp;
547  if((yp-ny) < 0) ny = yp;
548 
549 
550 
551  double xd = 0.0, yd = 0.0;
552  double dist = 0.0;
553 
554  for(int xi = -nx; xi <= nx; ++xi) {
555  xd = fabs(xi);
556  for(int yi = -ny; yi <= ny; ++yi) {
557 
558  yd = fabs(yi);
559  dist = util::pythag(xd,yd);
560  if(dist == 0.0) dist = 0.707107;
561 
562  BinValue = (fHspaceW->GetBinContent(xp+xi,yp+yi))/dist;
563  peaktotal += BinValue;
564  rho += BinValue*fHspaceW->GetBinCenterX(xp+xi);
565  theta += BinValue*fHspaceW->GetBinCenterY(yp+yi);
566 
567  }
568  }
569 
570  rho = rho/peaktotal;
571  theta = theta/peaktotal;
572 
573  }
574 
575 
576 
577  //......................................................................
578 
580  {
581  //
582  // This function does several things:
583  // 1. It sets the weights of all hits within 6 cm of the line to zero.
584  // 2. Then it resets the weights of the end points of the line to one.
585  // 3. Then it applies more scrubbing to remove stray hits that have no neighbors.
586  //
587 
588  // These variable were used in an old scheme in which I used a Fermi distribution function
589  // to determine the new hit weights based on their distance to the line. I left it in for
590  // potential future use.
591  // double dc = 4.0; // distance cutoff
592  // double ds = 1.0; // distance spread
593 
594 
595  // loop over all points, determine which are on that line and set their weights to zero
596  double m = 0.0;
597  double b = 0.0;
598  double s = sin(theta);
599  double c = cos(theta);
600 
601  // calculate slope/intercept in detector coordinates
602  if (s == 0.0) s = 1.0E-9;
603  m = -c/s;
604  b = r/s - m*fH.fZ0 + fH.fXY0;
605 
606  double d = 0.0; // perpendicular distance of point to line
607 
608  // numbers needed for determining which points are end points
609  double zmin = 1.0e9, zmax = -1.0e9;
610  unsigned int amax = 999999;
611  unsigned int amin = 999999;
612  // amax2 & amin2 are used for the case that there are 2 end points (when there are two points
613  // in the same plane within fRmDist (cm) of the line.)
614  // NOTE: There could be three points that fit this, but I only allow for keeping two.
615  unsigned int amax2 = 999999;
616  unsigned int amin2 = 999999;
617 
618  for(unsigned int a = 0; a < h.size(); ++a) {
619 
620  // if a point is on the line, set its weight to zero
621  if(h[a].fView == fView) {
622 
623  d = geo::DsqrToLine(h[a].fZ, h[a].fXY, 0.0, b, 100.0, 100.0*m + b);
624  d = sqrt(d);
625  // h[a].fW *= 1.0 - 1.0/(exp((d - dc)/ds) + 1.0); // The Fermi weight function
626  if(d < fRmDst) {
627  h[a].fW = 0.0;
628 
629  // determine end points
630  if(h[a].fZ == zmax) { amax2 = a; }
631  if(h[a].fZ == zmin) { amin2 = a; }
632  if(h[a].fZ > zmax) { zmax = h[a].fZ; amax = a; }
633  if(h[a].fZ < zmin) { zmin = h[a].fZ; amin = a; }
634  }
635 
636  } // end if fView
637  } // end for a
638 
639  // Reset the weights of the line end points to one.
640  if(amax != 999999) h[amax].fW = 1.0;
641  if(amin != 999999) h[amin].fW = 1.0;
642  if(amax2 != 999999) h[amax2].fW = 1.0;
643  if(amin2 != 999999) h[amin2].fW = 1.0;
644 
645 
646 
647  // Rescrub the hit list to remove points with no neighbors within sqrt(ScrubDistSq) (cm).
648  double dminsq = 0.0, dsq = 0.0;
649  double ScrubDistSq = 900.0;
650 
651  // We will keep track of which points are within sqrt(ScrubDistSq) (cm) of
652  // another point (i.e. - which points have a buddy) to make this more efficient.
653  std::vector<int> BuddyList;
654  for(unsigned int a = 0; a < h.size(); ++a) BuddyList.push_back(0);
655 
656  // For each point, loop over all others to determine if it has a "buddy" within sqrt(ScrubDistSq).
657  // If it has no "buddies," scrub it! (Set the weight to zero.)
658  for(unsigned int i = 0; i < h.size(); ++i) {
659  if(BuddyList[i] > 0) continue;
660  if(h[i].fView != fView) continue;
661  if(h[i].fW == 0.0) continue;
662  dminsq = 1.0e9;
663  for(unsigned int j = 0; j < h.size(); ++j) {
664  if(i == j) continue;
665  if(h[j].fView != fView) continue;
666  if(h[j].fW == 0.0) continue;
667  dsq = (h[i].fZ - h[j].fZ)*(h[i].fZ - h[j].fZ) + (h[i].fXY - h[j].fXY)*(h[i].fXY - h[j].fXY);
668  if(dsq < ScrubDistSq) {
669  BuddyList[i] += 1;
670  BuddyList[j] += 1;
671  }
672  if(dsq < dminsq) dminsq = dsq;
673  }
674  if(dminsq >= ScrubDistSq) h[i].fW = 0.0;
675  }
676 
677  }
678 
679 
680 
681 } // end namespace hough
682 ////////////////////////////////////////////////////////////////////////
double fXY0
X/Y offset applied to hits.
Definition: HoughResult.h:59
int FindBinY(float y) const
double fRhoSz
Size of rho bins (cm)
Definition: MultiHough2P.h:79
double fZ0
Z offset applied to hits.
Definition: HoughResult.h:60
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Perform a "2 point" Hough transform on a collection of hits.
Definition: Hough2P.cxx:20
Float_t y1[n_points_granero]
Definition: compare.C:5
LiteTH2 * fHspaceW
Hough transform.
Definition: MultiHough2P.h:90
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Float_t x1[n_points_granero]
Definition: compare.C:5
double fRsqr
Distance^2 cut on points (cm^2)
Definition: MultiHough2P.h:84
std::vector< HoughPeak > fPeak
List of peaks found in Hough space.
Definition: HoughResult.h:61
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
int fXwinSz
Smoothing size (bins) in x.
Definition: MultiHough2P.h:81
int fYwinSz
Smoothing size (bins) in y.
Definition: MultiHough2P.h:82
void abs(TH1 *hist)
double fSigma
Assumed location error on points (cm)
Definition: MultiHough2P.h:83
float GetBinCenterX(int ix) const
TH2F * MapToTH2F() const
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
#define M_PI
Definition: SbMath.h:34
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
double fThetaSz
Size of theta bins (degrees)
Definition: MultiHough2P.h:80
double dist
Definition: runWimpSim.h:113
MultiHough2P(geo::View_t v, double rhosz, double thetasz, int rhosm, int thetasm, double rsqr, double pco, int loops, double rmdist)
Construct the 2-point Hough transform.
void MultiMap(rb::HitList h)
int FindBinX(float x) const
const XML_Char * s
Definition: expat.h:262
void SortPeaks()
Sort Hough peaks from best (highest) to worst.
Definition: HoughResult.cxx:83
void ReweightHits(double rho, double theta, rb::HitList &h)
Just the essential information required for track fitting.
Definition: HitList.h:40
rb::HoughResult fH
Results of Hough transformation;.
Definition: MultiHough2P.h:91
Data for a single peak in Hough space.
Definition: HoughResult.h:18
double fPco
Number of sigma abouve average peak height to use as cutoff.
Definition: MultiHough2P.h:85
const double a
double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
In two dimensions, return the perpendicular distance from a point (x0,y0) to a line defined by end po...
Definition: Geo.cxx:338
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
int fLoops
Max number of Multi-Hough loops.
Definition: MultiHough2P.h:86
void Map(const rb::HitList &h)
Collect Geo headers and supply basic geometry functions.
Float_t d
Definition: plot.C:236
caf::StandardRecord * sr
const double j
Definition: BetheBloch.cxx:29
void SetBinContent(int ix, int iy, float val)
double fRmDst
Distance cut for removing points from a Hough line.
Definition: MultiHough2P.h:87
static Double_t Gaus(const Double_t x, const Double_t sigma)
float GetBinContent(int ix, int iy) const
float bin[41]
Definition: plottest35.C:14
LiteTH2(int nx, float x0, float x1, int ny, float y0, float y1)
double sigma(TH1F *hist, double percentile)
float GetBinCenterY(int iy) const
double fTNShi
Upper edge of time slice transformed.
Definition: HoughResult.h:58
LiteTH2(const LiteTH2 &h)
int GetNbinsX() const
int GetNbinsY() const
TH2F * ToTH2F() const
T sin(T number)
Definition: d0nt_math.hpp:132
double fTNSlo
Low edge of time slice transformed.
Definition: HoughResult.h:57
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
void BuildMap(double xlo, double xhi, double ylo, double yhi)
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
void RefinePeak(double &rho, double &theta, int xp, int yp)
static const double nb
Definition: Units.h:89
void RhoTheta(double x1, double y1, double x2, double y2, double *rho, double *theta, double *sigmarho, double *sigmatheta)
Double_t sum
Definition: plot.C:31
void AddBinContent(int ix, int iy, float val)
geo::View_t fView
Which detector view?
Definition: MultiHough2P.h:78
Float_t w
Definition: plot.C:20
T atan2(T number)
Definition: d0nt_math.hpp:72