Added drawing capabilities
[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowResolution.cxx
CommitLineData
39eefe19 1/** @file
2 @brief Implementation of an Resolution class */
3#include "flow/AliFMDFlowResolution.h"
4#include "flow/AliFMDFlowUtil.h"
5#include "flow/AliFMDFlowBessel.h"
824e71c1 6#include <TGraphErrors.h>
7#include <TH1.h>
8#include <TH2.h>
39eefe19 9#include <iostream>
10//#include <cmath>
11
12//====================================================================
13void
14AliFMDFlowResolution::Add(Double_t psiA, Double_t psiB)
15{
16 Double_t diff = NormalizeAngle(fOrder * (psiA - psiB));
17 Double_t contrib = cos(diff);
18 AliFMDFlowStat::Add(contrib);
19}
20
21//____________________________________________________________________
22Double_t
23AliFMDFlowResolution::Correction(UShort_t k) const
24{
25 Double_t e;
26 return Correction(k, e);
27}
28
29//____________________________________________________________________
30Double_t
31AliFMDFlowResolution::Correction(UShort_t, Double_t& e2) const
32{
33 e2 = fSqVar / fN;
34 return sqrt(2) * sqrt(fabs(fAverage));
35}
36
824e71c1 37//____________________________________________________________________
38void
39AliFMDFlowResolution::Draw(Option_t* option)
40{
41 TGraph* g = new TGraph(100);
42 for (UShort_t i = 0; i < g->GetN(); i++) {
43 Double_t x = -1. + 2. / 100 * i;
44 Double_t y = sqrt(2) * sqrt(fabs(x));
45 g->SetPoint(i, x, y);
46 }
47 g->SetName("naive_res");
48 g->SetTitle("Naive Resolution Function");
49 g->GetHistogram()->SetXTitle("<cos(n(|#Psi_{A}-#Psi_{B}|))>");
50 g->GetHistogram()->SetYTitle("<cos(n(#Psi-#Psi_{R}))>");
51 g->GetHistogram()->SetDirectory(0);
52 g->Draw(Form("lh %s", option));
53}
54
55
39eefe19 56//====================================================================
57Double_t
58AliFMDFlowResolutionStar::Correction(UShort_t k, Double_t& e2) const
59{
60 if (k > 4) return 0;
61 Double_t delta = 0;
62 Double_t chi = Chi(fAverage, k, delta);
63 Double_t dr = 0;
64 Double_t res = Res(sqrt(2) * chi, k, dr);
824e71c1 65 e2 = pow(dr * delta,2);
39eefe19 66 return res;
67}
68//____________________________________________________________________
69Double_t
70AliFMDFlowResolutionStar::Correction(UShort_t k) const
71{
72 Double_t e;
73 return Correction(k, e);
74}
75//____________________________________________________________________
76Double_t
77AliFMDFlowResolutionStar::Chi(Double_t res, UShort_t k,
78 Double_t& delta) const
79{
824e71c1 80 delta = 1;
39eefe19 81 Double_t chi = 2;
82 Double_t dr = 0;
83 for (UInt_t i = 0; i < 15; i++) {
84 if (Res(chi, k, dr) < res) chi += delta;
85 else chi -= delta;
86 delta /= 2;
87 }
88 return chi;
89}
90//____________________________________________________________________
91Double_t
824e71c1 92AliFMDFlowResolutionStar::Res(Double_t chi, UShort_t k, Double_t& dr) const
39eefe19 93{
94 // The resolution function is
95 //
96 // sqrt(pi) x exp(-x/4) (f1(x^2/4) + f2(x^2/4))
97 // r(x) = --------------------------------------------
98 // 2 sqrt(2)
99 //
100 //
101 // = c x (f1(y) - f2(y))
102 //
103 // where f1 is the modified Bessel function first kind I_{(k-1)/2},
104 // and f2 is the modified Bessel function of the first kind
105 // I_{(k+1)/2}, and
106 //
107 // sqrt(pi) exp(-x^2/4)
108 // c = --------------------, y = x^2/4
109 // 2 sqrt(2)
110 //
111 // The derivative of the resolution function is
112 //
113 // c
114 // r'(y) = - (4 sqrt(y) (f1'(y) - f2'(y)) - (4 y - 2)(f1(y) - f2(y)))
115 // 2
116 //
117 // c r(y)
118 // = - (4 sqrt(y) (f1'(y) - f2'(y))) + --------- - sqrt(y) r(y)
119 // 2 2 sqrt(y)
120 //
121 // Since dI_n(x)/dx = I_(n-1)(x) - n / x I_n(x), and substituting
122 // f3(y) = I_((k-3)/2)(y)
123 //
124 // c
125 // r'(y) = - ((4 - 2 k) f1(y) - (4 y + 2 k) f2(y) + 4 y f3(y))
126 // 2
127 //
128 Double_t y = chi * chi / 4;
129 Double_t c = sqrt(M_PI) * exp(-y) / (2 * sqrt(2));
130
131 // i[0] = I[(k-1)/2-1], i[1] = I[(k-1)/2], i[2] = I[(k-1)/2+1]
132 Double_t i[3], di[3];
133 AliFMDFlowBessel::Inu(Double_t(k-3)/2, Double_t(k+1)/2, y, i, di);
134
135 Double_t r = c * chi * (i[2] + i[1]);
136 dr = c / 2 * ((4 - 2*k)*i[1] - (4*y + 2*k)*i[2] + 4*y*i[0]);
137
138 return r;
139}
140
824e71c1 141//____________________________________________________________________
142void
143AliFMDFlowResolutionStar::Draw(Option_t* option)
144{
145 TString opt(option);
146 opt.ToLower();
147 Bool_t chi = opt.Contains("chi");
148
149 TH2* h = new TH2D(Form("star_%s_frame", (chi ? "chi" : "res")),
150 Form("STAR %s Function for k=(1,2,4)",
151 (chi ? "Chi" : "Resolution")),
152 100, 0, 1, 100, 0, (chi ? 3 : 1.5));
153 h->SetXTitle("<cos(n(|#Psi_{A}-#Psi_{B}|))>");
154 h->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
155 h->SetStats(0);
156 h->SetDirectory(0);
157 h->Draw();
158
159 for (UShort_t k = 1; k <= 4; k++) {
160 if (k == 3) continue;
161
162 TGraphErrors* g = new TGraphErrors(100);
163 for (UShort_t i = 0; i < g->GetN(); i++) {
164 Double_t e2 = 0;
165 Double_t x = 1. / 100 * i;
166 Double_t y = 0;
167 Double_t c = Chi(x, k, e2);
168 if (chi) y = c;
169 else {
170 Double_t dr = 0;
171 y = Res(sqrt(2) * c, k, dr);
172 e2 = pow(dr * e2,2);
173 }
174 g->SetPoint(i, x, y);
175 g->SetPointError(i, 0, sqrt(e2));
176 }
177 g->SetLineColor(k);
178 g->SetName(Form("star_%s_k%d", (chi ? "chi" : "res"), k));
179 g->SetTitle(Form("STAR %s Function for k=%d",
180 (chi ? "Chi" : "Resolution"), k));
181 g->GetHistogram()->SetXTitle("<cos(n(|#Psi_{A}-#Psi_{B}|))>");
182 g->GetHistogram()->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
183 g->GetHistogram()->SetDirectory(0);
184 g->Draw(Form("l %s same", option));
185 }
186}
187
39eefe19 188//====================================================================
189void
824e71c1 190AliFMDFlowResolutionTDR::Clear(Option_t*)
39eefe19 191{
192 fN = 0;
193 fLarge = 0;
194}
195//____________________________________________________________________
196void
197AliFMDFlowResolutionTDR::Add(Double_t psi_a, Double_t psi_b)
198{
199 Double_t a = fabs(psi_a - psi_b);
200 if (a >= .5 * M_PI) fLarge++;
201 fN++;
202}
203//____________________________________________________________________
204Double_t
205AliFMDFlowResolutionTDR::Correction(UShort_t k, Double_t& e2) const
206{
207 // From nucl-ex/9711003v2
208 //
209 //
210 // <cos n Delta phi> =
211 //
212 // sqrt(pi)
213 // -------- chi exp(-z) (I_((n-1)/2)(z) + I_((n+1)/2)(z))
214 // 2
215 //
216 // where z = chi^2 / 2
217 //
824e71c1 218 if (fLarge == 0) {
219 std::cerr << "TDR: K = 0" << std::endl;
220 return -1;
221 }
222 if (fN == 0) {
223 std::cerr << "TDR: N = 0" << std::endl;
224 return -1;
225 }
226 Double_t r = Double_t(fLarge) / fN;
39eefe19 227 Double_t echi2 = 0;
824e71c1 228 Double_t y = Chi2Over2(r, echi2);
229 return Res(k, y, echi2, e2);
230}
231
232
233
234//____________________________________________________________________
235Double_t
236AliFMDFlowResolutionTDR::Res(UShort_t k, Double_t y, Double_t echi2,
237 Double_t& e2) const
238{
239 // y = chi^2 / 2
39eefe19 240 Double_t chi = sqrt(2 * y);
241 Double_t c = sqrt(M_PI) * exp(-y) / 2;
824e71c1 242
39eefe19 243 // i[0] = I[(k-1)/2-1], i[1] = I[(k-1)/2], i[2] = I[(k-1)/2+1]
244 Double_t i[3], di[3];
245 AliFMDFlowBessel::Inu(Double_t(k-3)/2, Double_t(k+1)/2, y, i, di);
246
247 Double_t r = c * chi * (i[2] + i[1]);
248 Double_t dr = c / 2 * ((4 - 2*k)*i[1] - (4*y + 2*k)*i[2] + 4*y*i[0]);
249 e2 = dr * dr * echi2;
250 return r;
251}
824e71c1 252
39eefe19 253//____________________________________________________________________
254Double_t
255AliFMDFlowResolutionTDR::Correction(UShort_t k) const
256{
257 Double_t e;
258 return Correction(k, e);
259}
260//____________________________________________________________________
261Double_t
824e71c1 262AliFMDFlowResolutionTDR::Chi2Over2(Double_t r, Double_t& e2) const
39eefe19 263{
264 // From nucl-ex/9711003v2
265 //
266 // N equal to the total number of events, and k equal to the
267 // number of events where |Psi_A - Psi_B| > pi/2, we get
268 //
269 // k exp(-chi^2/2) k
270 // - = ------------- => chi^2 / 2 = - log 2 -
271 // N 2 N
272 //
273 // N
274 // => chi^2 / 2 = log -- (*)
275 // k
276 //
277 // 2 k
278 // => chi^2 = - 2 log ---
279 // N
280 //
281 // => chi = -/+ sgrt (- 2 log (2k/N))
282 //
283 // (*) this is what is returned.
284 //
285 // The error on chi is given by
286 //
287 // dchi 1
288 // d^2chi = (------)^2 d^2(k/N) = - -------------------- d^2(k/N)
289 // d(k/N) 4 (k/N)^2 log(2 k/N)
290 //
291 // where
292 // k k^2
293 // (1 - 2 -) dk^2 + --- dN^2
294 // N N^2
295 // d^2(k/N) = --------------------------
296 // N^2
297 //
298 //
299 // with dk = sqrt(k) and dN = sqrt(N), this reduces to
300 //
301 // k k^2 N k k^2 k^2
302 // (1 - 2 -) k + ----- N --- - 2 --- + ---
303 // N N^2 N N N
304 // d^2(k/N) = --------------------- = --------------------------
305 // N^2 N^2
306 //
307 // k N - k^2 k (N - k) k (1 - k / N)
308 // = --------- = - ------- = - -----------
309 // N^3 N N^2 N N
310 //
311 // Identifying r = k/N, we get
312 //
313 //
314 // d(k/N) = sqrt(r (1 - r) / N)
315 //
316 // which is the well-known result for Binomial errors.
317 // Alternatively, one could compute these error using an Baysian
318 // confidence level of say 68.3% (see for example
319 // http://home.fnal.gov/~paterno/images/effic.pdf).
320 //
321 // Our final expression for the error on chi then becomes
322 //
323 // 1
324 // d^2chi = - -------------------- r (1 - r) / N
325 // 4 r^2 log(2 r)
326 //
327 // 1 - r r - 1
328 // = - -------------- = ------------
329 // 4 r N log(2 r) 4 k log(2 r)
824e71c1 330 if (r == 0) {
331 std::cerr << "TDR: Large/All = " << r << " <= 0!" << std::endl;
332 return 0;
39eefe19 333 }
824e71c1 334 Double_t ratio = 1. / (2*r); // Double_t(fN) / (2 * fLarge);
39eefe19 335 if (ratio <= 0) {
336 std::cerr << "TDR: Large/All = " << ratio << " <= 0!" << std::endl;
337 return -1;
338 }
339 Double_t chi2over2 = log(ratio);
340 if (chi2over2 < 0) {
341 std::cerr << "TDR: log(" << ratio << ") = " << chi2over2
342 << " < 0" << std::endl;
343 return -1;
344 }
824e71c1 345 if (fLarge != 0) e2 = (r - 1) / (4 * fLarge * log(2 * r));
346 else e2 = (r - 1) / (4 * r * log(2 * r));
39eefe19 347 return chi2over2;
348}
349
824e71c1 350//____________________________________________________________________
351void
352AliFMDFlowResolutionTDR::Draw(Option_t* option)
353{
354 TString opt(option);
355 opt.ToLower();
356 Bool_t chi = opt.Contains("chi");
357
358 TH2* h = new TH2D(Form("tdr_%s_frame", (chi ? "chi" : "res")),
359 Form("TDR %s Function for k=(1,2,4)",
360 (chi ? "Chi" : "Resolution")),
361 100, 0, 1, 100, 0, (chi ? 3 : 1.5));
362 h->SetXTitle("K/N");
363 h->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
364 h->SetStats(0);
365 h->Draw();
366 h->SetDirectory(0);
367
368 for (UShort_t k = 1; k <= 4; k++) {
369 if (k == 3) continue;
370
371 TGraphErrors* g = new TGraphErrors;
372 Int_t i = 0;
373 for (Double_t x = 0.02; x < 0.5; x += 0.01) {
374 Double_t e2 = 0;
375 Double_t y = 0;
376 Double_t c = Chi2Over2(x, e2);
377 if (chi) y = sqrt(2 * c);
378 else {
379 Double_t dr = 0;
380 y = Res(k, c, e2, dr);
381 e2 = dr;
382 }
383 g->SetPoint(i, x, y);
384 g->SetPointError(i, 0, sqrt(fabs(e2)));
385 i++;
386 }
387 g->SetLineColor(k);
388 g->SetName(Form("tdr_%s_k%d", (chi ? "chi" : "res"), k));
389 g->SetTitle(Form("TDR %s Function for k=%d",
390 (chi ? "Chi" : "Resolution"), k));
391 g->GetHistogram()->SetXTitle("K/N");
392 g->GetHistogram()->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
393 g->GetHistogram()->SetDirectory(0);
394 g->Draw(Form("l %s %s", (k == 0 ? "h" : "same"), option));
395 if (chi) break;
396 }
397}
39eefe19 398
399//____________________________________________________________________
400//
401// EOF
402//