]>
Commit | Line | Data |
---|---|---|
97e94238 | 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 | */ | |
39eefe19 | 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 | |
97e94238 | 129 | Double_t tooHigh = 1; // upper edge estimate |
130 | Double_t tooLow = low; | |
131 | Double_t test = 0; | |
39eefe19 | 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++) { | |
97e94238 | 136 | test = 0.5 * (tooLow + tooHigh); |
39eefe19 | 137 | integral = BetaAB(low, test, k, n); |
97e94238 | 138 | if (integral > c) tooHigh = test; |
139 | else tooLow = test; | |
39eefe19 | 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 | |
97e94238 | 150 | Double_t tooLow = 0; // lower edge estimate |
151 | Double_t tooHigh = high; | |
39eefe19 | 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++) { | |
97e94238 | 157 | test = 0.5 * (tooHigh + tooLow); |
39eefe19 | 158 | integral = BetaAB(test, high, k, n); |
97e94238 | 159 | if (integral > c) tooLow = test; |
160 | else tooHigh = test; | |
39eefe19 | 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 | { | |
97e94238 | 171 | const Int_t kItMax = 100; |
172 | const Double_t kGold = 0.3819660; | |
173 | const Double_t kEps = 1e-10; | |
39eefe19 | 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 | ||
97e94238 | 187 | for (iter = 1; iter <= kItMax; iter++) { |
39eefe19 | 188 | Double_t xm = .5 * (a + b); |
97e94238 | 189 | Double_t tol1 = tol * fabs(x) + kEps; |
39eefe19 | 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); | |
97e94238 | 208 | d = kGold * e; |
39eefe19 | 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); | |
97e94238 | 219 | d = kGold * e; |
39eefe19 | 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 |