c -----------------------------------------------------------------
c POILIM.F
c
c calculates Poisson upper limits given number of events
c observed, expected background, and desired confidence level,
c folding in errors on background and signal acceptance
c
c                                    John Conway  October 1994
c                                                 January 1998
c -----------------------------------------------------------------

      character*10 resp
      
      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 expected number of '',
     .         ''background events '',
     .         ''(real): '')'
      accept '(f10.4)',sbg

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

      type '($,'' desired accuracy: '')'
      accept '(f10.4)',prec

      call poilim(nev,sac,xbg,sbg,conf,prec,plim,perr)
 
      print '(//,'' limit: less than'',
     .         f10.3,'' +/-'',f6.3,'' signal events '',//)',plim,perr

      end


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

      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)

c record of guesses

      integer nguess,nsuc,nexp,nx,ns
      real    xpoi(1000),xtst(1000),stst(1000)

c ntuple storage
      real xtu(100)

c counters, etc.

      integer i,j,k,n,nmax

      real test,sest,dtst,dcdp,dpo,poi,p,plo,phi

      logical done,finished

c ---------------------------------------------------------------------------

      call ntuset(1,3,' ','poilim.ntu',1)

c idiot called the routine

      plim=0.
      perr = 0.
      if (nev.lt.0) return

c initial, generic guess

      dpo = 3.0
      poi = 5.0
      nmax = 1000000

c Phase I - get close to result

      print '(//,'' Phase I - get close '')'

      print '(//,''      N     delta(N)    CL      sigma(CL) '',
     .           ''   dev    trials '')'

      finished = .false.

      do while(.not.finished)

        done = .false.

        nexp = 0
        nsuc = 0

        do while(.not.done)

c perform 10k experiments
          nx = 10000
          call poiexp(nev,sac,xbg,sbg,poi,nx,ns)

          nexp = nexp + nx
          nsuc = nsuc + ns

c calculate CL and error

          test = float(nsuc)/float(nexp)
          sest = sqrt(test*(1.0-test)/float(nexp))
          if(sest.eq.0.) sest = 1.108/float(nexp)

          dtst = test-cl

c if we are more than 4 sigma away, make new guess

          if(sest.ne.0.) then
            if(abs(dtst/sest).gt.4.0) done = .true.
          endif

c too many trials - move to new point

          if(nexp.ge.nmax) done = .true.

c print current status

          print '($,5f10.6,i8,a)',poi,dpo,test,sest,dtst,
     .                              nexp,char(13)
c          call flush(6)

        enddo

        print '(5f10.6,i8)',poi,dpo,test,sest,dtst,nexp

c we exited the loop, either because we were 4 sigma
c away from the CL or made too many trials
 
        if(nexp.lt.nmax) then
          if(test.gt.cl) then
            dpo=dpo/1.9
            poi=poi-dpo
          else
            dpo=dpo*1.1
            poi=poi+dpo
          endif
        endif

        if(nexp.ge.nmax) finished = .true.

      enddo

c Phase II - map out region near desired value

      print '(//,'' Phase II - scan region and perform fit '',/)'

      nx = 100000
      nguess = 0

      do while(dpo.gt.prec/2.)

        plo = 0.01 * nint((poi-0.15)*100.)
        phi = 0.01 * nint((poi+0.15)*100.)
        nx = nx * 2.

        do p = plo,phi,0.02

          call poiexp(nev,sac,xbg,sbg,p,nx,ns)

c calculate CL and error

          test = float(ns)/float(nx)
          sest = sqrt(test*(1.0-test)/float(nx))
          if(sest.eq.0.) sest = 1.108/float(nx)

c record values
          nguess = nguess + 1

          xpoi(nguess) = p
          xtst(nguess) = test
          stst(nguess) = sest

        enddo

        call qfit(nguess,xpoi,xtst,stst,cl,poi,dpo,dcdp)

        plim = poi
        perr = dpo

        print '('' upper limit: '',f10.4,'' +/- '',f6.4,f10.5)',
     .     plim,perr,dcdp

      enddo

      print '(//)'

      return
      end
      

      subroutine poiexp(nev,sac,xbg,sbg,poi,nexp,nsuc)

      implicit none

