C
C     Estimating small signals by using maximum likelihood
C     and Poisson statistics
C    
C     M. D. Hannam and W. J. Thompson
C     Department of Physics and Astronomy
C     University of North Carolina
C     Chapel Hill, NC 27599-3255, USA
C
C     Submitted to Nuclear Instruments and Methods A 
C    
C     Documentation follows equations in this paper.   
C    
C     Enquiries to e-mail: wmjthompson@msn.com
C    
C     Version 1.0, November 1998
C
C       
      integer counts(4096), startch, finch, center, range,
     *first_attempt, nchoice
      real bi(4096), si(4096), st_dev, loglike, sig, back, SS, BB,
     *log_likelihood, maxbl, maxsl,
     *sigma2, lower_sig_S, upper_sig_S, lower_sig_B,
     *upper_sig_B, errl, oversigS, oversigB, sigSa, sigBa, pi,
     *Ba, Sa, aa, bbb, cc, RR, areafactor,
     *first_norm_log, first_norm, first_norm_llS, first_norm_llB,
     *first_norm_likelihoodB, first_norm_likelihoodS,req_confS_percent 
      real norm_logB, norm_logS, start_normB, start_normS,
     *Bll1, Bll2, Sll1, Sll2, likelihoodB1, likelihoodB2,
     *likelihoodS1, likelihoodS2, deltaB, deltaS, exptB1, exptB2,
     *exptS1, exptS2, bratio, sratio, newsigS, newsigB, mean, det,
     *C11, C12, C22, center_ll, intres,
     *ri, term, muSquared,S0, S1, S2, S3, Fs, Fb, norm2B, norm2S,
     *confB, confS, Fbr, Fsr, Fb2, Fs2, KK, bstart, bres, bfin, sstart,
     *sres, sfin, rho, req_sig_S, req_sig_B, req_confS, req_confB,
     *mBG,bBG
      integer i,ii,channel,b,s,diff2,brange,srange,iosize,nope
      common bi,si,counts, startch, finch
      character*30 inputfile
C
      write(6,*) ' ***THE INPUT SPECTRUM MUST START WITH CHANNEL 1***'
      write(6,*)
      write(6,*)' Enter input filename:'
      read(5,*) inputfile
      open(UNIT=2,FILE=inputfile,STATUS='OLD')
      write(6,*)
      write(6,*) ' Enter startch, finch, center, st_dev'     
      read(5,*) startch, finch, center, st_dev      
      write(6,*)
      write(6,*) ' Constant (0) or linear (1) background?'
      read(5,*) nchoice
      
      if(nchoice.eq.1) then
	   write(6,*) ' Enter m and b in mi+b:'
           read(5,*) mBG,bBG
      endif
      
      OPEN (UNIT=10,FILE='poissonMLE.out',STATUS='unknown')
c
      write(*,100)
  100	format(/" PoissonMLE: ", 
     *"Estimating small signals by MLE for Poisson statistics",
     */," Version 1.0" )
C
      pi = 4.0*atan(1.0)
      nope = 1
      first_attempt = 1
C
      write(10,356) inputfile
  356   format(/," Spectrum: ",a30)
      
      write(*,110) startch, finch
      write(10,110) startch, finch
  110	format(/," Start and finish channels : ",I4," and ",I4)
      write(*,120) center, st_dev
      write(10,120) center, st_dev
  120	format(" Signal centered on channel: ",I4,
     */,       " with standard deviation   : ",F8.3)

      write(10,*)     
      if(nchoice.eq.0) then
          write(10,*) ' Constant BG b=1'
	else
	  write(10,*) ' Linear BG mx+b'
	  write(10,*) '  m=',mBG
	  write(10,*) '  b=',bBG
      endif
     
C     Spectrum input: channel number and counts 
      do  i = 1, finch
        read(2,*) channel, counts(i)
      end do
      range = finch - startch
C
C     Test signal shape s(i) is a Gaussian centered on 'center',
C     standard deviation 'st_dev' and unit height at center 
      sigma2 = 2*st_dev*st_dev
      do i = startch, finch
        diff2 = (i - center)*(i - center)
        si(i) = exp( - diff2 / sigma2)
C       Test background shape b(i) 
      
        if(nchoice.eq.0) then
           bi(i) = 1
         else
	   bi(i)=mBG*i+bBG
	endif
	
      end do
      areafactor = st_dev*sqrt(2*pi)
C
C     The rest of PoissonMLE analyses any spectrum <= 4096 channels 
C
C     Analytical estimates of background, signal,
C     and standard deviations; Section 4.3 
      S0 = 0.0
      S1 = 0.0
      S2 = 0.0
      S3 = 0.0
      Fs = 0.0
      Fb = 0.0
      do i = startch, finch
C       Sums according to Eq.(27) 
        ri = si(i) / bi(i)    
        term = counts(i)
        S0 = S0 + term   
        term = term * ri
        S1 = S1 + term   
        term = term * ri
        S2 = S2 + term   
        S3 = S3 + term*ri
C       Sums over b(i) and s(i) from Eqs.(17) and (18)   
        Fb = Fb + bi(i)  
        Fs = Fs + si(i)
      end do
