[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowEfficiency.cxx

/* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1 of
 * the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA
 */
#include <flow/AliFMDFlowEfficiency.h>
#include <cmath>
#include <iostream>
template <typename T>
T sign(const T& x, const T& y) 
{
  return (y >= 0 ? fabs(x) : -fabs(x));
}

//____________________________________________________________________
Double_t
AliFMDFlowEfficiency::LnGamma(Double_t z)
{
  if (z <= 0) return 0;
  Double_t c[] = { 2.5066282746310005, 76.18009172947146, 
		   -86.50532032941677, 24.01409824083091,  
		   -1.231739572450155, 0.1208650973866179e-2, 
		   -0.5395239384953e-5};
  Double_t x = z;
  Double_t y = x;
  Double_t t = x + 5.5;
  t        = (x+.5 * log(t) - t);
  Double_t s = 1.000000000190015;
  for (UInt_t i = 1; i < 7; i++) { 
    y  += 1;
    s  += c[i] / y;
  }
  Double_t v = t + log(c[0]*s/x);
  return v;
}

//____________________________________________________________________
Double_t
AliFMDFlowEfficiency::BetaCf(Double_t x, Double_t a, Double_t b)
{
  Int_t    itmax = 500;   // Maximum # of iterations
  Double_t eps   = 3e-14; // Precision. 
  Double_t fpmin = 1e-30;
  
  Double_t qab = a + b;
  Double_t qap = a + 1;
  Double_t qam = a - 1;
  Double_t c   = 1;
  Double_t d   = 1 - qab * x  / qap;
  if (fabs(d) < fpmin) d = fpmin;
  d          = 1/d;
  Double_t h   = d;
  Int_t m;
  for (m = 1; m <= itmax; m++) { 
    Int_t    m2 =  m * 2;
    Double_t aa =  m * (b - m) * x / ((qam + m2) * (a + m2));
    d         =  1 + aa * d;
    if (fabs(d) < fpmin) d = fpmin;
    c         =  1 + aa / c;
    if (fabs(c) < fpmin) d = fpmin;
    d         =  1 / d;
    h         *= d * c;
    aa        =  -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
    d         = 1 + aa * d;
    if (fabs(d) < fpmin) d = fpmin;
    c         =  1 + aa / c;
    if (fabs(c) < fpmin) d = fpmin;
    d         =  1 / d;
    Double_t dd = d * c;
    h         *= dd;
    if (fabs(dd - 1) <= eps) break;
  }
  if (m > itmax) { 
    std::cerr << "a or b too big, or max # iterations too small"
	      << std::endl;
  }
  return h;
}

  

//____________________________________________________________________
Double_t
AliFMDFlowEfficiency::IBetaI(Double_t a, Double_t b, Double_t x)
{
  Double_t bt = 0;
  if (x < 0 || x > 1) { 
    std::cerr << "x not in [0,1]" << std::endl;
    return 0;
  }
  if (x == 0 || x == 1) bt = 0;
  else 
    bt = (exp(LnGamma(a + b) - LnGamma(a) - LnGamma(b) 
	      + a * log(x) + b * log(1 - x)));
  if (x < (a + 1) / (a + b + 2)) 
    return bt * BetaCf(x, a, b) / a;
  return 1 - bt * BetaCf(1-x, b, a)/ b;
}

//____________________________________________________________________
Double_t
AliFMDFlowEfficiency::BetaAB(Double_t a, Double_t b, Int_t k, Int_t n)
{
  if (a == b) return 0;
  Int_t c1 = k + 1;
  Int_t c2 = n - k + 1;
  return IBetaI(c1, c2, b) - IBetaI(c1, c2, a);
}

//____________________________________________________________________
Double_t
AliFMDFlowEfficiency::SearchUpper(Double_t low, Int_t k, Int_t n, Double_t c)
{
  Double_t integral = BetaAB(low, 1, k, n);
  if (integral == c) return 1; // lucky -- this is the solution
  if (integral <  c) return -1;// no solution exists
  Double_t tooHigh = 1;         // upper edge estimate
  Double_t tooLow  = low;
  Double_t test    = 0;

  // Use a bracket-and-bisect search.  Loop 20 times, to end up with a
  // root guaranteed accurate to better than 0.1%.  
  for (UInt_t i = 0; i < 20; i++) { 
    test     = 0.5 * (tooLow + tooHigh);
    integral = BetaAB(low, test, k, n);
    if (integral > c) tooHigh = test;
    else              tooLow  = test;
  }
  return test;
}
//____________________________________________________________________
Double_t
AliFMDFlowEfficiency::SearchLower(Double_t high, Int_t k, Int_t n, Double_t c)
{
  Double_t integral = BetaAB(0, high, k, n);
  if (integral == c) return 0; // lucky -- this is the solution
  if (integral <  c) return -1;// no solution exists
  Double_t tooLow  = 0;         // lower edge estimate
  Double_t tooHigh = high;
  Double_t test     = 0;

  // Use a bracket-and-bisect search.  Loop 20 times, to end up with a
  // root guaranteed accurate to better than 0.1%.  
  for (UInt_t i = 0; i < 20; i++) { 
    test     = 0.5 * (tooHigh + tooLow);
    integral = BetaAB(test, high, k, n);
    if (integral > c) tooLow  = test;
    else              tooHigh = test;
  }
  return test;
}

//____________________________________________________________________
Double_t 
AliFMDFlowEfficiency::Brent(Double_t ax, Double_t bx, Double_t cx, 
			    Double_t& xmin, Int_t k, Int_t n, Double_t c, 
			    Double_t tol)
{
  const Int_t    kItMax = 100;
  const Double_t kGold  = 0.3819660;
  const Double_t kEps   = 1e-10;
  
  Int_t iter;
  Double_t a  = (ax < cx ? ax : cx);
  Double_t b  = (ax > cx ? ax : cx);
  Double_t d  = 0;
  Double_t e  = 0;
  Double_t x  = bx;
  Double_t w  = bx;
  Double_t v  = bx;
  Double_t fw = Length(x, k, n, c);
  Double_t fv = fw;
  Double_t fx = fv;
  
  for (iter = 1; iter <= kItMax; iter++) { 
    Double_t xm   = .5 * (a + b);
    Double_t tol1 = tol * fabs(x) + kEps;
    Double_t tol2 = 2 * tol1;
    if (fabs(x - xm) < (tol2 - .5 * (b-a))) {
      xmin = x;
      return fx;
    }
    if (fabs(e) > tol1) { 
      Double_t r = (x - w) * (fx - fv);
      Double_t q = (x - v) * (fx - fw);
      Double_t p = (x - v) * q - (x - w) * r;
      q        = 2 * (q - r);
      if (q > 0) p *= -1;
      q        = fabs(q);
      Double_t t = e;
      e        = d;
      if (fabs(p) >= fabs(.5 * q * t) || 
	  p <= q * (a - x)            || 
	  p >= q * (b - x)) {
	e = (x > xm ? a - x : b - x);
	d = kGold * e;
      }
      else {
	d        = p / q;
	Double_t u = x + d;
	if (u - a < tol2 || b - u < tol2) 
	  d = sign(tol1, xm - x);
      }
    }
    else { 
      e = (x >= xm ? a - x : b - x);
      d = kGold * e;
    }
    Double_t u  = (fabs(d) >= tol1 ? x + d : x + sign(tol1, d));
    Double_t fu = Length(u, k, n, c);
    if (fu <= fx) { 
      if (u >= x) a = x;
      else        b = x;
      v  = w;
      w  = x;
      x  = u;
      fv = fw;
      fw = fx;
      fx = fu;
    }
    else { 
      if (u < x) a = u;
      else       b = u;
      if (fu <= fw || w == x) { 
	v  = w;
	w  = u;
	fv = fw;
	fw = fu;
      }
      else if (fu <= fv || v == x || v == w) { 
	v  = u;
	fv = fu;
      }
    }
  }
  std::cerr << "Brett: too many iterations" << std::endl;
  xmin = x;
  return fx;
}
  
//____________________________________________________________________
Double_t 
AliFMDFlowEfficiency::Length(Double_t l, Int_t k, Int_t n, Double_t c)
{
  Double_t h = SearchUpper(l, k, n, c);
  if (h == -1) return 2;
  return (h-l);
}

//____________________________________________________________________
Double_t 
AliFMDFlowEfficiency::Interval(Int_t k, Int_t n, 
			       Double_t& low, Double_t& high, 
			       Double_t c)
{
  if (n == 0) { 
    low  = 0;
    high = 1;
    return .5;
  }
  
  Double_t e = Double_t(k) / n;
  Double_t l, h;
  
  if (k == 0) { 
    l = 0;
    h = SearchUpper(l, k, n, c);
  }
  else if (k == n) { 
    h = 1;
    l = SearchLower(h, k, n, c);
  }
  else { 
    Brent(0, .5, 1, l, n, k, 1e-9, c);
    h = l + Length(l, n, k, c);
  }
  low  = l;
  high = h;
  return e;
}

//____________________________________________________________________
//
// EOF
//
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;
97e94238	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
97e94238	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;
97e94238	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