CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions
DTLinearFit Class Reference

#include <DTLinearFit.h>

Public Member Functions

 DTLinearFit ()
 Constructor. More...
 
void fit (const std::vector< float > &x, const std::vector< float > &y, int ndat, const std::vector< float > &sigy, float &slope, float &intercept, double &chi2, float &covss, float &covii, float &covsi) const
 
void fit3par (const std::vector< float > &xfit, const std::vector< float > &yfit, const std::vector< int > &lfit, const int nptfit, const std::vector< float > &sigy, float &aminf, float &bminf, float &cminf, double &chi2fit, const bool debug) const
 
void fit4Var (const std::vector< float > &xfit, const std::vector< float > &yfit, const std::vector< int > &lfit, const std::vector< double > &tfit, const int nptfit, float &aminf, float &bminf, float &cminf, float &vminf, double &chi2fit, const bool vdrift_4parfit, const bool debug) const
 
void fitNpar (const int npar, const std::vector< float > &xfit, const std::vector< float > &yfit, const std::vector< int > &lfit, const std::vector< double > &tfit, const std::vector< float > &sigy, float &aminf, float &bminf, float &cminf, float &vminf, double &chi2fit, const bool debug) const
 
 ~DTLinearFit ()
 Destructor. More...
 

Detailed Description

Description:

detailed description

Author
Stefano Lacaprara - INFN Legnaro stefa.nosp@m.no.l.nosp@m.acapr.nosp@m.ara@.nosp@m.pd.in.nosp@m.fn.i.nosp@m.t

Definition at line 26 of file DTLinearFit.h.

Constructor & Destructor Documentation

DTLinearFit::DTLinearFit ( )

Constructor.

Definition at line 18 of file DTLinearFit.cc.

18  {
19 
20 }
DTLinearFit::~DTLinearFit ( )

Destructor.

Definition at line 23 of file DTLinearFit.cc.

23  {
24 }

Member Function Documentation

void DTLinearFit::fit ( const std::vector< float > &  x,
const std::vector< float > &  y,
int  ndat,
const std::vector< float > &  sigy,
float &  slope,
float &  intercept,
double &  chi2,
float &  covss,
float &  covii,
float &  covsi 
) const

Definition at line 27 of file DTLinearFit.cc.

References delta, i, and j.

Referenced by DTSegmentUpdator::fit(), and trackingPlots.Iteration::modules().

37 {
38  chi2=0;
39  float sy = 0, sx = 0;
40  float s11 = 0, sxy = 0, sxx = 0;
41  for (int i = 0; i != nptfit; i++) {
42  float sy2 = sigy[i]*sigy[i];
43  sy += y[i] / sy2;
44  sxy += x[i]*y[i] / sy2;
45  s11 += 1. / sy2;
46  sx += x[i] / sy2;
47  sxx += x[i]*x[i] / sy2;
48  }
49 
50  float delta = s11*sxx - sx*sx;
51 
52  if (delta==0) return;
53 
54  intercept = (sy*sxx - sxy*sx) / delta;
55  slope = (sxy*s11 - sy*sx) / delta;
56 
57  covii = sxx / delta;
58  covss = s11 / delta;
59  covsi = -sx / delta;
60 
61  for (int j=0; j<nptfit; j++){
62  const double ypred = intercept + slope*x[j];
63  const double dy = (y[j] - ypred)/sigy[j];
64  chi2 += dy*dy;
65  }
66 
67 }
delta
dbl * delta
Definition: mlp_gen.cc:36
i
int i
Definition: DBlmapReader.cc:9
slope
static const double slope[3]
Definition: CastorTimeSlew.cc:6
j
int j
Definition: DBlmapReader.cc:9
beam_dqm_sourceclient-live_cfg.chi2
tuple chi2
Definition: beam_dqm_sourceclient-live_cfg.py:337
void DTLinearFit::fit3par ( const std::vector< float > &  xfit,
const std::vector< float > &  yfit,
const std::vector< int > &  lfit,
const int  nptfit,
const std::vector< float > &  sigy,
float &  aminf,
float &  bminf,
float &  cminf,
double &  chi2fit,
const bool  debug = 0 
) const