c perform specified number of pseudoexperiments

      integer nev     ! number of events
      real    sac     ! error on acceptance
      real    xbg     ! expected background
      real    sbg     ! error on background
      real    poi     ! guess for upper limit
      integer nexp    ! number of pseudoexperiments
      integer nsuc    ! number of successes (returned)

      integer n,ntry,ierror

      integer nbkg,nsig

      real a,b,xsig,xbkg

      logical done

      n    = 0
      nsuc = 0

      do while(n.lt.nexp)

c get mean expected signal and background for this experiment

        xsig = -1.
        xbkg = -1.
        do while(xsig.lt.0..or.xbkg.lt.0.)
          call rannor(a,b)
          xsig = poi*(1.0 + sac*a)
          xbkg = xbg + sbg*b
        enddo

c generate poisson number of signal and background

        if(xbkg.gt.0.) then
          call poissn(xbkg,nbkg,ierror)
        else
          nbkg = 0
        endif

        if(ierror.ne.0) then
          print '('' error generating Poisson number: '',f10.3,2i10)',
     .      xbkg,nbkg,ierror
        endif

        if(xsig.gt.0.) then
          call poissn(xsig,nsig,ierror)
        else
          nsig = 0
        endif

        if(ierror.ne.0) then
          print '('' error generating Poisson number: '',f10.3,2i10)',
     .      xsig,nsig,ierror
        endif

c calculate fraction with total number of events exceeds nev, but
c still has fewer background than nev

        ntry = nbkg + nsig
        if(nbkg.le.nev) then
          n = n + 1
          if(ntry.gt.nev) nsuc = nsuc + 1
        endif

      enddo

      return
      end




      subroutine pfit(nguess,xpoi,xtst,stst,prob,poi,dpo,dcdp)

c fit line to data from CL versus N

      integer nguess
      real xpoi(*),xtst(*),stst(*),poi,dpo,dcdp,dmax

      sum   = 0.
      sumx  = 0.
      sumxx = 0.
      sumy  = 0.
      sumyy = 0.
      sumxy = 0.

      do i=1,nguess
        sum   = sum   + 1.0/stst(i)**2
        sumx  = sumx  + xpoi(i)/stst(i)**2
        sumxx = sumxx + xpoi(i)*xpoi(i)/stst(i)**2
        sumy  = sumy  + xtst(i)/stst(i)**2
        sumyy = sumyy + xtst(i)*xtst(i)/stst(i)**2
        sumxy = sumxy + xtst(i)*xpoi(i)/stst(i)**2
      enddo

      det = sum*sumxx - sumx**2

      if(det.ne.0.) then

c determine linear parameters, extrapolate to desired CL

        a = (sumy*sumxx-sumxy*sumx)/det
        b = (sumxy*sum-sumy*sumx)/det

        if(sum.gt.4.and.b.ne.0.) then

          poi = (prob-a)/b
          dcdp = b
      
          dc = 1.0/sum + sum/det*(poi-sumx/sum)**2

          dc = sqrt(dc)

          dpo = dc/dcdp

        else 

c here we handle those rare cases where we have zero slope

          print '('' error in pfit: zero slope! '')'

        endif

      else

c here we handle those rare cases where the determinant vanishes

        print '('' error in pfit: zero determinant! '')'

      endif

      if(poi.lt.0) poi = 0.

      return
      end


      subroutine qfit(nguess,xpoi,xtst,stst,prob,poi,dpo,dcdp)

      implicit none

c fit quadratic to data from CL versus N

      integer nguess
      real xpoi(*),xtst(*),stst(*),prob,poi,dpo,dcdp,dmax

      integer nord
      parameter(nord=3)

      real*8 a(nord),cov(nord,nord),v(100),ppoi,dc2
      real chisq,chipr,poi1,poi2,arg
      integer ndof

      call fitlsq(nguess,xpoi,xtst,stst,
     .            nord,a,cov,chisq,ndof,chipr)

      arg = a(2)**2 - 4.0*(a(1)-prob)*a(3)

      if(arg.gt.0.) then
        poi1 = (-a(2) + sqrt(arg))/2.0/a(3)
        poi2 = (-a(2) - sqrt(arg))/2.0/a(3)
      else
        print '('' negative discriminant! '')'
      endif

