HoughVertexAlg.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Find event vertex using a Hough transform
3 /// \author Christopher Backhouse - bckhouse@caltech.edu
4 ////////////////////////////////////////////////////////////////////////
6 #include "HoughVertex/LiteTH2.h"
7 #include "Utilities/func/MathUtil.h"
8 
9 #include "TH2.h"
10 
11 #include <cassert>
12 
13 #include "fhiclcpp/ParameterSet.h"
14 
15 namespace hv
16 {
17 
18  // Constructor
20  {
21  fPadV = pset.get<unsigned int>("PadV");
22  fPadZmin = pset.get<unsigned int>("PadZmin");
23  fPadZmax = pset.get<unsigned int>("PadZmax");
24  fNbinsVZ = pset.get<int>("NbinsVZ");
25  fZoomPadSize = pset.get<unsigned int>("ZoomPadSize");
26  fSigmaZ = pset.get<float>("SigmaZ");
27  }
28 
29  // Destructor
31 
32  //......................................................................
33  TVector3 HoughVertexAlg::FindVertex(const std::vector<Vector2D>& xhits,
34  const std::vector<Vector2D>& yhits,
35  TH2F** xmap, TH2F** ymap) const
36  {
37  // Find the event "bounding box"
38  float minx = FLT_MAX, maxx = -FLT_MAX;
39  float miny = FLT_MAX, maxy = -FLT_MAX;
40  float minz = FLT_MAX, maxz = -FLT_MAX;
41  for(const Vector2D& hit: xhits){
42  minx = std::min(hit.v, minx);
43  maxx = std::max(hit.v, maxx);
44  minz = std::min(hit.z, minz);
45  maxz = std::max(hit.z, maxz);
46  }
47  for(const Vector2D& hit: yhits){
48  miny = std::min(hit.v, miny);
49  maxy = std::max(hit.v, maxy);
50  minz = std::min(hit.z, minz);
51  maxz = std::max(hit.z, maxz);
52  }
53 
54  // All these constants were determined by plotting mc.nu.vtx vs slc.boxmin
55  // and slc.boxmax from CAFs (for true vertices inside the detector).
56 
57  // 50cm padding in x and y
58  minx -= fPadV;
59  maxx += fPadV;
60  miny -= fPadV;
61  maxy += fPadV;
62 
63  maxz = std::min(minz+fPadZmax, maxz);
64  minz -= fPadZmin;
65 
66  // cm-sized bins
67  const int Nx = int(maxx-minx)+1;
68  const int Ny = int(maxy-miny)+1;
69  const int Nz = int(maxz-minz)+1;
70 
71  LiteTH2 houghx(Nz, minz, maxz, Nx, minx, maxx);
72  LiteTH2 houghy(Nz, minz, maxz, Ny, miny, maxy);
73 
74  FillMaps(houghx, houghy, xhits, yhits);
75 
76  if(xmap) *xmap = houghx.ToTH2F();
77  if(ymap) *ymap = houghy.ToTH2F();
78 
79  TVector3 vtx = FindMaximum(houghx, houghy);
80 
81  // Zoom into vertex with mm precision
82  LiteTH2 houghxzoom(fNbinsVZ, vtx.Z()-fZoomPadSize, vtx.Z()+fZoomPadSize, fNbinsVZ, vtx.X()-fZoomPadSize, vtx.X()+fZoomPadSize);
83  LiteTH2 houghyzoom(fNbinsVZ, vtx.Z()-fZoomPadSize, vtx.Z()+fZoomPadSize, fNbinsVZ, vtx.Y()-fZoomPadSize, vtx.Y()+fZoomPadSize);
84 
85  FillMaps(houghxzoom, houghyzoom, xhits, yhits);
86 
87  return FindMaximum(houghxzoom, houghyzoom);
88  }
89 
90  //......................................................................
91  void HoughVertexAlg::FillMaps(LiteTH2& houghx, LiteTH2& houghy,
92  const std::vector<Vector2D>& xhits,
93  const std::vector<Vector2D>& yhits) const
94  {
95 
96  for(int view: {0, 1}){
97  const std::vector<Vector2D>& hits = (view == 0) ? xhits : yhits;
98  LiteTH2& hough = (view == 0) ? houghx : houghy;
99 
100  const float dy = hough.GetBinCenterY(1)-hough.GetBinCenterY(0);
101 
102  for(unsigned int i = 0; i < hits.size()-1; ++i){
103  const Vector2D hi = hits[i];
104 
105  for(unsigned int j = i+1; j < hits.size(); ++j){
106  const Vector2D hj = hits[j];
107 
108  // Due to the slight tilt of planes this basically never happens
109  if(hi.z == hj.z) continue;
110 
111  // Solution for y=mx+c
112  const float m = (hj.v-hi.v)/(hj.z-hi.z);
113  const float c = hi.v - m*hi.z;
114 
115  for(int ix = 0; ix < hough.GetNbinsX()+2; ++ix){
116 
117  const float x = hough.GetBinCenterX(ix);
118 
119  double wz = 1.;
120 
121  const float minz = std::min(hi.z, hj.z);
122 
123  const float sigmaz = fSigmaZ; // cm
124 
125  // Weights downstream of the center of the
126  // most-upstream of the two cells decrease with
127  // positive half of a gaussian to smooth out the z-resolution
128  // The gaussian sigma = sqrt(sigmaz) for clarity
129  if(x > minz) wz = exp(-util::sqr((x-minz)/(2*sigmaz)));
130 
131  const float y0 = m*x+c;
132 
133  // Error propagation
134  const float sigma_y_sq =
135  util::sqr( hj.dv * (x-hi.z)/(hj.z-hi.z) ) +
136  util::sqr( hi.dv * (hj.z-x)/(hj.z-hi.z) ) +
137  util::sqr( hj.dz * (hj.v-hi.v)*(hi.z-x)/util::sqr(hj.z-hi.z) ) +
138  util::sqr( hi.dz * (hj.v-hi.v)*(hj.z-x)/util::sqr(hj.z-hi.z) );
139 
140  const float sigma_y = sqrt(sigma_y_sq);
141 
142  const unsigned int Ny = hough.GetNbinsY();
143  const double iy0 = hough.FindBinY(y0-3*sigma_y);
144  double y = hough.GetBinCenterY(iy0)-dy;
145  const double denom = 1/(9*sigma_y_sq);
146  const double inv_sigma = 1/sigma_y;
147  for(unsigned int iy = iy0; iy < Ny+2; ++iy){
148  y += dy;
149 
150  if(y < y0-3*sigma_y) continue;
151  if(y > y0+3*sigma_y) break;
152 
153  // exp() call in computation of gaussian shows up in profiles
154  // const float w = exp(-util::sqr(y-y0)/(2*sigma_y_sq)) / sigma_y;
155 
156  // Triweight looks the most like normal, and has the nice
157  // property that it goes to zero at 3 when scaled to the same SD
158  // as the normal.
159  const double w = util::cube(1-util::sqr(y-y0)*denom) * inv_sigma;
160 
161  hough.AddBinContent(ix, iy, w*wz);
162  } // end for iy
163  } // end for ix
164  } // end for j
165  } // end for i
166  } // end for view
167  } // end for FillMaps
168 
169  //......................................................................
170  TVector3 HoughVertexAlg::FindMaximum(const LiteTH2& houghx,
171  const LiteTH2& houghy) const
172  {
173  float bestw = -1;
174  float bestz = 0, bestx = 0, besty = 0;
175 
176  // OK, both maps have same binning on x (z) axis.
177  for(int iz = 0; iz < houghx.GetNbinsX()+2; ++iz){
178  const float z = houghx.GetBinCenterX(iz);
179 
180  float local_bestx = 0, local_besty = 0;
181  float local_bestwx = -1;
182  float local_bestwy = -1;
183 
184  for(int ix = 0; ix < houghx.GetNbinsY()+2; ++ix){
185  const float x = houghx.GetBinCenterY(ix);
186 
187  if(houghx.GetBinContent(iz, ix) > local_bestwx){
188  local_bestwx = houghx.GetBinContent(iz, ix);
189  local_bestx = x;
190  }
191  } // end for ix
192 
193  for(int iy = 0; iy < houghy.GetNbinsY()+2; ++iy){
194  const float y = houghy.GetBinCenterY(iy);
195  if(houghy.GetBinContent(iz, iy) > local_bestwy){
196  local_bestwy = houghy.GetBinContent(iz, iy);
197  local_besty = y;
198  }
199  } // end for iy
200 
201  if(local_bestwx + local_bestwy > bestw){
202  bestw = local_bestwx+local_bestwy;
203  bestz = z;
204  bestx = local_bestx;
205  besty = local_besty;
206  }
207  } // end for iz
208 
209  return TVector3(bestx, besty, bestz);
210  }
211 } // namespace
Accumulation space for HoughVertex.
Definition: LiteTH2.h:20
T max(const caf::Proxy< T > &a, T b)
void AddBinContent(int ix, int iy, float val)
Definition: LiteTH2.h:90
void FillMaps(LiteTH2 &houghx, LiteTH2 &houghy, const std::vector< Vector2D > &xhits, const std::vector< Vector2D > &yhits) const
Fill the Hough space.
unsigned int fPadZmax
Padding extent, minimum of zmin + PadZmax and zmax.
Perform a "2 point" Hough transform on a collection of hits.
Definition: Hough2P.cxx:20
double maxy
T sqrt(T number)
Definition: d0nt_math.hpp:156
float GetBinCenterX(int ix) const
Definition: LiteTH2.h:48
TH2F * ToTH2F() const
Definition: LiteTH2.h:97
T cube(T x)
More efficient cube function than pow(x,3)
Definition: MathUtil.h:26
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
int GetNbinsY() const
Definition: LiteTH2.h:46
float fSigmaZ
Sigma (Width) [cm] of Gaussian smear in Z is sqrt() of this value.
TVector3 FindMaximum(const LiteTH2 &houghx, const LiteTH2 &houghy) const
Search the Hough space for the highest peak.
TSpline3 hi("hi", xhi, yhi, 18,"0")
void hits()
Definition: readHits.C:15
double dy[NP][NC]
float GetBinCenterY(int iy) const
Definition: LiteTH2.h:54
int fNbinsVZ
Number of bins in V/Z for cm-precision scan.
T get(std::string const &key) const
Definition: ParameterSet.h:231
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
const double j
Definition: BetheBloch.cxx:29
z
Definition: test.py:28
HoughVertexAlg(fhicl::ParameterSet const &pset)
Find event vertex using a Hough transform.
Definition: event.h:1
unsigned int fPadZmin
Padding extent upstream in Z.
int GetNbinsX() const
Definition: LiteTH2.h:45
T min(const caf::Proxy< T > &a, T b)
TVector3 FindVertex(const std::vector< Vector2D > &xhits, const std::vector< Vector2D > &yhits, TH2F **xmap=0, TH2F **ymap=0) const
unsigned int fPadV
Accumulation space dimensions; bounding box of hits in slice.
Float_t w
Definition: plot.C:20
unsigned int fZoomPadSize
Size of pad XvX [cm] for mm-precision (zoom) scan.
float GetBinContent(int ix, int iy) const
Definition: LiteTH2.h:60
int FindBinY(float y) const
Definition: LiteTH2.h:75