Definition at line 220 of file DTLinearFit.cc.

References debug.

229  {
230 
231  float vminf;
232  vector <double> tfit( nptfit, 0.);
233 
234  fitNpar(3,xfit,yfit,lfit,tfit,sigy,aminf,bminf,cminf,vminf,chi2fit,debug);
235 }
DTLinearFit::fitNpar
void fitNpar(const int npar, const std::vector< float > &xfit, const std::vector< float > &yfit, const std::vector< int > &lfit, const std::vector< double > &tfit, const std::vector< float > &sigy, float &aminf, float &bminf, float &cminf, float &vminf, double &chi2fit, const bool debug) const
Definition: DTLinearFit.cc:70
debug
#define debug
Definition: HDRShower.cc:19
void DTLinearFit::fit4Var ( const std::vector< float > &  xfit,
const std::vector< float > &  yfit,
const std::vector< int > &  lfit,
const std::vector< double > &  tfit,
const int  nptfit,
float &  aminf,
float &  bminf,
float &  cminf,
float &  vminf,
double &  chi2fit,
const bool  vdrift_4parfit = 0,
const bool  debug = 0 
) const

Definition at line 238 of file DTLinearFit.cc.

References a, b, gather_cfg::cout, debug, and j.

Referenced by DTSegmentUpdator::calculateT0corr().