c     print '('' prob: '',3f10.3)',chipr,poi1,poi2

      if(abs(poi-poi1).lt.abs(poi-poi2)) then
        poi = poi1
      else
        poi = poi2
      endif

      dcdp = a(2) + 2.0*a(3)*poi

      ppoi = poi
      dc2 = cov(1,1) + 2.0*cov(1,2)*ppoi + 2.0*cov(1,3)*ppoi**2
     .      + cov(2,2)*ppoi**2 + 2.0*cov(2,3)*ppoi**3 + cov(3,3)*ppoi**4

      dpo = sqrt(dc2)/dcdp

      return
      end


      subroutine fitlsq(npoint,x,y,s,
     .                  nord,a,cov,chisq,ndof,chipr)

c linear least squares polynomial fit to data points

      logical scantu

      integer i,j,k,npoint
      real x(*),y(*),s(*)

      real*8 xx,yy,ss

      real*8 a(nord),cov(nord,nord),v(100)
      real chisq,chipr
      integer ndof

c accumulate vectors/matrices
      call vzero(cov,nord*nord*2)
      call vzero(a,nord*2)
      call vzero(v,nord*2)

      do k=1,npoint

        xx = x(k)
        yy = y(k)
        ss = s(k)

        if(yy.ne.0..and.ss.ne.0..and.
     .     (xlof.eq.xhif.or.(xlof.le.xx.and.xx.le.xhif))) then
          do i=1,nord
            if(i.ne.1) then
              v(i) = v(i) + yy*(xx**(i-1))/(ss**2)
            else
              v(i) = v(i) + yy/(ss**2)
            endif
            do j=1,nord
              if(i.ne.1.or.j.ne.1) then
                cov(i,j) = cov(i,j) + (xx**(i+j-2))/(ss**2)
              else
                cov(i,j) = cov(i,j) + 1./(ss**2)
              endif
            enddo
          enddo
        endif

      enddo

c invert variance matrix
      call invert(nord,cov,cov)

c calculate a
      do i = 1,nord
        a(i) = 0.
        do j=1,nord
          a(i) = a(i) + cov(i,j)*v(j)
        enddo
      enddo

c calculate chi squared
      chisq = 0.
      ndof = 0
      do k=1,npoint
        xx = x(k)
        yy = y(k)
        ss = s(k)
        if(yy.ne.0..and.ss.ne.0..and.
     .     (xlof.eq.xhif.or.(xlof.le.xx.and.xx.le.xhif))) then
          yf = a(1)
          do i=2,nord
            yf = yf + a(i)*(xx**(i-1))
          enddo
          chisq = chisq + (yy-yf)**2/ss**2
          ndof = ndof + 1
        endif
      enddo
      ndof = ndof - nord

      chipr = chiprb(chisq,ndof)

      return
      end



      function chiprb(x,ndf)  
      chiprb=0.0    
      test=x/ndf    
      if(test.ge.20.0) return 
      a=0.5*x       
      ibigx=0       
      if(a.gt.90.0) ibigx=1   
      y=0.0         
      if(x.ge.0.0.and.ndf.ge.1) go to 1 
      chiprb=-1.0   
      print 500     
  500 format(' chiprb error -- x < 0 or ndf < 1') 
      return        
    1 xd=x
      ieven=0       
      ntst=ndf/2    
      ntst=ntst*2   
      if(ntst.eq.ndf) ieven=1 
      if(ieven.eq.1) y=exp(-a)
      if(ndf.gt.2.and.ibigx.eq.0) y=exp(-a)       
      s=y 
      if(ieven.eq.0) s=erfc(sqrt(a))    
      if(ndf.gt.2) go to 10   
      chiprb=s      
      return        
   10 xd=0.5*(ndf-1)
      z=1.0         
      if(ieven.eq.0) z=0.5    
      if(ibigx.eq.0) go to 200
      e=0.0         
      if(ieven.eq.0) e=0.5723649429     
      c=alog(a)     
  100 if(z.gt.xd) go to 150   
      e=alog(z)+e   
      s=exp(c*z-a-e)+s        
      z=z+1.0       
      go to 100     
  150 chiprb=s      
      return        
  200 e=1.0         
      c=0.0         
      if(a.eq.0.0) go to 400  
      if(z.gt.xd) go to 400   
      if(ieven.eq.0) e=0.5641895435/sqrt(a)       
  300 e=e*a/z       
      c=c+e         
      z=z+1.0       
      if(z.le.xd) go to 300   
  400 chiprb=c*y+s  
      return        
      end 


      FUNCTION ERFXX(X)       
