bug fix
[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowResolution.cxx
CommitLineData
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 39AliFMDFlowResolution::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 48AliFMDFlowResolution::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//____________________________________________________________________
60AliFMDFlowResolution&
61AliFMDFlowResolution::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 76void
77AliFMDFlowResolution::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//____________________________________________________________________
90Double_t
91AliFMDFlowResolution::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//____________________________________________________________________
101Double_t
102AliFMDFlowResolution::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;
2a4540e6 110 return sqrt(Double_t(2)) * sqrt(fabs(fAverage));
39eefe19 111}
112
824e71c1 113//____________________________________________________________________
114void
115AliFMDFlowResolution::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;
2a4540e6 123 Double_t y = sqrt(Double_t(2)) * sqrt(fabs(x));
824e71c1 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//____________________________________________________________________
135void
136AliFMDFlowResolution::Browse(TBrowser* b)
137{
138 b->Add(&fContrib);
139}
824e71c1 140
39eefe19 141//====================================================================
142Double_t
143AliFMDFlowResolutionStar::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;
2a4540e6 154 Double_t res = Res(sqrt(Double_t(2)) * chi, k, dr);
824e71c1 155 e2 = pow(dr * delta,2);
39eefe19 156 return res;
157}
158//____________________________________________________________________
159Double_t
160AliFMDFlowResolutionStar::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//____________________________________________________________________
169Double_t
170AliFMDFlowResolutionStar::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//____________________________________________________________________
191Double_t
824e71c1 192AliFMDFlowResolutionStar::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;
2a4540e6 229 Double_t c = sqrt(M_PI) * exp(-y) / (2 * sqrt(Double_t(2)));
39eefe19 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//____________________________________________________________________
242void
243AliFMDFlowResolutionStar::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;
2a4540e6 276 y = Res(sqrt(Double_t(2)) * c, k, dr);
824e71c1 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 294AliFMDFlowResolutionTDR::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 306AliFMDFlowResolutionTDR::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//____________________________________________________________________
318AliFMDFlowResolutionTDR&
319AliFMDFlowResolutionTDR::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 332void
824e71c1 333AliFMDFlowResolutionTDR::Clear(Option_t*)
39eefe19 334{
97e94238 335 // Clear internal variables.
336 // Parameters:
337 // options Ignored
39eefe19 338 fN = 0;
339 fLarge = 0;
340}
341//____________________________________________________________________
342void
97e94238 343AliFMDFlowResolutionTDR::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//____________________________________________________________________
360Double_t
361AliFMDFlowResolutionTDR::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//____________________________________________________________________
391Double_t
392AliFMDFlowResolutionTDR::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//____________________________________________________________________
444Double_t
445AliFMDFlowResolutionTDR::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//____________________________________________________________________
455Double_t
824e71c1 456AliFMDFlowResolutionTDR::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//____________________________________________________________________
545void
546AliFMDFlowResolutionTDR::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
598//____________________________________________________________________
9b98d361 599void
600AliFMDFlowResolutionTDR::Browse(TBrowser* b)
601{
602 AliFMDFlowResolution::Browse(b);
603 b->Add(&fNK);
604}
605
606//____________________________________________________________________
39eefe19 607//
608// EOF
609//