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