C         
C     ENTRY POINTS ... ERF, ERFC, FREQ. 
C         
C     THESE THREE FUNCTIONS ARE COMPUTED FROM THE RATIONAL APPROXIMAT-
C     IONS OF W.J.CODY, MATHEMATICS OF COMPUTATION, VOLUME 22 (1969), 
C     PAGES 631-637.
C         
C     FOR ABS(X) LE 0.47 THE BASIC FUNCTION IS ERF.  FOR ABS(X) GT 0.47         
C     THE BASIC FUNCTION IS ERFC.  THE FINAL RESULT IS OBTAINED IN TERM         
C     OF THE BASIC FUNCTION AS SHOWN IN THE FOLLOWING TABLE.  A=ABS(X)
C     AND C=SQRT(2).
C       FUNCTION     A.LE.(0.47)        
C       REQUIRED       (ALL X)          (NEGATIVE X)       (POSITIVE X)         
C         
C       ERF(X)         ERF(X)             ERFC(A)-1          1-ERFC(X)
C         
C       ERFC(X)       1-ERF(X)            2-ERFC(A)           ERFC(X) 
C         
C       FREQ(X)   0.5+0.5*ERF(X/C)      0.5*ERFC(A/C)    1-0.5*ERFC(X/C)        
C         
C         
      DIMENSION P1(4),Q1(4), P2(8),Q2(8), P3(5),Q3(5)       
      DATA P1/2.42667 95523 0532 E2,    
     1        2.19792 61618 2942 E1,    
     2        6.99638 34886 1914 E0,    
     3       -3.56098 43701 8154 E-2/   
      DATA Q1/2.15058 87586 9861 E2,    
     1        9.11649 05404 5149 E1,    
     2        1.50827 97630 4078 E1,    
     3        1.00000 00000 0000 E0/    
      DATA P2/3.00459 26102 0162 E2,    
     1        4.51918 95371 1873 E2,    
     2        3.39320 81673 4344 E2,    
     3        1.52989 28504 6940 E2,    
     4        4.31622 27222 0567 E1,    
     5        7.21175 82508 8309 E0,    
     6        5.64195 51747 8974 E-1,   
     7       -1.36864 85738 2717 E-7/   
      DATA Q2/3.00459 26095 6983 E2,    
     1        7.90950 92532 7898 E2,    
     2        9.31354 09485 0610 E2,    
     3        6.38980 26446 5631 E2,    
     4        2.77585 44474 3988 E2,    
     5        7.70001 52935 2295 E1,    
     6        1.27827 27319 6294 E1,    
     7        1.00000 00000 0000 E0/    
      DATA P3/-2.99610 70770 3542 E-3,  
     1        -4.94730 91062 3251 E-2,  
     2        -2.26956 59353 9687 E-1,  
     3        -2.78661 30860 9648 E-1,  
     4        -2.23192 45973 4185 E-2/  
      DATA Q3/1.06209 23052 8468 E-2,   
     1        1.91308 92610 7830 E-1,   
     2        1.05167 51070 6793 E0,    
     3        1.98733 20181 7135 E0,    
     4        1.00000 00000 0000 E0/    
      DATA CONST1/ 0.70710 67811 86548 /
C        ( CONST1= SQRT(1/2) )
      DATA CONST2/0.56418 95835 47756 / 
C        ( CONST2= SQRT(1/PI) )         
C         
C         
C  ENTRY POINTS. SET IENTRY (=1 FOR ERF, =2 FOR ERFC, =3 FOR FREQ).   
C         
      ENTRY ERF(X)  
      IENTRY=1      
      T=X 
      A=ABS(T)      
      IF(A.LE.6.0) GO TO 11   
      ERFXX=SIGN(1.0,X)       
      GO TO 90      
