bug fix
[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowEfficiency.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 #include <flow/AliFMDFlowEfficiency.h>
19 #include <cmath>
20 #include <iostream>
21 template <typename T>
22 T sign(const T& x, const T& y) 
23 {
24   return (y >= 0 ? fabs(x) : -fabs(x));
25 }
26
27 //____________________________________________________________________
28 Double_t
29 AliFMDFlowEfficiency::LnGamma(Double_t z)
30 {
31   if (z <= 0) return 0;
32   Double_t c[] = { 2.5066282746310005, 76.18009172947146, 
33                    -86.50532032941677, 24.01409824083091,  
34                    -1.231739572450155, 0.1208650973866179e-2, 
35                    -0.5395239384953e-5};
36   Double_t x = z;
37   Double_t y = x;
38   Double_t t = x + 5.5;
39   t        = (x+.5 * log(t) - t);
40   Double_t s = 1.000000000190015;
41   for (UInt_t i = 1; i < 7; i++) { 
42     y  += 1;
43     s  += c[i] / y;
44   }
45   Double_t v = t + log(c[0]*s/x);
46   return v;
47 }
48
49 //____________________________________________________________________
50 Double_t
51 AliFMDFlowEfficiency::BetaCf(Double_t x, Double_t a, Double_t b)
52 {
53   Int_t    itmax = 500;   // Maximum # of iterations
54   Double_t eps   = 3e-14; // Precision. 
55   Double_t fpmin = 1e-30;
56   
57   Double_t qab = a + b;
58   Double_t qap = a + 1;
59   Double_t qam = a - 1;
60   Double_t c   = 1;
61   Double_t d   = 1 - qab * x  / qap;
62   if (fabs(d) < fpmin) d = fpmin;
63   d          = 1/d;
64   Double_t h   = d;
65   Int_t m;
66   for (m = 1; m <= itmax; m++) { 
67     Int_t    m2 =  m * 2;
68     Double_t aa =  m * (b - m) * x / ((qam + m2) * (a + m2));
69     d         =  1 + aa * d;
70     if (fabs(d) < fpmin) d = fpmin;
71     c         =  1 + aa / c;
72     if (fabs(c) < fpmin) d = fpmin;
73     d         =  1 / d;
74     h         *= d * c;
75     aa        =  -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
76     d         = 1 + aa * d;
77     if (fabs(d) < fpmin) d = fpmin;
78     c         =  1 + aa / c;
79     if (fabs(c) < fpmin) d = fpmin;
80     d         =  1 / d;
81     Double_t dd = d * c;
82     h         *= dd;
83     if (fabs(dd - 1) <= eps) break;
84   }
85   if (m > itmax) { 
86     std::cerr << "a or b too big, or max # iterations too small"
87               << std::endl;
88   }
89   return h;
90 }
91
92   
93
94 //____________________________________________________________________
95 Double_t
96 AliFMDFlowEfficiency::IBetaI(Double_t a, Double_t b, Double_t x)
97 {
98   Double_t bt = 0;
99   if (x < 0 || x > 1) { 
100     std::cerr << "x not in [0,1]" << std::endl;
101     return 0;
102   }
103   if (x == 0 || x == 1) bt = 0;
104   else 
105     bt = (exp(LnGamma(a + b) - LnGamma(a) - LnGamma(b) 
106               + a * log(x) + b * log(1 - x)));
107   if (x < (a + 1) / (a + b + 2)) 
108     return bt * BetaCf(x, a, b) / a;
109   return 1 - bt * BetaCf(1-x, b, a)/ b;
110 }
111
112 //____________________________________________________________________
113 Double_t
114 AliFMDFlowEfficiency::BetaAB(Double_t a, Double_t b, Int_t k, Int_t n)
115 {
116   if (a == b) return 0;
117   Int_t c1 = k + 1;
118   Int_t c2 = n - k + 1;
119   return IBetaI(c1, c2, b) - IBetaI(c1, c2, a);
120 }
121
122 //____________________________________________________________________
123 Double_t
124 AliFMDFlowEfficiency::SearchUpper(Double_t low, Int_t k, Int_t n, Double_t c)
125 {
126   Double_t integral = BetaAB(low, 1, k, n);
127   if (integral == c) return 1; // lucky -- this is the solution
128   if (integral <  c) return -1;// no solution exists
129   Double_t tooHigh = 1;         // upper edge estimate
130   Double_t tooLow  = low;
131   Double_t test    = 0;
132
133   // Use a bracket-and-bisect search.  Loop 20 times, to end up with a
134   // root guaranteed accurate to better than 0.1%.  
135   for (UInt_t i = 0; i < 20; i++) { 
136     test     = 0.5 * (tooLow + tooHigh);
137     integral = BetaAB(low, test, k, n);
138     if (integral > c) tooHigh = test;
139     else              tooLow  = test;
140   }
141   return test;
142 }
143 //____________________________________________________________________
144 Double_t
145 AliFMDFlowEfficiency::SearchLower(Double_t high, Int_t k, Int_t n, Double_t c)
146 {
147   Double_t integral = BetaAB(0, high, k, n);
148   if (integral == c) return 0; // lucky -- this is the solution
149   if (integral <  c) return -1;// no solution exists
150   Double_t tooLow  = 0;         // lower edge estimate
151   Double_t tooHigh = high;
152   Double_t test     = 0;
153
154   // Use a bracket-and-bisect search.  Loop 20 times, to end up with a
155   // root guaranteed accurate to better than 0.1%.  
156   for (UInt_t i = 0; i < 20; i++) { 
157     test     = 0.5 * (tooHigh + tooLow);
158     integral = BetaAB(test, high, k, n);
159     if (integral > c) tooLow  = test;
160     else              tooHigh = test;
161   }
162   return test;
163 }
164
165 //____________________________________________________________________
166 Double_t 
167 AliFMDFlowEfficiency::Brent(Double_t ax, Double_t bx, Double_t cx, 
168                             Double_t& xmin, Int_t k, Int_t n, Double_t c, 
169                             Double_t tol)
170 {
171   const Int_t    kItMax = 100;
172   const Double_t kGold  = 0.3819660;
173   const Double_t kEps   = 1e-10;
174   
175   Int_t iter;
176   Double_t a  = (ax < cx ? ax : cx);
177   Double_t b  = (ax > cx ? ax : cx);
178   Double_t d  = 0;
179   Double_t e  = 0;
180   Double_t x  = bx;
181   Double_t w  = bx;
182   Double_t v  = bx;
183   Double_t fw = Length(x, k, n, c);
184   Double_t fv = fw;
185   Double_t fx = fv;
186   
187   for (iter = 1; iter <= kItMax; iter++) { 
188     Double_t xm   = .5 * (a + b);
189     Double_t tol1 = tol * fabs(x) + kEps;
190     Double_t tol2 = 2 * tol1;
191     if (fabs(x - xm) < (tol2 - .5 * (b-a))) {
192       xmin = x;
193       return fx;
194     }
195     if (fabs(e) > tol1) { 
196       Double_t r = (x - w) * (fx - fv);
197       Double_t q = (x - v) * (fx - fw);
198       Double_t p = (x - v) * q - (x - w) * r;
199       q        = 2 * (q - r);
200       if (q > 0) p *= -1;
201       q        = fabs(q);
202       Double_t t = e;
203       e        = d;
204       if (fabs(p) >= fabs(.5 * q * t) || 
205           p <= q * (a - x)            || 
206           p >= q * (b - x)) {
207         e = (x > xm ? a - x : b - x);
208         d = kGold * e;
209       }
210       else {
211         d        = p / q;
212         Double_t u = x + d;
213         if (u - a < tol2 || b - u < tol2) 
214           d = sign(tol1, xm - x);
215       }
216     }
217     else { 
218       e = (x >= xm ? a - x : b - x);
219       d = kGold * e;
220     }
221     Double_t u  = (fabs(d) >= tol1 ? x + d : x + sign(tol1, d));
222     Double_t fu = Length(u, k, n, c);
223     if (fu <= fx) { 
224       if (u >= x) a = x;
225       else        b = x;
226       v  = w;
227       w  = x;
228       x  = u;
229       fv = fw;
230       fw = fx;
231       fx = fu;
232     }
233     else { 
234       if (u < x) a = u;
235       else       b = u;
236       if (fu <= fw || w == x) { 
237         v  = w;
238         w  = u;
239         fv = fw;
240         fw = fu;
241       }
242       else if (fu <= fv || v == x || v == w) { 
243         v  = u;
244         fv = fu;
245       }
246     }
247   }
248   std::cerr << "Brett: too many iterations" << std::endl;
249   xmin = x;
250   return fx;
251 }
252   
253 //____________________________________________________________________
254 Double_t 
255 AliFMDFlowEfficiency::Length(Double_t l, Int_t k, Int_t n, Double_t c)
256 {
257   Double_t h = SearchUpper(l, k, n, c);
258   if (h == -1) return 2;
259   return (h-l);
260 }
261
262 //____________________________________________________________________
263 Double_t 
264 AliFMDFlowEfficiency::Interval(Int_t k, Int_t n, 
265                                Double_t& low, Double_t& high, 
266                                Double_t c)
267 {
268   if (n == 0) { 
269     low  = 0;
270     high = 1;
271     return .5;
272   }
273   
274   Double_t e = Double_t(k) / n;
275   Double_t l, h;
276   
277   if (k == 0) { 
278     l = 0;
279     h = SearchUpper(l, k, n, c);
280   }
281   else if (k == n) { 
282     h = 1;
283     l = SearchLower(h, k, n, c);
284   }
285   else { 
286     Brent(0, .5, 1, l, n, k, 1e-9, c);
287     h = l + Length(l, n, k, c);
288   }
289   low  = l;
290   high = h;
291   return e;
292 }
293
294 //____________________________________________________________________
295 //
296 // EOF
297 //
298
299     
300   
301