19 #include "Utilities/func/MathUtil.h" 29 int ny,
float y0,
float y1)
33 fBins =
new float[(nx+2)*(ny+2)];
34 memset(
fBins, 0, (nx+2)*(ny+2)*
sizeof(
float));
61 const int bin = iy*(
fNx+2)+ix;
62 if(bin < 0 || bin > (
fNx+2)*(
fNy+2))
return 0;
69 if(ret >
fNx+1)
return fNx+1;
76 if(ret >
fNy+1)
return fNy+1;
81 const int bin = iy*(
fNx+2)+ix;
87 const int bin = iy*(
fNx+2)+ix;
94 for(
int ix = 0; ix <
fNx+2; ++ix)
95 for(
int iy = 0; iy <
fNy+2; ++iy)
123 fSigma(3.0/
sqrt(12.)),
135 double x2,
double y2,
136 double* rho,
double*
theta,
137 double* sigmarho,
double* sigmatheta)
139 *theta =
atan2(x1-x2,y2-y1);
140 *rho = 0.5*(
cos(*theta)*(x1+
x2)+
sin(*theta)*(y1+
y2));
148 while (*theta>=
M_PI) {
156 *sigmatheta = M_SQRT2*
fSigma/
d;
191 nth,-0.5*
M_PI/180.0,179.5*
M_PI/180.0);
200 ret->SetTitle(
";rho [cm];theta [radians]");
217 const Double_t arg = x/
sigma;
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);
237 double tlo=9E9, thi=-9e9;
240 for (i=0; i<h.size(); ++
i) {
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;
262 for (i=0; i<h.size(); ++
i) {
264 for (j=i+1; j<h.size(); ++
j) {
266 const bool hasweight = h[
i].fW != 0.0 && h[
j].fW != 0.0;
267 if(!hasweight)
continue;
272 if (rsqr>
fRsqr)
continue;
278 if (
abs(h[i].fXY - h[j].fXY) < 3.9 && rsqr <
fRsqr/4.0 )
continue;
281 double sigmarho, sigmatheta;
285 &sigmarho, &sigmatheta);
294 if (ixlo<1) ixlo = 1;
295 if (iylo<1) iylo = 1;
303 for (iy=iylo; iy<=iyhi; ++iy) {
305 yprob[iy] =
Gaus(y-theta, sigmatheta);
308 for (ix=ixlo; ix<=ixhi; ++ix) {
310 const double xprob =
Gaus(x-rho, sigmarho);
311 for (iy=iylo; iy<=iyhi; ++iy) {
315 w = 2.0*h[
i].fW*h[
j].fW/(h[
i].fW+h[
j].fW) *
359 Ps =
sqrt(Ps/nb - Pave*Pave);
364 if(th == 0.0)
return;
368 int xpeak = 0, ypeak = 0;
406 if(peak < th) flag = 1;
407 if(mh >=
fLoops) flag = 1;
412 if(mh >=
fLoops) flag = 1;
453 for(
unsigned int i = 0;
i <
fH.
fPeak.size(); ++
i) {
459 if(peak < th) flag = 1;
489 for(
int a = 0;
a < wx;
a++) {
490 for(
int b = 0;
b < wy;
b++) {
501 for (i=1; i<=nx; ++
i) {
502 for(j=1; j<=ny; ++
j) {
514 if (iy>ny) iy = iy%ny;
515 if (iy<1) iy = ny+iy;
536 double peaktotal = 0.0;
537 double BinValue = 0.0;
543 if((xp+nx) >= nXbins) nx = (nXbins-1) - xp;
544 if((xp-nx) < 0) nx = xp;
546 if((yp+ny) >= nYbins) ny = (nYbins-1) - yp;
547 if((yp-ny) < 0) ny = yp;
551 double xd = 0.0, yd = 0.0;
554 for(
int xi = -nx; xi <= nx; ++xi) {
556 for(
int yi = -ny; yi <= ny; ++yi) {
560 if(dist == 0.0) dist = 0.707107;
563 peaktotal += BinValue;
571 theta = theta/peaktotal;
598 double s =
sin(theta);
599 double c =
cos(theta);
602 if (s == 0.0) s = 1.0E-9;
609 double zmin = 1.0e9,
zmax = -1.0e9;
610 unsigned int amax = 999999;
611 unsigned int amin = 999999;
615 unsigned int amax2 = 999999;
616 unsigned int amin2 = 999999;
618 for(
unsigned int a = 0;
a < h.size(); ++
a) {
630 if(h[
a].fZ ==
zmax) { amax2 =
a; }
631 if(h[
a].fZ == zmin) { amin2 =
a; }
633 if(h[
a].fZ < zmin) { zmin = h[
a].fZ; amin =
a; }
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;
648 double dminsq = 0.0, dsq = 0.0;
649 double ScrubDistSq = 900.0;
653 std::vector<int> BuddyList;
654 for(
unsigned int a = 0;
a < h.size(); ++
a) BuddyList.push_back(0);
658 for(
unsigned int i = 0;
i < h.size(); ++
i) {
659 if(BuddyList[
i] > 0)
continue;
661 if(h[
i].fW == 0.0)
continue;
663 for(
unsigned int j = 0;
j < h.size(); ++
j) {
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) {
672 if(dsq < dminsq) dminsq = dsq;
674 if(dminsq >= ScrubDistSq) h[
i].fW = 0.0;
double fXY0
X/Y offset applied to hits.
int FindBinY(float y) const
double fRhoSz
Size of rho bins (cm)
double fZ0
Z offset applied to hits.
fvar< T > fabs(const fvar< T > &x)
Perform a "2 point" Hough transform on a collection of hits.
Float_t y1[n_points_granero]
LiteTH2 * fHspaceW
Hough transform.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Float_t x1[n_points_granero]
double fRsqr
Distance^2 cut on points (cm^2)
std::vector< HoughPeak > fPeak
List of peaks found in Hough space.
Vertical planes which measure X.
int fXwinSz
Smoothing size (bins) in x.
int fYwinSz
Smoothing size (bins) in y.
double fSigma
Assumed location error on points (cm)
float GetBinCenterX(int ix) const
T sqr(T x)
More efficient square function than pow(x,2)
Horizontal planes which measure Y.
double fThetaSz
Size of theta bins (degrees)
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
void SortPeaks()
Sort Hough peaks from best (highest) to worst.
void ReweightHits(double rho, double theta, rb::HitList &h)
Just the essential information required for track fitting.
rb::HoughResult fH
Results of Hough transformation;.
Data for a single peak in Hough space.
double fPco
Number of sigma abouve average peak height to use as cutoff.
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...
fvar< T > exp(const fvar< T > &x)
int fLoops
Max number of Multi-Hough loops.
void Map(const rb::HitList &h)
Collect Geo headers and supply basic geometry functions.
void SetBinContent(int ix, int iy, float val)
double fRmDst
Distance cut for removing points from a Hough line.
static Double_t Gaus(const Double_t x, const Double_t sigma)
float GetBinContent(int ix, int iy) const
LiteTH2(int nx, float x0, float x1, int ny, float y0, float y1)
float GetBinCenterY(int iy) const
double fTNShi
Upper edge of time slice transformed.
LiteTH2(const LiteTH2 &h)
double fTNSlo
Low edge of time slice transformed.
double pythag(double x, double y)
2D Euclidean distance
void BuildMap(double xlo, double xhi, double ylo, double yhi)
assert(nhit_max >=nhit_nbins)
void RefinePeak(double &rho, double &theta, int xp, int yp)
void RhoTheta(double x1, double y1, double x2, double y2, double *rho, double *theta, double *sigmarho, double *sigmatheta)
void AddBinContent(int ix, int iy, float val)
geo::View_t fView
Which detector view?