C         
      ENTRY ERFC(X) 
      IENTRY=2      
      T=X 
      IF(T.GE.-6.0) GO TO 10  
      ERFXX=2.0     
      GO TO 90      
C         
      ENTRY FREQ(X) 
      IENTRY=3      
      T=CONST1*X    
      IF(T.LE.6.0) GO TO 10   
      ERFXX=1.0     
      GO TO 90      
C         
C  SELECT BASIC FUNCTION. SET IBASIC (=1 FOR ERF, =2 FOR ERFC).       
C         
   10 A=ABS(T)      
   11 S=T**2        
      IF(A.GT.0.47) GO TO 20  
C         
C  IBASIC=1.  SET Y=ERF(T).   
C         
      IBASIC=1      
      Y=T*(P1(1)+S*(P1(2)+S*(P1(3)+S*P1(4) )))    
     1   /(Q1(1)+S*(Q1(2)+S*(Q1(3)+S*Q1(4) )))    
      GO TO 30      
C         
C  IBASIC=2. SET Y=ERFC(A).   
C         
   20 IBASIC=2      
      IF(A.GT.4.0) GO TO 21   
      Y=EXP(-S)*(P2(1)+A*(P2(2)+A*(P2(3)+A*(P2(4)+A*(P2(5)+ 
     1        A*(P2(6)+A*(P2(7)+A*P2(8) )))))))   
     2         /(Q2(1)+A*(Q2(2)+A*(Q2(3)+A*(Q2(4)+A*(Q2(5)+ 
     3        A*(Q2(6)+A*(Q2(7)+A*Q2(8) )))))))   
      GO TO 30      
C         
   21 Y=0.0         
      IF(A.GT.26.0) GO TO 30  
      R=1.0/A       
      U=R**2        
      Y=R*EXP(-S)*( CONST2 +  
     1     U*(P3(1)+U*(P3(2)+U*(P3(3)+U*(P3(4)+U*P3(5) )))) 
     2      /(Q3(1)+U*(Q3(2)+U*(Q3(3)+U*(Q3(4)+U*Q3(5) )))) )         
C         
C  EXPRESS FINAL RESULT IN TERMS OF Y.  
C         
   30 IF(IENTRY.NE.1) GO TO 40
      IF(IBASIC.EQ.2) GO TO 31
      ERFXX=Y       
      GO TO 90      
   31 ERFXX=1.0-Y   
      IF(X.LT.0.0) ERFXX=-ERFXX         
      GO TO 90      
C         
   40 IF(IENTRY.NE.2) GO TO 50
      IF(IBASIC.EQ.2) GO TO 41
      ERFXX=1.0-Y   
      GO TO 90      
   41 ERFXX=Y       
      IF(X.LT.0.0) ERFXX=2.0-Y
      GO TO 90      
C         
   50 R=0.5*Y       
      IF(IBASIC.EQ.2) GO TO 51
      ERFXX=0.5+R   
      GO TO 90      
   51 ERFXX=1.0-R   
      IF(X.LT.0.0) ERFXX=R    
   90 RETURN        
      END 


      subroutine invert(n,a,b)
c
c  routine to invert an n x n matrix a, n < 16,
c  putting the result in b
c
      real*8 a(n,n),b(n,n),q(15,15),e
      do 1 i=1,15
      do 1 j=1,15
    1 q(i,j)=0.000000000000000
      do 5 i=1,n
      do 5 j=1,n
    5 b(i,j)=a(i,j)
      do 10 i=1,n
   10 if(b(i,i) .ne. 0.) q(i,i)=1.
      do 40 k=1,n
      if(b(k,k) .eq. 0.) go to 40
      do 30 i=1,n
      if(i .eq. k) go to 30
      e=b(i,k)
      do 20 j=1,n
      b(i,j)=b(i,j) - e*(b(k,j)/b(k,k))
      if(j .eq. k) b(i,j)=0.
      q(i,j)=q(i,j) - e*(q(k,j)/b(k,k))
   20 continue
   30 continue
   40 continue
      do 50 i=1,n
      e=b(i,i)
      if(e .eq. 0.) e=1.
      do 50 j=1,n
   50 b(i,j)=q(i,j)/e
      return
      end
