]>
Commit | Line | Data |
---|---|---|
97e94238 | 1 | /* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk> |
2 | * | |
3 | * This library is free software; you can redistribute it and/or | |
4 | * modify it under the terms of the GNU Lesser General Public License | |
5 | * as published by the Free Software Foundation; either version 2.1 of | |
6 | * the License, or (at your option) any later version. | |
7 | * | |
8 | * This library is distributed in the hope that it will be useful, but | |
9 | * WITHOUT ANY WARRANTY; without even the implied warranty of | |
10 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
11 | * Lesser General Public License for more details. | |
12 | * | |
13 | * You should have received a copy of the GNU Lesser General Public | |
14 | * License along with this library; if not, write to the Free Software | |
15 | * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 | |
16 | * USA | |
17 | */ | |
39eefe19 | 18 | /** @file |
19 | @brief Implementation of an Resolution class */ | |
97e94238 | 20 | //____________________________________________________________________ |
21 | // | |
22 | // Calculate the event plane resolution. | |
23 | // Input is the observed phis. | |
24 | // There's a number of implementations of this. | |
25 | // One is based on a naive interpretation of the Voloshin paper | |
26 | // Another is taken from aliroot/PWG2/FLOW | |
27 | // An finally one is based on Ollitraut's article. | |
39eefe19 | 28 | #include "flow/AliFMDFlowResolution.h" |
29 | #include "flow/AliFMDFlowUtil.h" | |
30 | #include "flow/AliFMDFlowBessel.h" | |
824e71c1 | 31 | #include <TGraphErrors.h> |
32 | #include <TH1.h> | |
33 | #include <TH2.h> | |
39eefe19 | 34 | #include <iostream> |
9b98d361 | 35 | #include <TBrowser.h> |
39eefe19 | 36 | //#include <cmath> |
37 | ||
38 | //==================================================================== | |
9b98d361 | 39 | AliFMDFlowResolution::AliFMDFlowResolution(UShort_t n) |
40 | : fOrder(n), | |
41 | fContrib("contrib", Form("cos(%d(#Psi_{A} - #Psi_{B}))", fOrder), 100,-1,1) | |
42 | { | |
43 | fContrib.SetDirectory(0); | |
44 | fContrib.SetXTitle(Form("cos(%d(#Psi_{A} - #Psi_{B}))", fOrder)); | |
45 | } | |
46 | ||
47 | //____________________________________________________________________ | |
97e94238 | 48 | AliFMDFlowResolution::AliFMDFlowResolution(const AliFMDFlowResolution& o) |
49 | : AliFMDFlowStat(o), | |
9b98d361 | 50 | fOrder(o.fOrder), |
51 | fContrib(o.fContrib) | |
97e94238 | 52 | { |
53 | // Copy constructor | |
54 | // Parameters: | |
55 | // o Object to copy from. | |
9b98d361 | 56 | fContrib.SetDirectory(0); |
57 | fContrib.SetXTitle(Form("cos(%d(#Psi_{A} - #Psi_{B}", fOrder)); | |
97e94238 | 58 | } |
59 | //____________________________________________________________________ | |
60 | AliFMDFlowResolution& | |
61 | AliFMDFlowResolution::operator=(const AliFMDFlowResolution& o) | |
62 | { | |
63 | // Assignment operator | |
64 | // Parameter: | |
65 | // o Object to assign from | |
66 | // | |
67 | AliFMDFlowStat::operator=(o); | |
68 | fOrder = o.fOrder; | |
9b98d361 | 69 | |
70 | fContrib.Reset(); | |
71 | fContrib.Add(&o.fContrib); | |
72 | ||
97e94238 | 73 | return *this; |
74 | } | |
75 | //____________________________________________________________________ | |
39eefe19 | 76 | void |
77 | AliFMDFlowResolution::Add(Double_t psiA, Double_t psiB) | |
78 | { | |
97e94238 | 79 | // add data point |
80 | // Parameters: | |
81 | // psiA A sub-event plane angle Psi_A in [0,2pi] | |
82 | // psiB B sub-event plane angle Psi_B in [0,2pi] | |
39eefe19 | 83 | Double_t diff = NormalizeAngle(fOrder * (psiA - psiB)); |
84 | Double_t contrib = cos(diff); | |
85 | AliFMDFlowStat::Add(contrib); | |
9b98d361 | 86 | fContrib.Fill(contrib); |
39eefe19 | 87 | } |
88 | ||
89 | //____________________________________________________________________ | |
90 | Double_t | |
91 | AliFMDFlowResolution::Correction(UShort_t k) const | |
92 | { | |
97e94238 | 93 | // Get the correction for harmonic strength of order @a k |
94 | // k The harminic strenght order to get the correction for | |
95 | // Returns <cos(n(psi - psi_R))> | |
39eefe19 | 96 | Double_t e; |
97 | return Correction(k, e); | |
98 | } | |
99 | ||
100 | //____________________________________________________________________ | |
101 | Double_t | |
102 | AliFMDFlowResolution::Correction(UShort_t, Double_t& e2) const | |
103 | { | |
97e94238 | 104 | // Get the correction for harmonic strength of order k |
105 | // Parameters: | |
106 | // k The harminic strenght order to get the correction for | |
107 | // e2 The square error on the correction | |
108 | // Returns <cos(n(psi - psi_R))> | |
39eefe19 | 109 | e2 = fSqVar / fN; |
110 | return sqrt(2) * sqrt(fabs(fAverage)); | |
111 | } | |
112 | ||
824e71c1 | 113 | //____________________________________________________________________ |
114 | void | |
115 | AliFMDFlowResolution::Draw(Option_t* option) | |
116 | { | |
97e94238 | 117 | // Draw this corrrection function |
118 | // Parameters: | |
119 | // option String of options. Passed to TGraph::Draw | |
824e71c1 | 120 | TGraph* g = new TGraph(100); |
121 | for (UShort_t i = 0; i < g->GetN(); i++) { | |
122 | Double_t x = -1. + 2. / 100 * i; | |
123 | Double_t y = sqrt(2) * sqrt(fabs(x)); | |
124 | g->SetPoint(i, x, y); | |
125 | } | |
126 | g->SetName("naive_res"); | |
127 | g->SetTitle("Naive Resolution Function"); | |
128 | g->GetHistogram()->SetXTitle("<cos(n(|#Psi_{A}-#Psi_{B}|))>"); | |
129 | g->GetHistogram()->SetYTitle("<cos(n(#Psi-#Psi_{R}))>"); | |
130 | g->GetHistogram()->SetDirectory(0); | |
131 | g->Draw(Form("lh %s", option)); | |
132 | } | |
133 | ||
9b98d361 | 134 | //____________________________________________________________________ |
135 | void | |
136 | AliFMDFlowResolution::Browse(TBrowser* b) | |
137 | { | |
138 | b->Add(&fContrib); | |
139 | } | |
824e71c1 | 140 | |
39eefe19 | 141 | //==================================================================== |
142 | Double_t | |
143 | AliFMDFlowResolutionStar::Correction(UShort_t k, Double_t& e2) const | |
144 | { | |
97e94238 | 145 | // Get the correction for harmonic strength of order k |
146 | // Parameters: | |
147 | // k The harminic strenght order to get the correction for | |
148 | // e2 The square error on the correction | |
149 | // Returns <cos(n(psi - psi_R))> | |
39eefe19 | 150 | if (k > 4) return 0; |
151 | Double_t delta = 0; | |
152 | Double_t chi = Chi(fAverage, k, delta); | |
153 | Double_t dr = 0; | |
154 | Double_t res = Res(sqrt(2) * chi, k, dr); | |
824e71c1 | 155 | e2 = pow(dr * delta,2); |
39eefe19 | 156 | return res; |
157 | } | |
158 | //____________________________________________________________________ | |
159 | Double_t | |
160 | AliFMDFlowResolutionStar::Correction(UShort_t k) const | |
161 | { | |
97e94238 | 162 | // Get the correction for harmonic strength of order @a k |
163 | // k The harminic strenght order to get the correction for | |
164 | // Returns <cos(n(psi - psi_R))> | |
39eefe19 | 165 | Double_t e; |
166 | return Correction(k, e); | |
167 | } | |
168 | //____________________________________________________________________ | |
169 | Double_t | |
170 | AliFMDFlowResolutionStar::Chi(Double_t res, UShort_t k, | |
171 | Double_t& delta) const | |
172 | { | |
97e94238 | 173 | // Get chi |
174 | // Parameters: | |
175 | // res First shot at the resolution. | |
176 | // k Order | |
177 | // delta On return, the last step size in \chi - | |
178 | // which is taken to be delta chi | |
179 | // Returns chi | |
824e71c1 | 180 | delta = 1; |
39eefe19 | 181 | Double_t chi = 2; |
182 | Double_t dr = 0; | |
183 | for (UInt_t i = 0; i < 15; i++) { | |
184 | if (Res(chi, k, dr) < res) chi += delta; | |
185 | else chi -= delta; | |
186 | delta /= 2; | |
187 | } | |
188 | return chi; | |
189 | } | |
190 | //____________________________________________________________________ | |
191 | Double_t | |
824e71c1 | 192 | AliFMDFlowResolutionStar::Res(Double_t chi, UShort_t k, Double_t& dr) const |
39eefe19 | 193 | { |
194 | // The resolution function is | |
195 | // | |
196 | // sqrt(pi) x exp(-x/4) (f1(x^2/4) + f2(x^2/4)) | |
197 | // r(x) = -------------------------------------------- | |
198 | // 2 sqrt(2) | |
199 | // | |
200 | // | |
201 | // = c x (f1(y) - f2(y)) | |
202 | // | |
203 | // where f1 is the modified Bessel function first kind I_{(k-1)/2}, | |
204 | // and f2 is the modified Bessel function of the first kind | |
205 | // I_{(k+1)/2}, and | |
206 | // | |
207 | // sqrt(pi) exp(-x^2/4) | |
208 | // c = --------------------, y = x^2/4 | |
209 | // 2 sqrt(2) | |
210 | // | |
211 | // The derivative of the resolution function is | |
212 | // | |
213 | // c | |
214 | // r'(y) = - (4 sqrt(y) (f1'(y) - f2'(y)) - (4 y - 2)(f1(y) - f2(y))) | |
215 | // 2 | |
216 | // | |
217 | // c r(y) | |
218 | // = - (4 sqrt(y) (f1'(y) - f2'(y))) + --------- - sqrt(y) r(y) | |
219 | // 2 2 sqrt(y) | |
220 | // | |
221 | // Since dI_n(x)/dx = I_(n-1)(x) - n / x I_n(x), and substituting | |
222 | // f3(y) = I_((k-3)/2)(y) | |
223 | // | |
224 | // c | |
225 | // r'(y) = - ((4 - 2 k) f1(y) - (4 y + 2 k) f2(y) + 4 y f3(y)) | |
226 | // 2 | |
227 | // | |
228 | Double_t y = chi * chi / 4; | |
229 | Double_t c = sqrt(M_PI) * exp(-y) / (2 * sqrt(2)); | |
230 | ||
231 | // i[0] = I[(k-1)/2-1], i[1] = I[(k-1)/2], i[2] = I[(k-1)/2+1] | |
232 | Double_t i[3], di[3]; | |
233 | AliFMDFlowBessel::Inu(Double_t(k-3)/2, Double_t(k+1)/2, y, i, di); | |
234 | ||
235 | Double_t r = c * chi * (i[2] + i[1]); | |
236 | dr = c / 2 * ((4 - 2*k)*i[1] - (4*y + 2*k)*i[2] + 4*y*i[0]); | |
237 | ||
238 | return r; | |
239 | } | |
240 | ||
824e71c1 | 241 | //____________________________________________________________________ |
242 | void | |
243 | AliFMDFlowResolutionStar::Draw(Option_t* option) | |
244 | { | |
97e94238 | 245 | // Draw this corrrection function |
246 | // Parameters: | |
247 | // option String of options. Passed to TGraph::Draw | |
248 | // Options: | |
249 | // chi Draw chi rather than the resolution. | |
824e71c1 | 250 | TString opt(option); |
251 | opt.ToLower(); | |
252 | Bool_t chi = opt.Contains("chi"); | |
253 | ||
254 | TH2* h = new TH2D(Form("star_%s_frame", (chi ? "chi" : "res")), | |
255 | Form("STAR %s Function for k=(1,2,4)", | |
256 | (chi ? "Chi" : "Resolution")), | |
257 | 100, 0, 1, 100, 0, (chi ? 3 : 1.5)); | |
258 | h->SetXTitle("<cos(n(|#Psi_{A}-#Psi_{B}|))>"); | |
259 | h->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>")); | |
260 | h->SetStats(0); | |
261 | h->SetDirectory(0); | |
262 | h->Draw(); | |
263 | ||
264 | for (UShort_t k = 1; k <= 4; k++) { | |
265 | if (k == 3) continue; | |
266 | ||
267 | TGraphErrors* g = new TGraphErrors(100); | |
268 | for (UShort_t i = 0; i < g->GetN(); i++) { | |
269 | Double_t e2 = 0; | |
270 | Double_t x = 1. / 100 * i; | |
271 | Double_t y = 0; | |
272 | Double_t c = Chi(x, k, e2); | |
273 | if (chi) y = c; | |
274 | else { | |
275 | Double_t dr = 0; | |
276 | y = Res(sqrt(2) * c, k, dr); | |
277 | e2 = pow(dr * e2,2); | |
278 | } | |
279 | g->SetPoint(i, x, y); | |
280 | g->SetPointError(i, 0, sqrt(e2)); | |
281 | } | |
282 | g->SetLineColor(k); | |
283 | g->SetName(Form("star_%s_k%d", (chi ? "chi" : "res"), k)); | |
284 | g->SetTitle(Form("STAR %s Function for k=%d", | |
285 | (chi ? "Chi" : "Resolution"), k)); | |
286 | g->GetHistogram()->SetXTitle("<cos(n(|#Psi_{A}-#Psi_{B}|))>"); | |
287 | g->GetHistogram()->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>")); | |
288 | g->GetHistogram()->SetDirectory(0); | |
289 | g->Draw(Form("l %s same", option)); | |
290 | } | |
291 | } | |
292 | ||
39eefe19 | 293 | //==================================================================== |
9b98d361 | 294 | AliFMDFlowResolutionTDR::AliFMDFlowResolutionTDR(UShort_t n) |
295 | : AliFMDFlowResolution(n), | |
296 | fLarge(0), | |
297 | fNK("nk", "Number of events with (#Psi_{}A-#Psi_{B})>90^{o}", 2, 0, 2) | |
298 | { | |
299 | fNK.SetDirectory(0); | |
300 | fNK.SetYTitle("# events"); | |
301 | fNK.GetXaxis()->SetBinLabel(1, "k"); | |
302 | fNK.GetXaxis()->SetBinLabel(2, "N"); | |
303 | } | |
304 | ||
305 | //____________________________________________________________________ | |
97e94238 | 306 | AliFMDFlowResolutionTDR::AliFMDFlowResolutionTDR(const |
307 | AliFMDFlowResolutionTDR& o) | |
308 | : AliFMDFlowResolution(o), | |
9b98d361 | 309 | fLarge(o.fLarge), |
310 | fNK(o.fNK) | |
311 | { | |
312 | fNK.SetDirectory(0); | |
313 | fNK.SetYTitle("# events"); | |
314 | fNK.GetXaxis()->SetBinLabel(1, "k"); | |
315 | fNK.GetXaxis()->SetBinLabel(2, "N"); | |
316 | } | |
97e94238 | 317 | //____________________________________________________________________ |
318 | AliFMDFlowResolutionTDR& | |
319 | AliFMDFlowResolutionTDR::operator=(const AliFMDFlowResolutionTDR& o) | |
320 | { | |
321 | // Assignment operator | |
322 | // Parameter: | |
323 | // o Object to assign from | |
324 | // | |
325 | AliFMDFlowResolution::operator=(o); | |
326 | fLarge = o.fLarge; | |
9b98d361 | 327 | fNK.Reset(); |
328 | fNK.Add(&fNK); | |
97e94238 | 329 | return *this; |
330 | } | |
331 | //____________________________________________________________________ | |
39eefe19 | 332 | void |
824e71c1 | 333 | AliFMDFlowResolutionTDR::Clear(Option_t*) |
39eefe19 | 334 | { |
97e94238 | 335 | // Clear internal variables. |
336 | // Parameters: | |
337 | // options Ignored | |
39eefe19 | 338 | fN = 0; |
339 | fLarge = 0; | |
340 | } | |
341 | //____________________________________________________________________ | |
342 | void | |
97e94238 | 343 | AliFMDFlowResolutionTDR::Add(Double_t psiA, Double_t psiB) |
39eefe19 | 344 | { |
97e94238 | 345 | // add data point. If |psi_a - psi_b| >= pi/2 increase large |
346 | // counter. | |
347 | // Parameters: | |
348 | // psiA A sub-event plane angle Psi_A in [0,2pi] | |
349 | // psiB B sub-event plane angle Psi_B in [0,2pi] | |
9b98d361 | 350 | AliFMDFlowResolution::Add(psiA, psiB); |
97e94238 | 351 | Double_t a = fabs(psiA - psiB); |
9b98d361 | 352 | if (a >= .5 * M_PI) { |
353 | fNK.Fill(.5); | |
354 | fLarge++; | |
355 | } | |
356 | fNK.Fill(1.5); | |
39eefe19 | 357 | fN++; |
358 | } | |
359 | //____________________________________________________________________ | |
360 | Double_t | |
361 | AliFMDFlowResolutionTDR::Correction(UShort_t k, Double_t& e2) const | |
362 | { | |
363 | // From nucl-ex/9711003v2 | |
364 | // | |
365 | // | |
366 | // <cos n Delta phi> = | |
367 | // | |
368 | // sqrt(pi) | |
369 | // -------- chi exp(-z) (I_((n-1)/2)(z) + I_((n+1)/2)(z)) | |
370 | // 2 | |
371 | // | |
372 | // where z = chi^2 / 2 | |
373 | // | |
824e71c1 | 374 | if (fLarge == 0) { |
9b98d361 | 375 | // std::cerr << "TDR: K = 0" << std::endl; |
824e71c1 | 376 | return -1; |
377 | } | |
378 | if (fN == 0) { | |
9b98d361 | 379 | // std::cerr << "TDR: N = 0" << std::endl; |
824e71c1 | 380 | return -1; |
381 | } | |
382 | Double_t r = Double_t(fLarge) / fN; | |
39eefe19 | 383 | Double_t echi2 = 0; |
824e71c1 | 384 | Double_t y = Chi2Over2(r, echi2); |
385 | return Res(k, y, echi2, e2); | |
386 | } | |
387 | ||
388 | ||
389 | ||
390 | //____________________________________________________________________ | |
391 | Double_t | |
392 | AliFMDFlowResolutionTDR::Res(UShort_t k, Double_t y, Double_t echi2, | |
393 | Double_t& e2) const | |
394 | { | |
97e94238 | 395 | // The resolution function is |
396 | // | |
397 | // sqrt(pi) x exp(-x/4) (f1(x^2/4) + f2(x^2/4)) | |
398 | // r(x) = -------------------------------------------- | |
399 | // 2 sqrt(2) | |
400 | // | |
401 | // | |
402 | // = c x (f1(y) - f2(y)) | |
403 | // | |
404 | // where f1 is the modified Bessel function first kind I_{(k-1)/2}, | |
405 | // and f2 is the modified Bessel function of the first kind | |
406 | // I_{(k+1)/2}, and | |
407 | // | |
408 | // sqrt(pi) exp(-x^2/4) | |
409 | // c = --------------------, y = x^2/4 | |
410 | // 2 sqrt(2) | |
411 | // | |
412 | // The derivative of the resolution function is | |
413 | // | |
414 | // c | |
415 | // r'(y) = - (4 sqrt(y) (f1'(y) - f2'(y)) - (4 y - 2)(f1(y) - f2(y))) | |
416 | // 2 | |
417 | // | |
418 | // c r(y) | |
419 | // = - (4 sqrt(y) (f1'(y) - f2'(y))) + --------- - sqrt(y) r(y) | |
420 | // 2 2 sqrt(y) | |
421 | // | |
422 | // Since dI_n(x)/dx = I_(n-1)(x) - n / x I_n(x), and substituting | |
423 | // f3(y) = I_((k-3)/2)(y) | |
424 | // | |
425 | // c | |
426 | // r'(y) = - ((4 - 2 k) f1(y) - (4 y + 2 k) f2(y) + 4 y f3(y)) | |
427 | // 2 | |
428 | // | |
824e71c1 | 429 | // y = chi^2 / 2 |
39eefe19 | 430 | Double_t chi = sqrt(2 * y); |
431 | Double_t c = sqrt(M_PI) * exp(-y) / 2; | |
824e71c1 | 432 | |
39eefe19 | 433 | // i[0] = I[(k-1)/2-1], i[1] = I[(k-1)/2], i[2] = I[(k-1)/2+1] |
434 | Double_t i[3], di[3]; | |
435 | AliFMDFlowBessel::Inu(Double_t(k-3)/2, Double_t(k+1)/2, y, i, di); | |
436 | ||
437 | Double_t r = c * chi * (i[2] + i[1]); | |
438 | Double_t dr = c / 2 * ((4 - 2*k)*i[1] - (4*y + 2*k)*i[2] + 4*y*i[0]); | |
439 | e2 = dr * dr * echi2; | |
440 | return r; | |
441 | } | |
824e71c1 | 442 | |
39eefe19 | 443 | //____________________________________________________________________ |
444 | Double_t | |
445 | AliFMDFlowResolutionTDR::Correction(UShort_t k) const | |
446 | { | |
97e94238 | 447 | // Get the correction for harmonic strength of order k |
448 | // Parameters: | |
449 | // k The harminic strenght order to get the correction for | |
450 | // Returns <cos(n(psi - psi_R))> | |
39eefe19 | 451 | Double_t e; |
452 | return Correction(k, e); | |
453 | } | |
454 | //____________________________________________________________________ | |
455 | Double_t | |
824e71c1 | 456 | AliFMDFlowResolutionTDR::Chi2Over2(Double_t r, Double_t& e2) const |
39eefe19 | 457 | { |
458 | // From nucl-ex/9711003v2 | |
459 | // | |
460 | // N equal to the total number of events, and k equal to the | |
461 | // number of events where |Psi_A - Psi_B| > pi/2, we get | |
462 | // | |
463 | // k exp(-chi^2/2) k | |
464 | // - = ------------- => chi^2 / 2 = - log 2 - | |
465 | // N 2 N | |
466 | // | |
467 | // N | |
468 | // => chi^2 / 2 = log -- (*) | |
469 | // k | |
470 | // | |
471 | // 2 k | |
472 | // => chi^2 = - 2 log --- | |
473 | // N | |
474 | // | |
475 | // => chi = -/+ sgrt (- 2 log (2k/N)) | |
476 | // | |
477 | // (*) this is what is returned. | |
478 | // | |
479 | // The error on chi is given by | |
480 | // | |
481 | // dchi 1 | |
482 | // d^2chi = (------)^2 d^2(k/N) = - -------------------- d^2(k/N) | |
483 | // d(k/N) 4 (k/N)^2 log(2 k/N) | |
484 | // | |
485 | // where | |
486 | // k k^2 | |
487 | // (1 - 2 -) dk^2 + --- dN^2 | |
488 | // N N^2 | |
489 | // d^2(k/N) = -------------------------- | |
490 | // N^2 | |
491 | // | |
492 | // | |
493 | // with dk = sqrt(k) and dN = sqrt(N), this reduces to | |
494 | // | |
495 | // k k^2 N k k^2 k^2 | |
496 | // (1 - 2 -) k + ----- N --- - 2 --- + --- | |
497 | // N N^2 N N N | |
498 | // d^2(k/N) = --------------------- = -------------------------- | |
499 | // N^2 N^2 | |
500 | // | |
501 | // k N - k^2 k (N - k) k (1 - k / N) | |
502 | // = --------- = - ------- = - ----------- | |
503 | // N^3 N N^2 N N | |
504 | // | |
505 | // Identifying r = k/N, we get | |
506 | // | |
507 | // | |
508 | // d(k/N) = sqrt(r (1 - r) / N) | |
509 | // | |
510 | // which is the well-known result for Binomial errors. | |
511 | // Alternatively, one could compute these error using an Baysian | |
512 | // confidence level of say 68.3% (see for example | |
513 | // http://home.fnal.gov/~paterno/images/effic.pdf). | |
514 | // | |
515 | // Our final expression for the error on chi then becomes | |
516 | // | |
517 | // 1 | |
518 | // d^2chi = - -------------------- r (1 - r) / N | |
519 | // 4 r^2 log(2 r) | |
520 | // | |
521 | // 1 - r r - 1 | |
522 | // = - -------------- = ------------ | |
523 | // 4 r N log(2 r) 4 k log(2 r) | |
824e71c1 | 524 | if (r == 0) { |
525 | std::cerr << "TDR: Large/All = " << r << " <= 0!" << std::endl; | |
526 | return 0; | |
39eefe19 | 527 | } |
824e71c1 | 528 | Double_t ratio = 1. / (2*r); // Double_t(fN) / (2 * fLarge); |
39eefe19 | 529 | if (ratio <= 0) { |
530 | std::cerr << "TDR: Large/All = " << ratio << " <= 0!" << std::endl; | |
531 | return -1; | |
532 | } | |
533 | Double_t chi2over2 = log(ratio); | |
534 | if (chi2over2 < 0) { | |
535 | std::cerr << "TDR: log(" << ratio << ") = " << chi2over2 | |
536 | << " < 0" << std::endl; | |
537 | return -1; | |
538 | } | |
824e71c1 | 539 | if (fLarge != 0) e2 = (r - 1) / (4 * fLarge * log(2 * r)); |
540 | else e2 = (r - 1) / (4 * r * log(2 * r)); | |
39eefe19 | 541 | return chi2over2; |
542 | } | |
543 | ||
824e71c1 | 544 | //____________________________________________________________________ |
545 | void | |
546 | AliFMDFlowResolutionTDR::Draw(Option_t* option) | |
547 | { | |
97e94238 | 548 | // Draw this corrrection function |
549 | // Parameters: | |
550 | // option String of options. Passed to TGraph::Draw | |
551 | // Options: | |
552 | // chi Draw chi rather than the resolution. | |
824e71c1 | 553 | TString opt(option); |
554 | opt.ToLower(); | |
555 | Bool_t chi = opt.Contains("chi"); | |
556 | ||
557 | TH2* h = new TH2D(Form("tdr_%s_frame", (chi ? "chi" : "res")), | |
558 | Form("TDR %s Function for k=(1,2,4)", | |
559 | (chi ? "Chi" : "Resolution")), | |
560 | 100, 0, 1, 100, 0, (chi ? 3 : 1.5)); | |
561 | h->SetXTitle("K/N"); | |
562 | h->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>")); | |
563 | h->SetStats(0); | |
564 | h->Draw(); | |
565 | h->SetDirectory(0); | |
566 | ||
567 | for (UShort_t k = 1; k <= 4; k++) { | |
568 | if (k == 3) continue; | |
569 | ||
570 | TGraphErrors* g = new TGraphErrors; | |
571 | Int_t i = 0; | |
572 | for (Double_t x = 0.02; x < 0.5; x += 0.01) { | |
573 | Double_t e2 = 0; | |
574 | Double_t y = 0; | |
575 | Double_t c = Chi2Over2(x, e2); | |
576 | if (chi) y = sqrt(2 * c); | |
577 | else { | |
578 | Double_t dr = 0; | |
579 | y = Res(k, c, e2, dr); | |
580 | e2 = dr; | |
581 | } | |
582 | g->SetPoint(i, x, y); | |
583 | g->SetPointError(i, 0, sqrt(fabs(e2))); | |
584 | i++; | |
585 | } | |
586 | g->SetLineColor(k); | |
587 | g->SetName(Form("tdr_%s_k%d", (chi ? "chi" : "res"), k)); | |
588 | g->SetTitle(Form("TDR %s Function for k=%d", | |
589 | (chi ? "Chi" : "Resolution"), k)); | |
590 | g->GetHistogram()->SetXTitle("K/N"); | |
591 | g->GetHistogram()->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>")); | |
592 | g->GetHistogram()->SetDirectory(0); | |
593 | g->Draw(Form("l %s %s", (k == 0 ? "h" : "same"), option)); | |
594 | if (chi) break; | |
595 | } | |
596 | } | |
39eefe19 | 597 | |
9b98d361 | 598 | //____________________________________________________________________ |
599 | void | |
600 | AliFMDFlowResolutionTDR::Browse(TBrowser* b) | |
601 | { | |
602 | AliFMDFlowResolution::Browse(b); | |
603 | b->Add(&fNK); | |
604 | } | |
605 | ||
39eefe19 | 606 | //____________________________________________________________________ |
607 | // | |
608 | // EOF | |
609 | // |