57fdc0e2d5c556e34ee9c8563a07cedd2c29f8df
[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 <iostream>
7 //#include <cmath>
8
9 //====================================================================
10 void 
11 AliFMDFlowResolution::Add(Double_t psiA, Double_t psiB) 
12
13   Double_t diff    = NormalizeAngle(fOrder * (psiA - psiB));
14   Double_t contrib = cos(diff);
15   AliFMDFlowStat::Add(contrib);
16 }
17
18 //____________________________________________________________________
19 Double_t 
20 AliFMDFlowResolution::Correction(UShort_t k) const 
21
22   Double_t e;
23   return Correction(k, e); 
24 }
25
26 //____________________________________________________________________
27 Double_t 
28 AliFMDFlowResolution::Correction(UShort_t, Double_t& e2) const 
29
30   e2 = fSqVar / fN;
31   return sqrt(2) * sqrt(fabs(fAverage));
32 }
33
34 //====================================================================
35 Double_t 
36 AliFMDFlowResolutionStar::Correction(UShort_t k, Double_t& e2) const 
37
38   if (k > 4) return 0;
39   Double_t delta = 0;
40   Double_t chi   = Chi(fAverage, k, delta);
41   Double_t dr    = 0;
42   Double_t res   = Res(sqrt(2) * chi, k, dr);
43   e2           = pow(dr * delta,2);
44   return res;
45 }
46 //____________________________________________________________________
47 Double_t 
48 AliFMDFlowResolutionStar::Correction(UShort_t k) const 
49
50   Double_t e;
51   return Correction(k, e); 
52 }
53 //____________________________________________________________________
54 Double_t 
55 AliFMDFlowResolutionStar::Chi(Double_t res, UShort_t k, 
56                               Double_t& delta) const 
57 {
58   delta        = 1;
59   Double_t chi   = 2;
60   Double_t dr    = 0;
61   for (UInt_t i = 0; i < 15; i++) { 
62     if (Res(chi, k, dr) < res) chi += delta;
63     else                       chi -= delta;
64     delta /= 2;
65   }
66   return chi;
67 }
68 //____________________________________________________________________
69 Double_t 
70 AliFMDFlowResolutionStar::Res(Double_t chi, UShort_t k, 
71                                     Double_t& dr) const 
72
73   // The resolution function is 
74   // 
75   //          sqrt(pi) x exp(-x/4) (f1(x^2/4) + f2(x^2/4))
76   //   r(x) = --------------------------------------------
77   //                          2 sqrt(2) 
78   // 
79   //        
80   //        = c x (f1(y) - f2(y))
81   //
82   // where f1 is the modified Bessel function first kind I_{(k-1)/2}, 
83   // and f2 is the modified Bessel function of the first kind
84   // I_{(k+1)/2}, and 
85   // 
86   //          sqrt(pi) exp(-x^2/4) 
87   //      c = --------------------,   y = x^2/4
88   //              2 sqrt(2)
89   // 
90   // The derivative of the resolution function is 
91   //
92   //            c 
93   //    r'(y) = - (4 sqrt(y) (f1'(y) - f2'(y)) - (4 y - 2)(f1(y) - f2(y)))
94   //            2
95   // 
96   //            c                                    r(y)   
97   //          = - (4 sqrt(y) (f1'(y) - f2'(y))) + --------- - sqrt(y) r(y)
98   //            2                                 2 sqrt(y)       
99   // 
100   // Since dI_n(x)/dx = I_(n-1)(x) - n / x I_n(x), and substituting 
101   // f3(y) = I_((k-3)/2)(y) 
102   // 
103   //            c  
104   //    r'(y) = - ((4 - 2 k) f1(y) - (4 y + 2 k) f2(y) + 4 y f3(y))
105   //            2   
106   // 
107   Double_t y  = chi * chi / 4;
108   Double_t c  = sqrt(M_PI) * exp(-y) / (2 * sqrt(2));   
109   
110   // i[0] = I[(k-1)/2-1], i[1] = I[(k-1)/2], i[2] = I[(k-1)/2+1]
111   Double_t i[3], di[3];
112   AliFMDFlowBessel::Inu(Double_t(k-3)/2, Double_t(k+1)/2, y, i, di);
113   
114   Double_t r  = c * chi * (i[2] + i[1]);
115   dr        = c / 2 * ((4 - 2*k)*i[1] - (4*y + 2*k)*i[2] + 4*y*i[0]);
116
117   return r;  
118 }
119
120 //====================================================================
121 void 
122 AliFMDFlowResolutionTDR::Clear() 
123 {
124   fN = 0;
125   fLarge = 0;
126 }
127 //____________________________________________________________________
128 void 
129 AliFMDFlowResolutionTDR::Add(Double_t psi_a, Double_t psi_b)
130
131   Double_t a = fabs(psi_a - psi_b);
132   if (a >= .5 * M_PI) fLarge++;
133   fN++;
134 }
135 //____________________________________________________________________
136 Double_t 
137 AliFMDFlowResolutionTDR::Correction(UShort_t k, Double_t& e2) const 
138
139   // From nucl-ex/9711003v2 
140   //
141   //
142   //   <cos n Delta phi> = 
143   //
144   //         sqrt(pi)
145   //         -------- chi exp(-z) (I_((n-1)/2)(z) + I_((n+1)/2)(z))
146   //            2
147   // 
148   // where z = chi^2 / 2
149   //
150   Double_t echi2 = 0;
151   Double_t y     = Chi2Over2(echi2);
152   Double_t chi   = sqrt(2 * y);
153   Double_t c     = sqrt(M_PI) * exp(-y) / 2;
154   
155   // i[0] = I[(k-1)/2-1], i[1] = I[(k-1)/2], i[2] = I[(k-1)/2+1]
156   Double_t i[3], di[3];
157   AliFMDFlowBessel::Inu(Double_t(k-3)/2, Double_t(k+1)/2, y, i, di);
158   
159   Double_t r  = c * chi * (i[2] + i[1]);
160   Double_t dr = c / 2 * ((4 - 2*k)*i[1] - (4*y + 2*k)*i[2] + 4*y*i[0]);
161   e2        = dr * dr * echi2;
162   return r;
163 }
164 //____________________________________________________________________
165 Double_t 
166 AliFMDFlowResolutionTDR::Correction(UShort_t k) const 
167
168   Double_t e;
169   return Correction(k, e); 
170 }
171 //____________________________________________________________________
172 Double_t 
173 AliFMDFlowResolutionTDR::Chi2Over2(Double_t& e2) const 
174 {
175   // From nucl-ex/9711003v2 
176   //
177   // N equal to the total number of events, and k equal to the
178   // number of events where |Psi_A - Psi_B| > pi/2, we get 
179   //
180   //   k   exp(-chi^2/2)                               k
181   //   - = -------------     =>    chi^2 / 2 = - log 2 -
182   //   N          2                                    N
183   //
184   //                                               N
185   //                         =>    chi^2 / 2 = log --   (*)
186   //                                               k
187   //            
188   //                                                   2 k
189   //                         =>    chi^2     = - 2 log ---
190   //                                                    N
191   //                                                            
192   //                         =>    chi       = -/+ sgrt (- 2 log (2k/N))
193   //
194   // (*) this is what is returned.  
195   //
196   //  The error on chi is given by 
197   //
198   //              dchi                           1
199   //    d^2chi = (------)^2 d^2(k/N) = - -------------------- d^2(k/N)
200   //              d(k/N)                 4 (k/N)^2 log(2 k/N)
201   // 
202   // where 
203   //                      k         k^2
204   //               (1 - 2 -) dk^2 + --- dN^2
205   //                      N         N^2
206   //   d^2(k/N) = --------------------------
207   //                           N^2 
208   // 
209   // 
210   // with dk = sqrt(k) and dN = sqrt(N), this reduces to 
211   //
212   //                     k      k^2        N k     k^2   k^2
213   //              (1 - 2 -) k + ----- N    --- - 2 --- + ---
214   //                     N      N^2         N       N     N
215   //   d^2(k/N) = --------------------- = --------------------------
216   //                           N^2                 N^2 
217   //
218   //              k N - k^2   k (N - k)   k (1 - k / N)
219   //            = --------- = - ------- = - -----------
220   //                 N^3      N   N^2     N      N
221   //
222   // Identifying r = k/N, we get 
223   //
224   //            
225   //   d(k/N)   = sqrt(r (1 - r) / N) 
226   //
227   // which is the well-known result for Binomial errors.
228   // Alternatively, one could compute these error using an Baysian
229   // confidence level of say 68.3% (see for example 
230   // http://home.fnal.gov/~paterno/images/effic.pdf). 
231   // 
232   // Our final expression for the error on chi then becomes 
233   //
234   //                      1          
235   //   d^2chi = - -------------------- r (1 - r) / N
236   //              4 r^2 log(2 r)
237   //
238   //                 1 - r            r - 1
239   //          = - -------------- = ------------
240   //              4 r N log(2 r)   4 k log(2 r)  
241
242   if (fLarge == 0) { 
243     std::cerr << "TDR: Large = 0" << std::endl;
244     return -1;
245   }
246   Double_t ratio = Double_t(fN) / (2 * fLarge);
247   if (ratio <= 0) {
248     std::cerr << "TDR: Large/All = " << ratio << " <= 0!" << std::endl;
249     return -1;
250   }
251   Double_t chi2over2  = log(ratio);
252   if (chi2over2 < 0) { 
253     std::cerr << "TDR: log(" << ratio << ") = " << chi2over2 
254               << " < 0" << std::endl; 
255     return -1;
256   }
257   Double_t r = Double_t(fLarge) / fN;
258   e2 = (r - 1) / (4 * fLarge * log(2 * r));
259   return chi2over2;
260 }
261
262
263 //____________________________________________________________________
264 //
265 // EOF
266 //