
c BAYES - program to check equivalency of Bayes and frequentist 
c         upper limit calculations

      implicit none

      integer nev  ! number of observed events
      real*4  sac  ! relative uncertainty in acceptance
      real*4  xbg  ! expected background
      real*4  sbg  ! uncertainty in background
      real*4  cl   ! desired CL
      real*4  prec ! desired precision
      real*4  plim ! upper limit on Poisson process
      real*4  perr ! error on upper limit on Poisson process

      type '($,'' enter number of observed events: '')'
      accept '(i)',nev

      type '($,'' enter relative error on acceptance: '')'
      accept '(f10.4)',sac

      type '($,'' enter expected number of background events '',
     .        ''(real): '')'
      accept '(f10.4)',xbg

      type '($,'' enter error on number of background events '',
     .        ''(real): '')'
      accept '(f10.4)',sbg

      type '($,'' enter confidence level: '')'
      accept '(f10.4)',cl

      type '($,'' integration step size (0.02 or less): '')'
      accept '(f10.4)',prec

      call baylim(nev,sac,xbg,sbg,cl,prec,plim,perr)

       print '(//,'' limit: less than'',
     .         f10.3,'' +/-'',f6.3,'' signal events '',//)',plim,perr

      end


      subroutine baylim(nev,sac,xbg,sbg,cl,prec,plim,perr)

c calculate Bayesian upper limit on Poisson process

      implicit none

c arguments

      integer nev     ! number of events
      real    sac     ! error on acceptance
      real    xbg     ! expected background
      real    sbg     ! error on background
      real    cl      ! desired CL
      real    prec    ! desired precision
      real    plim    ! Poisson upper limit  (returned)
      real    perr    ! error on upper limit (returned)

      real    xev,xevmax,dxev

      logical done

      real*8 blike,xlike

      integer i,nlist,icl
      real blist(10000),xlist(10000),bint,bsum,bcl

      real xtu(100)

c      call ntuset(1,2,' ','bayes.ntu',1)

c scan likelihood L(nev|xev,...)

      xevmax = 10.0*nev
      if(xevmax.lt.10.) xevmax = 10.
      dxev = prec
      nlist = 0
      bsum = 0.

      done = .false.

      xev = dxev/2.

c      do xev = dxev/2.,xevmax,dxev

      do while(.not.done)

        xlike = blike(nev,sac,xbg,sbg,xev)

        nlist = nlist + 1
        xlist(nlist) = xev
        blist(nlist) = xlike

        bsum = bsum + xlike

        print '($,i5,2g12.5,a)',nev,xev,xlike,char(13)

        xtu(1) = xev
        xtu(2) = xlike

c        call ntufil(1,xtu)

        xev = xev + dxev
        if(blist(nlist)/blist(1).lt.1.0e-6) done = .true.
        if(xev.gt.xevmax) done = .true.

      enddo

c find place just below threshold
      bint = 0.
      icl = 0
      do i=1,nlist
        if(bint         .le.cl*bsum.and.
     .     bint+blist(i).gt.cl*bsum) then
          icl = i
          bcl = bint
        endif
        bint = bint + blist(i)
      enddo

c interpolate linearly
      plim = xlist(icl) + (xlist(icl+1)-xlist(icl))*
     .                     (cl*bsum-bcl)/blist(icl+1)

      perr = 0.

      return
      end

      real*8 function blike(nev,sac,xbg,sbg,xev)

c return likelihood to observe nev events given expected

      implicit none

      integer nev     ! number of events observed
      real    xev     ! number of events expected
      real    sac     ! error on acceptance
      real    xbg     ! expected background
      real    sbg     ! error on background

      integer i,nmax
      parameter(nmax=100000)

      real*4  a,b
      real*8  yyev,ssev,yyex,xxev,ssac,xxbg,ssbg,yybg,xint,xxx
      real*8  stirling

      logical done


c perform double gaussian integral by Monte Carlo

      xxev = nev
      xxbg = xbg
      ssbg = sbg
      ssac = sac
      xxev = xev

      ssev = xxev*ssac

      xint = 0.

      done = .false.

      do i=1,nmax

c pick expected background and signal from Gaussian

        yybg = -1.
        yyev = -1.

        do while(yybg.lt.0..or.
c    .           yybg.gt.nev.or.
     .           yyev.lt.0.)
          call rannor(a,b)
          yybg = xxbg + a*ssbg
          yyev = xxev*(1.0+b*ssac)
        enddo

c total expected         
        yyex = yybg + yyev

c value of integrand

        xxx = yyex**nev * dexp(-yyex) / stirling(nev)  ! Poisson term

c     .      dexp(-(xxbg-yybg)**2/2.0/ssbg**2)        ! gaussian for bkgd
c     .      dexp(-(xxev-yyev)**2/2.0/ssev**2)        ! gaussian for bkgd
c     .      /(2.0*pi)/ssbg/ssev

        xint = xint + xxx

      enddo

      blike = xint/float(nmax)

      return
      end

      real*8 function stirling(n)

c Stirling's approximation for n! with first correction term

      real*8 x,fact(0:9)
      parameter (root2pi=2.506628)
      data fact/1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./

      if(n.ge.0.and.n.le.9) then
        stirling=fact(n)
        return
      endif
      if(n.gt.0) then
        x=n
        stirling=root2pi*sqrt(x)*x**n*exp(-x)*(1.+1./12./x)
      else
        stirling=1.
      endif
      return
      end