C     Equation (31) 
      aa = Fb*S3 - Fs*S2
      bbb = Fb*S2 - Fs*S1 
      cc = Fb*S1 - Fs*S0
C     Equation (30) 
      RR = (cc/bbb) * (1 + aa*cc/(bbb*bbb))
C
C     Background and signal estimates 
      Ba = (S0 - RR*S1 + (RR**2)*S2) / Fb    
      Sa = RR*Ba
C
C     Large-N-limit sigmas and elements of covariance matrix
C     Section 4.3 
      oversigB = 0.0     
      oversigS = 0.0 
      K = 0.0
      do i = startch, finch
        muSquared = (Sa*si(i)+Ba*bi(i))**2
        oversigB = oversigB + counts(i)*bi(i)**2/muSquared
        oversigS = oversigS + counts(i)*si(i)**2/muSquared
C       Eq.(39)        
        K = K + counts(i)*bi(i)*si(i)/muSquared
      end do
      sigBa = 1.0/oversigB
      sigSa = 1.0/oversigS
C
C     Analytical error estimates Section 4.3  
C     Summations in Eq. (36) 
      Fbr = 0.0 
      Fsr = 0.0 
      Fb2 = 0.0 
      Fs2 = 0.0
      do i = startch, finch
        mean = Ba*bi(i) + Sa*si(i)
        bratio = bi(i)/mean   
        sratio = si(i)/mean
        Fb2 = Fb2 + counts(i)*bratio**2
        Fbr = Fbr + bi(i)*(counts(i) - mean)/mean
        Fs2 = Fs2 + counts(i)*sratio**2
        Fsr = Fsr + si(i)*(counts(i) - mean)/mean
      end do
      newsigB = (sqrt(Fbr**2 + Fb2) + Fbr)/Fb2
      newsigS = (sqrt(Fsr**2 + Fs2) + Fsr)/Fs2
      write(*,130) Ba, newsigB, Sa, newsigS
      write(10,130) Ba, newsigB, Sa, newsigS
  130	format(/," Analytical background level: ",F8.3," +- ",F7.3,
     */," Analytical signal level    : ",F8.3," +- ",F7.3)
C
C     Grid search for B and S based on analytical estimates  
      do while ( nope == 1 )
        if ( first_attempt > 0 ) then
          bstart = Ba - 2.0*sigBa
          bfin = Ba + 2.0*sigBa
          sstart = Sa - 2.0*sigSa  
          sfin = Sa + 2.0*sigSa
          bres = (bfin - bstart) / 10.0
          sres = (sfin - sstart) / 10.0
          iosize = 0.0
C         Start with Gaussian 1-sigma C.L. in S
          req_confS = 0.683 
        else
         write(*,140) 
         write(10,140) 
  140    format(/," Enter background level (B) search values as",/,
     *	" start_B  end_B  B_step")
         write(*,*) "[ Enter B_step as zero to quit PoissonMLE ]"
         read(*,*) bstart, bfin, bres
	 write(10,*) bstart, bfin, bres
         if ( bres == 0.0 ) then
C         Graceful exit from PoissonMLE
          close(10)
          call exit
         end if
         write(*,150) 
	 write(10,150) 
  150    format(/," Enter signal level (S) search values as",/,
     *	" start_S  end_S  S_step")
         read(*,*) sstart, sfin, sres
	 write(10,*) sstart, sfin, sres
         write(*,*) "Enter required S confidence in percent: "
         read(*,*) req_confS_percent
         req_confS = 0.01*req_confS_percent
         req_confB = req_confS
        end if
C
      first_attempt = 0
      brange = (bfin - bstart)/bres + 1.1     
      srange = (sfin - sstart)/sres	+ 1.1
      maxbl = -3000000.0
      maxsl = -3000000.0
C     Loop over background level range 
      back = bstart
      do b = 1, brange
C       Loop over signal level range 
        sig = sstart
        do s = 1, srange
          log_likelihood = loglike(sig, back)
          if (log_likelihood > maxsl) then
           maxsl = log_likelihood
           SS = sig
          end if
          sig = sig + sres 
        end do 
        if (maxsl > maxbl) then
          maxbl = maxsl
          BB = back
        end if
        back = back + bres
      end do

C     Standard deviation estimates by lnL - 1/2 method Section 4.3 
      sig = SS 
      back = BB
      errl = 100000.0
      do while ( errl > maxbl - 0.5 )
        errl = loglike(sig, back)
C       Decrement of signal level  
        sig = sig - 0.1
      end do
      lower_sig_S = SS - sig
      sig = SS 
      back = BB
      errl = 100000.0
      do while ( errl > maxbl - 0.5 )
        errl = loglike(sig, back)
C       Increment of signal level  
        sig = sig + 0.1
      end do
      upper_sig_S = -SS + sig
      sig = SS 
      back = BB
      errl = 100000.0
      do while ( errl > maxbl - 0.5 )
        errl = loglike(sig, back)
