coverity 15108 fixed
[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowResolution.h
1 // -*- mode: C++ -*-
2 /* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk>
3  *
4  * This library is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU Lesser General Public License
6  * as published by the Free Software Foundation; either version 2.1 of
7  * the License, or (at your option) any later version.
8  *
9  * This library is distributed in the hope that it will be useful, but
10  * WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12  * Lesser General Public License for more details.
13  *
14  * You should have received a copy of the GNU Lesser General Public
15  * License along with this library; if not, write to the Free Software
16  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
17  * USA
18  */
19 /** @file 
20     @brief Declaration of an Resolution class */
21 //____________________________________________________________________
22 //
23 // Calculate the event plane resolution. 
24 // Input is the observed phis. 
25 // There's a number of implementations of this. 
26 #ifndef ALIFMDFLOWRESOLUTION_H
27 #define ALIFMDFLOWRESOLUTION_H
28 #include <flow/AliFMDFlowStat.h>
29 #include <TH1D.h>
30 #include <TH2D.h>
31 class TBrowser;
32
33 //______________________________________________________
34 /** @class AliFMDFlowResolution flow/AliFMDFlowResolution.h <flow/AliFMDFlowResolution.h>
35     @brief Class to calculate the event plane resolution based on two
36     sub-events 
37     @ingroup a_basic
38     
39     This class calculates the event plane resolution based on the
40     basic formulas given in Phys. Rev. @b C58, 1671.   That is, the
41     resolution is given by 
42     @f[ 
43     R_{k} = \langle\cos(km(\Psi_m-\Psi_R))\rangle 
44     @f] 
45     where @f$ \Psi_R@f$ is the unknown @e true event plane angle, and
46     @f$ n = km@f$
47
48     Using two random sub-events, @f$ A, B@f$ we get that 
49     @f[
50     \langle\cos(n(\Psi^A_m-\Psi^B_m))\rangle = 
51     \langle\cos(n(\Psi^A_m-\Psi_R))\rangle 
52     \langle\cos(n(\Psi^B_m-\Psi_R))\rangle 
53     @f]
54     If the sub-events are of equal size, and randomly chosen, then we
55     get that 
56     @f[ 
57     \langle\cos(n(\Psi^A_m-\Psi_R))\rangle = 
58     \sqrt{\langle\cos(n(\Psi^A_m-\Psi^B_m))}
59     @f]
60     and it follows that 
61     @f{eqnarray*}
62     \langle\cos(km(\Psi_m-\Psi_R))\rangle & = & 
63     \sqrt{2}\langle\cos(n(\Psi^A_m-\Psi_R))\rangle\\
64     \sqrt{2}\sqrt{\langle\cos(n(\Psi^A_m-\Psi^B_m))}
65     @f}
66
67     Hence, the event-plane resolution is simply the square root of the
68     scaled average distance between the two sub-events, multiplied by
69     @f$ \sqrt{s}@f$ 
70
71     The error is therefor @f$ \sqrt{s}@f$ times the variance of the
72     cosine of the distance between the two sub-events. 
73 */
74 class AliFMDFlowResolution : public AliFMDFlowStat
75 {
76 public:
77   /** Constructor
78       @param n Harmonic order */
79   AliFMDFlowResolution(UShort_t n=0);
80   /** Destructor */
81   virtual ~AliFMDFlowResolution() {}
82   /** Copy constructor 
83       @param o Object to copy from */ 
84   AliFMDFlowResolution(const AliFMDFlowResolution& o);
85   /** Assignment operator
86       @param o Object to copy from */ 
87   AliFMDFlowResolution& operator=(const AliFMDFlowResolution& o);
88
89   /** @{ 
90       @name Processing */
91   /** add data point 
92       @param psiA A sub-event plane angle @f$ \Psi_A \in[0,2\pi]@f$
93       @param psiB B sub-event plane angle @f$ \Psi_B \in[0,2\pi]@f$ */
94   virtual void Add(Double_t psiA, Double_t psiB);
95   /** @} */
96
97   /** @{ 
98       @name Information */
99   /** Get the correction for harmonic strength of order @a k 
100       @param k  The harminic strenght order to get the correction for
101       @param e2 The square error on the correction 
102       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
103   virtual Double_t Correction(UShort_t k, Double_t& e2) const;
104   /** Get the correction for harmonic strength of order @a k 
105       @param k The harminic strenght order to get the correction for
106       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
107   virtual Double_t Correction(UShort_t k) const;
108   /** Get the harmnic order */
109   UShort_t Order() const { return fOrder; }
110   /** @} */
111   
112   /** @{ 
113       @name Utility */
114   /** Draw this corrrection function 
115       @param option Options passed to drawing */
116   virtual void Draw(Option_t* option=""); //*MENU*
117   /** This is a folder */ 
118   Bool_t IsFolder() const { return kTRUE; }
119   /** Look through this object */ 
120   virtual void Browse(TBrowser* b); 
121   /** @} */
122
123   /** @{ 
124       @name histograms */ 
125   /** @return Contrib histogram */
126   const TH1& ContribHistogram() const { return fContrib; }
127   /** @} */
128 protected:
129   /** Order */
130   UShort_t fOrder; // Order 
131   /** Contributions */
132   TH1D fContrib;
133   /** Define for ROOT I/O */
134   ClassDef(AliFMDFlowResolution,1);
135 };
136
137 //______________________________________________________
138 /** 
139     @ingroup a_basic
140     The event plane angle resolution function is 
141     @f[
142     R_k(\chi) = \frac{\pi} \chi e^{-\chi^2/4}\left( 
143     I_{\frac{k-1}2}(\chi^2/4) + I_{\frac{k+1}2}(\chi^2/4)\right)
144     @f]
145     Where @f$ I_n(x)@f$ is the modified Bessel function of the first
146     kind.   Identifying 
147     @f[
148     y = \chi^2/4\quad C=\frac{\sqrt{\pi}e^{-y}}{2\sqrt{2}}
149     @f]
150     and the short hands @f$ f2(y) = I_{\frac{k-1}2}@f$ and @f$ f3(y) =
151     I_{\frac{k+1}2}(y)@f$ we can write this more compact as 
152     @f[ 
153     R_k(y) = C y (f2(y) + f3(y))
154     @f] 
155     The derivative of the resolution function is 
156     @f[ 
157     R_k'(y) = \frac{C}{2}\left[4\sqrt{y}\left(f2'(y)-f3'(y)\right) 
158     - (4 y - 2)\left(f2(y) + f3(y)\right)\right]
159     @f]
160     Since 
161     @f[
162     I_\nu'(x) = I_{\nu-1}(x) - \frac{\nu}{x} I_\nu(x)\quad,
163     @f]
164     and setting @f$ f1(y) = I_{\frac{k-3}2}(y)@f$, we get 
165     @f[ 
166     R_k'(y) = \frac{C}{2}\left[4y f1(y) + (4-2k) f2(y) - (4y+2k)
167     f3(y)\right]
168     @f]
169
170     In this class, the argument @f$ \chi@f$ is estimated by finding
171     the minima of @f$ R_k(\chi)@f$ near the average of
172     @f$\cos(n(\Psi_A-\Psi_B))@f$.  The error @f$ \delta\chi@f$ is
173     estimated as the largest step size in the minimisation. 
174
175     The total error on the correction is then 
176     @f[ 
177     \delta R_k = R_k'(\chi) \delta\chi
178     @f] 
179 */
180
181 class AliFMDFlowResolutionStar : public AliFMDFlowResolution
182 {
183 public:
184   /** Constructor
185       @param n Harmonic order */
186   AliFMDFlowResolutionStar(UShort_t n=0) : AliFMDFlowResolution(n) {}
187   /** Destructor */
188   virtual ~AliFMDFlowResolutionStar() {}
189
190   /** @{ 
191       @name Information */
192   /** Get the correction for harmonic strength of order @a k 
193       @param k The harminic strenght order to get the correction for
194       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
195   virtual Double_t Correction(UShort_t k) const;
196   /** Get the correction for harmonic strength of order @a k 
197       @param k  The harminic strenght order to get the correction for
198       @param e2 The square error on the correction 
199       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
200   virtual Double_t Correction(UShort_t k, Double_t& e2) const;
201   /** Get @f$ \chi@f$ 
202       @param res  First shot at the resolution. 
203       @param k    Order 
204       @param delta On return, the last step size in @f$ \chi@f$ -
205       which is taken to be @f$ \delta\chi@f$  
206       @return @f$\chi@f$ */
207   virtual Double_t Chi(Double_t res, UShort_t k, Double_t& delta) const;
208   /** @} */
209
210   /** @{ 
211       @name Utilities */
212   /** Draw this corrrection function 
213       @param option Options passed to drawing */
214   virtual void Draw(Option_t* option=""); //*MENU*
215   /** @} */
216 protected:
217   /** Calculate resolution 
218       @param chi @f$ \chi@f$
219       @param k   Order factor 
220       @param dr  On return, the derivative of @f$ R(\chi)@f$
221       @return 
222       @f[ 
223       \frac{\sqrt{\pi/2}}{2}\chi e^{-\chi^2/4}
224       (I_{\frac{(k-1)}{2}}(\chi^2/4)+ I_{\frac{(k+1)}{2}}(\chi^2/4))
225       @f] */
226   Double_t Res(Double_t chi, UShort_t k, Double_t& dr) const;
227   /** Define for ROOT I/O */
228   ClassDef(AliFMDFlowResolutionStar,1);
229 };
230
231 //______________________________________________________
232 /** 
233     @ingroup a_basic
234
235     For more on the event plane angle resolution function, please
236     refer to the class description of ResolutionStar. 
237
238     In this class @f$ \chi@f$ is calculated from the ratio @f$ k/N@f$
239     of events with @f$ |\Psi_A - \Psi_B| > \pi/2@f$ to the total
240     number of events. 
241
242     The pre-print @c nucl-ex/9711003v2 gives the formula 
243     @f[ 
244     \frac{k}{N} = \frac{e^{-\chi^2/2}}{2}
245     @f] 
246     for @f$ \chi@f$.  Note, that this differs from the @f$ \chi@f$
247     used in ResolutionStar by a factor of @f$ 1/sqrt{2}@f$. 
248     
249     We can isolate @f$ \chi = \mp\sqrt{-2\log(2k/N)}@f$ from the
250     equation above.  
251
252     Since @f$ r=k/N@f$ is obviously a efficiency-like ratio, we get
253     that error @f$ \delta r@f$ is given by Binomial errors 
254     @f[ 
255     \delta r =  \sqrt{r\frac{1 - r}{N}}\quad.
256     @f] 
257     The total error @f$ \delta\chi@f$ then becomes 
258     @f[ 
259     \delta^2\chi = \left(\frac{d\chi}{dr}\right)^2 \delta^2r 
260     = \frac{r - 1}{4 k log(2 r)}
261     @f]
262
263     The total error on the correction is, as in
264     ResolutionStar, then given by 
265     @f[ 
266     \delta R_k = R_k'(\chi) \delta\chi
267     @f]
268 */
269 class AliFMDFlowResolutionTDR : public AliFMDFlowResolution
270 {
271 public:
272   /** Constructor
273       @param n Harmonic order */
274   AliFMDFlowResolutionTDR(UShort_t n=0); 
275   /** DTOR */
276   ~AliFMDFlowResolutionTDR() {}
277   /** Copy constructor 
278       @param o  Object to copy from */ 
279   AliFMDFlowResolutionTDR(const AliFMDFlowResolutionTDR& o);
280   /** Assignment operator 
281       @param o  Object to assign from 
282       @return reference to this */ 
283   AliFMDFlowResolutionTDR& operator=(const AliFMDFlowResolutionTDR& o);
284
285   /** @{ 
286       @name Processing */
287   virtual void Clear(Option_t* option="");
288   /** add a data  point */
289   virtual void Add(Double_t psiA, Double_t psiB);
290   /** @} */
291
292   /** @{ 
293       @name Information */
294   /** Get the correction for harmonic strength of order @a k 
295       @param k The harminic strenght order to get the correction for
296       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
297   virtual Double_t Correction(UShort_t k) const;
298   /** Get the correction for harmonic strength of order @a k 
299       @param e2 The square error on the correction 
300       @param k  The harminic strenght order to get the correction for
301       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
302   virtual Double_t Correction(UShort_t k, Double_t& e2) const;
303   /** Get @f$ \chi^2/2@f$ 
304       @param e2 The square error on the correction 
305       @return @f$ \chi^2/2@f$ */
306   virtual Double_t Chi2Over2(Double_t r, Double_t& e2) const;
307   /** @} */
308
309   /** @{ 
310       @name Utilities */
311   /** Draw this corrrection function 
312       @param option Options passed to drawing */
313   virtual void Draw(Option_t* option=""); //*MENU*
314   /** Look through this object */ 
315   virtual void Browse(TBrowser* b); 
316   /** @} */
317
318   /** @{ 
319       @name histograms */ 
320   /** @return NK histogram */
321   const TH1& NKHistogram() const { return fNK; }
322   /** @} */
323 protected:
324   /** Calculate resolution 
325       @param k    Order factor 
326       @param y    @f$ \chi^2/2@f$
327       @param echi @f$\delta\chi@f$ 
328       @param dr   On return, the derivative of @f$ R(\chi)@f$
329       @return 
330       @f[ 
331       \frac{\sqrt{\pi/2}}{2}\chi e^{-\chi^2/2}
332       (I_{\frac{(k-1)}{2}}(\chi^2/2)+ I_{\frac{(k+1)}{2}}(\chi^2/2))
333       @f] */
334   Double_t Res(UShort_t k, Double_t y, Double_t echi2, Double_t& e2) const;
335   /** Number of events with large diviation */
336   ULong_t fLarge; // Number of events with large angle 
337   /** Histogram of large and all angles */ 
338   TH1D fNK;
339   /** Define for ROOT I/O */
340   ClassDef(AliFMDFlowResolutionTDR,1);
341 };
342
343
344 #endif
345 //
346 // EOF
347 //