#include "MyBayesUL.h"
#if !defined(__CINT__) || defined(__MAKECINT__)
  
ClassImp(MyBayesUL)
   ; 

#include "TMath.h"
#include "Riostream.h"

#endif

MyBayesUL::MyBayesUL()
{
}

MyBayesUL::~MyBayesUL()
{
}

Double_t
MyBayesUL::GetClassicalUL(Int_t n_obs, Double_t n_exp,
			  Double_t cl)
{
   Double_t n_ul = 0;

   if (n_obs == 0)
      n_ul = -TMath::Log(1-cl) - n_exp;
   else {
      n_ul = 0.5*chisin(cl, 2*(n_obs+1)) - n_exp;
   }

   return n_ul;
}

Double_t
MyBayesUL::GetUL(Int_t n_obs, Double_t n_exp,
		 Double_t cl)
{
   Double_t n_ul = 0;

   if (n_obs == 0)
      n_ul = -TMath::Log(1-cl);
   else {
      Double_t beta = 1 - cl;
#if 0
      Double_t p = 0;
      for (Int_t n=0; n<=n_obs; ++n) {
	 p += TMath::Poisson(n, n_exp);
      }
#else
      Double_t p = 1 - TMath::Gamma(n_obs+1, n_exp);
#endif
      Double_t beta2 = beta * p;
      n_ul = 0.5*chisin(1-beta2, 2*(n_obs+1)) - n_exp;
      if (n_ul < -TMath::Log(beta))
	 n_ul = -TMath::Log(beta);
   }

   return n_ul;
}

/*
* Mathlib gen
*
*
C     This corresponds to CHISIN,IF=DOUBLE and CHISIN64,IF=-DOUBLE
      FUNCTION CHISIN(Q,N)

C     Computes the inverse of the
C     chi**2 distribution with n degrees of freedom.
C     Based on R.B. Goldstein,  ALGORITHM 451 Chi-Square Quantiles,
C     Collected Algorithms from CACM
C     Note the complementary definition of the integral!
*/
Double_t
MyBayesUL::chisin(Double_t Q, Int_t N)
{
   const Double_t Z1 = 1;
   const Double_t HF = Z1/2;

   const Double_t C1[] =
      { 1.565326E-3, 1.060438E-3,-6.950356E-3,-1.323293E-2, 2.277679E-2,
	-8.986007E-3,-1.513904E-2};
   const Double_t C2[] =
      {-1.450117E-3, 2.530010E-3, 5.169654E-3,-1.153761E-2, 1.128186E-2,
       2.607083E-2,-2.237368E-1};
   const Double_t C3[] =
     {9.780499E-5,-8.426812E-4, 3.125580E-3,-8.553069E-3, 1.348028E-4,
      4.713941E-1, 1.0000886};
   const Double_t A[] =
      { 1.264616E-2,-1.425296E-2, 1.400483E-2,-5.886090E-3,-1.091214E-2,
	-2.304527E-2, 3.135411E-3,-2.728484E-4,-9.699681E-3, 1.316872E-2,
	2.618914E-2,-2.222222E-1, 5.406674E-5, 3.483789E-5,-7.274761E-4,
	3.292181E-3,-8.729713E-3, 4.714045E-1, 1};

   Double_t H = 0;

   if ((Q < 0) || (Q >= 1)) {
      H=0;
      cout << "CHISIN>> error, Q = " << Q << endl;
   }
   else if (N <= 0) {
      H=0;
      cout << "CHISIN>> error, N = " << Q << endl;
   }
   else if (Q == 0) {
      H=0;
   }
   else if (N == 1) {
      H = TMath::Power(gausin(HF*(1-Q)),2);
   }
   else if (N == 2) {
      H = -2*TMath::Log(1-Q);
   }
   else {
      Double_t F1=Z1/N;
      Double_t T=gausin(Q);
      Double_t F2=TMath::Sqrt(F1)*T;
      if (N < (2*Int_t(4*TMath::Abs(T)))) {
	 Double_t S1=C1[1-1];
	 Double_t S2=C2[1-1];
	 Double_t S3=C3[1-1];
	 for (Int_t i=2; i<=7; ++i) {
	    S1 = C1[i-1]+S1*F2;
	    S2 = C2[i-1]+S2*F2;
	    S3 = C3[i-1]+S3*F2;
	 }
	 H = (S1*F1+S2)*F1+S3;
      }
      else {
	 H =
	    (
	       ((A[1-1]+A[2-1]*F2)*F1
		+
		(((A[3-1]+A[4-1]*F2)*F2+A[5-1])*F2+A[6-1])
		  )*F1
	       +(((((A[7-1]+A[8-1]*F2)*F2+A[9-1])*F2+A[10-1])*F2
		  +A[11-1])*F2+A[12-1])
	       )*F1
	    +
	    (((((A[13-1]*F2+A[14-1])*F2+A[15-1])*F2+A[16-1])*F2+A[17-1])*F2*F2
	      +A[18-1])*F2
	    +A[19-1];

      }

      H = N * TMath::Power(H,3);
   }

   return H;
}

