Added drawing capabilities
[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowResolution.cxx
1 /** @file 
2     @brief Implementation of an Resolution class */
3 #include "flow/AliFMDFlowResolution.h"
4 #include "flow/AliFMDFlowUtil.h"
5 #include "flow/AliFMDFlowBessel.h"
6 #include <TGraphErrors.h>
7 #include <TH1.h>
8 #include <TH2.h>
9 #include <iostream>
10 //#include <cmath>
11
12 //====================================================================
13 void 
14 AliFMDFlowResolution::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 //____________________________________________________________________
22 Double_t 
23 AliFMDFlowResolution::Correction(UShort_t k) const 
24
25   Double_t e;
26   return Correction(k, e); 
27 }
28
29 //____________________________________________________________________
30 Double_t 
31 AliFMDFlowResolution::Correction(UShort_t, Double_t& e2) const 
32
33   e2 = fSqVar / fN;
34   return sqrt(2) * sqrt(fabs(fAverage));
35 }
36
37 //____________________________________________________________________
38 void
39 AliFMDFlowResolution::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
56 //====================================================================
57 Double_t 
58 AliFMDFlowResolutionStar::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);
65   e2             = pow(dr * delta,2);
66   return res;
67 }
68 //____________________________________________________________________
69 Double_t 
70 AliFMDFlowResolutionStar::Correction(UShort_t k) const 
71
72   Double_t e;
73   return Correction(k, e); 
74 }
75 //____________________________________________________________________
76 Double_t 
77 AliFMDFlowResolutionStar::Chi(Double_t res, UShort_t k, 
78                               Double_t& delta) const 
79 {
80   delta          = 1;
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 //____________________________________________________________________
91 Double_t 
92 AliFMDFlowResolutionStar::Res(Double_t chi, UShort_t k, Double_t& dr) const 
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
141 //____________________________________________________________________
142 void
143 AliFMDFlowResolutionStar::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
188 //====================================================================
189 void 
190 AliFMDFlowResolutionTDR::Clear(Option_t*) 
191 {
192   fN = 0;
193   fLarge = 0;
194 }
195 //____________________________________________________________________
196 void 
197 AliFMDFlowResolutionTDR::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 //____________________________________________________________________
204 Double_t 
205 AliFMDFlowResolutionTDR::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   //
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;
227   Double_t echi2 = 0;
228   Double_t y     = Chi2Over2(r, echi2);
229   return Res(k, y, echi2, e2);
230 }
231  
232   
233
234 //____________________________________________________________________
235 Double_t 
236 AliFMDFlowResolutionTDR::Res(UShort_t k, Double_t y, Double_t echi2, 
237                              Double_t& e2) const
238 {
239   // y = chi^2 / 2
240   Double_t chi   = sqrt(2 * y);
241   Double_t c     = sqrt(M_PI) * exp(-y) / 2;
242
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 }
252
253 //____________________________________________________________________
254 Double_t 
255 AliFMDFlowResolutionTDR::Correction(UShort_t k) const 
256
257   Double_t e;
258   return Correction(k, e); 
259 }
260 //____________________________________________________________________
261 Double_t 
262 AliFMDFlowResolutionTDR::Chi2Over2(Double_t r, Double_t& e2) const 
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)  
330   if (r == 0) { 
331     std::cerr << "TDR: Large/All = " << r << " <= 0!" << std::endl;
332     return 0;
333   }
334   Double_t ratio = 1. / (2*r); // Double_t(fN) / (2 * fLarge);
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   }
345   if (fLarge != 0) e2 = (r - 1) / (4 * fLarge * log(2 * r));
346   else             e2 = (r - 1) / (4 * r * log(2 * r));
347   return chi2over2;
348 }
349
350 //____________________________________________________________________
351 void
352 AliFMDFlowResolutionTDR::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 }
398
399 //____________________________________________________________________
400 //
401 // EOF
402 //