Added drawing capabilities
[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowResolution.h
1 // -*- mode: C++ -*-
2 /** @file 
3     @brief Declaration of an Resolution class */
4 #ifndef FLOW_RESOLUTION_H
5 #define FLOW_RESOLUTION_H
6 #include <flow/AliFMDFlowStat.h>
7
8 //______________________________________________________
9 /** @class AliFMDFlowResolution flow/AliFMDFlowResolution.h <flow/AliFMDFlowResolution.h>
10     @brief Class to calculate the event plane resolution based on two
11     sub-events 
12     @ingroup a_basic
13     
14     This class calculates the event plane resolution based on the
15     basic formulas given in Phys. Rev. @b C58, 1671.   That is, the
16     resolution is given by 
17     @f[ 
18     R_{k} = \langle\cos(km(\Psi_m-\Psi_R))\rangle 
19     @f] 
20     where @f$ \Psi_R@f$ is the unknown @e true event plane angle, and
21     @f$ n = km@f$
22
23     Using two random sub-events, @f$ A, B@f$ we get that 
24     @f[
25     \langle\cos(n(\Psi^A_m-\Psi^B_m))\rangle = 
26     \langle\cos(n(\Psi^A_m-\Psi_R))\rangle 
27     \langle\cos(n(\Psi^B_m-\Psi_R))\rangle 
28     @f]
29     If the sub-events are of equal size, and randomly chosen, then we
30     get that 
31     @f[ 
32     \langle\cos(n(\Psi^A_m-\Psi_R))\rangle = 
33     \sqrt{\langle\cos(n(\Psi^A_m-\Psi^B_m))}
34     @f]
35     and it follows that 
36     @f{eqnarray*}
37     \langle\cos(km(\Psi_m-\Psi_R))\rangle & = & 
38     \sqrt{2}\langle\cos(n(\Psi^A_m-\Psi_R))\rangle\\
39     \sqrt{2}\sqrt{\langle\cos(n(\Psi^A_m-\Psi^B_m))}
40     @f}
41
42     Hence, the event-plane resolution is simply the square root of the
43     scaled average distance between the two sub-events, multiplied by
44     @f$ \sqrt{s}@f$ 
45
46     The error is therefor @f$ \sqrt{s}@f$ times the variance of the
47     cosine of the distance between the two sub-events. 
48 */
49 class AliFMDFlowResolution : public AliFMDFlowStat
50 {
51 public:
52   /** Constructor
53       @param n Harmonic order */
54   AliFMDFlowResolution(UShort_t n) : fOrder(n) {}
55   /** Destructor */
56   virtual ~AliFMDFlowResolution() {}
57   /** add data point 
58       @param psiA A sub-event plane angle @f$ \Psi_A \in[0,2\pi]@f$
59       @param psiB B sub-event plane angle @f$ \Psi_B \in[0,2\pi]@f$ */
60   virtual void Add(Double_t psiA, Double_t psiB);
61   /** Get the correction for harmonic strength of order @a k 
62       @param k  The harminic strenght order to get the correction for
63       @param e2 The square error on the correction 
64       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
65   virtual Double_t Correction(UShort_t k, Double_t& e2) const;
66   /** Get the correction for harmonic strength of order @a k 
67       @param k The harminic strenght order to get the correction for
68       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
69   virtual Double_t Correction(UShort_t k) const;
70   /** Get the harmnic order */
71   UShort_t Order() const { return fOrder; }
72   /** Draw this corrrection function 
73       @param option Options passed to drawing */
74   virtual void Draw(Option_t* option=""); //*MENU*
75 protected:
76   /** Order */
77   UShort_t fOrder;
78   /** Define for ROOT I/O */
79   ClassDef(AliFMDFlowResolution,1);
80 };
81
82 //______________________________________________________
83 /** 
84     @ingroup a_basic
85     The event plane angle resolution function is 
86     @f[
87     R_k(\chi) = \frac{\pi} \chi e^{-\chi^2/4}\left( 
88     I_{\frac{k-1}2}(\chi^2/4) + I_{\frac{k+1}2}(\chi^2/4)\right)
89     @f]
90     Where @f$ I_n(x)@f$ is the modified Bessel function of the first
91     kind.   Identifying 
92     @f[
93     y = \chi^2/4\quad C=\frac{\sqrt{\pi}e^{-y}}{2\sqrt{2}}
94     @f]
95     and the short hands @f$ f2(y) = I_{\frac{k-1}2}@f$ and @f$ f3(y) =
96     I_{\frac{k+1}2}(y)@f$ we can write this more compact as 
97     @f[ 
98     R_k(y) = C y (f2(y) + f3(y))
99     @f] 
100     The derivative of the resolution function is 
101     @f[ 
102     R_k'(y) = \frac{C}{2}\left[4\sqrt{y}\left(f2'(y)-f3'(y)\right) 
103     - (4 y - 2)\left(f2(y) + f3(y)\right)\right]
104     @f]
105     Since 
106     @f[
107     I_\nu'(x) = I_{\nu-1}(x) - \frac{\nu}{x} I_\nu(x)\quad,
108     @f]
109     and setting @f$ f1(y) = I_{\frac{k-3}2}(y)@f$, we get 
110     @f[ 
111     R_k'(y) = \frac{C}{2}\left[4y f1(y) + (4-2k) f2(y) - (4y+2k)
112     f3(y)\right]
113     @f]
114
115     In this class, the argument @f$ \chi@f$ is estimated by finding
116     the minima of @f$ R_k(\chi)@f$ near the average of
117     @f$\cos(n(\Psi_A-\Psi_B))@f$.  The error @f$ \delta\chi@f$ is
118     estimated as the largest step size in the minimisation. 
119
120     The total error on the correction is then 
121     @f[ 
122     \delta R_k = R_k'(\chi) \delta\chi
123     @f] 
124 */
125
126 class AliFMDFlowResolutionStar : public AliFMDFlowResolution
127 {
128 public:
129   /** Constructor
130       @param n Harmonic order */
131   AliFMDFlowResolutionStar(UShort_t n) 
132     : AliFMDFlowResolution(n) {}
133   /** Destructor */
134   ~AliFMDFlowResolutionStar() {}
135   /** Get the correction for harmonic strength of order @a k 
136       @param k The harminic strenght order to get the correction for
137       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
138   virtual Double_t Correction(UShort_t k) const;
139
140   /** Get the correction for harmonic strength of order @a k 
141       @param k  The harminic strenght order to get the correction for
142       @param e2 The square error on the correction 
143       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
144   virtual Double_t Correction(UShort_t k, Double_t& e2) const;
145   /** Get @f$ \chi@f$ 
146       @param res  First shot at the resolution. 
147       @param k    Order 
148       @param delta On return, the last step size in @f$ \chi@f$ -
149       which is taken to be @f$ \delta\chi@f$  
150       @return @f$\chi@f$ */
151   virtual Double_t Chi(Double_t res, UShort_t k, Double_t& delta) const;
152   /** Draw this corrrection function 
153       @param option Options passed to drawing */
154   virtual void Draw(Option_t* option=""); //*MENU*
155 protected:
156   /** Calculate resolution 
157       @param chi @f$ \chi@f$
158       @param k   Order factor 
159       @param dr  On return, the derivative of @f$ R(\chi)@f$
160       @return 
161       @f[ 
162       \frac{\sqrt{\pi/2}}{2}\chi e^{-\chi^2/4}
163       (I_{\frac{(k-1)}{2}}(\chi^2/4)+ I_{\frac{(k+1)}{2}}(\chi^2/4))
164       @f] */
165   Double_t Res(Double_t chi, UShort_t k, Double_t& dr) const;
166   /** Define for ROOT I/O */
167   ClassDef(AliFMDFlowResolutionStar,1);
168 };
169
170 //______________________________________________________
171 /** 
172     @ingroup a_basic
173
174     For more on the event plane angle resolution function, please
175     refer to the class description of ResolutionStar. 
176
177     In this class @f$ \chi@f$ is calculated from the ratio @f$ k/N@f$
178     of events with @f$ |\Psi_A - \Psi_B| > \pi/2@f$ to the total
179     number of events. 
180
181     The pre-print @c nucl-ex/9711003v2 gives the formula 
182     @f[ 
183     \frac{k}{N} = \frac{e^{-\chi^2/2}}{2}
184     @f] 
185     for @f$ \chi@f$.  Note, that this differs from the @f$ \chi@f$
186     used in ResolutionStar by a factor of @f$ 1/sqrt{2}@f$. 
187     
188     We can isolate @f$ \chi = \mp\sqrt{-2\log(2k/N)}@f$ from the
189     equation above.  
190
191     Since @f$ r=k/N@f$ is obviously a efficiency-like ratio, we get
192     that error @f$ \delta r@f$ is given by Binomial errors 
193     @f[ 
194     \delta r =  \sqrt{r\frac{1 - r}{N}}\quad.
195     @f] 
196     The total error @f$ \delta\chi@f$ then becomes 
197     @f[ 
198     \delta^2\chi = \left(\frac{d\chi}{dr}\right)^2 \delta^2r 
199     = \frac{r - 1}{4 k log(2 r)}
200     @f]
201
202     The total error on the correction is, as in
203     ResolutionStar, then given by 
204     @f[ 
205     \delta R_k = R_k'(\chi) \delta\chi
206     @f]
207 */
208 class AliFMDFlowResolutionTDR : public AliFMDFlowResolution
209 {
210 public:
211   /** Constructor
212       @param n Harmonic order */
213   AliFMDFlowResolutionTDR(UShort_t n) 
214     : AliFMDFlowResolution(n), fLarge(0) {}
215   /** DTOR */
216   ~AliFMDFlowResolutionTDR() {}
217   virtual void Clear(Option_t* option="");
218   /** add a data  point */
219   virtual void Add(Double_t psi_a, Double_t psi_b);
220   /** Get the correction for harmonic strength of order @a k 
221       @param k The harminic strenght order to get the correction for
222       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
223   virtual Double_t Correction(UShort_t k) const;
224   /** Get the correction for harmonic strength of order @a k 
225       @param e2 The square error on the correction 
226       @param k  The harminic strenght order to get the correction for
227       @return @f$ \langle\cos(n(\psi_n - \psi_R))\rangle@f$ */ 
228   virtual Double_t Correction(UShort_t k, Double_t& e2) const;
229   /** Get @f$ \chi^2/2@f$ 
230       @param e2 The square error on the correction 
231       @return @f$ \chi^2/2@f$ */
232   virtual Double_t Chi2Over2(Double_t r, Double_t& e2) const;
233   /** Draw this corrrection function 
234       @param option Options passed to drawing */
235   virtual void Draw(Option_t* option=""); //*MENU*
236 protected:
237   /** Calculate resolution 
238       @param k    Order factor 
239       @param y    @f$ \chi^2/2@f$
240       @param echi @f$\delta\chi@f$ 
241       @param dr   On return, the derivative of @f$ R(\chi)@f$
242       @return 
243       @f[ 
244       \frac{\sqrt{\pi/2}}{2}\chi e^{-\chi^2/2}
245       (I_{\frac{(k-1)}{2}}(\chi^2/2)+ I_{\frac{(k+1)}{2}}(\chi^2/2))
246       @f] */
247   Double_t Res(UShort_t k, Double_t y, Double_t echi2, Double_t& e2) const;
248   /** Number of events with large diviation */
249   ULong_t fLarge;
250   /** Define for ROOT I/O */
251   ClassDef(AliFMDFlowResolutionTDR,1);
252 };
253
254
255 #endif
256 //
257 // EOF
258 //