/*
*
* $Id: gausin64.F,v 1.1.1.1 1996/04/01 15:02:43 mclareni Exp $
*
* $Log: gausin64.F,v $
* Revision 1.1.1.1  1996/04/01 15:02:43  mclareni
* Mathlib gen
*
*
      FUNCTION  GAUSIN(P)
      FUNCTION DGAUSN(P)
C     Computes a "Normal Deviate"
C     Based on G.W. Hill & A.W. Davis, Algorithm 442 Normal Deviate
C     Collected Algorithms from CACM
*/
Double_t
MyBayesUL::gausin(Double_t P)
{
   const Double_t C = 2.50662827463100050;
   const Double_t Z1 = 1;
   const Double_t HF = Z1/2;
   const Double_t C1 = 3*Z1/4;
   const Double_t C2 = 7*Z1/8;
   const Double_t C3 = Z1/3;

   Double_t H = 0;

   if ((P <= 0) || (P >= 1)) {
      H = 0;
      cout << "GAUSIN>> error, P = " << P << endl;      
   }
   else if (P == HF) {
      H=0;
   }
   else {
      Double_t X = P;
      if (P > HF) X=1-P;
      X = TMath::Sqrt(-2*TMath::Log(X));
      X = X-((7.47395*X+494.877)*X+1637.720)/
	 (((X+117.9407)*X+908.401)*X+659.935);
      if(P < HF) X=-X;
      Double_t S=X*X;;
      Double_t Z=C*(P-freq(X))*TMath::Exp(HF*S);
      H=(((((C1*S+C2)*Z+X)*X+HF)*C3*Z+HF*X)*Z+1)*Z+X;
   }

   return H;
}

