#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<assert.h>

#include "genlimit.h"


/* Joel Heinrich  13 April 2005

This program, example.c, computes cross section upper limit for
example described in CDF note 7587.  See

http://www-cdf.fnal.gov/physics/statistics/statistics_software.html

and

http://www-cdf.fnal.gov/publications/cdf7587_genlimit.pdf


*/





double rlx();  /* uniform generator, calls CERNLIB's RANLUX */
void rlxrandomize(); /* initializes RANLUX with random seed */

/* generates gamma distribution deviates with mean a and variance a */
double gamma_gen(double a,double (*ru)());


#define NCHAN 4

int main(){
  const double beta=0.95;
  const int nchan=NCHAN, nens=150000;
  const PRIOR prior=flat;
  EB* ens = (EB*)malloc(nchan*nens*sizeof(EB));
  int nobs[NCHAN];
  int i,k;
  double xgl[300], lwgl[300];
  const int nglmax = sizeof(xgl)/sizeof(xgl[0]);
  int ngl=0;
  double limit=0, uncert=0;

  rlxrandomize();


  nobs[0] = 6; /*  6 events observed in channel 0 */
  nobs[1] = 1; /*  1 event  observed in channel 1 */
  nobs[2] = 1; /*  1 event  observed in channel 2 */
  nobs[3] = 8; /*  8 events observed in channel 3 */

  /* generate prior ensemble */

  k=0;
  for(i=0;i<nens;++i) {
    const double gam25=gamma_gen(25,rlx); /* for chan 2 and 3 background */

    /* channel 0 */
    ens[k].e = 0.04*gamma_gen(25,rlx);    /* mean = 1   , sigma = 0.2  */
    ens[k++].b = gamma_gen(4,rlx);        /* mean = 4   , sigma = 2    */
    
    /* channel 1 */
    ens[k].e = 0.01*gamma_gen(100,rlx);   /* mean = 1   , sigma = 0.1  */
    ens[k++].b = 0.1*gamma_gen(16,rlx);   /* mean = 1.6 , sigma = 0.4  */
    
    /* channel 2 */
    ens[k].e = 1e-4*gamma_gen(1e4,rlx);   /* mean = 1   , sigma = 0.01 */
    ens[k++].b = 0.01*gam25;              /* mean = 0.25, sigma = 0.05 */
    
    /* channel 3 */
    ens[k].e = 0.04*gamma_gen(25,rlx);    /* mean = 1   , sigma = 0.2  */
    ens[k++].b = 0.2*gam25;               /* mean = 5   , sigma = 1    */

    /* expected backgrounds of channels 2 and 3 are completely correlated */
    
  }
  
  /* print out the first 5 lines of the prior ensemble */

  for(k=0;k<nchan;++k)
    printf("  e%d     b%d     ",k,k);

  printf("\n");

  for(i=0;i<5;++i) {
    for(k=0;k<nchan;++k) {
      printf("%.4f %.4f   ",ens[i*nchan+k].e,ens[i*nchan+k].b);
    }
    printf("\n");
  }
  fflush(stdout);



  /* calculate and print out limit */

  ngl = 0;  /* make sure set to zero before first call */
  limit = cslimit(beta,nchan,nens,nobs,ens,&ngl,xgl,lwgl,prior,&uncert);

  /* assertion to catch error in dimension of arrays for Gauss-Laguerre
     abscissas and weights */
  assert(ngl <= nglmax );

  printf("upper limit at %.1f%% credibility level:   ",beta*100);
  
  printf("%.4f +/- %.4f\n",limit,uncert);
  fflush(stdout);
  


  free(ens);
  return 0;
}