GridSearch.cxx
Go to the documentation of this file.
1 ///
2 /// \file GridSearch.cxx
3 /// \brief Produce a list of vertex and prong seeds for ElasticArms
4 /// \version $Id: GridSearch.cxx,v 1.7 2012-10-04 21:23:23 messier Exp $
5 /// \author messier@indiana.edu
6 ///
9 
10 #include <iostream>
11 #include <cmath>
12 #include <algorithm>
13 
14 #include "RecoBase/HoughResult.h"
15 
16 namespace earms {
17 
18  //......................................................................
19 
20  GridSearch::GridSearch(double xmn, double xmx,
21  double ymn, double ymx,
22  double zmn, double zmx,
23  double colinear,
24  unsigned int nseeddir) :
25  fXmin(xmn),
26  fXmax(xmx),
27  fYmin(ymn),
28  fYmax(ymx),
29  fZmin(zmn),
30  fZmax(zmx),
31  fColinear(colinear),
32  fNseedDir(nseeddir)
33  { }
34 
35  //......................................................................
36 
37  void GridSearch::AddVertex(double x, double y, double z)
38  {
39  if (x<fXmin) return;
40  if (y<fYmin) return;
41  if (z<fZmin) return;
42  if (x>fXmax) return;
43  if (y>fYmax) return;
44  if (z>fZmax) return;
45  fX.push_back(x);
46  fY.push_back(y);
47  fZ.push_back(z);
48  }
49 
50  //......................................................................
51 
53  const std::vector<double>& mx,
54  const std::vector<double>& bx,
55  const std::vector<double>& my,
56  const std::vector<double>& by,
57  const std::vector<double>& f)
58  {
59  unsigned int i, j;
60  std::vector< std::pair<double, double> > xz;
61  std::vector< std::pair<double, double> > yz;
62 
63  //
64  // Sort the xz and yz hits by z
65  //
66  for (i=0; i<arms.fN; ++i) {
67  if (arms.fView[i]==0) {
68  xz.push_back(std::pair<double,double>(arms.fZ[i], arms.fXorY[i]));
69  }
70  else if (arms.fView[i]==1) {
71  yz.push_back(std::pair<double,double>(arms.fZ[i], arms.fXorY[i]));
72  }
73  else abort();
74  }
75  sort(xz.begin(), xz.end());
76  sort(yz.begin(), yz.end());
77 
78  //
79  // Grab points from the specified locations in the hit lists
80  //
81  unsigned int ix;
82  unsigned int iy;
83  double xseed, zxseed;
84  double yseed, zyseed;
85  for (i=0; i<f.size(); ++i) {
86  ix = fabs(f[i])*xz.size()/100;
87  iy = fabs(f[i])*yz.size()/100;
88 
89 
90  //
91  // In forward direction, pick the hits
92  //
93  if (f[i]>0.0) {
94  zxseed = xz[ix].first;
95  xseed = xz[ix].second;
96  zyseed = yz[iy].first;
97  yseed = yz[iy].second;
98  }
99  else {
100  //
101  // In backward direction, pick mirror images of the hits
102  // reflected across the boundary of the first hit
103  //
104  zxseed = 2.0*xz[0].first -xz[ix].first;
105  xseed = 2.0*xz[0].second-xz[ix].second;
106  zyseed = 2.0*yz[0].first -yz[iy].first;
107  yseed = 2.0*yz[0].second-yz[iy].second;
108  }
109 
110  //
111  // A point on the line in the yz view
112  //
113  for (j=0; j<my.size(); ++j) {
114  this->AddVertex(xseed, my[j]*zxseed+by[j], zxseed);
115  }
116  //
117  // A point on the line in the xz view
118  //
119  for (j=0; j<mx.size(); ++j) {
120  this->AddVertex(mx[j]*zyseed+bx[j], yseed, zyseed);
121  }
122  //
123  // A point that takes information just from the hits
124  //
125  this->AddVertex(xseed, yseed, 0.5*(zxseed+zyseed));
126  }
127  }
128 
129  //......................................................................
130 
131  void GridSearch::AddDirection(double theta, double phi)
132  {
133  fTheta.push_back(theta);
134  fPhi. push_back(phi);
135  }
136 
137  //......................................................................
138 
140  const rb::HoughResult& hry,
141  unsigned int nmx,
142  unsigned int nmatch,
143  bool flipped)
144  {
145  //
146  // For there to be any intersections we need at minimum of two lines
147  // in x and one in y
148  //
149  if ((hrx.fPeak.size()<2) && (hry.fPeak.size()<1)){
150  //try flipping before we exit
151  if (flipped==false) {
152  this->AddHoughIntersections(hry,hrx,nmx,nmatch,true);
153  }
154  return;
155  }
156 
157  //
158  // Only use the top set of results
159  //
160  unsigned int nx = hrx.fPeak.size();
161  if (nx>nmx) nx = nmx;
162 
163  unsigned int i, j, k;
164  //
165  // Loop over all the intersections in the x view. Grab y location
166  // from the dominant line in that view.
167  //
168  double x, y, z;
169  double my, by;
170  double dot;
171 
172  unsigned int ny = nmatch;
173  if (hry.fPeak.size()<ny) ny = hry.fPeak.size();
174 
175  for (i=0; i<nx; ++i) {
176  double mi, bi;
177  hrx.SlopeIcept(i, &mi, &bi);
178  for (j=i+1; j<nx; ++j) {
179  double mj, bj;
180  hrx.SlopeIcept(j, &mj, &bj);
181 
182  // Skip lines paired with themselves
183  if (mi==mj) continue;
184 
185  // Also skip lines that are nearly colinear
186  dot = fabs((mi*mj+1.0)/sqrt(mi*mi+1.0)/sqrt(mj*mj+1.0));
187  if (dot>fColinear) continue;
188 
189  z = (bj-bi)/(mi-mj);
190  x = 0.5*((mi+mj)*z+(bi+bj));
191 
192 
193  //
194  // What y value to assign to this x-z pair? Check the top
195  // nmatch hough lines in the other view.
196  //
197  for (k=0; k<ny; ++k) {
198  hry.SlopeIcept(k, &my, &by);
199 
200  y = my*z + by;
201 
202 
203  if (flipped) this->AddVertex(y,x,z);
204  else this->AddVertex(x,y,z);
205  }
206 
207  } // loop on j hough lines in x view
208  } // loop on i hough lines x view
209 
210  //
211  // Repeat exchanging the meaning of x and y
212  //
213  if (flipped==false) {
214  this->AddHoughIntersections(hry,hrx,nmx,nmatch,true);
215  }
216  }
217 
218  //......................................................................
219 
220  ///
221  /// Set a standard set of directions uniformly distributed on a
222  /// sphere using the "Golden Section spiral"
223  /// (http://www.xsi-blog.com/archives/115)
224  ///
226  {
227  static const double nd = fNseedDir;
228  static const double inc = M_PI*(3.0-sqrt(5.0));
229  static const double off = 2.0/nd;
230 
231  unsigned int i;
232  double x, y, z;
233  double r, phi;
234  for (i=0; i<fNseedDir; ++i) {
235  y = (i+0.5)*off-1.0;
236  r = sqrt(1.0-y*y);
237  phi = i*inc;
238  x = r*cos(phi);
239  z = r*sin(phi);
240 
241  fTheta.push_back(acos(z/sqrt(x*x+y*y+z*z)));
242  fPhi. push_back(atan2(y,x));
243  }
244  }
245 
246  //......................................................................
247 
248  double GridSearch::ChooseBestProngDir(ElasticArms& arms, unsigned int a)
249  {
250  unsigned int b;
251  unsigned int j, jbest=0;
252  double E = 99E99;
253  double Ebest = 99E99;
254 
255  //
256  // Loop over all directions in the table
257  //
258  for (j=0; j<fTheta.size(); ++j) {
259  arms.SetArm(a, fTheta[j], fPhi[j]);
260  //
261  // Check this direction against arms which have already been set.
262  // If it's too close to an existing arm, skip it.
263  //
264  double dot;
265  double maxdot = -1.0;
266  double mindot = 1.0;
267  for (b=0; b<a; ++b) {
268  dot =
269  arms.fdXds[a]*arms.fdXds[b]+
270  arms.fdYds[a]*arms.fdYds[b]+
271  arms.fdZds[a]*arms.fdZds[b];
272  if (dot>maxdot) maxdot = dot;
273  if (dot<mindot) mindot = dot;
274  }
275  if (maxdot> fColinear) continue;
276  if (mindot<-fColinear) continue;
277 
278  //
279  // Test the goodness of this prong. Remember it if its the best so
280  // far.
281  //
282  arms.UpdateMia();
283  arms.UpdateVia();
284  E = arms.E();
285  if (E<Ebest) {
286  jbest = j;
287  Ebest = E;
288  }
289  }
290  arms.SetArm(a, fTheta[jbest], fPhi[jbest]);
291  return Ebest;
292  }
293 
294  //......................................................................
295 
297  {
298  unsigned int i, j;
299  double E = 99E99;
300  double Ebest = 99E99;
301  double xbest=0;
302  double ybest=0;
303  double zbest=0;
304  for (i=0; i<fX.size(); ++i) {
305  arms.SetVertex(fX[i], fY[i], fZ[i]);
306  for (j=0; j<arms.fM; ++j) {
307  arms.SetArm(j, 0.5*M_PI, 2.0*M_PI*(float)j/(float)arms.fM);
308  }
309  for (j=0; j<arms.fM; ++j){
310  E = this->ChooseBestProngDir(arms, j);
311  }
312  if (E<Ebest) {
313  Ebest = E;
314  xbest = fX[i];
315  ybest = fY[i];
316  zbest = fZ[i];
317  }
318  }
319  //
320  // Leave the vertex and arm directions in their best configuration
321  //
322  arms.SetVertex(xbest, ybest, zbest);
323  for (j=0; j<arms.fM; ++j) {
324  arms.SetArm(j, 0.5*M_PI, 2.0*M_PI*(float)j/(float)arms.fM);
325  }
326  for (j=0; j<arms.fM; ++j) E = this->ChooseBestProngDir(arms, j);
327  return E;
328  }
329 } // end namespace earms
330 ////////////////////////////////////////////////////////////////////////
double fYmin
Only consider vertex locations inside a box [cm].
Definition: GridSearch.h:108
std::vector< double > fPhi
List of prong directions.
Definition: GridSearch.h:119
void AddDirection(double theta, double phi)
Definition: GridSearch.cxx:131
void SetStandardDirections()
Definition: GridSearch.cxx:225
std::vector< HoughPeak > fPeak
List of peaks found in Hough space.
Definition: HoughResult.h:61
double fColinear
|dot product|>fColinear means lines are colinear
Definition: GridSearch.h:112
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::vector< double > fZ
List of vertex z coorfinates.
Definition: GridSearch.h:117
std::vector< double > fTheta
List of prong directions.
Definition: GridSearch.h:118
T acos(T number)
Definition: d0nt_math.hpp:54
void AddVtxPoints(const ElasticArms &arms, const std::vector< double > &mx, const std::vector< double > &bx, const std::vector< double > &my, const std::vector< double > &by, const std::vector< double > &f)
Definition: GridSearch.cxx:52
#define M_PI
Definition: SbMath.h:34
Elastic Arms vertex finding algorithm.
void AddHoughIntersections(const rb::HoughResult &hrx, const rb::HoughResult &hry, unsigned int nmx, unsigned int nmh, bool flipped=false)
Definition: GridSearch.cxx:139
unsigned int fNseedDir
How many default seed directoins to check.
Definition: GridSearch.h:113
double fZmax
Only consider vertex locations inside a box [cm].
Definition: GridSearch.h:111
unsigned int fM
Number of arms.
Definition: ElasticArms.h:129
std::vector< double > fXorY
Transverse hit positions (cm)
Definition: ElasticArms.h:109
std::vector< double > fdZds
Track dz/ds.
Definition: ElasticArms.h:132
void SlopeIcept(unsigned int i, double *m, double *b) const
Slope and intercepts of Hough lines.
Definition: HoughResult.cxx:62
std::vector< double > fY
List of vertex y coordinates.
Definition: GridSearch.h:116
Float_t E
Definition: plot.C:20
const double a
const double j
Definition: BetheBloch.cxx:29
unsigned int fN
Data to fit.
Definition: ElasticArms.h:108
double fXmax
Only consider vertex locations inside a box [cm].
Definition: GridSearch.h:107
z
Definition: test.py:28
void AddVertex(double x, double y, double z)
Definition: GridSearch.cxx:37
double FindBestVtx(ElasticArms &arms)
Definition: GridSearch.cxx:296
double ChooseBestProngDir(ElasticArms &arms, unsigned int i)
Definition: GridSearch.cxx:248
std::vector< double > fZ
Z locations of hits (cm)
Definition: ElasticArms.h:111
T sin(T number)
Definition: d0nt_math.hpp:132
double fXmin
Only consider vertex locations inside a box [cm].
Definition: GridSearch.h:106
void SetArm(unsigned int i, double theta, double phi)
double fYmax
Only consider vertex locations inside a box [cm].
Definition: GridSearch.h:109
Data resulting from a Hough transform on the cell hit positions.
std::vector< double > fdXds
Track dx/ds.
Definition: ElasticArms.h:130
std::vector< double > fX
List of vertex x coordinates.
Definition: GridSearch.h:115
GridSearch(double xmn, double xmx, double ymn, double ymx, double zmn, double zmx, double colinear, unsigned int nseeddir)
Definition: GridSearch.cxx:20
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
TRandom3 r(0)
std::vector< int > fView
Which view is the hit in?
Definition: ElasticArms.h:110
double fZmin
Only consider vertex locations inside a box [cm].
Definition: GridSearch.h:110
std::vector< double > fdYds
Track dy/ds.
Definition: ElasticArms.h:131
Summary from a Hough transform applied to a group of cell hits.
Definition: HoughResult.h:35
void SetVertex(double x, double y, double z)
T atan2(T number)
Definition: d0nt_math.hpp:72