/*
*
* $Id: freq.F,v 1.1.1.1 1996/04/01 15:01:53 mclareni Exp $
*
* $Log: freq.F,v $
* Revision 1.1.1.1  1996/04/01 15:01:53  mclareni
* Mathlib gen
*
*
      FUNCTION FREQ(X)
*/
Double_t
MyBayesUL::freq(Double_t X)
{
   const Double_t Z1 = 1;
   const Double_t HF = Z1/2;
   const Double_t C1 = 0.56418958354775629;
   const Double_t W2 = 1.41421356237309505;
   const Double_t RW2 = 1/W2;
   const Double_t P1[] = {
      2.4266795523053175E+2,
      2.1979261618294152E+1,
      6.9963834886191355E+0,
      -3.5609843701815385E-2};
   const Double_t Q1[] = {
      +2.1505887586986120E+2,
      +9.1164905404514901E+1,
      +1.5082797630407787E+1,
      +1};
   const Double_t P2[] = {
      +3.00459261020161601E+2,
      +4.51918953711872942E+2,
      +3.39320816734343687E+2,
      +1.52989285046940404E+2,
      +4.31622272220567353E+1,
      +7.21175825088309366E+0,
      +5.64195517478973971E-1,
      -1.36864857382716707E-7};
   const Double_t Q2[] = {
      +3.00459260956983293E+2,
      +7.90950925327898027E+2,
      +9.31354094850609621E+2,
      +6.38980264465631167E+2,
      +2.77585444743987643E+2,
      +7.70001529352294730E+1,
      +1.27827273196294235E+1,
      +1};
   const Double_t P3[] = {
      -2.99610707703542174E-3,
      -4.94730910623250734E-2,
      -2.26956593539686930E-1,
      -2.78661308609647788E-1,
      -2.23192459734184686E-2};
   const Double_t Q3[] = {
      +1.06209230528467918E-2,
      +1.91308926107829841E-1,
      +1.05167510706793207E+0,
      +1.98733201817135256E+0,
      +1};

   Double_t H = 0;
   Double_t HC = 0;

   Double_t V = RW2*TMath::Abs(X);
   if (V < HF) {
      Double_t Y=V*V;
      Double_t AP=P1[3];
      Double_t AQ=Q1[3];
      for (Int_t i=2; i>=0; --i) {
	 AP=P1[i]+Y*AP;
	 AQ=Q1[i]+Y*AQ;
      }

      H=V*AP/AQ;
      HC=1-H;
   }
   else if (V < 4) {
      Double_t AP=P2[7];
      Double_t AQ=Q2[7];
      for (Int_t i=6; i>=0; --i) {
	 AP=P2[i]+V*AP;
	 AQ=Q2[i]+V*AQ;
      }
      HC=TMath::Exp(-V*V)*AP/AQ;
      H=1-HC;
   }
   else {
      Double_t Y=1/(V*V);
      Double_t AP=P3[4];
      Double_t AQ=Q3[4];
      for (Int_t I=3; I>=0; --I) {
	 AP=P3[I]+Y*AP;
	 AQ=Q3[I]+Y*AQ;
      }
      HC=TMath::Exp(-V*V)*(C1+Y*AP/AQ)/V;
      H=1-HC;
   }

   return
      (X > 0) ? (HF+HF*H) : (HF*HC);
}


/*
 TRolke impl. ???
*/
/*
classical
 */
Double_t
MyBayesUL::GetUL_TRolke(Double_t n_obs,
			Double_t cl)
{
   Double_t n_ul = 0;

   if (n_obs == 0) {
      n_ul = -TMath::Log(1-cl);
   }
   else {
      n_ul = 0.5*Chi2Percentile(2*(n_obs+1), 1-cl);
   }

   return n_ul;
}


//______________________________________________________________________
/*
  stolen from TRolke class
*/
Double_t //TRolke
MyBayesUL::Chi2Percentile(Double_t df,Double_t CL1)
{
 //
 // Finds the Chi-square argument x such that the integral
 // from x to infinity of the Chi-square density is equal
 // to the given cumulative probability y.
 //
 // This is accomplished using the inverse gamma integral
 // function and the relation
 //    x/2 = igami( df/2, y );
 //
 //   Author Jan Conrad (CERN) Jan.Conrad@cern.ch
 //   Copyright by Stephen L. Moshier (Cephes Math Library)

  return  2*InverseIncompleteGamma(df/2,CL1);
}

