RLFit.cxx
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////
2 /// \file RLFit.cxx
3 /// \brief Make a linear fit to a collection of raw hits
4 /// \version $Id: RLFit.cxx,v 1.9 2011-11-23 20:57:37 bckhouse Exp $
5 /// \author messier@indiana.edu
6 /////////////////////////////////////////////////////////////////////
7 #include <cassert>
8 #include <cmath>
9 
11 
13 #include "Geometry/Geometry.h"
14 #include "GeometryObjects/Geo.h"
16 #include "RawData/RawDigit.h"
17 #include "TrackFit/RLFit.h"
18 #include "Utilities/func/MathUtil.h"
19 
20 namespace trk
21 {
22  //.................................................................
24  unsigned int i1,
25  unsigned int i2) :
26  fR(8),
27  fT(8)
28  {
29  //
30  // The schedule of the anealing. Units are cm^2. First iteration is
31  // not used as weights come from the seeds.
32  //
33  fR[0] = 0.0; fT[0] = 0.0;
34  fR[1] = 36.0; fT[1] = 1000.0;
35  fR[2] = 36.0; fT[2] = 500.0;
36  fR[3] = 36.0; fT[3] = 250.0;
37  fR[4] = 36.0; fT[4] = 100.0;
38  fR[5] = 36.0; fT[5] = 50.0;
39  fR[6] = 36.0; fT[6] = 25.0;
40  fR[7] = 36.0; fT[7] = 10.0;
41  //
42  // The geometry box
43  //
45  fBoxXlo = -geo->DetHalfWidth();
46  fBoxXhi = geo->DetHalfWidth();
47  fBoxYlo = -geo->DetHalfHeight();
48  fBoxYhi = geo->DetHalfHeight();
49  fBoxZlo = 0.0;
50  fBoxZhi = geo->DetLength();
51 
52  //
53  // Pull out the hit information we need for a fit
54  //
56  for (unsigned int i=i1; i<=i2; ++i) {
57  int plane = cmap->GetPlane(d[i].get());
58  int cell = cmap->GetCell(d[i].get());
59  int view = geo->Plane(plane)->View();
60 
61  double xyz[3];
62  geo->Plane(plane)->Cell(cell)->GetCenter(xyz);
63 
64  if (view==geo::kX) {
65  fXhit.push_back(d[i].get());
66  fXx.push_back(xyz[0]);
67  fZx.push_back(xyz[2]);
68  fWx.push_back(0);
69  }
70  else if (view==geo::kY) {
71  fYhit.push_back(d[i].get());
72  fYy.push_back(xyz[1]);
73  fZy.push_back(xyz[2]);
74  fWy.push_back(0);
75  }
76  }
77  this->SeedWeights(fXx, fZx, fYy, fZy, fWx);
78  this->SeedWeights(fYy, fZy, fXx, fZx, fWy);
79  }
80 
81  //......................................................................
82  void RLFit::SeedWeights(const std::vector<double>& x,
83  const std::vector<double>& zx,
84  const std::vector<double>& y,
85  const std::vector<double>& zy,
86  std::vector<double>& w)
87  {
88  for (unsigned int i=0; i<zx.size(); ++i) {
89  double d;
90  w[i] = 1.0E-9;
91  //
92  // Find closest hit in same view
93  //
94  for (unsigned int j=0; j<zx.size(); ++j) {
95  if (i==j) continue;
96  d = util::pythag(zx[i]-zx[j], x[i]-x[j]);
97  if (d<30.0) w[i] += 1.0;
98  }
99  //
100  // Find closest hit in orthogonal view
101  //
102  for (unsigned int j=0; j<zy.size(); ++j) {
103  d = fabs(zx[i]-zy[j]);
104  if (d<30.0) w[i] += 0.1;
105  }
106  }
107  }
108 
109  //......................................................................
110  void RLFit::Fit()
111  {
112  //
113  // Make a line fit to the hits. To improve preformance in a noisy
114  // environment, use a deterministic annealing algoritm
115  //
116  double x1 = 0;
117  double x2 = 0;
118  double zx1 = 0;
119  double zx2 = 0;
120  double y1 = 0;
121  double y2 = 0;
122  double zy1 = 0;
123  double zy2 = 0;
124  for (unsigned int i=0; i<fR.size(); ++i) {
125  static const double eps = 1.0E-6;
126  //
127  // Calculate weights for this iteration. First iteration uses seed
128  // weights
129  //
130  if (i>0) {
131  for (unsigned int j=0; j<fWx.size(); ++j) {
132  double d = geo::DsqrToLine(fZx[j], fXx[j], zx1, x1, zx2, x2);
133 
134  fWx[j] = (1.0+exp(-fR[i]/fT[i]))/(1.0+exp((d-fR[i])/fT[i]));
135  //
136  // Penalty for hits outside the detector in the other view
137  //
138  double y = fMzy*fZx[j] + fBzy;
139  if (y<fBoxYlo) fWx[j] *= exp(-(fBoxYlo-y)/fT[i]);
140  if (y>fBoxYhi) fWx[j] *= exp(-(y-fBoxYhi)/fT[i]);
141 
142  if (fWx[j]<eps) fWx[j] = eps;
143  }
144  for (unsigned int j=0; j<fWy.size(); ++j) {
145  double d = geo::DsqrToLine(fZy[j], fYy[j], zy1, y1, zy2, y2);
146 
147  fWy[j] = (1.0+exp(-fR[i]/fT[i]))/(1.0+exp((d-fR[i])/fT[i]));
148  //
149  // Penalty for hits outside the detector in the other view
150  //
151  double x = fMzx*fZy[j] + fBzx;
152  if (x<fBoxXlo) fWy[j] *= exp(-(fBoxXlo-x)/fT[i]);
153  if (x>fBoxXhi) fWy[j] *= exp(-(x-fBoxXhi)/fT[i]);
154 
155  if (fWy[j]<eps) fWy[j] = eps;
156  }
157  }
158  //
159  // Perform the fits to the two views
160  //
161  geo::LinFitMinDperp(fZx, fXx, fWx, &zx1, &x1, &zx2, &x2);
162  geo::LinFitMinDperp(fZy, fYy, fWy, &zy1, &y1, &zy2, &y2);
163  //
164  // Compute slopes and intercepts
165  //
166  double dzx = zx2-zx1;
167  if (dzx<=0.0) dzx = 1.0E-9;
168  fMzx = (x2-x1)/dzx;
169  fBzx = (x1*zx2-x2*zx1)/dzx;
170 
171  double dzy = zy2-zy1;
172  if (dzy<=0.0) dzy = 1.0E-9;
173  fMzy = (y2-y1)/dzy;
174  fBzy = (y1*zy2-y2*zy1)/dzy;
175  }
176  fChi2 = 0.0;
177  fDOF = 0.0;
178  for (unsigned int j=0; j<fWx.size(); ++j) {
179  double d = geo::DsqrToLine(fZx[j], fXx[j], zx1, x1, zx2, x2);
180  fChi2 += fWx[j]*d;
181  fDOF += fWx[j];
182  }
183  for (unsigned int j=0; j<fWy.size(); ++j) {
184  double d = geo::DsqrToLine(fZy[j], fYy[j], zy1, y1, zy2, y2);
185  fChi2 += fWy[j]*d;
186  fDOF += fWy[j];
187  }
188 
189  // Force north-going assumption on both views...
190  if (zx1>zx2) { std::swap(x1, x2); std::swap(zx1, zx2); }
191  if (zy1>zy2) { std::swap(y1, y2); std::swap(zy1, zy2); }
192 
193  //
194  // Compute slopes and intercepts
195  //
196  double dzx = zx2-zx1;
197  if (dzx<=0.0) dzx = 1.0E-9;
198  fMzx = (x2-x1)/dzx;
199  fBzx = (x1*zx2-x2*zx1)/dzx;
200 
201  double dzy = zy2-zy1;
202  if (dzy<=0.0) dzy = 1.0E-9;
203  fMzy = (y2-y1)/dzy;
204  fBzy = (y1*zy2-y2*zy1)/dzy;
205 
206  //
207  // Adjust z ranges to completely contain hits on the track. "On
208  // track" is defined as having significant weight
209  //
210  fZ1 = fZx[0];
211  fZ2 = fZx[0]+0.01;
212  for (unsigned int i=0; i<fWx.size(); ++i) {
213  if (fWx[i]>0.001) {
214  if (fZx[i]<fZ1) fZ1 = fZx[i];
215  if (fZx[i]>fZ2) fZ2 = fZx[i];
216  }
217  }
218  for (unsigned int i=0; i<fWy.size(); ++i) {
219  if (fWy[i]>0.001) {
220  if (fZy[i]<fZ1) fZ1 = fZy[i];
221  if (fZy[i]>fZ2) fZ2 = fZy[i];
222  }
223  }
224  //
225  // Use line fits to compute x and y locations at the start and end z
226  // positions
227  //
228  fX1 = fMzx*fZ1 + fBzx;
229  fX2 = fMzx*fZ2 + fBzx;
230 
231  fY1 = fMzy*fZ1 + fBzy;
232  fY2 = fMzy*fZ2 + fBzy;
233  }
234 
235  //......................................................................
237  {
238  if (fZ1>fZ2) {
239  std::swap(fX1, fX2);
240  std::swap(fY1, fY2);
241  std::swap(fZ1, fZ2);
242  this->UpdateSlopes();
243  }
244  }
245 
246  //......................................................................
248  {
249  if (fY2>fY1) {
250  std::swap(fX1, fX2);
251  std::swap(fY1, fY2);
252  std::swap(fZ1, fZ2);
253  this->UpdateSlopes();
254  }
255  }
256 
257  //......................................................................
259  {
260  fMzx = (fX2-fX1)/(fZ2-fZ1);
261  fBzx = (fX1*fZ2-fX2*fZ1)/(fZ2-fZ1);
262 
263  fMzy = (fY2-fY1)/(fZ2-fZ1);
264  fBzy = (fY1*fZ2-fY2*fZ1)/(fZ2-fZ1);
265  }
266 
267  //......................................................................
269  {
270  this->AssumeNorthGoing();
271 
272  double dz = (fZ2-fZ1);
273  double dx = (fX2-fX1);
274  double dy = (fY2-fY1);
275 
276  return dz/sqrt(dx*dx+dy*dy+dz*dz);
277  }
278 
279  //......................................................................
280  double RLFit::PhiBeam()
281  {
282  this->AssumeNorthGoing();
283 
284  double dx = (fX2-fX1);
285  double dy = (fY2-fY1);
286 
287  return atan2(dy,dx);
288  }
289 
290  //......................................................................
292  {
293  this->AssumeDownGoing();
294 
295  double dz = (fZ2-fZ1);
296  double dx = (fX2-fX1);
297  double dy = (fY2-fY1);
298 
299  return -dy/sqrt(dx*dx+dy*dy+dz*dz);
300  }
301 
302  //......................................................................
304  {
305  this->AssumeDownGoing();
306 
307  double dx = (fX2-fX1);
308  double dz = (fZ2-fZ1);
309 
310  return atan2(-dx,-dz);
311  }
312 
313  //......................................................................
314  void RLFit::SetRTSchedule(const std::vector<double>& r,
315  const std::vector<double>& t)
316  {
317  assert(r.size()==r.size());
318  fR.resize(r.size());
319  fT.resize(t.size());
320  for (unsigned int i=0; i<r.size(); ++i) {
321  fR[i] = r[i];
322  fT[i] = t[i];
323  }
324  }
325 
326 } // end namespace trk
327 ////////////////////////////////////////////////////////////////////////
double PhiBeam()
Definition: RLFit.cxx:280
double fDOF
Sum of weights * deltas.
Definition: RLFit.h:120
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
std::vector< double > fXx
x positions (cm)
Definition: RLFit.h:93
double fY2
End y positon (cm)
Definition: RLFit.h:109
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double fBoxYlo
Detector box Y edge.
Definition: RLFit.h:114
double fMzy
Slope in y-z view.
Definition: RLFit.h:102
Float_t y1[n_points_granero]
Definition: compare.C:5
std::vector< const rawdata::RawDigit * > fYhit
Y view hits.
Definition: RLFit.h:92
Float_t x1[n_points_granero]
Definition: compare.C:5
std::vector< double > fZx
z positions (cm)
Definition: RLFit.h:94
double fBoxZhi
Detector box Z edge.
Definition: RLFit.h:117
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::vector< double > fWy
y view weights
Definition: RLFit.h:98
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::vector< const rawdata::RawDigit * > fXhit
X view hits.
Definition: RLFit.h:91
double DetLength() const
std::vector< double > fWx
x view weights
Definition: RLFit.h:95
const PlaneGeo * Plane(unsigned int i) const
double PhiCosmic()
Definition: RLFit.cxx:303
void UpdateSlopes()
Definition: RLFit.cxx:258
double fMzx
Slope in x-z view.
Definition: RLFit.h:100
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
double fX1
Start x position (cm)
Definition: RLFit.h:105
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
double fY1
Start y position (cm)
Definition: RLFit.h:106
Track finder for cosmic rays.
std::vector< double > fYy
y positions (cm)
Definition: RLFit.h:96
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
void AssumeDownGoing()
Definition: RLFit.cxx:247
double fZ1
Start z position (cm)
Definition: RLFit.h:107
double fBoxXhi
Detector box X edge.
Definition: RLFit.h:113
double dy[NP][NC]
double dx[NP][NC]
double CosThetaCosmic()
Definition: RLFit.cxx:291
double dz[NP][NC]
void AssumeNorthGoing()
Definition: RLFit.cxx:236
double LinFitMinDperp(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared perpendicular dista...
Definition: Geo.cxx:449
Make a linear fit to a collection of raw hits.
void SeedWeights(const std::vector< double > &x, const std::vector< double > &zx, const std::vector< double > &y, const std::vector< double > &zy, std::vector< double > &w)
Definition: RLFit.cxx:82
void SetRTSchedule(const std::vector< double > &r, const std::vector< double > &t)
Definition: RLFit.cxx:314
void Fit()
Definition: RLFit.cxx:110
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
Collect Geo headers and supply basic geometry functions.
Float_t d
Definition: plot.C:236
unsigned short GetPlane(const rawdata::RawDigit *dig)
Definition: CMap.cxx:285
std::vector< double > fR
Range of weight function.
Definition: RLFit.h:87
const double j
Definition: BetheBloch.cxx:29
double fBoxZlo
Detector box Z edge.
Definition: RLFit.h:116
double DetHalfHeight() const
std::vector< double > fZy
z positions (cm)
Definition: RLFit.h:97
std::vector< double > fT
Transition width of weight function.
Definition: RLFit.h:88
double CosThetaBeam()
Definition: RLFit.cxx:268
double fBzy
Intercept in y-z view.
Definition: RLFit.h:103
double DetHalfWidth() const
double fZ2
End z position (cm)
Definition: RLFit.h:110
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
double fChi2
Definition: RLFit.h:119
double fBzx
Intercept in x-z view.
Definition: RLFit.h:101
RLFit(const std::vector< art::Ptr< rawdata::RawDigit > > &d, unsigned int i1, unsigned int i2)
Definition: RLFit.cxx:23
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
double fBoxXlo
Detector box X edge.
Definition: RLFit.h:112
unsigned short GetCell(const rawdata::RawDigit *dig)
Definition: CMap.cxx:327
Helper for AttenCurve.
Definition: Path.h:10
Float_t w
Definition: plot.C:20
Definition: fwd.h:28
Encapsulate the geometry of one entire detector (near, far, ndos)
double fX2
End x position (cm)
Definition: RLFit.h:108
double fBoxYhi
Detector box Y edge.
Definition: RLFit.h:115
T atan2(T number)
Definition: d0nt_math.hpp:72