FMD/flow/AliFMDFlowEfficiency.cxx

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