//______________________________________________________________________
Double_t //TRolke
MyBayesUL::InverseIncompleteGamma(Double_t df,Double_t CL1)
{
  // calculates the inverse of the incomplete (complemented)gamma integral
  // for df degrees of freedom and fraction CL1
  // * Given p, the function finds x such that
  // Starting with the approximate value
  //          3
  //  x = df t
  //
  // where
  //
  //  t = 1 - d - InverseNormal(p) sqrt(df)
  //
  // and
  //
  // d = 1/9df,
  //
  // the routine performs up to 10 Newton iterations to find the
  // root of 1- TMath::Gamma(df,CL1) - p = 0.
  // Author Jan Conrad (CERN) Jan.Conrad@cern.ch
  // Copyright by Stephen L. Moshier (Cephes Math Library)


  const Double_t MACHEP = 1.11022302462515654042E-12;   /* 2**-53 */
  const Double_t MAXLOG = 8.8029691931113054295988E1;
  const Double_t MAXNUM = 1.701411834604692317316873e38;
  const Double_t dithresh = 5.0*MACHEP;

  // bound the solution
  Double_t x0 = MAXNUM;
  Double_t y1 = 0;
  Double_t x1 = 0;
  Double_t yh = 1.0;

  // approximation to the inverse function
  Double_t d = 1/(9*df);
  Double_t y = ( 1 - d - InverseNormal(CL1) * TMath::Sqrt(d) );
  Double_t x = df * y * y * y;
  Double_t lgm = TMath::LnGamma(df);

  Int_t i,dir;

  for (i = 0; i < 10; i++) {
     if (x > x0 || x < x1) goto ihalve;
     y = 1 - TMath::Gamma(df,x);
     if (y < y1 || y > yh) goto ihalve;
     if (y < CL1) {
	x0 = x;
	y1 = y;
     } else {
	x1 =x;
        yh =y;
     }

     // compute derivative of the function at this point
     d = (df  - 1)*TMath::Log(x) -x - TMath::LnGamma(df);
     if (d < -MAXLOG ) goto ihalve;
     d = -TMath::Exp(d);

     // compute the step to the next approximatioin of x
     d = (y -CL1)/d;
     if (TMath::Abs(d/x) < MACHEP) goto done;
     x = x -d;

     // Resort to interval halving if Newton iteration did not converge

     ihalve:
	d = 0.0625;
        if ( x0 == MAXNUM ) {
           if( x <= 0 ) x = 1;
           while( x0 == MAXNUM ) {
	      x = (1 + d) * x;
	      y = 1 - TMath::Gamma(df,x);
	      if (y < CL1 ) {
	         x0 = x;
	         y1 = y;
	         break;
	      }
	      d = d +d;
           }
        }
        d = 0.5;
  }

  dir = 0;

  for (i=0; i < 400; i++) {
      x = x1 + d * (x0-x1);
      y = 1 - TMath::Gamma(df,x);

      lgm = (x0-x1)/(x1+x0);
      if (TMath::Abs(lgm) < dithresh ) break;

      lgm = (y-CL1)/CL1;
      if (TMath::Abs(lgm) < dithresh ) break;
      if (y >= CL1) {
	  x1 = x;
	  yh = y;
	  if (dir < 0) {
	      dir = 0;
	      d   = 0.5;
	  }
	  else if(dir > 0)
	    d = 0.5*d + 0.5;
	  else
	    d = (CL1-y1)/(yh-y1);
	  dir += 1;
      } else {
	  x0 = x;
	  CL1 = y;
	  if (dir > 0 ) {
	      dir = 0;
	      d   = 0.5;
	  }
	  else if(dir < -1)
	    d = 0.5* d;
	  else
	    d = (CL1-y1)/(yh-y1);
	  dir -= 1;
      }

   }

   if ( x == 0) {
      cout << "ERROR" << endl;
   }

   done:
   return x;
}


