        PROGRAM POIS_BKGD
C--
C  This program calculates the combined 90 and 95% confidence level
C  upper limits for a signal where we assume a the number of observed
C  events is a result of a background (with some uncertainty associated
C  with its mean rate) and a potential signal.  This programme uses the 
C  PDG method to estimate a conservative upper limit on the signal rate. 
C
C--
C  The program prompts the user for the necessary information.
C--
C       Pekka K. Sinervo
C       University of Toronto and CDF
C
C       June 7, 1993
C--
        IMPLICIT NONE
C
C
        INTEGER         DUMP
        INTEGER         LUPRT, LUTERM
        PARAMETER      (LUPRT = 6)
        PARAMETER      (LUTERM= 5)
        INTEGER         I
        INTEGER         N_CL
C # limits to calculate
        PARAMETER      (N_CL = 2)               
C CL limits to calculate
        REAL            CL_PROB(N_CL)           
C Stores preliminary CL
        REAL            PRELIMINARY_CL(N_CL)    
C
C # observed events
        INTEGER         N_OBSERVED              
C background rate & sigma
        REAL            MU_BKGD, SIGMA_BKGD     
        REAL*8          SUM_BKGD, SUM_SIG_BKGD
        REAL*8          NEXT_BKGD_TERM, NEXT_SIG_BKGD_TERM
        INTEGER         I_EVENTS
        INTEGER         I_CL
        REAL            EVENTS_PRODUCED,        CL_FUZZ
        REAL            PREVIOUS_EVENTS_PRODUCED
        REAL            TOTAL_ACCEPTANCE
        REAL            TEMP
        REAL            DELTA_CL
        REAL            NEW_PROB,       PREVIOUS_PROB
        REAL            STEP
        INTEGER         N_TRIES,        MAX_TRIES
        INTEGER         N_EXPERIMENTS,  OLD_N_EXPERIMENTS
C
        INTEGER         MODE
C--
C  External functions referenced
C--
C Routine to perform MC experiments
        REAL            MC_POISSON_EXP  
C--
C  Initialize...
C--
        CALL RANINT
C--
C  Get the information from the user.
C--
C--
C  Do we want to calculate a background significance or an confidence
C  level interval...
C--
        WRITE(LUPRT, '('' Choices are: 1)  confidence level interval'',
     >          /,     ''              2)  background significance'')')
        READ(LUTERM, *) MODE
C--
C  Go do the background calculation elsewhere...
C--
        IF (MODE.EQ.2) THEN
          CALL BKGPRB
          STOP
        ENDIF
C-
C  Get the info for a confidence level calculation
C--
        WRITE(LUPRT, '('' What is the observed number of events?'')')
        READ(LUTERM, *) N_OBSERVED
        IF (N_OBSERVED.LT.0) THEN
          WRITE(LUPRT, '('' Illegal # of events---run aborted'')')
          STOP 4
        ENDIF
C--
        WRITE(LUPRT, '('' What is the estimated background rate?'')')
        READ(LUTERM, *) MU_BKGD
        IF (MU_BKGD.LT.0.0) THEN
          WRITE(LUPRT, '('' Illegal background rate---run aborted'')')
          STOP 4
        ENDIF
C--
        WRITE(LUPRT, '('' What is the absolute uncertainty (in events)''
     &,
     >  '' of the background rate?'')')
        READ(LUTERM, *) SIGMA_BKGD
        IF (SIGMA_BKGD.LT.0.0) THEN
          WRITE(LUPRT, '('' Illegal uncertainty---run aborted'')')
          STOP 4
        ENDIF
