coverity 15108 fixed
[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowResolution.cxx
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  */
18 /** @file 
19     @brief Implementation of an Resolution class */
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.
28 #include "flow/AliFMDFlowResolution.h"
29 #include "flow/AliFMDFlowUtil.h"
30 #include "flow/AliFMDFlowBessel.h"
31 #include <TGraphErrors.h>
32 #include <TH1.h>
33 #include <TH2.h>
34 #include <iostream>
35 #include <TBrowser.h>
36 //#include <cmath>
37
38 //====================================================================
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 //____________________________________________________________________
48 AliFMDFlowResolution::AliFMDFlowResolution(const AliFMDFlowResolution& o)
49   : AliFMDFlowStat(o), 
50     fOrder(o.fOrder),
51     fContrib(o.fContrib)
52 {
53   // Copy constructor 
54   // Parameters: 
55   //   o   Object to copy from. 
56   fContrib.SetDirectory(0);
57   fContrib.SetXTitle(Form("cos(%d(#Psi_{A} - #Psi_{B}", fOrder));
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;
69
70   fContrib.Reset();
71   fContrib.Add(&o.fContrib);
72
73   return *this;
74 }
75 //____________________________________________________________________
76 void 
77 AliFMDFlowResolution::Add(Double_t psiA, Double_t psiB) 
78
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]
83   Double_t diff    = NormalizeAngle(fOrder * (psiA - psiB));
84   Double_t contrib = cos(diff);
85   AliFMDFlowStat::Add(contrib);
86   fContrib.Fill(contrib);
87 }
88
89 //____________________________________________________________________
90 Double_t 
91 AliFMDFlowResolution::Correction(UShort_t k) const 
92
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))>
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
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))>
109   e2 = fSqVar / fN;
110   return sqrt(Double_t(2)) * sqrt(fabs(fAverage));
111 }
112
113 //____________________________________________________________________
114 void
115 AliFMDFlowResolution::Draw(Option_t* option) 
116 {
117   // Draw this corrrection function 
118   // Parameters: 
119   //   option   String of options. Passed to TGraph::Draw
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(Double_t(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
134 //____________________________________________________________________
135 void
136 AliFMDFlowResolution::Browse(TBrowser* b)
137 {
138   b->Add(&fContrib);
139 }
140
141 //====================================================================
142 Double_t 
143 AliFMDFlowResolutionStar::Correction(UShort_t k, Double_t& e2) const 
144
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))>
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(Double_t(2)) * chi, k, dr);
155   e2             = pow(dr * delta,2);
156   return res;
157 }
158 //____________________________________________________________________
159 Double_t 
160 AliFMDFlowResolutionStar::Correction(UShort_t k) const 
161
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))>
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 {
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
180   delta          = 1;
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 
192 AliFMDFlowResolutionStar::Res(Double_t chi, UShort_t k, Double_t& dr) const 
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(Double_t(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
241 //____________________________________________________________________
242 void
243 AliFMDFlowResolutionStar::Draw(Option_t* option) 
244 {
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. 
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(Double_t(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
293 //====================================================================
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 //____________________________________________________________________
306 AliFMDFlowResolutionTDR::AliFMDFlowResolutionTDR(const 
307                                                  AliFMDFlowResolutionTDR& o)
308   : AliFMDFlowResolution(o), 
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 }
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;
327   fNK.Reset();
328   fNK.Add(&fNK);
329   return *this;
330 }
331 //____________________________________________________________________
332 void 
333 AliFMDFlowResolutionTDR::Clear(Option_t*) 
334 {
335   // Clear internal variables. 
336   // Parameters: 
337   //   options   Ignored
338   fN = 0;
339   fLarge = 0;
340 }
341 //____________________________________________________________________
342 void 
343 AliFMDFlowResolutionTDR::Add(Double_t psiA, Double_t psiB)
344
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]
350   AliFMDFlowResolution::Add(psiA, psiB);
351   Double_t a = fabs(psiA - psiB);
352   if (a >= .5 * M_PI) { 
353     fNK.Fill(.5);
354     fLarge++;
355   }
356   fNK.Fill(1.5);
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   //
374   if (fLarge == 0) { 
375     // std::cerr << "TDR: K = 0" << std::endl;
376     return -1;
377   }
378   if (fN == 0) { 
379     // std::cerr << "TDR: N = 0" << std::endl;
380     return -1;
381   }
382   Double_t r     = Double_t(fLarge) / fN;
383   Double_t echi2 = 0;
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 {
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   // 
429   // y = chi^2 / 2
430   Double_t chi   = sqrt(2 * y);
431   Double_t c     = sqrt(M_PI) * exp(-y) / 2;
432
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 }
442
443 //____________________________________________________________________
444 Double_t 
445 AliFMDFlowResolutionTDR::Correction(UShort_t k) const 
446
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))>
451   Double_t e;
452   return Correction(k, e); 
453 }
454 //____________________________________________________________________
455 Double_t 
456 AliFMDFlowResolutionTDR::Chi2Over2(Double_t r, Double_t& e2) const 
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)  
524   if (r == 0) { 
525     std::cerr << "TDR: Large/All = " << r << " <= 0!" << std::endl;
526     return 0;
527   }
528   Double_t ratio = 1. / (2*r); // Double_t(fN) / (2 * fLarge);
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   }
539   if (fLarge != 0) e2 = (r - 1) / (4 * fLarge * log(2 * r));
540   else             e2 = (r - 1) / (4 * r * log(2 * r));
541   return chi2over2;
542 }
543
544 //____________________________________________________________________
545 void
546 AliFMDFlowResolutionTDR::Draw(Option_t* option) 
547 {
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. 
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 }
597
598 //____________________________________________________________________
599 void
600 AliFMDFlowResolutionTDR::Browse(TBrowser* b)
601 {
602   AliFMDFlowResolution::Browse(b);
603   b->Add(&fNK);
604 }
605
606 //____________________________________________________________________
607 //
608 // EOF
609 //