#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 cscutlimit(double beta,double smax,
		  int nchan,int nens,const int nobs[],const EB* ens,
		  int* ngl,double xgl[],double lwgl[],
		  PRIOR prior,double* uncertainty) {
  const double norm = csint(0,nchan,nens,nobs,ens,ngl,xgl,lwgl,prior,NULL);
  const double tail = csint(smax,nchan,nens,nobs,ens,ngl,xgl,lwgl,prior,NULL);
  const double eps=1.0e-6;
  double limit = galim(beta,nchan,nens,nobs,ens), rpdf=0;
  double dl=limit, lo=0, hi=smax;

  if(beta<=0) return 0;
  if(beta>=1) return smax;

  if(limit>smax || limit<0) dl = limit = 0.5*smax;

  while(fabs(dl)>1.0e-10*limit) {
    double pbeta =
      1-(csint(limit,nchan,nens,nobs,ens,ngl,xgl,lwgl,prior,NULL)-tail)/
      (norm-tail);
    rpdf = 1/cspdf(limit,norm-tail,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;

    csint2cut(0,limit,smax,nchan,nens,nobs,ens,ngl,xgl,lwgl,prior,
	      &i1,&i2,&v11,&v12,&v22);
    *uncertainty = rpdf*sqrt(v22 + v11*c*c - 2*v12*c)/(norm-tail)   ;
  }


  return limit;
}