249  {
250 
251  if (debug) cout << "Entering Fit4Var" << endl;
252 
253  const double sigma = 0.0295;// FIXME errors can be inserted .just load them/that is the usual TB resolution value for DT chambers
254  vector<float> sigy;
255 
256  aminf = 0.;
257  bminf = 0.;
258  cminf = -999.;
259  vminf = 0.;
260 
261  int nppar = 0;
262  double chi2fitN2 = -1. ;
263  double chi2fit3 = -1.;
264  double chi2fitN3 = -1. ;
265  double chi2fitN4 = -1.;
266  float bminf3 = bminf;
267  float aminf3 = aminf;
268  float cminf3 = cminf;
269  int nppar2 = 0;
270  int nppar3 = 0;
271  int nppar4 = 0;
272 
273  cminf = -999.;
274  vminf = 0.;
275 
276  for (int j=0; j<nptfit; j++)
277  sigy.push_back(sigma);
278 
279  float a = 0.;
280  float b = 0.;
281  float covss, covii, covsi;
282 
283  fit(xfit,yfit,nptfit,sigy,b,a,chi2fit,covss,covii,covsi);
284 
285  bminf = b;
286  aminf = a;
287  nppar = 2;
288  nppar2 = nppar;
289  chi2fitN2 = chi2fit/(nptfit-2);
290 
291  // cout << "dt0 = 0 chi2fit = " << chi2fit << " slope = "<<b<<endl;
292 
293  if (nptfit >= 3) {
294 
295  fitNpar(3,xfit,yfit,lfit,tfit,sigy,bminf,aminf,cminf,vminf,chi2fit,debug);
296  chi2fit3 = chi2fit;
297 
298  if (cminf!=-999.) {
299  nppar = 3;
300  if (nptfit>3)
301  chi2fitN3 = chi2fit /(nptfit-3);
302  } else {
303  bminf = b;
304  aminf = a;
305  chi2fitN3 = chi2fit /(nptfit-2);
306  }
307 
308  bminf3 = bminf;
309  aminf3 = aminf;
310  cminf3 = cminf;
311  nppar3 = nppar;
312 
313  if (debug) {
314  cout << "dt0= 0 : slope 2 = " << b << " pos in = " << a << " chi2fitN2 = " << chi2fitN2
315  << " nppar = " << nppar2 << " nptfit = " << nptfit << endl;
316  cout << "dt0 = 0 : slope 3 = " << bminf << " pos out = " << aminf << " chi2fitN3 = "
317  << chi2fitN3 << " nppar = " << nppar3 << " T0_ev ns = " << cminf/0.00543 << endl;
318  cout << " vdrift_4parfit "<< vdrift_4parfit<<endl;
319  }
320 
321 
322  if (nptfit>=5) {
323  fitNpar(4,xfit,yfit,lfit,tfit,sigy,bminf,aminf,cminf,vminf,chi2fit,debug);
324 
325  if (vminf!=0) {
326  nppar = 4;
327  if (nptfit<=nppar){
328  chi2fitN4=-1;
329  } else{
330  chi2fitN4= chi2fit / (nptfit-nppar);
331  }
332  } else {
333  if (nptfit <= nppar) chi2fitN4=-1;
334  else chi2fitN4 = chi2fit / (nptfit-nppar);
335  }
336 
337  if (fabs(vminf) >= 0.29) {
338  // for safety and for code construction..dont accept correction on dv/vdrift greater then 0.09
339  vminf = 0.;
340  cminf = cminf3;
341  aminf = aminf3;
342  bminf = bminf3;
343  nppar = 3;
344  chi2fit = chi2fit3;
345  }
346 
347  }
348 
349  if (!vdrift_4parfit){ //if not required explicitly leave the t0 and track step as at step 3
350  // just update vdrift value vmin for storing in the segments for monitoring
351  cminf = cminf3;
352  aminf = aminf3;
353  bminf = bminf3;
354  nppar = 3;
355  chi2fit = chi2fit3;
356  }
357 
358  nppar4 = nppar;
359 
360  } //end nptfit >=3
361 
362  if (debug) {
363  cout << " dt0= 0 : slope 4 = " << bminf << " pos out = " << aminf <<" chi2fitN4 = " << chi2fitN4
364  << " nppar= " << nppar4 << " T0_ev ns= " << cminf/0.00543 <<" delta v = " << vminf <<endl;
365  cout << nptfit << " nptfit " << " end chi2fit = " << chi2fit/ (nptfit-nppar ) << " T0_ev ns= " << cminf/0.00543 << " delta v = " << vminf <<endl;
366  }
367 
368  if ( fabs(vminf) >= 0.09 && debug ) { //checks only vdrift less then 10 % accepted
369 // cout << "vminf gt 0.09 det= " << endl;
370 // cout << "dt0= 0 : slope 4 = "<< bminf << " pos out = " << aminf << " chi2fitN4 = " << chi2fitN4
371 // << " T0_ev ns = " << cminf/0.00543 << " delta v = "<< vminf << endl;
372 // cout << "dt0 = 0 : slope 2 = "<< b << " pos in = " << a <<" chi2fitN2 = " << chi2fitN2
373 // << " nppar = " << nppar-1 << " nptfit = " << nptfit <<endl;
374 // cout << "dt0 = 0 : slope 3 = " << bminf << " pos out = " << aminf << " chi2fitN3 = "
375 // << chi2fitN3 << " T0_ev ns = " << cminf/0.00543 << endl;
376  cout << nptfit <<" nptfit "<< " end chi2fit = " << chi2fit << "T0_ev ns= " << cminf/0.00543 << "delta v = "<< vminf <<endl;
377  }
378 
379  if (nptfit != nppar) chi2fit = chi2fit / (nptfit-nppar);
380 }
DTLinearFit::fit
void fit(const std::vector< float > &x, const std::vector< float > &y, int ndat, const std::vector< float > &sigy, float &slope, float &intercept, double &chi2, float &covss, float &covii, float &covsi) const
Definition: DTLinearFit.cc:27
j
int j
Definition: DBlmapReader.cc:9
DTLinearFit::fitNpar
void fitNpar(const int npar, const std::vector< float > &xfit, const std::vector< float > &yfit, const std::vector< int > &lfit, const std::vector< double > &tfit, const std::vector< float > &sigy, float &aminf, float &bminf, float &cminf, float &vminf, double &chi2fit, const bool debug) const
Definition: DTLinearFit.cc:70
debug
#define debug
Definition: HDRShower.cc:19
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:145
void DTLinearFit::fitNpar ( const int  npar,
const std::vector< float > &  xfit,
const std::vector< float > &  yfit,
const std::vector< int > &  lfit,
const std::vector< double > &  tfit,
const std::vector< float > &  sigy,
float &  aminf,
float &  bminf,
float &  cminf,
float &  vminf,
double &  chi2fit,
const bool  debug = 0 
) const

Definition at line 70 of file DTLinearFit.cc.

References alignmentValidation::c1, j, and margin.