C--
C  Now we go the data--dump it back to the user...
C--
        WRITE(LUPRT, '( //,'' Calculate Poisson CL with background'',
     >                                             '' channels'',
     >          //,'' Number of observed events = '',I5,
     >          //,'' Background event rate     = '',F7.1,
     >          //,'' Uncertainty on background = '',F7.1)')
     >          N_OBSERVED, MU_BKGD, SIGMA_BKGD
C--
C  Do the calculation first without including the systematic uncertainty.
C--
        DO I_CL = 1, N_CL
C--
C  We make an initial guess at the rate to get the right scale
C  for the EVENTS_PRODUCED.
C--
          EVENTS_PRODUCED = MU_BKGD
C--
C  If we have no observed events, use a rough guess at the overall
C  CL limit.
C--
          IF (EVENTS_PRODUCED.LE.0.0) THEN
            EVENTS_PRODUCED = 3.0
          ENDIF
C--
C  Now verify that with the user.
C--
          WRITE(LUPRT,'(//'' Initial guess at limit is '',G10.3,
     >            '':  Enter your guess:'')')
     >       EVENTS_PRODUCED
          READ(LUTERM, *) EVENTS_PRODUCED
C--
C  Initialize other bookkeeping.
C--
C Prior guess...
          PREVIOUS_EVENTS_PRODUCED = 0.0  
C Calculated probability
          NEW_PROB = 0.0                  
C ...and previous step.
          PREVIOUS_PROB = 0.0             
          N_EXPERIMENTS = 0
          DELTA_CL = 9999999.
          N_TRIES  = 0
C Avoid run-away search....
          MAX_TRIES = 30        
C Don't dump distribution during search.
          DUMP = -6             
C--
C  Iterate till we converge on an answer.
C--
          WRITE(LUPRT, '(/,'' *** Calculate CL Using PDG Formula'',
     >     '' (does not account for bkgd uncertainty) ***'')')
C Required max change in probability
          CL_FUZZ = 0.0005        
          DO WHILE ((ABS(DELTA_CL).GT.CL_FUZZ .OR.
     >                         (N_TRIES.LT.3) .OR.
     >                         (STEP.GT.5.0)) .AND.
     >                   N_TRIES.LT.MAX_TRIES)
C Store prob from last step
            PREVIOUS_PROB = NEW_PROB    
C and calculate new one.
C
C--
C  Use the PDG formula to determine the CL this corresponds to.
C--
            SUM_BKGD = 1.0
            SUM_SIG_BKGD = EXP(-EVENTS_PRODUCED)
            NEXT_BKGD_TERM = 1.0
            NEXT_SIG_BKGD_TERM = EXP(-EVENTS_PRODUCED)
C
            DO I_EVENTS = 1, N_OBSERVED
              NEXT_BKGD_TERM = NEXT_BKGD_TERM * MU_BKGD / FLOAT(I_EVENTS
     &)
              SUM_BKGD = SUM_BKGD + NEXT_BKGD_TERM
C
              NEXT_SIG_BKGD_TERM = NEXT_SIG_BKGD_TERM * 
     >          (MU_BKGD + EVENTS_PRODUCED) / FLOAT(I_EVENTS)
              SUM_SIG_BKGD = SUM_SIG_BKGD + NEXT_SIG_BKGD_TERM
            ENDDO
C--
C  Calculate the CL interval this corresponds to.
C--
            NEW_PROB = SUM_SIG_BKGD / SUM_BKGD
C--
C  Figure out where to go in terms of events produced.  Calculate a
C  Newton-Ralphson step and keep it within reason.
C  On first go around, have to iterate to get two measured points.
C--
            IF (N_TRIES.EQ.0) THEN
              IF (NEW_PROB.LT.CL_PROB(I_CL)) THEN
                STEP =  -EVENTS_PRODUCED * 0.3
              ELSE
                STEP =   EVENTS_PRODUCED * 0.3
              ENDIF
            ELSE
              IF (PREVIOUS_PROB.EQ.NEW_PROB) THEN
                WRITE (LUPRT, '(/,'' ALERT...not converging'',
     >              ''...enter next step value: '')')
                READ(LUTERM,*) STEP
              ELSE
                STEP =
     >           ((PREVIOUS_EVENTS_PRODUCED - EVENTS_PRODUCED) /
     >            (PREVIOUS_PROB            - NEW_PROB) ) *
     >            (CL_PROB(I_CL) - NEW_PROB)
              ENDIF
            ENDIF
C--
C  Keep the step reasonable.
C--
            IF (ABS(STEP).GT.EVENTS_PRODUCED*0.2) THEN
              STEP = SIGN(EVENTS_PRODUCED*0.2, STEP)
            ENDIF
            IF (ABS(STEP).LT.1.0E-2) THEN
              STEP = SIGN(1.0E-2, STEP)
            ENDIF
            TEMP = EVENTS_PRODUCED + STEP
C--
C  Update book-keeping.
C--
            PREVIOUS_EVENTS_PRODUCED = EVENTS_PRODUCED
            EVENTS_PRODUCED          = TEMP
            PREVIOUS_PROB            = NEW_PROB
C--
C  See how close we are to final answer...
C--
            DELTA_CL = NEW_PROB - CL_PROB(I_CL)
            N_TRIES  = N_TRIES + 1
C--
            WRITE(LUPRT, '(/,'' Iteration Number '', I3,
     >          /,'' Signal rate  '',F7.1,''  Prob is '',G10.2,
     >       ''  Updated Signal rate '',F7.1)')
     >          N_TRIES, PREVIOUS_EVENTS_PRODUCED, PREVIOUS_PROB,
     >          EVENTS_PRODUCED
          ENDDO
C--
C  Store preliminary result away.
C--
          PRELIMINARY_CL(I_CL) = EVENTS_PRODUCED
C
          WRITE(LUPRT,'(//,'' +'',65(''-''),''+'')')
          WRITE(LUPRT,'('' |    Result of '',F5.2,'' CL calculation '',
     >      30X,''|'')') CL_PROB(I_CL)
          WRITE(LUPRT,'('' |    Limit is  '',G10.3,11X,
     >      30X,''|'')') EVENTS_PRODUCED
          WRITE(LUPRT,'('' +'',65(''-''),''+'')')
C
        ENDDO
C--
C  Now start the MC search for the appropriate confidence levels.
C--
        WRITE(LUPRT,'(//,'' *** MC calculation of CL '',
     >     '' (does include bkgd uncertainty) ***'')')
C
        DO I_CL = 1, N_CL
C--
C  We make an initial guess at the rate to get the right scale
C  for the EVENTS_PRODUCED.
C--
          EVENTS_PRODUCED = PRELIMINARY_CL(I_CL)
C--
C  If we have no observed events, use a rough guess at the overall
C  CL limit.
C--
          IF (EVENTS_PRODUCED.LE.0.0) THEN
            EVENTS_PRODUCED = MU_BKGD
          ENDIF
C--
C  Now verify that with the user.
C--
          WRITE(LUPRT,'(//,'' Initial guess at limit is '',G10.3,
     >            '':  Enter your guess:'')')
     >       EVENTS_PRODUCED
          READ(LUTERM, *) EVENTS_PRODUCED
C--
C  Initialize other bookkeeping.
C--
C Prior guess...
          PREVIOUS_EVENTS_PRODUCED = 0.0  
C Calculated probability
          NEW_PROB = 0.0                  
C ...and previous step.
          PREVIOUS_PROB = 0.0             
          N_EXPERIMENTS = 0
          DELTA_CL = 9999999.
          N_TRIES  = 0
C Avoid run-away search....
          MAX_TRIES = 30        
C Don't dump distribution during search.
          DUMP = -6             
C--
C  Iterate till we converge on an answer.
C--
C Required max change in probability
          CL_FUZZ = 0.0005        
          DO WHILE ((ABS(DELTA_CL).GT.CL_FUZZ .OR.
     >                         (N_TRIES.LT.3) .OR.
     >                         (STEP.GT.5.0)) .AND.
     >                   N_TRIES.LT.MAX_TRIES)
C Store prob from last step
            PREVIOUS_PROB = NEW_PROB    
C and calculate new one.
C
C--
C  Set the number of experiments to do based on how close we are
C  getting...
C--
            OLD_N_EXPERIMENTS = N_EXPERIMENTS
            IF (N_TRIES.LT.2) THEN
              N_EXPERIMENTS = 20000
            ELSEIF(DELTA_CL.LT.0.001) THEN
              N_EXPERIMENTS = 500000
            ELSEIF(DELTA_CL.LT.0.01) THEN
              N_EXPERIMENTS = 200000
            ELSEIF(DELTA_CL.LT.0.02) THEN
              N_EXPERIMENTS = 100000
            ELSE
              N_EXPERIMENTS = 50000
            ENDIF
            N_EXPERIMENTS = MAX(N_EXPERIMENTS, OLD_N_EXPERIMENTS)
C--
C  Now perform the MC calculation, assuming we perform the requisite number
C  of experiments.
C--
            NEW_PROB = MC_POISSON_EXP(N_EXPERIMENTS, N_OBSERVED, 
     >        MU_BKGD, SIGMA_BKGD, EVENTS_PRODUCED,
     >        DUMP)
C--
C  Figure out where to go in terms of events produced.  Calculate a
C  Newton-Ralphson step and keep it within reason.
C  On first go around, have to iterate to get two measured points.
C--
            IF (N_TRIES.EQ.0) THEN
              IF (NEW_PROB.LT.CL_PROB(I_CL)) THEN
                STEP =  -EVENTS_PRODUCED * 0.3
              ELSE
                STEP =   EVENTS_PRODUCED * 0.3
              ENDIF
            ELSE
              IF (PREVIOUS_PROB.EQ.NEW_PROB) THEN
                WRITE (LUPRT, '(/,'' ALERT...not converging'',
     >              ''...enter next step value: '')')
                READ(LUTERM,*) STEP
              ELSE
                STEP =
     >           ((PREVIOUS_EVENTS_PRODUCED - EVENTS_PRODUCED) /
     >            (PREVIOUS_PROB            - NEW_PROB) ) *
     >            (CL_PROB(I_CL) - NEW_PROB)
              ENDIF
            ENDIF
C--
C  Keep the step reasonable.
C--
            IF (ABS(STEP).GT.EVENTS_PRODUCED*0.2) THEN
              STEP = SIGN(EVENTS_PRODUCED*0.2, STEP)
            ENDIF
            IF (ABS(STEP).LT.1.0E-2) THEN
              STEP = SIGN(1.0E-2, STEP)
            ENDIF
            TEMP = EVENTS_PRODUCED + STEP
C--
C  Update book-keeping.
C--
            PREVIOUS_EVENTS_PRODUCED = EVENTS_PRODUCED
            EVENTS_PRODUCED          = TEMP
            PREVIOUS_PROB            = NEW_PROB
C--
C  See how close we are to final answer...
C--
            DELTA_CL = NEW_PROB - CL_PROB(I_CL)
            N_TRIES  = N_TRIES + 1
C--
            WRITE(LUPRT, '(/,'' Iteration Number '', I3,
     >          /,'' Signal rate  '',F7.1,''  Prob is '',G10.2,
     >       ''  Updated Signal rate '',F7.1)')
     >          N_TRIES, PREVIOUS_EVENTS_PRODUCED, PREVIOUS_PROB,
     >          EVENTS_PRODUCED
          ENDDO
C--
C  Report results...
C--
          WRITE(LUPRT,'(//,'' +'',65(''-''),''+'')')
          WRITE(LUPRT,'('' |    Result of '',F5.2,'' CL calculation '',
     >      30X,''|'')') CL_PROB(I_CL)
          WRITE(LUPRT,'('' |    Limit is  '',G10.3,11X,
     >      30X,''|'')') EVENTS_PRODUCED
          WRITE(LUPRT,'('' +'',65(''-''),''+'')')
C
        ENDDO
C--
C  OK, we are done...
C--
        STOP
        END
C--
        SUBROUTINE BKGPRB
C--
C  This function calculates the significance of a background
C  fluctiuating to produce an observed signal.
C--
        IMPLICIT NONE
C
        INTEGER         LUPRT, LUTERM
        PARAMETER      (LUPRT = 6)
        PARAMETER      (LUTERM= 5)
C
        INTEGER I
        INTEGER J, IERR, STAT
        INTEGER I_EXP, N_EXP, N_DUMP, I_FAKE
        INTEGER         I_EVT
        INTEGER         N_TAGGABLE, N_TAG_1, N_TAG_2, N_TAG_12
        REAL            SEL_EFF, SIGMA, TAG_1_EFF, TAG_2_EFF
        INTEGER         MAX_EXP
        PARAMETER      (MAX_EXP = 100000)
        REAL            SIGMA_MEASURED(MAX_EXP)
        REAL            SIGMA_COUNTING(MAX_EXP)
        LOGICAL         TAG_1, TAG_2
        LOGICAL DUMP
        REAL            SUM_1, SUM_2
        REAL            MEAN_COUNTING, STDEV_COUNTING
        REAL            MEAN_MEASURED, STDEV_MEASURED
C--
C  External functions
C--
        REAL GAUPRB
        REAL RCDF
C
C  How many experiments
C--
        WRITE(LUPRT, '('' How many experiments & how many dumped?'')')
        READ(LUTERM, *) N_EXP, N_DUMP
C--
C  Set the dump flag on
C--
        DUMP = .TRUE.
C--
C  Read in the starting parameters
C--
        WRITE(LUPRT, '('' Enter produced events, sel eff, tag effs'')')
        READ(LUTERM,*) SIGMA, SEL_EFF, TAG_1_EFF, TAG_2_EFF
C
        DO I_EXP=1, N_EXP
C
          DUMP = (I_EXP.LE.N_DUMP)
C--
C  Now first calculate the number of taggable events.
C--
          CALL POISSN(SIGMA*SEL_EFF, N_TAGGABLE, IERR)
C--
C  Now generate the number of tagged events...
C--
          N_TAG_12 = 0
          N_TAG_1 = 0
          N_TAG_2 = 0
          DO I_EVT = 1, N_TAGGABLE
            TAG_1 = (RCDF(1).LT.TAG_1_EFF) 
            TAG_2 = (RCDF(1).LT.TAG_2_EFF)
            IF (TAG_1.AND.TAG_2) THEN
              N_TAG_12 = N_TAG_12 + 1
            ENDIF
            IF (TAG_1) THEN
              N_TAG_1 = N_TAG_1 + 1
            ENDIF
            IF (TAG_2) THEN
              N_TAG_2 = N_TAG_2 + 1
            ENDIF
          ENDDO
C--
C  Now calculate the expected cross section
C--
          SIGMA_MEASURED(I_EXP) = 
     >          (N_TAG_1 + N_TAG_2)/(TAG_1_EFF + TAG_2_EFF)
     >          / SEL_EFF
C--
C  Do it another way...
C--
          SIGMA_COUNTING(I_EXP) =
     >          (N_TAG_1 + N_TAG_2 - N_TAG_12)/
     >          (TAG_1_EFF + TAG_2_EFF - TAG_1_EFF*TAG_2_EFF) 
     >          / SEL_EFF
C--
C  Write out the result...
C--
          IF (DUMP) THEN
            WRITE(LUPRT, '('' # taggable, tagged (1,2,12), sigma_obs '',
     >          4I4, 2F10.2)') N_TAGGABLE, N_TAG_1, N_TAG_2, N_TAG_12,
     >                       SIGMA_MEASURED(I_EXP),
     >                       SIGMA_COUNTING(I_EXP)
          ENDIF
C
        ENDDO
C--
C  Now analyse the results.
C--
        SUM_1 = 0.0
        SUM_2 = 0.0
        DO I_EXP=1, N_EXP
          SUM_1 = SUM_1 + SIGMA_COUNTING(I_EXP)
          SUM_2 = SUM_2 + SIGMA_COUNTING(I_EXP)**2
        ENDDO
        MEAN_COUNTING = SUM_1/N_EXP
        STDEV_COUNTING = SQRT(SUM_2/N_EXP - MEAN_COUNTING**2)
C
        SUM_1 = 0.0
        SUM_2 = 0.0
        DO I_EXP=1, N_EXP
          SUM_1 = SUM_1 + SIGMA_MEASURED(I_EXP)
          SUM_2 = SUM_2 + SIGMA_MEASURED(I_EXP)**2
        ENDDO
        MEAN_MEASURED = SUM_1/N_EXP
        STDEV_MEASURED = SQRT(SUM_2/N_EXP - MEAN_MEASURED**2)
C
        WRITE(LUPRT,'(/,'' Mean and Sigma (using tags) = '',2G14.3)')
     >          MEAN_MEASURED, STDEV_MEASURED
        WRITE(LUPRT,'(/,'' Mean and Sigma (using evts) = '',2G14.3)')
     >          MEAN_COUNTING, STDEV_COUNTING
        RETURN
        END
        INTEGER FUNCTION MC_POISSON_EXP()
        MC_POISSON_EXP = 0
        RETURN
        END