C       Decrement of background level
        back = back - 0.1
      end do
      lower_sig_B = BB - back 
      confB = 0.0 
      confS = 0.0 
      norm2B = 0.0 
      norm2S = 0.0
      sig= SS 
      back= BB
      errl= 100000.0
      do while ( errl > maxbl - 0.5 )
        errl = loglike(sig, back)
C       Increment of background level
        back = back + 0.1
      end do
      upper_sig_B = -BB + back
C
C     Likelihood normalization and integration for confidence limits
C     Section 4.2 
      first_norm_log = loglike(SS, BB)
      first_norm = exp(first_norm_log / range)
      intres = 25.0
      deltaB = upper_sig_B / intres
      deltaS = upper_sig_S / intres
      start_normB = BB - 4.0*intres*deltaB
      start_normS = SS - 4.0*intres*deltaS
C     Normalization integrals Eqs.(21) to (24) 
      do i = 0, 8.0*intres
        exptB1 = start_normB + i*deltaB
        exptS1 = start_normS + i*deltaS
        first_norm_llB = loglike(SS,exptB1)-first_norm_log
        first_norm_llS = loglike(exptS1,BB)-first_norm_log
        first_norm_likelihoodB = exp(first_norm_llB)
        first_norm_likelihoodS = exp(first_norm_llS)
        norm2B =  norm2B + first_norm_likelihoodB*deltaB
        norm2S =  norm2S + first_norm_likelihoodS*deltaS
      end do
      norm_logB = first_norm_log + log(norm2B)
      norm_logS = first_norm_log + log(norm2S)
C
      ii = 0
      exptB2 = BB   
      exptS1 = SS
      Sll1 = loglike(exptS1, BB) - norm_logS
      Bll2 = loglike(SS, exptB2) - norm_logB
      likelihoodS1 = exp(Sll1)
      likelihoodB2 = exp(Bll2)
      center_ll = likelihoodS1
      confB = confB + likelihoodB2*deltaB    
      confS = confS + likelihoodS1*deltaS
C     Integrate out to desired confidence interval in S 
      do while ( confS < req_confS )
        exptB1 = BB - float(ii)*deltaB
        exptB2 = BB + float(ii)*deltaB
        exptS1 = SS - float(ii)*deltaS
        exptS2 = SS + float(ii)*deltaS
        Bll1 = loglike(SS, exptB1) - norm_logB
        Sll1 = loglike(exptS1, BB) - norm_logS
        Bll2 = loglike(SS, exptB2) - norm_logB
        Sll2 = loglike(exptS2, BB) - norm_logS
        likelihoodB1 = exp(Bll1)
        likelihoodS1 = exp(Sll1)
        likelihoodB2 = exp(Bll2)
        likelihoodS2 = exp(Sll2)
        confB = confB + (likelihoodB1 + likelihoodB2)*deltaB
        confS = confS + (likelihoodS1 + likelihoodS2)*deltaS
        ii = ii + 1
      end do
      req_sig_B = float(ii)*deltaB         
      req_sig_S = float(ii)*deltaS
      write(*,160) BB,req_sig_B,100*confB,SS,req_sig_S,100*confS
      write(10,160) BB,req_sig_B,100*confB,SS,req_sig_S,100*confS
  160 format(/," Background level (B): ",F8.3," +- ",F7.3,
     *" with C.L. ",F5.1," %",
     */," Signal level (S)    : ",F8.3," +- ",F7.3,
     *" with C.L. ",F5.1," %")
      write(*,170) SS*areafactor,req_sig_S*areafactor,100*confS
      write(10,170) SS*areafactor,req_sig_S*areafactor,100*confS
  170	format(" Total signal (T)    : ",F8.3," +- ",F7.3,
     *" with C.L. ",F5.1," %")

C     Covariance matrix calculation 
      oversigB = 0.0 
      oversigS = 0.0 
      KK = 0.0
      do i = startch, finch
        muSquared = counts(i)/(SS*si(i)+BB*bi(i))**2
        oversigB = oversigB + muSquared*bi(i)**2
        oversigS = oversigS + muSquared*si(i)**2
C       Eq.(39)
        KK = KK + muSquared*bi(i)*si(i)
      end do
      sigBa = 1.0/sqrt(oversigB)   
      sigSa = 1.0/sqrt(oversigS)
      det = 1.0 - (KK*sigBa*sigSa)**2
C     Covariance matrix elements   
      C11 = sigBa/(det*sigSa)
      C22 = sigSa/(det*sigBa)
      C12 = - KK*sigSa*sigBa/det
      rho = C12/sqrt(C11*C22)
      write(*,180) rho
      write(10,180) rho
      write(*,*) "----------------------------------------------------"
      write(10,*) "----------------------------------------------------"
  180	format(/," Correlation coefficient: ",F7.3)
      end do

      end
C
C
      function loglike(Signal, Background)
C     Log likelihood for Poisson statistics Eqs.(4), (13) 
C
      real bi(4096), si(4096), Signal, Background, mu, loglike
      integer counts(4096), j, startch, finch
      common bi,si,counts, startch, finch
C
      loglike = 0.0
      do j = startch, finch
C     Eq.(13)
        mu = Background*bi(j) + Signal*si(j)
        loglike = loglike + log(mu)*counts(j) - mu
      end do
      end