Lutz.cxx
Go to the documentation of this file.
1 ///
2 /// \file Lutz.cxx
3 /// \brief "Breakpoint" track fitter, G. Lutz, NIM A273 (1988) 349-361
4 /// \author messier@indiana.edu
5 ///
6 
7 #include "Lutz.h"
9 
10 #include <iostream>
11 
12 namespace bpfit {
13  Lutz::Lutz(unsigned int n, unsigned int N) :
14  fn(n),
15  fN(N),
16  fz(n),
17  fXI(n),
18  fwx(n),
19  fZ(N),
20  fwSJ(N),
21  fZeta(n),
22  fD(N+2,n),
23  fAinv(N+2),
24  fa(0),
25  fb(0),
26  fAlpha(fN),
27  fx(n)
28  {
29  if(N >= n)
30  throw BPException(__FILE__, __LINE__, kNn);
31 
32  for(unsigned int i = 0; i < n; i++)
33  fZeta[i].resize(N);
34  }
35 
36  //......................................................................
37 
38  void Lutz::SetMeasurement(unsigned int i, double z, double x, double sigx)
39  {
40  if(i >= fn)
41  throw BPException(__FILE__, __LINE__, kINDEX);
42  if(z != z || x != x || sigx != sigx)
43  throw BPException(__FILE__, __LINE__, kNAN);
44  if(sigx <= 0.0)
45  throw BPException(__FILE__, __LINE__, kSIGX, sigx);
46 
47  fz[i] = z;
48  fXI[i] = x;
49  fwx[i] = 1.0/(sigx*sigx);
50  }
51 
52  //......................................................................
53 
54  void Lutz::SetScatteringPlane(unsigned int j, double Z, double sigSJ)
55  {
56  if(j >= fN)
57  throw BPException(__FILE__, __LINE__, kINDEX);
58  if(Z != Z || sigSJ != sigSJ)
59  throw BPException(__FILE__, __LINE__, kNAN);
60  if(sigSJ <= 0.0)
61  throw BPException(__FILE__, __LINE__, kSIGSJ, sigSJ);
62 
63  fZ[j] = Z;
64  fwSJ[j] = 1.0/(sigSJ*sigSJ);
65  }
66 
67  //......................................................................
68 
70  {
71  unsigned int i, j;
72  for (i=0; i<fn; ++i) {
73  for (j=0; j<fN; ++j) {
74  fZeta[i][j] = fz[i]-fZ[j];
75  if (fZeta[i][j]<0.0) fZeta[i][j] = 0.0;
76  }
77  }
78  }
79 
80  //......................................................................
81  //
82  // Equations 31.
83  //
84  void Lutz::CalcD()
85  {
86  unsigned int i, j;
87  for (i=0; i<fn; ++i) {
88  fD[0][i] = fwx[i];
89  fD[1][i] = fwx[i]*fz[i];
90  for (j=0; j<fN; ++j) {
91  fD[j+2][i] = fwx[i]*fZeta[i][j];
92  }
93  }
94  }
95 
96  //......................................................................
97  //
98  // Equations 29.
99  //
101  {
102  unsigned int i, J, K;
103 
104  double sum00 = 0;
105  double sum01 = 0;
106  double sum11 = 0;
107  double fwxfz;
108  for (i=0; i<fn; ++i) {
109  fwxfz = fwx[i]*fz[i];
110  sum00 += fwx[i];
111  sum01 += fwxfz;
112  sum11 += fwxfz*fz[i];
113  }
114  fAinv[0][0] = sum00;
115  fAinv[0][1] = sum01;
116  fAinv[1][0] = sum01;
117  fAinv[1][1] = sum11;
118 
119  double sum0J2;
120  double sum1J2;
121  double fwxZeta;
122  for (J=0; J<fN; ++J) {
123  sum0J2 = 0;
124  sum1J2 = 0;
125  for (i=0; i<fn; ++i) {
126  fwxZeta = fwx[i]*fZeta[i][J];
127  sum0J2 += fwxZeta;
128  sum1J2 += fz[i]*fwxZeta;
129  }
130  fAinv[0][J+2] = sum0J2;
131  fAinv[J+2][0] = sum0J2;
132  fAinv[1][J+2] = sum1J2;
133  fAinv[J+2][1] = sum1J2;
134  }
135 
136  double sum;
137  for (J=0; J<fN; ++J) {
138  for (K=J; K<fN; ++K) {
139  sum = 0;
140  for (i=0; i<fn; ++i) sum += fwx[i]*fZeta[i][J]*fZeta[i][K];
141  if (J==K) sum += fwSJ[J];
142  fAinv[J+2][K+2] = sum;
143  fAinv[K+2][J+2] = sum;
144  }
145  }
146 
147  double det = 0.0;
148 
149  /// Inverting the matrix
150  /// Sometimes this operation fails, producing error: matrix is singular
151  /// Thus we need to put it into try/catch block to avoid the hard exit.
152  try {
153  fAinv.InvertFast(&det);
154  }
155  catch(...) {
156  std::cerr << __FILE__ << ":" << __LINE__
157  << " Matrix inversion failed" << std::endl;
158  return false;
159  } // end of inverting the matrix
160 
161  if (det<=0.0) {
162  std::cerr << __FILE__ << ":" << __LINE__
163  << " Matrix inversion failed" << std::endl;
164  return false;
165  }
166  return true;
167  }
168 
169  //......................................................................
170  //
171  // Equation 35
172  //
173  void Lutz::Calcq()
174  {
175  unsigned int j;
176  //
177  // Solve according to equation 35. Don't bother with the second term
178  // which is by construction zero.
179  //
181  q *= fD;
182  q *= fAinv;
183 
184  //
185  // Unpack the resutls
186  //
187  fa = q[0];
188  fb = q[1];
189  for (j=0; j<fN; ++j) fAlpha[j] = q[2+j];
190  }
191 
192  //......................................................................
193 
194  void Lutz::Calcx()
195  {
196  unsigned int i, j;
197  for (i=0; i<fn; ++i) {
198  fx[i] = fa + fb*fz[i];
199  for (j=0; j<fN; ++j) {
200  if (fZeta[i][j]>0.0) fx[i] += fAlpha[j]*fZeta[i][j];
201  }
202  }
203  }
204 
205  //......................................................................
206 
207  double Lutz::X(double z) const
208  {
209  unsigned int j;
210  double zeta;
211  double x = fa + fb*z;
212  for (j=0; j<fN; ++j) {
213  zeta = z-fZ[j];
214  if (zeta>0.0) x += fAlpha[j]*zeta;
215  }
216  return x;
217  }
218 
219  //......................................................................
220 
221  void Lutz::Fit(double *a, double* b, double* alpha, double* chi2)
222  {
223  bool ok;
224  unsigned int j;
225  //
226  // Make the fit
227  //
228  this->CalcZeta();
229  this->CalcD();
230  ok = this->CalcAinv();
231  if (ok) {
232  this->Calcq();
233  this->Calcx();
234  }
235  else this->NullResult();
236  //
237  // Unpack the results
238  //
239  if (a!=0) *a = fa;
240  if (b!=0) *b = fb;
241  if (alpha!=0) {
242  for (j=0; j<fN; ++j) alpha[j] = fAlpha[j];
243  }
244  if (chi2!=0) {
245  *chi2 = this->Chi2XI() + this->Chi2Beta();
246  }
247  }
248 
249  //......................................................................
250 
251  double Lutz::Chi2XIi(unsigned int i) const
252  {
253  double chi;
254 
255  chi = fXI[i]-fx[i];
256  chi *= chi;
257  chi *= fwx[i];
258 
259  return chi;
260  }
261 
262  //......................................................................
263 
264  double Lutz::Chi2BetaJ(unsigned int j) const
265  {
266  double chi;
267 
268  chi = 0.0-fAlpha[j];
269  chi *= chi;
270  chi *= fwSJ[j];
271 
272  return chi;
273  }
274 
275  //......................................................................
276 
277  double Lutz::Chi2XI() const
278  {
279  unsigned int i;
280  double sum = 0.0;
281  for (i=0; i<fn; ++i) sum += this->Chi2XIi(i);
282  return sum;
283  }
284 
285  //......................................................................
286 
287  double Lutz::Chi2Beta() const
288  {
289  unsigned int j;
290  double sum = 0.0;
291  for (j=0; j<fN; ++j) sum += this->Chi2BetaJ(j);
292  return sum;
293  }
294 
295  //....................................................................
296 
298  {
299  unsigned int i, j;
300  fa = 0.0;
301  fb = 0.0;
302  for (i=0; i<fn; ++i) fx[i] = 0.0;
303  for (j=0; j<fN; ++j) fAlpha[j] = 0.0;
304  }
305 }
306 
307 ////////////////////////////////////////////////////////////////////////
void Calcx()
Definition: Lutz.cxx:194
TMatrixTSym< double > fAinv
System of equations, Eqn. 29.
Definition: Lutz.h:99
std::vector< double > fz
Measurement plane longitudinal locations.
Definition: Lutz.h:92
const unsigned int kSIGX
non-positive value of sigx
Definition: BPException.h:23
double Chi2XIi(unsigned int i) const
After fit, contribution to chi^2 of the ith measurement.
Definition: Lutz.cxx:251
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
void SetScatteringPlane(unsigned int j, double Z, double sigSJ)
Set Z location of a scattering plane j.
Definition: Lutz.cxx:54
void SetMeasurement(unsigned int i, double z, double x, double sigx)
Add measurements.
Definition: Lutz.cxx:38
Lutz(unsigned int n, unsigned int N)
Construct a fitter for n measurements and N scattering planes.
Definition: Lutz.cxx:13
const unsigned int kNAN
nan encountered
Definition: BPException.h:22
std::vector< double > fZ
Location of scattering planes.
Definition: Lutz.h:95
OStream cerr
Definition: OStream.cxx:7
"Break-point" track fitter
unsigned int fn
Number of measurements.
Definition: Lutz.h:90
double fa
Fitted intercept.
Definition: Lutz.h:100
double chi2()
double X(double z) const
After fit, best-fit x location of track at location z.
Definition: Lutz.cxx:207
Double_t K
const unsigned int kNn
number of scat surf >= number of measurements
Definition: BPException.h:20
Float_t Z
Definition: plot.C:38
void Fit(double *a, double *b, double *alpha, double *chi2)
Do the fit.
Definition: Lutz.cxx:221
TMatrixT< double > fD
Eqn. 31.
Definition: Lutz.h:98
TVectorT< double > fXI
Measurements in transverse location.
Definition: Lutz.h:93
const double a
void resize(T &x, std::vector< int > dims)
Definition: resize.hpp:41
double fb
Fitted slope.
Definition: Lutz.h:101
std::vector< double > fwx
Measurement weight (1/sigx^2)
Definition: Lutz.h:94
double Chi2XI() const
Definition: Lutz.cxx:277
std::vector< double > fwSJ
Scattering angle weight (1/sigSJ^2)
Definition: Lutz.h:96
const double j
Definition: BetheBloch.cxx:29
double Chi2BetaJ(unsigned int j) const
After fit, contribution to chi^2 of the jth scattering plane.
Definition: Lutz.cxx:264
z
Definition: test.py:28
const unsigned int kINDEX
index problem
Definition: BPException.h:21
std::vector< double > fAlpha
Fitted scatters.
Definition: Lutz.h:102
const unsigned int kSIGSJ
non-positive value of sigSJ
Definition: BPException.h:24
std::vector< std::vector< double > > fZeta
measurement-scattering plane distances
Definition: Lutz.h:97
unsigned int fN
Number of scattering planes.
Definition: Lutz.h:91
void Calcq()
Definition: Lutz.cxx:173
const hit & b
Definition: hits.cxx:21
void CalcZeta()
Definition: Lutz.cxx:69
bool CalcAinv()
Definition: Lutz.cxx:100
Double_t sum
Definition: plot.C:31
double Chi2Beta() const
Definition: Lutz.cxx:287
void CalcD()
Definition: Lutz.cxx:84
void NullResult()
Definition: Lutz.cxx:297
std::vector< double > fx
Fitted track locations.
Definition: Lutz.h:103
TFile fa("Li7.root")