]>
Commit | Line | Data |
---|---|---|
a2bea7ce | 1 | /************************************************************************** |
2 | * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | ||
b0f5e3fc | 18 | ////////////////////////////////////////////////////////////////////////// |
a2bea7ce | 19 | // Alice ITS class to help keep statistical information. Can also be // |
20 | // used to fit data to lines and other 2 dimentional sytistical // | |
21 | // operations. // | |
b0f5e3fc | 22 | // // |
23 | // version: 0.0.0 Draft. // | |
24 | // Date: April 18 1999 // | |
25 | // By: Bjorn S. Nilsen // | |
a2bea7ce | 26 | // Updated: 1.0.0, Date: September 6 2007, By: Bjorn S. Nilsen // |
b0f5e3fc | 27 | // // |
28 | ////////////////////////////////////////////////////////////////////////// | |
a2bea7ce | 29 | #include <stdio.h> // ios::fmtflags fmt used in PrintAscii |
30 | #include "Riostream.h" // IO functions. | |
31 | #include "TMath.h" // TMath::Sqrt() function used. | |
32 | #include "AliITSstatistics2.h" // Also defined TObject {base class} | |
b0f5e3fc | 33 | |
34 | ClassImp(AliITSstatistics2) | |
35 | ||
36 | // | |
492b8715 | 37 | AliITSstatistics2::AliITSstatistics2() : |
38 | TObject(), // Base Class | |
39 | fN(-1), // number of enetries -1 => Uninitilized | |
40 | fOrder(0), // maximum moment of distributions (^n) | |
a2bea7ce | 41 | fX(0), //[fOrder] array of sums of x^n |
42 | fYx(0), //[fOrder] array of sums of (xy)^n | |
43 | fY(0), //[fOrder] array of sums of y^n | |
44 | fW(0) //[fOrder] array of sums of w^n (weights) | |
45 | //,fDig(5) // The number of significant digits to keep | |
46 | //,fOver(0) //! In case of numerical precistion problems | |
47 | { | |
492b8715 | 48 | // default constructor |
49 | // Inputs: | |
50 | // none. | |
51 | // Outputs: | |
52 | // none. | |
53 | // Return: | |
54 | // A default constructed AliITSstatistics class | |
55 | ||
b0f5e3fc | 56 | return; |
57 | } | |
492b8715 | 58 | //______________________________________________________________________ |
59 | AliITSstatistics2::AliITSstatistics2(Int_t order) : | |
60 | TObject(), // Base Class | |
61 | fN(0), // number of enetries -1 => Uninitilized | |
62 | fOrder(order), // maximum moment of distributions (^n) | |
63 | fX(new Double_t[order]), //[fOrder] array of sums of x^n | |
64 | fYx(new Double_t[order]), //[fOrder] array of sums of (xy)^n | |
65 | fY(new Double_t[order]), //[fOrder] array of sums of y^n | |
a2bea7ce | 66 | fW(new Double_t[order]) //[fOrder] array of sums of w^n (weights) |
67 | //,fDig(5) // The number of significant digits to keep | |
68 | //,fOver(0) //! In case of numeerical precistion problems | |
69 | { // constructor to maximum moment/order order | |
492b8715 | 70 | // Inputs: |
71 | // Int_t order The maximum moment of distributions {for example x^n} | |
72 | Int_t i; | |
b0f5e3fc | 73 | |
492b8715 | 74 | for(i=0;i<order;i++) { |
75 | fX[i] = 0.0; | |
76 | fY[i] = 0.0; | |
77 | fYx[i] = 0.0; | |
78 | fW[i] = 0.0; | |
79 | } // end for i | |
b0f5e3fc | 80 | fN = 0; |
81 | return; | |
82 | } | |
492b8715 | 83 | //______________________________________________________________________ |
b0f5e3fc | 84 | AliITSstatistics2::~AliITSstatistics2(){ |
492b8715 | 85 | // destructor |
86 | // Inputs: | |
87 | // none. | |
88 | // Outputs: | |
89 | // none. | |
90 | // Return: | |
91 | // none. | |
92 | ||
07025d6d | 93 | if(fX!=0) delete[] fX; |
94 | if(fY!=0) delete[] fY; | |
95 | if(fYx!=0) delete[] fYx; | |
96 | if(fW!=0) delete[] fW; | |
a2bea7ce | 97 | // if(fOver!=0) delete fOver; fOver=0; |
98 | // fDig=0; | |
07025d6d | 99 | fX = 0; |
100 | fY = 0; | |
101 | fYx = 0; | |
102 | fW = 0; | |
b0f5e3fc | 103 | fN = 0; |
104 | fOrder = 0; | |
105 | } | |
b0f5e3fc | 106 | //_______________________________________________________________ |
107 | AliITSstatistics2& AliITSstatistics2::operator=(AliITSstatistics2 &source){ | |
492b8715 | 108 | // operator = |
109 | // Inputs: | |
110 | // AliITSstaticstics2 &source The source of this copy. | |
111 | // Outputs: | |
112 | // none. | |
113 | // Return: | |
114 | // A copy of the source class | |
115 | ||
116 | if(this==&source) return *this; | |
117 | TObject::operator=(source); | |
118 | Reset(source.GetOrder()); | |
119 | fN = source.GetN(); | |
120 | fOrder=source.GetOrder(); | |
121 | for(Int_t i=0;i<source.fOrder;i++){ | |
122 | this->fX[i] = source.fX[i]; | |
123 | this->fYx[i] = source.fYx[i]; | |
124 | this->fY[i] = source.fY[i]; | |
125 | this->fW[i] = source.fW[i]; | |
126 | } // end for i | |
a2bea7ce | 127 | // this->fDig = source.fDig; |
128 | // if(fOver!=0) this->fOver = new AliITSstatistics2(*(source.fOver)); | |
129 | // else fOver=0; | |
492b8715 | 130 | return *this; |
b0f5e3fc | 131 | } |
132 | //_______________________________________________________________ | |
492b8715 | 133 | AliITSstatistics2::AliITSstatistics2(AliITSstatistics2 &source): |
134 | TObject(source), // Base Class | |
135 | fN(source.GetN()), // number of enetries -1 => Uninitilized | |
136 | fOrder(source.GetOrder()),// maximum moment of distributions (^n) | |
137 | fX(new Double_t[source.GetOrder()]),//[fOrder] array of sums of x^n | |
138 | fYx(new Double_t[source.GetOrder()]),//[fOrder] array of sums of (xy)^n | |
139 | fY(new Double_t[source.GetOrder()]),//[fOrder] array of sums of y^n | |
a2bea7ce | 140 | fW(new Double_t[source.GetOrder()]) //[fOrder] array of sums of w^n (weights) |
141 | //,fDig(source.fDig) // The number of significant digits to keep | |
142 | //,fOver(0) //! In case of numerical precistion problems | |
143 | { | |
492b8715 | 144 | // Copy constructor |
145 | // Inputs: | |
146 | // AliITSstatistics2 & source the source of this copy | |
147 | // Outputs: | |
148 | // none. | |
149 | // Return: | |
150 | // A copy of the source. | |
151 | ||
492b8715 | 152 | for(Int_t i=0;i<source.fOrder;i++){ |
153 | this->fX[i] = source.fX[i]; | |
154 | this->fYx[i] = source.fYx[i]; | |
155 | this->fY[i] = source.fY[i]; | |
156 | this->fW[i] = source.fW[i]; | |
157 | } // end for i | |
a2bea7ce | 158 | //if(fOver!=0) this->fOver = new AliITSstatistics2(*(source.fOver)); |
159 | return; | |
b0f5e3fc | 160 | } |
492b8715 | 161 | //______________________________________________________________________ |
162 | void AliITSstatistics2::Reset(Int_t order){ | |
163 | // Reset/zero all statistics variables statistics | |
164 | // Inputs: | |
165 | // none. | |
166 | // Outputs: | |
167 | // none. | |
168 | // Return: | |
169 | // none. | |
170 | Int_t i; | |
b0f5e3fc | 171 | |
492b8715 | 172 | for(i=0;i<fOrder;i++) { |
173 | fX[i] = 0.0; | |
174 | fY[i] = 0.0; | |
175 | fYx[i] = 0.0; | |
176 | fW[i] = 0.0; | |
177 | } // end for i | |
b0f5e3fc | 178 | fN = 0; |
492b8715 | 179 | if(order<0) return; // just zero |
180 | if(fX!=0) delete[] fX; | |
181 | if(fY!=0) delete[] fY; | |
182 | if(fYx!=0) delete[] fYx; | |
183 | if(fW!=0) delete[] fW; | |
184 | fX = 0; | |
185 | fY = 0; | |
186 | fYx = 0; | |
187 | fW = 0; | |
188 | fN = 0; | |
189 | fOrder = 0; | |
190 | if(order==0) return; | |
191 | fOrder = order; | |
192 | fX = new Double_t[fOrder]; | |
193 | fY = new Double_t[fOrder]; | |
194 | fYx = new Double_t[fOrder]; | |
195 | fW = new Double_t[fOrder]; | |
a2bea7ce | 196 | //if(fOver!=0) delete fOver; fOver = 0; |
b0f5e3fc | 197 | return; |
198 | } | |
a2bea7ce | 199 | //---------------------------------------------------------------------- |
200 | /* | |
201 | void SetSignificantDigits(Int_t d){ | |
202 | // Sets the number of significant digits. If adding a value to | |
203 | // one of this class' arrays looses significance at the fDig | |
204 | // level, a new instance of this class is created to keep | |
205 | // signigicance at or better than fDig level. if fDig<0, then | |
206 | // this feature it disabled and significance can be lost. | |
207 | // Inputs: | |
208 | // Int_t d The new significance level | |
209 | // Outputs: | |
210 | // none. | |
211 | // Return: | |
212 | // none. | |
213 | ||
214 | fDig = d; | |
215 | } | |
216 | */ | |
492b8715 | 217 | //______________________________________________________________________ |
b0f5e3fc | 218 | void AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0){ |
492b8715 | 219 | // add next x,y pair to statistics |
220 | // Inputs: | |
221 | // Double_t y y value of pair | |
222 | // Double_t x x value of pair | |
223 | // Double_t w weight of pair | |
224 | // Outputs: | |
225 | // none. | |
226 | // Return: | |
227 | // none. | |
b0f5e3fc | 228 | Double_t xs=1.0,ys=1.0,yxs=1.0,ws=1.0; |
229 | Int_t i; | |
e8189707 | 230 | const Double_t kBig=1.0e+38; |
231 | ||
232 | if(y>kBig || x>kBig || w>kBig) return; | |
a2bea7ce | 233 | /* If problem with precision, then creat/fill fOver |
234 | as a partical sum to be added to "this" later. | |
235 | if(????fDig){ | |
236 | if(fOver==0){ | |
237 | fOver = new AliITSstatistics2(fOrder); | |
238 | } // end if fOver==0 | |
239 | fOver->AddValue(y,x,w); | |
240 | return; | |
241 | } // end if(???) | |
242 | */ | |
b0f5e3fc | 243 | fN++; |
a2bea7ce | 244 | for(i=0;i<GetOrder();i++){ |
b0f5e3fc | 245 | xs *= x; |
246 | ys *= y; | |
247 | yxs *= x*y; | |
248 | ws *= w; | |
07025d6d | 249 | fX[i] += xs*w; |
250 | fY[i] += ys*w; | |
251 | fYx[i] += yxs*w; | |
252 | fW[i] += ws; | |
b0f5e3fc | 253 | } // end for i |
254 | } | |
492b8715 | 255 | //______________________________________________________________________ |
256 | Double_t AliITSstatistics2::GetXNth(Int_t order)const{ | |
257 | // This give the unbiased estimator for the RMS. | |
258 | // Inputs: | |
259 | // Int_t order the order of x^n value to be returned | |
260 | // Output: | |
261 | // none. | |
262 | // Return: | |
263 | // The value sum{x^n}. | |
b0f5e3fc | 264 | Double_t s; |
265 | ||
a2bea7ce | 266 | if(GetWN(1)!=0.0 && order<=GetOrder()) s = GetXN(order)/GetWN(1); |
b0f5e3fc | 267 | else { |
268 | s = 0.0; | |
492b8715 | 269 | Error("GetXNth","error fOrder=%d fN=%d fW[0]=%f\n", |
a2bea7ce | 270 | GetOrder(),GetN(),GetWN(1)); |
b0f5e3fc | 271 | } // end else |
272 | return s; | |
273 | } | |
492b8715 | 274 | //______________________________________________________________________ |
275 | Double_t AliITSstatistics2::GetYNth(Int_t order)const{ | |
276 | // This give the unbiased estimator for the RMS. | |
277 | // Inputs: | |
278 | // Int_t order the order of y^n value to be returned | |
279 | // Outputs: | |
280 | // none. | |
281 | // Return: | |
282 | // The value sum{y^n} | |
b0f5e3fc | 283 | Double_t s; |
284 | ||
a2bea7ce | 285 | if(GetWN(1)!=0.0&&order<=GetOrder()) s = GetYN(order)/GetWN(1); |
b0f5e3fc | 286 | else { |
287 | s = 0.0; | |
492b8715 | 288 | Error("GetYNth","fOrder=%d fN=%d fW[0]=%f\n", |
a2bea7ce | 289 | GetOrder(),GetN(),GetWN(1)); |
b0f5e3fc | 290 | } // end else |
291 | return s; | |
292 | } | |
492b8715 | 293 | //______________________________________________________________________ |
294 | Double_t AliITSstatistics2::GetYXNth(Int_t order)const{ | |
295 | // This give the unbiased estimator for the RMS. | |
296 | // Inputs: | |
297 | // Int_t order the order of (xy)^n value to be returned | |
298 | // Outputs: | |
299 | // none. | |
300 | // Return: | |
301 | // The value sum{(xy)^n} | |
b0f5e3fc | 302 | Double_t s; |
303 | ||
a2bea7ce | 304 | if(GetWN(1)!=0.0&&order<=GetOrder()) s = GetYXN(order)/GetWN(1); |
b0f5e3fc | 305 | else { |
306 | s = 0.0; | |
492b8715 | 307 | Error("GetYXNth","fOrder=%d fN=%d fW[0]=%f\n", |
a2bea7ce | 308 | GetOrder(),GetN(),GetWN(1)); |
b0f5e3fc | 309 | } // end else |
310 | return s; | |
311 | } | |
492b8715 | 312 | //______________________________________________________________________ |
313 | Double_t AliITSstatistics2::GetRMSX()const{ | |
314 | // This give the unbiased estimator for the RMS. | |
315 | // Inputs: | |
316 | // none. | |
317 | // Outputs: | |
318 | // none. | |
319 | // Return: | |
320 | // The rms value | |
b0f5e3fc | 321 | Double_t x,x2,w,ww,s; |
322 | ||
323 | x = GetMeanX(); // first order | |
324 | x2 = GetXNth(2); // second order | |
a2bea7ce | 325 | w = GetWN(1); // first order |
326 | ww = GetWN(2); // second order | |
b0f5e3fc | 327 | |
328 | if(w*w==ww) return (-1.0); | |
329 | s = (x2-x*x)*w*w/(w*w-ww); | |
330 | return TMath::Sqrt(s); | |
331 | } | |
492b8715 | 332 | //______________________________________________________________________ |
333 | Double_t AliITSstatistics2::GetRMSY()const{ | |
334 | // This give the unbiased estimator for the RMS. | |
335 | // Inputs: | |
336 | // none. | |
337 | // Outputs: | |
338 | // none. | |
339 | // Return: | |
340 | // The rms value | |
b0f5e3fc | 341 | Double_t x,x2,w,ww,s; |
342 | ||
343 | x = GetMeanY(); // first order | |
344 | x2 = GetYNth(2); // second order | |
a2bea7ce | 345 | w = GetWN(1); // first order |
346 | ww = GetWN(2); // second order | |
b0f5e3fc | 347 | |
348 | if(w*w==ww) return (-1.0); | |
349 | s = (x2-x*x)*w*w/(w*w-ww); | |
350 | return TMath::Sqrt(s); | |
351 | } | |
492b8715 | 352 | //______________________________________________________________________ |
353 | Double_t AliITSstatistics2::GetRMSYX()const{ | |
354 | // This give the unbiased estimator for the RMS. | |
355 | // Inputs: | |
356 | // none. | |
357 | // Outputs: | |
358 | // none. | |
359 | // Return: | |
360 | // The rms value | |
b0f5e3fc | 361 | Double_t x,x2,w,ww,s; |
362 | ||
363 | x = GetMeanYX(); // first order | |
364 | x2 = GetYXNth(2); // second order | |
a2bea7ce | 365 | w = GetWN(1); // first order |
366 | ww = GetWN(2); // second order | |
b0f5e3fc | 367 | |
368 | if(w*w==ww) return (-1.0); | |
369 | s = (x2-x*x)*w*w/(w*w-ww); | |
370 | return TMath::Sqrt(s); | |
371 | } | |
492b8715 | 372 | //______________________________________________________________________ |
373 | Double_t AliITSstatistics2::GetErrorMeanY()const{ | |
374 | //This is the error in the mean or the square root of the | |
375 | // variance of the mean. | |
376 | // Inputs: | |
377 | // none. | |
378 | // Outputs: | |
379 | // none. | |
380 | // Return: | |
381 | // The error on the mean | |
b0f5e3fc | 382 | Double_t rms,w,ww,s; |
383 | ||
384 | rms = GetRMSY(); | |
a2bea7ce | 385 | w = GetWN(1); |
386 | ww = GetWN(2); | |
b0f5e3fc | 387 | s = rms*rms*ww/(w*w); |
388 | return TMath::Sqrt(s); | |
389 | } | |
492b8715 | 390 | //______________________________________________________________________ |
391 | Double_t AliITSstatistics2::GetErrorMeanX()const{ | |
392 | //This is the error in the mean or the square root of the | |
393 | // variance of the mean. | |
394 | // Inputs: | |
395 | // none. | |
396 | // Outputs: | |
397 | // none. | |
398 | // Return: | |
399 | // The error on the mean | |
b0f5e3fc | 400 | Double_t rms,w,ww,s; |
401 | ||
402 | rms = GetRMSX(); | |
a2bea7ce | 403 | w = GetWN(1); |
404 | ww = GetWN(2); | |
b0f5e3fc | 405 | s = rms*rms*ww/(w*w); |
406 | return TMath::Sqrt(s); | |
407 | } | |
492b8715 | 408 | //______________________________________________________________________ |
409 | Double_t AliITSstatistics2::GetErrorMeanYX()const{ | |
410 | //This is the error in the mean or the square root of the | |
411 | // variance of the mean. | |
412 | // Inputs: | |
413 | // none. | |
414 | // Outputs: | |
415 | // none. | |
416 | // Return: | |
417 | // The error on the mean | |
b0f5e3fc | 418 | Double_t rms,w,ww,s; |
419 | ||
420 | rms = GetRMSYX(); | |
a2bea7ce | 421 | w = GetWN(1); |
422 | ww = GetWN(2); | |
b0f5e3fc | 423 | s = rms*rms*ww/(w*w); |
424 | return TMath::Sqrt(s); | |
425 | } | |
492b8715 | 426 | //______________________________________________________________________ |
427 | Double_t AliITSstatistics2::GetErrorRMSY()const{ | |
428 | // This is the error in the mean or the square root of the variance | |
429 | // of the mean. at this moment this routine is only defined for | |
430 | // weights=1. | |
431 | // Inputs: | |
432 | // none. | |
433 | // Outputs: | |
434 | // none. | |
435 | // Return: | |
436 | // The error on the rms | |
b0f5e3fc | 437 | Double_t x,x2,x3,x4,w,ww,m2,m4,n,s; |
438 | ||
a2bea7ce | 439 | if(GetWN(1)!=(Double_t)GetN()||GetN()<4) return (-1.); |
b0f5e3fc | 440 | x = GetMeanY(); // first order |
441 | x2 = GetYNth(2); // second order | |
a2bea7ce | 442 | w = GetWN(1); // first order |
443 | ww = GetWN(2); // second order | |
b0f5e3fc | 444 | if(w*w==ww) return (-1.0); |
445 | s = (x2-x*x)*w*w/(w*w-ww); | |
446 | ||
447 | m2 = s; | |
448 | n = (Double_t) GetN(); | |
449 | x3 = GetYNth(3); | |
450 | x4 = GetYNth(4); | |
492b8715 | 451 | // This equation assumes that all of the weights are equal to 1. |
b0f5e3fc | 452 | m4 = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x); |
453 | s = (m4-(n-3.)*m2*m2/(n-1.))/n; | |
454 | return TMath::Sqrt(s); | |
455 | } | |
492b8715 | 456 | //______________________________________________________________________ |
457 | Double_t AliITSstatistics2::GetErrorRMSX()const{ | |
458 | // This is the error in the mean or the square root of the variance | |
459 | // of the mean. at this moment this routine is only defined for | |
460 | // weights=1. | |
461 | // Inputs: | |
462 | // none. | |
463 | // Outputs: | |
464 | // none. | |
465 | // Return: | |
466 | // The error on the rms | |
b0f5e3fc | 467 | Double_t x,x2,x3,x4,w,ww,m2,m4,n,s; |
468 | ||
a2bea7ce | 469 | if(GetWN(1)!=(Double_t)GetN()||GetN()<4) return (-1.); |
b0f5e3fc | 470 | x = GetMeanX(); // first order |
471 | x2 = GetXNth(2); // second order | |
a2bea7ce | 472 | w = GetWN(1); // first order |
473 | ww = GetWN(2); // second order | |
b0f5e3fc | 474 | if(w*w==ww) return (-1.0); |
475 | s = (x2-x*x)*w*w/(w*w-ww); | |
476 | ||
477 | m2 = s; | |
478 | n = (Double_t) GetN(); | |
479 | x3 = GetXNth(3); | |
480 | x4 = GetXNth(4); | |
481 | // This equation assumes that all of the weights are equal to 1. | |
482 | m4 = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x); | |
483 | s = (m4-(n-3.)*m2*m2/(n-1.))/n; | |
484 | return TMath::Sqrt(s); | |
485 | } | |
492b8715 | 486 | //______________________________________________________________________ |
487 | Double_t AliITSstatistics2::GetErrorRMSYX()const{ | |
488 | // This is the error in the mean or the square root of the variance | |
489 | // of the mean. at this moment this routine is only defined for | |
490 | // weights=1. | |
491 | // Inputs: | |
492 | // none. | |
493 | // Outputs: | |
494 | // none. | |
495 | // Return: | |
496 | // The error on the rms | |
b0f5e3fc | 497 | Double_t x,x2,x3,x4,w,ww,m2,m4,n,s; |
498 | ||
a2bea7ce | 499 | if(GetWN(1)!=(Double_t)GetN()||GetN()<4) return (-1.); |
b0f5e3fc | 500 | x = GetMeanYX(); // first order |
501 | x2 = GetYXNth(2); // second order | |
a2bea7ce | 502 | w = GetWN(1); // first order |
503 | ww = GetWN(2); // second order | |
b0f5e3fc | 504 | if(w*w==ww) return (-1.0); |
505 | s = (x2-x*x)*w*w/(w*w-ww); | |
506 | ||
507 | m2 = s; | |
508 | n = (Double_t) GetN(); | |
509 | x3 = GetYXNth(3); | |
510 | x4 = GetYXNth(4); | |
492b8715 | 511 | // This equation assumes that all of the weights are equal to 1. |
b0f5e3fc | 512 | m4 = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x); |
513 | s = (m4-(n-3.)*m2*m2/(n-1.))/n; | |
514 | return TMath::Sqrt(s); | |
515 | } | |
516 | //_______________________________________________________________________ | |
a2bea7ce | 517 | Double_t AliITSstatistics2::FitToLine(Double_t &a,Double_t &ea, |
518 | Double_t &b,Double_t &eb)const{ | |
519 | // fit to y = ax+b returns Chi squared or -1.0 if an error. | |
520 | // The fitting is done by analitically minimizing | |
521 | /* | |
522 | Begin_Latex | |
523 | \begin{equation*} | |
524 | \Chi^{2}=\sum_{i} (y_{i}-a x_{i} -b)^{2} w_{i} | |
525 | \end{equation*} | |
526 | Where if the weight used in | |
527 | AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0) | |
528 | is of the form | |
529 | \begin{equation*} | |
530 | w_{i}=\frac{1}{\delta y^{2}}. | |
531 | \end{equation*} | |
532 | Then we get the typicall chi square minimization. | |
533 | End_Latex | |
534 | */ | |
492b8715 | 535 | // Inputs: |
536 | // none. | |
537 | // Outputs: | |
538 | // Double_t a The slope parameter | |
a2bea7ce | 539 | // Double_t ea Error on fitted slope parameter |
492b8715 | 540 | // Double_t b The intercept paramter |
a2bea7ce | 541 | // Double_t eb Error on fitted intercept parameter |
492b8715 | 542 | // Return: |
543 | // The Chi^2 of the fit | |
b0f5e3fc | 544 | Double_t c,d,e,f,g,h; |
545 | ||
a2bea7ce | 546 | a = ea = b = eb = 0.0; |
547 | if(GetOrder()<2 || GetN()<3){ | |
548 | Error("FitToLine","Order=%d<2 or N=%d<3",GetOrder(),GetN()); | |
492b8715 | 549 | return -1.0; |
550 | } // end if | |
b0f5e3fc | 551 | c = GetWN(1); |
552 | d = GetYN(1); | |
553 | e = GetXN(1); | |
554 | f = GetYXN(1); | |
555 | g = GetXN(2); | |
556 | h = c*g-e*e; | |
492b8715 | 557 | b = d*g-f*e; |
558 | a = c*f-d*e; | |
a2bea7ce | 559 | if(h==0.0){ |
492b8715 | 560 | Error("FitToLine","vertical line: fOrder=%d fN=%d " |
a2bea7ce | 561 | "GetWN(1)=%g X GetXN(2)=%g - GetXN(1)=%g^2 = 0", |
562 | GetOrder(),GetN(),c,g,e); | |
b0f5e3fc | 563 | return -1.0; |
564 | } // end if h | |
a2bea7ce | 565 | a = a/h; |
566 | b = b/h; | |
567 | // Now for the errors. | |
568 | ea = c*c*g+(a*a-1.0)*c*e*e; | |
569 | ea = ea/(h*h); | |
570 | if(ea<0.0){ | |
571 | Error("FitToLine","ea=%g is less than zero",ea); | |
572 | return -2.0; | |
573 | } // end if ea<0 | |
574 | ea = TMath::Sqrt(ea); | |
575 | eb = c*g*g-2.0*d*e*g-2.0*(1.0-b)*c*e*e*g+2.0*(1.0-b)*d*e*e*e+ | |
576 | GetYN(2)*e*e+(1.0-b)*(1.0-b)*c*e*e*e*e; | |
577 | eb = eb/(h*h); | |
578 | if(eb<0.0){ | |
579 | Error("FitToLine","eb=%g is less than zero",eb); | |
580 | return -2.0; | |
581 | } // end if ea<0 | |
582 | eb = TMath::Sqrt(eb); | |
583 | c = GetChiSquared(a,b); | |
584 | if(c<=0.0){ // must be a numerical precision problem. | |
585 | } // end if | |
586 | return c; | |
492b8715 | 587 | } |
588 | //_______________________________________________________________________ | |
589 | Double_t AliITSstatistics2::GetChiSquared(Double_t a,Double_t b)const{ | |
a2bea7ce | 590 | // Returns Chi^2 value of data to line y=ax+b with given a,b. |
591 | /* | |
592 | Begin_Latex | |
593 | Note: The Chi^2 value is computed from the expression | |
594 | \begin{equation*} | |
595 | \chi^{2}=\sum_{i}{w_{i}y_{i}^{2}} + b^{2}\sum_{i}{w_{i}} | |
596 | -2b\sum_{i}{w_{i}y_{i}}-2a\sum_{i}{w_{i}y_{i}x_{i}} | |
597 | +2ab\sum_{i}{w_{i}x_{i}} | |
598 | +a^{2}\sum_{i}w_{i}x_{i}^{2} | |
599 | \end{equation*} | |
600 | and not form the expression | |
601 | \begin{equation*} | |
602 | \chi^{2}= \sum_{i}{(y_{i}-ax_{i}-b)^{2}w_{i}. | |
603 | \end{equation*} | |
604 | Consiquently, there are occations when numerically these | |
605 | two expressions will not agree. In fact the form code here | |
606 | can give negitive values. This happens when the numerical | |
607 | significance is larger than the $\chi^{2}$ value. This should | |
608 | not be confused the the error values which can be returned. | |
609 | At present there is no check on the numberical significance | |
610 | of any results. | |
611 | End_Latex | |
612 | */ | |
492b8715 | 613 | // Inputs: |
614 | // Double_t a The slope parameter | |
615 | // Double_t b The intercept paramter | |
616 | // Outputs:: | |
617 | // none. | |
618 | // Return: | |
619 | // The Chi^2 of the fit | |
620 | Double_t c2; | |
621 | ||
622 | c2 = GetYN(2)+b*b*GetWN(1)+ | |
623 | a*a*GetXN(2)-2.0*b*GetYN(1)-2.0*a*GetYXN(1)+2.0*b*a*GetXN(1); | |
a2bea7ce | 624 | c2 /= (Double_t)GetN() - 2.0; |
492b8715 | 625 | return c2; |
626 | } | |
627 | //______________________________________________________________________ | |
628 | void AliITSstatistics2::PrintAscii(ostream *os)const{ | |
629 | // Print out class data values in Ascii Form to output stream | |
630 | // Inputs: | |
631 | // ostream *os Output stream where Ascii data is to be writen | |
632 | // Outputs: | |
633 | // none. | |
634 | // Return: | |
635 | // none. | |
636 | Int_t i; | |
637 | #if defined __GNUC__ | |
638 | #if __GNUC__ > 2 | |
a2bea7ce | 639 | ios::fmtflags fmt; // Standard IO format object, required for output. |
492b8715 | 640 | #else |
641 | Int_t fmt; | |
642 | #endif | |
643 | #else | |
644 | #if defined __ICC || defined __ECC || defined __xlC__ | |
645 | ios::fmtflags fmt; | |
646 | #else | |
647 | Int_t fmt; | |
648 | #endif | |
649 | #endif | |
650 | ||
a2bea7ce | 651 | *os << fN <<" "<< GetOrder(); |
492b8715 | 652 | fmt = os->setf(ios::scientific); // set scientific floating point output |
a2bea7ce | 653 | for(i=0;i<GetOrder();i++) *os <<" "<< GetXN(i+1); |
654 | for(i=0;i<GetOrder();i++) *os <<" "<< GetYXN(i+1); | |
655 | for(i=0;i<GetOrder();i++) *os <<" "<< GetYN(i+1); | |
656 | for(i=0;i<GetOrder();i++) *os <<" "<< GetWN(i+1); | |
657 | //if(fOver!=0) { *os << " " << fDig; | |
658 | //*os << " " << fOver; | |
659 | // } else *os << " " << fDig; | |
492b8715 | 660 | os->flags(fmt); // reset back to old Formating. |
661 | return; | |
662 | } | |
663 | //______________________________________________________________________ | |
664 | void AliITSstatistics2::ReadAscii(istream *is){ | |
665 | // Read in class data values in Ascii Form to output stream | |
666 | // Inputs: | |
667 | // istream *is Input stream where Ascii data is to be read in from | |
668 | // Outputs: | |
669 | // none. | |
670 | // Return: | |
671 | // none. | |
672 | Int_t i; | |
673 | ||
674 | *is >> i >> fOrder; | |
675 | Reset(fOrder); | |
676 | fN = i; | |
677 | for(i=0;i<fOrder;i++) *is >> fX[i]; | |
678 | for(i=0;i<fOrder;i++) *is >> fYx[i]; | |
679 | for(i=0;i<fOrder;i++) *is >> fY[i]; | |
680 | for(i=0;i<fOrder;i++) *is >> fW[i]; | |
a2bea7ce | 681 | //*is >> fDig; |
682 | // if(fDig>0) *is >> fOver; | |
683 | // else fDig *= -1; | |
492b8715 | 684 | } |
685 | //______________________________________________________________________ | |
686 | ostream &operator<<(ostream &os,const AliITSstatistics2 &s){ | |
687 | // Standard output streaming function | |
688 | // Inputs: | |
689 | // ostream &os output steam | |
690 | // AliITSstatistics2 &s class to be streamed. | |
691 | // Output: | |
692 | // none. | |
693 | // Return: | |
694 | // ostream &os The stream pointer | |
695 | ||
696 | s.PrintAscii(&os); | |
697 | return os; | |
698 | } | |
699 | //______________________________________________________________________ | |
700 | istream &operator>>(istream &is,AliITSstatistics2 &s){ | |
701 | // Standard inputput streaming function | |
702 | // Inputs: | |
703 | // istream &is input steam | |
704 | // AliITSstatistics2 &s class to be streamed. | |
705 | // Output: | |
706 | // none. | |
707 | // Return: | |
708 | // ostream &os The stream pointer | |
709 | ||
710 | s.ReadAscii(&is); | |
711 | return is; | |
b0f5e3fc | 712 | } |
a2bea7ce | 713 |