]>
Commit | Line | Data |
---|---|---|
b0f5e3fc | 1 | ////////////////////////////////////////////////////////////////////////// |
2 | // Alice ITS class to help keep statistical information // | |
3 | // // | |
4 | // version: 0.0.0 Draft. // | |
5 | // Date: April 18 1999 // | |
6 | // By: Bjorn S. Nilsen // | |
7 | // // | |
8 | ////////////////////////////////////////////////////////////////////////// | |
9 | #include <stdio.h> | |
10 | #include <math.h> | |
11 | #include "TMath.h" | |
12 | #include "AliITSstatistics2.h" | |
13 | ||
14 | ClassImp(AliITSstatistics2) | |
15 | ||
16 | // | |
17 | AliITSstatistics2::AliITSstatistics2() : TObject(){ | |
18 | // | |
19 | // default constructor | |
20 | // | |
07025d6d | 21 | fX = 0; |
22 | fY = 0; | |
23 | fYx = 0; | |
24 | fW = 0; | |
b0f5e3fc | 25 | fN = 0; |
26 | fOrder = 0; | |
27 | return; | |
28 | } | |
29 | ||
30 | ||
31 | AliITSstatistics2::AliITSstatistics2(Int_t order) : TObject(){ | |
32 | // | |
33 | // constructor to maximum moment/order order | |
34 | // | |
b0f5e3fc | 35 | fOrder = order; |
07025d6d | 36 | fX = new Double_t[order]; |
37 | fY = new Double_t[order]; | |
38 | fYx = new Double_t[order]; | |
39 | fW = new Double_t[order]; | |
40 | for(Int_t i=0;i<order;i++) {fX[i] = 0.0;fY[i] = 0.0; | |
41 | fYx[i] = 0.0; fW[i] = 0.0;} | |
b0f5e3fc | 42 | fN = 0; |
43 | return; | |
44 | } | |
45 | ||
46 | AliITSstatistics2::~AliITSstatistics2(){ | |
47 | // | |
48 | // destructor | |
49 | // | |
07025d6d | 50 | if(fX!=0) delete[] fX; |
51 | if(fY!=0) delete[] fY; | |
52 | if(fYx!=0) delete[] fYx; | |
53 | if(fW!=0) delete[] fW; | |
54 | fX = 0; | |
55 | fY = 0; | |
56 | fYx = 0; | |
57 | fW = 0; | |
b0f5e3fc | 58 | fN = 0; |
59 | fOrder = 0; | |
60 | } | |
61 | ||
62 | //_______________________________________________________________ | |
63 | AliITSstatistics2& AliITSstatistics2::operator=(AliITSstatistics2 &source){ | |
64 | // operator = | |
65 | ||
e8189707 | 66 | if(this==&source) return *this; |
67 | if(source.fOrder!=0){ | |
68 | this->fOrder = source.fOrder; | |
69 | this->fN = source.fN; | |
07025d6d | 70 | this->fX = new Double_t[this->fOrder]; |
71 | this->fW = new Double_t[this->fOrder]; | |
e8189707 | 72 | for(Int_t i=0;i<source.fOrder;i++){ |
07025d6d | 73 | this->fX[i] = source.fX[i]; |
74 | this->fW[i] = source.fW[i]; | |
e8189707 | 75 | } // end for i |
76 | }else{ | |
07025d6d | 77 | this->fX = 0; |
78 | this->fW = 0; | |
e8189707 | 79 | this->fN = 0; |
80 | this->fOrder = 0; | |
81 | }// end if source.fOrder!=0 | |
82 | return *this; | |
b0f5e3fc | 83 | } |
84 | //_______________________________________________________________ | |
85 | AliITSstatistics2::AliITSstatistics2(AliITSstatistics2 &source){ | |
86 | // Copy constructor | |
87 | ||
e8189707 | 88 | if(this==&source) return; |
89 | if(source.fOrder!=0){ | |
90 | this->fOrder = source.fOrder; | |
91 | this->fN = source.fN; | |
07025d6d | 92 | this->fX = new Double_t[this->fOrder]; |
93 | this->fW = new Double_t[this->fOrder]; | |
e8189707 | 94 | for(Int_t i=0;i<source.fOrder;i++){ |
07025d6d | 95 | this->fX[i] = source.fX[i]; |
96 | this->fW[i] = source.fW[i]; | |
e8189707 | 97 | } // end for i |
98 | }else{ | |
07025d6d | 99 | this->fX = 0; |
100 | this->fW = 0; | |
e8189707 | 101 | this->fN = 0; |
102 | this->fOrder = 0; | |
103 | }// end if source.fOrder!=0 | |
b0f5e3fc | 104 | } |
105 | ||
106 | void AliITSstatistics2::Reset(){ | |
107 | // | |
108 | // Reset/zero statistics | |
109 | // | |
07025d6d | 110 | for(Int_t i=0;i<fOrder;i++) {fX[i] = 0.0;fY[i] = 0.0; |
111 | fYx[i] = 0.0; fW[i] = 0.0;} | |
b0f5e3fc | 112 | fN = 0; |
113 | return; | |
114 | } | |
115 | ||
116 | void AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0){ | |
117 | // | |
118 | // add next x,y pair to statistics | |
119 | // | |
120 | Double_t xs=1.0,ys=1.0,yxs=1.0,ws=1.0; | |
121 | Int_t i; | |
122 | ||
e8189707 | 123 | const Double_t kBig=1.0e+38; |
124 | ||
125 | if(y>kBig || x>kBig || w>kBig) return; | |
126 | ||
127 | ||
b0f5e3fc | 128 | fN++; |
129 | for(i=0;i<fOrder;i++){ | |
130 | xs *= x; | |
131 | ys *= y; | |
132 | yxs *= x*y; | |
133 | ws *= w; | |
07025d6d | 134 | fX[i] += xs*w; |
135 | fY[i] += ys*w; | |
136 | fYx[i] += yxs*w; | |
137 | fW[i] += ws; | |
b0f5e3fc | 138 | } // end for i |
139 | } | |
140 | ||
141 | Double_t AliITSstatistics2::GetXNth(Int_t order){ | |
142 | // | |
143 | // This give the unbiased estimator for the RMS. | |
144 | // | |
145 | ||
146 | Double_t s; | |
147 | ||
07025d6d | 148 | if(fW[0]!=0.0&&order<=fOrder) s = fX[order-1]/fW[0]; |
b0f5e3fc | 149 | else { |
150 | s = 0.0; | |
07025d6d | 151 | printf("AliITSstatistics2: error in GetNth: fOrder=%d fN=%d fW[0]=%f\n", |
152 | fOrder,fN,fW[0]); | |
b0f5e3fc | 153 | } // end else |
154 | return s; | |
155 | } | |
156 | Double_t AliITSstatistics2::GetYNth(Int_t order){ | |
157 | // | |
158 | // This give the unbiased estimator for the RMS. | |
159 | // | |
160 | Double_t s; | |
161 | ||
07025d6d | 162 | if(fW[0]!=0.0&&order<=fOrder) s = fY[order-1]/fW[0]; |
b0f5e3fc | 163 | else { |
164 | s = 0.0; | |
07025d6d | 165 | printf("AliITSstatistics2: error in GetNth: fOrder=%d fN=%d fW[0]=%f\n", |
166 | fOrder,fN,fW[0]); | |
b0f5e3fc | 167 | } // end else |
168 | return s; | |
169 | } | |
170 | Double_t AliITSstatistics2::GetYXNth(Int_t order){ | |
171 | // This give the unbiased estimator for the RMS. | |
172 | Double_t s; | |
173 | ||
07025d6d | 174 | if(fW[0]!=0.0&&order<=fOrder) s = fYx[order-1]/fW[0]; |
b0f5e3fc | 175 | else { |
176 | s = 0.0; | |
07025d6d | 177 | printf("AliITSstatistics2: error in GetNth: fOrder=%d fN=%d fW[0]=%f\n", |
178 | fOrder,fN,fW[0]); | |
b0f5e3fc | 179 | } // end else |
180 | return s; | |
181 | } | |
182 | Double_t AliITSstatistics2::GetRMSX(){ | |
183 | // This give the unbiased estimator for the RMS. | |
184 | Double_t x,x2,w,ww,s; | |
185 | ||
186 | x = GetMeanX(); // first order | |
187 | x2 = GetXNth(2); // second order | |
07025d6d | 188 | w = fW[0]; // first order - 1. |
189 | ww = fW[1]; // second order - 1. | |
b0f5e3fc | 190 | |
191 | if(w*w==ww) return (-1.0); | |
192 | s = (x2-x*x)*w*w/(w*w-ww); | |
193 | return TMath::Sqrt(s); | |
194 | } | |
195 | Double_t AliITSstatistics2::GetRMSY(){ | |
196 | // This give the unbiased estimator for the RMS. | |
197 | Double_t x,x2,w,ww,s; | |
198 | ||
199 | x = GetMeanY(); // first order | |
200 | x2 = GetYNth(2); // second order | |
07025d6d | 201 | w = fW[0]; // first order - 1. |
202 | ww = fW[1]; // second order - 1. | |
b0f5e3fc | 203 | |
204 | if(w*w==ww) return (-1.0); | |
205 | s = (x2-x*x)*w*w/(w*w-ww); | |
206 | return TMath::Sqrt(s); | |
207 | } | |
208 | Double_t AliITSstatistics2::GetRMSYX(){ | |
209 | // This give the unbiased estimator for the RMS. | |
210 | Double_t x,x2,w,ww,s; | |
211 | ||
212 | x = GetMeanYX(); // first order | |
213 | x2 = GetYXNth(2); // second order | |
07025d6d | 214 | w = fW[0]; // first order - 1. |
215 | ww = fW[1]; // second order - 1. | |
b0f5e3fc | 216 | |
217 | if(w*w==ww) return (-1.0); | |
218 | s = (x2-x*x)*w*w/(w*w-ww); | |
219 | return TMath::Sqrt(s); | |
220 | } | |
221 | Double_t AliITSstatistics2::GetErrorMeanY(){ | |
222 | //This is the error in the mean or the square root of the variance of the mean. | |
223 | Double_t rms,w,ww,s; | |
224 | ||
225 | rms = GetRMSY(); | |
07025d6d | 226 | w = fW[0]; |
227 | ww = fW[1]; | |
b0f5e3fc | 228 | s = rms*rms*ww/(w*w); |
229 | return TMath::Sqrt(s); | |
230 | } | |
231 | Double_t AliITSstatistics2::GetErrorMeanX(){ | |
232 | //This is the error in the mean or the square root of the variance of the mean. | |
233 | Double_t rms,w,ww,s; | |
234 | ||
235 | rms = GetRMSX(); | |
07025d6d | 236 | w = fW[0]; |
237 | ww = fW[1]; | |
b0f5e3fc | 238 | s = rms*rms*ww/(w*w); |
239 | return TMath::Sqrt(s); | |
240 | } | |
241 | Double_t AliITSstatistics2::GetErrorMeanYX(){ | |
242 | //This is the error in the mean or the square root of the variance of the mean. | |
243 | Double_t rms,w,ww,s; | |
244 | ||
245 | rms = GetRMSYX(); | |
07025d6d | 246 | w = fW[0]; |
247 | ww = fW[1]; | |
b0f5e3fc | 248 | s = rms*rms*ww/(w*w); |
249 | return TMath::Sqrt(s); | |
250 | } | |
251 | ||
252 | ||
253 | Double_t AliITSstatistics2::GetErrorRMSY(){ | |
254 | //This is the error in the mean or the square root of the variance of the mean. | |
255 | // at this moment this routine is only defined for weights=1. | |
256 | Double_t x,x2,x3,x4,w,ww,m2,m4,n,s; | |
257 | ||
07025d6d | 258 | if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.); |
b0f5e3fc | 259 | x = GetMeanY(); // first order |
260 | x2 = GetYNth(2); // second order | |
07025d6d | 261 | w = fW[0]; // first order - 1. |
262 | ww = fW[1]; // second order - 1. | |
b0f5e3fc | 263 | if(w*w==ww) return (-1.0); |
264 | s = (x2-x*x)*w*w/(w*w-ww); | |
265 | ||
266 | m2 = s; | |
267 | n = (Double_t) GetN(); | |
268 | x3 = GetYNth(3); | |
269 | x4 = GetYNth(4); | |
270 | // This equation assumes that all of the weights are equal to 1. | |
271 | m4 = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x); | |
272 | s = (m4-(n-3.)*m2*m2/(n-1.))/n; | |
273 | return TMath::Sqrt(s); | |
274 | } | |
275 | Double_t AliITSstatistics2::GetErrorRMSX(){ | |
276 | //This is the error in the mean or the square root of the variance of the mean. | |
277 | // at this moment this routine is only defined for weights=1. | |
278 | Double_t x,x2,x3,x4,w,ww,m2,m4,n,s; | |
279 | ||
07025d6d | 280 | if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.); |
b0f5e3fc | 281 | x = GetMeanX(); // first order |
282 | x2 = GetXNth(2); // second order | |
07025d6d | 283 | w = fW[0]; // first order - 1. |
284 | ww = fW[1]; // second order - 1. | |
b0f5e3fc | 285 | if(w*w==ww) return (-1.0); |
286 | s = (x2-x*x)*w*w/(w*w-ww); | |
287 | ||
288 | m2 = s; | |
289 | n = (Double_t) GetN(); | |
290 | x3 = GetXNth(3); | |
291 | x4 = GetXNth(4); | |
292 | // This equation assumes that all of the weights are equal to 1. | |
293 | m4 = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x); | |
294 | s = (m4-(n-3.)*m2*m2/(n-1.))/n; | |
295 | return TMath::Sqrt(s); | |
296 | } | |
297 | Double_t AliITSstatistics2::GetErrorRMSYX(){ | |
298 | //This is the error in the mean or the square root of the variance of the mean. | |
299 | // at this moment this routine is only defined for weights=1. | |
300 | Double_t x,x2,x3,x4,w,ww,m2,m4,n,s; | |
301 | ||
07025d6d | 302 | if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.); |
b0f5e3fc | 303 | x = GetMeanYX(); // first order |
304 | x2 = GetYXNth(2); // second order | |
07025d6d | 305 | w = fW[0]; // first order - 1. |
306 | ww = fW[1]; // second order - 1. | |
b0f5e3fc | 307 | if(w*w==ww) return (-1.0); |
308 | s = (x2-x*x)*w*w/(w*w-ww); | |
309 | ||
310 | m2 = s; | |
311 | n = (Double_t) GetN(); | |
312 | x3 = GetYXNth(3); | |
313 | x4 = GetYXNth(4); | |
314 | // This equation assumes that all of the weights are equal to 1. | |
315 | m4 = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x); | |
316 | s = (m4-(n-3.)*m2*m2/(n-1.))/n; | |
317 | return TMath::Sqrt(s); | |
318 | } | |
319 | //_______________________________________________________________________ | |
320 | Double_t AliITSstatistics2::FitToLine(Double_t &a,Double_t &b){ | |
321 | // fit to y = a+bx returns Chi squared or -1.0 if an error | |
322 | Double_t c,d,e,f,g,h; | |
323 | ||
324 | a = b = 0.0; | |
325 | if(fOrder<2 || fN<3) return -1.0; | |
326 | c = GetWN(1); | |
327 | d = GetYN(1); | |
328 | e = GetXN(1); | |
329 | f = GetYXN(1); | |
330 | g = GetXN(2); | |
331 | h = c*g-e*e; | |
332 | a = d*g-f*e; | |
333 | b = c*f-d*e; | |
334 | if(h!=0.0){ | |
335 | a = a/h; | |
336 | b = b/h; | |
337 | }else{ | |
338 | printf("AliITSstatistics2: Error in FitToLine vertical line\n"); | |
339 | return -1.0; | |
340 | } // end if h | |
341 | h = GetYN(2)+a*a*c+b*b*g-2.0*a*d-2.0*b*f+2.0*a*b*e; | |
342 | h /= (Double_t)fN - 2.0; | |
343 | return h; | |
344 | } |