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

#include "genlimit.h"

/*  Joel Heinrich  8 April 2005

Returns cross section posterior p.d.f. evaluated at s.

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

*/


double cspdf(double s,double norm,
	     int nchan,int nens,const int nobs[],const EB* ens,PRIOR prior) {
  int i,k;
  double sum = 0, lgp = 0;
  const EB* p = ens;

  assert(nens>0);
  assert(prior==flat || prior==corr);
  for(k=0;k<nchan;++k)
    lgp -= lgamma(nobs[k]+1);

  for(i=0;i<nens;++i) {
    double t = lgp, esum=0;
    for(k=0;k<nchan;++k) {
      const double mu = s*p->e + p->b;
      const int n = nobs[k];
      esum += p->e;
      t += ( (n>0) ? n*log(mu) : 0 ) - mu;
      ++p;
    }
    sum += (prior==flat) ? exp(t) : esum*exp(t);
  }
  return sum/(norm*nens);
}