#include <math.h>
#include <stddef.h>

#include "genlimit.h"


/*  Joel Heinrich  8 April 2005

Returns cross section upper limit.

*uncertainty (if not null pointer) returned with uncertainty of
limit due to Monte Carlo statistical fluctuations of the prior
ensemble.

See CDF note 7587 for details:
http://www-cdf.fnal.gov/publications/cdf7587_genlimit.pdf

*/




double cslimit(double beta,int nchan,int nens,const int nobs[],const EB* ens,
	       int* ngl,double xgl[],double lwgl[],
	       PRIOR prior,double* uncertainty) {
  const double eps=1.0e-6;
  double norm = csint(0,nchan,nens,nobs,ens,ngl,xgl,lwgl,prior,NULL);
  double limit = galim(beta,nchan,nens,nobs,ens);
  double dl=limit, rpdf=0;
  double lo=0, hi=1e200;

  while(fabs(dl)>1.0e-10*limit) {
    const double pbeta =
      1-csint(limit,nchan,nens,nobs,ens,ngl,xgl,lwgl,prior,NULL)/norm;
    rpdf = 1/cspdf(limit,norm,nchan,nens,nobs,ens,prior);
    if(pbeta>beta) {
      hi=limit*(1+eps);
    } else {
      lo=limit*(1-eps);
    }
    dl = (pbeta-beta)*rpdf;
    if(limit-dl>=lo && limit-dl<=hi) {
      limit -= dl;
    } else {
      dl = limit - 0.5*(hi+lo);
      limit = 0.5*(hi+lo);
    }
  }

  if (uncertainty) {
    double i1=0, i2=0, v11=0, v22=0, v12=0;
    const double c = 1-beta;
    csint2(0,limit,nchan,nens,nobs,ens,ngl,xgl,lwgl,prior,
	   &i1,&i2,&v11,&v12,&v22);
    *uncertainty = rpdf*sqrt(v22 + v11*c*c - 2*v12*c)/norm   ;
  }


  return limit;
}