//___________________________________________________________________
Double_t //TRolke
MyBayesUL::InverseNormal(Double_t CL1)
{

  // Returns the argument, x, for which the area under the
  // Gaussian probability density function (integrated from
  // minus infinity to x) is equal to y.
  //
  // For small arguments 0 < y < exp(-2), the program computes
  // z = sqrt( -2.0 * log(y) );  then the approximation is
  // x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
  // There are two rational functions P/Q, one for 0 < y < exp(-32)
  // and the other for y up to exp(-2).  For larger arguments,
  // w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
  //
  //   Author Jan Conrad (CERN) Jan.Conrad@cern.ch
  //   Copyright by Stephen L. Moshier (Cephes Math Library)

  const Double_t MAXNUM = 1.7976931348623158E308;
  const Double_t s2pi   = 2.50662827463100050242E0;


  const Int_t P0[20] = {5165,3678,64335,49229,60889,40609,32785,16472,56250,34822,22160,49228,49661,15319,56526,16427,51838,36598,54175,49139};

  const Int_t Q0[32] = {29945,53792,17813,16383,58806,33815,46210,16402,33236,13579,38670,16469,13919,22210,11982,49260,52223,43241,2131,16489,47065,59277,33377,49236,30051,49412,53165,16431,52693,64191,61148,49138};

  const Int_t P1[36] = {36571,39484,14535,16400,31106,37589,34413,16447,32703,17062,38103,16460,7081,4334,2638,16454,17982,37887,24242,16429,49057,31340,32313,16385,38519,17549,62441,49089,16839,61027,61717,49057,59279,27851,6366,48972};

  const Int_t Q1[32] = {54299,41275,36698,16431,10607,35310,45572,16454,17948,8844,43162,16452,55365,40327,5575,16430,47648,43068,2437,16388,562,52699,13068,49090,45853,17773,32590,49059,1566,31101,38079,48974};

  const Int_t P2[36] = {54644,59283,59112,16393,30731,49959,43313,16411,45228,62251,33454,16399,39436,6375,21532,16373,9809,62302,51781,16329,13645,36969,22096,16265,6162,59598,50098,16179,8756,19468,19497,16070,36857,12376,52396,15930};

  const Int_t Q2[32] = {59432,22155,6362,16408,14358,43730,28749,16397,25088,10927,2119,16374,15515,24534,44455,16331,51083,44470,31783,16267,10499,26111,32555,16181,52471,62454,17292,16072,27453,47222,10725,15933};


  Double_t x, y, z, y2, x0, x1;
  Int_t code;

  if( CL1 <= 0.0 ) return -MAXNUM;
  if( CL1 >= 1.0 ) return MAXNUM;
  code = 1;
  y = CL1;
  if (y > (1.0 - 0.13533528323661269189) ) { // 0.135... = exp(-2)
     y = 1.0 - y;
     code = 0;
  }

  if (y > 0.13533528323661269189 ) {
     y  = y - 0.5;
     y2 = y * y;
     x  = y + y * (y2 * EvalPolynomial( y2, P0, 4)/EvalMonomial( y2, Q0, 8 ));
     x  = x * s2pi;
     return x;
  }

  x  = TMath::Sqrt( -2.0 * TMath::Log(y) );
  x0 = x - TMath::Log(x)/x;

  z = 1.0/x;
  if ( x < 8.0 ) { // y > exp(-32) = 1.2664165549e-14
     x1 = z * EvalPolynomial( z, P1, 8 )/EvalMonomial( z, Q1, 8 );
  } else {
     x1 = z * EvalPolynomial( z, P2, 8 )/EvalMonomial( z, Q2, 8 );
  }
  x = x0 - x1;
  if ( code != 0 ) x = -x;
  return x;
}

//______________________________________________________________________
Double_t // TRolke
MyBayesUL::EvalPolynomial(Double_t x, const Int_t  coef[], Int_t N)
{
 // evaluate polynomial

  const Int_t   *p;
  p = coef;
  Double_t ans = *p++;
  Int_t i = N;

  do
     ans = ans * x  +  *p++;
  while( --i );

  return ans;
}

//______________________________________________________________________
Double_t //TRolke
MyBayesUL::EvalMonomial(Double_t x, const Int_t coef[], Int_t N)
{
   // evaluate mononomial

  Double_t ans;
  const Int_t   *p;

  p   = coef;
  ans = x + *p++;
  Int_t i = N-1;

  do
     ans = ans * x  + *p++;
  while( --i );

  return ans;
}