Referenced by DTSegmentUpdator::fit().

81  {
82 
83  int nptfit=xfit.size();
84  if (nptfit<npar) return;
85 
86  double sx = 0.;
87  double sxx = 0.;
88  double sy = 0.;
89  double sxy = 0.;
90  double sl = 0.;
91  double sl2 = 0.;
92  double sly = 0.;
93  double slx = 0.;
94  double st = 0.;
95  double st2 = 0.;
96  double slt = 0.;
97  double sltx = 0.;
98  double slty = 0.;
99 
100  for (int j=0; j<nptfit; j++){
101  sx += xfit[j];
102  sy += yfit[j];
103  sxx+= xfit[j]*xfit[j];
104  sxy+= xfit[j]*yfit[j];
105  sl += lfit[j];
106  sl2+= lfit[j]*lfit[j];
107  sly+= lfit[j]*yfit[j];
108  slx+= lfit[j]*xfit[j];
109  st += tfit[j];
110  st2+= tfit[j]*tfit[j];
111  slt+= lfit[j]*tfit[j];
112  sltx+= lfit[j]*tfit[j]*xfit[j];
113  slty+= lfit[j]*tfit[j]*yfit[j];
114  } //end loop
115 
116  const double d1 = sy;
117  const double d2 = sxy;
118  const double d3 = sly;
119  const double c1 = sl;
120  const double c2 = slx;
121  const double c3 = sl2;
122  const double b1 = sx;
123  const double b2 = sxx;
124  const double b3 = slx;
125  const double a1 = nptfit;
126  const double a2 = sx;
127  const double a3 = sl;
128 
129  const double b4 = b2*a1-b1*a2;
130  const double c4 = c2*a1-c1*a2;
131  const double d4 = d2*a1-d1*a2;
132  const double b5 = a1*b3-a3*b1;
133  const double c5 = a1*c3-a3*c1;
134  const double d5 = a1*d3-d1*a3;
135  const double a6 = slt;
136  const double b6 = sltx;
137  const double c6 = st;
138  const double v6 = st2;
139  const double d6 = slty;
140 
141  chi2fit = -1.;
142 
143  if (npar==3) {
144 
145  if (((c5*b4-c4*b5)*b4*a1)!=0) {
146  cminf = (d5*b4-d4*b5)/(c5*b4-c4*b5);
147  if (fabs(cminf)<0.000001) cminf=0;
148  aminf = d4/b4 -cminf *c4/b4;
149  bminf = (d1/a1 -cminf*c1/a1 -aminf*b1/a1);
150 
151  chi2fit = 0.;
152  for (int j=0; j<nptfit; j++){
153  const double ypred = aminf*xfit[j] + bminf;
154  const double ycorr = yfit[j]-cminf*lfit[j];
155  const double dy = (ycorr - ypred)/sigy[j];
156  chi2fit += dy*dy;
157 
158  if (tfit[j]==0.) continue;
159  double xwire = yfit[j]-tfit[j];
160 
161  // define a small safety margin
162  double margin = 0.1;
163 // cout << " pred: " << ypred << " hit: " << ycorr << "chi2: " << dy*dy << " t0: " << -cminf/0.00543 << endl;
164 
165  // Sanity checks - check that the hit after the fit is still withing the DT cell volume
166  if ((fabs(yfit[j]-xwire)>margin) && (fabs(ycorr-xwire)>margin) && ((yfit[j]-xwire)*(ycorr-xwire)<0)) {
167 // cout << " segmentpred: " << ypred << " hit: " << ycorr << " xwire: " << xwire << " yhit: " << yfit[j] << " t0: " << -cminf/0.00543 << endl;
168 // cout << " XXX hit moved across wire!!!" << endl;
169  chi2fit=-1.;
170  return;
171  }
172 
173  if (fabs(ypred-xwire)>2.1 + margin) {
174 // cout << " segmentpred: " << ypred << " hit: " << ycorr << " xwire: " << xwire << " yhit: " << yfit[j] << " t0: " << -cminf/0.00543 << endl;
175 // cout << " XXX segment outside chamber!!!" << endl;
176  chi2fit=-1.;
177  return;
178  }
179 
180  if (fabs(ycorr-xwire)>2.1 + margin) {
181 // cout << " segmentpred: " << ypred << " hit: " << ycorr << " xwire: " << xwire << " yhit: " << yfit[j] << " t0: " << -cminf/0.00543 << endl;
182 // cout << " XXX hit outside chamber!!!" << endl;
183  chi2fit=-1.;
184  return;
185  }
186  if (fabs(chi2fit<0.0001)) chi2fit=0;
187  } //end loop chi2
188  }
189  } else if (npar==4) {
190  const double det = (a1*a1*(b2*v6 - b6*b6) - a1*(a2*a2*v6 - 2*a2*a6*b6 + a6*a6*b2 + b2*c6*c6 + b3*(b3*v6 - 2*b6*c6))
191  + a2*a2*c6*c6 + 2*a2*(a3*(b3*v6 - b6*c6) - a6*b3*c6) + a3*a3*(b6*b6 - b2*v6)
192  + a6*(2*a3*(b2*c6 - b3*b6) + a6*b3*b3));
193 
194  if (det != 0) {
195  // computation of a, b, c and v
196  bminf = (a1*(a2*(b6*d6 - v6*d2) + a6*(b6*d2 - b2*d6) + d1*(b2*v6 - b6*b6)) - a2*(b3*(c6*d6 - v6*d3)
197  + c6*(b6*d3 - c6*d2)) + a3*(b2*(c6*d6 - v6*d3) + b3*(v6*d2 - b6*d6) + b6*(b6*d3 - c6*d2))
198  + a6*(b2*c6*d3 + b3*(b3*d6 - b6*d3 - c6*d2)) - d1*(b2*c6*c6 + b3*(b3*v6 - 2*b6*c6)))/det;
199  aminf = - (a1*a1*(b6*d6 - v6*d2) - a1*(a2*(a6*d6 - v6*d1) - a6*a6*d2 + a6*b6*d1 + b3*(c6*d6 - v6*d3)
200  + c6*(b6*d3 - c6*d2)) + a2*(a3*(c6*d6 - v6*d3) + c6*(a6*d3 - c6*d1)) + a3*a3*(v6*d2 - b6*d6)
201  + a3*(a6*(b3*d6 + b6*d3 - 2*c6*d2) - d1*(b3*v6 - b6*c6)) - a6*b3*(a6*d3 - c6*d1))/det;
202  cminf = -(a1*(b2*(c6*d6 - v6*d3) + b3*(v6*d2 - b6*d6) + b6*(b6*d3 - c6*d2)) + a2*a2*(v6*d3 - c6*d6)
203  + a2*(a3*(b6*d6 - v6*d2) + a6*(b3*d6 - 2*b6*d3 + c6*d2) - d1*(b3*v6 - b6*c6))
204  + a3*(d1*(b2*v6 - b6*b6) - a6*(b2*d6 - b6*d2)) + a6*(a6*(b2*d3 - b3*d2) - d1*(b2*c6 - b3*b6)))/det;
205  vminf = - (a1*a1*(b2*d6 - b6*d2) - a1*(a2*a2*d6 - a2*(a6*d2 + b6*d1) + a6*b2*d1 + b2*c6*d3
206  + b3*(b3*d6 - b6*d3 - c6*d2)) + a2*a2*c6*d3 + a2*(a3*(2*b3*d6 - b6*d3 - c6*d2) - b3*(a6*d3 + c6*d1))
207  + a3*a3*(b6*d2 - b2*d6) + a3*(a6*(b2*d3 - b3*d2) + d1*(b2*c6 - b3*b6)) + a6*b3*b3*d1)/det;
208 
209  chi2fit = 0.;
210  for (int j=0; j<nptfit; j++) {
211  const double ypred = aminf*xfit[j] + bminf;
212  const double dy = (yfit[j]+vminf*lfit[j]*tfit[j]-cminf*lfit[j] -ypred)/sigy[j];
213  chi2fit += dy*dy;
214  } //end loop chi2
215  }
216  }
217 }
margin
Double_t margin
Definition: trackSplitPlot.h:45
j
int j
Definition: DBlmapReader.cc:9
Generated for CMSSW Reference Manual by doxygen 1.8.5