        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--
C Confidence levels to calculate...
        CL_PROB(1) = 0.05       
        CL_PROB(2) = 0.10
        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
        REAL    MUWBBB, MUWCCB, MUWC, MUFAKE(3)
        REAL    SUWBBB, SUWCCB, SUWC, SUFAKE(3)
        REAL    EFWBBB(3), EFWCCB(3), EFWC(3)
        REAL    SFWBBB(3), SFWCCB(3), SFWC(3)
C
        REAL    RANNXT
        REAL    SMEARED_BKGD
        INTEGER NWBBB, NWCCB, NWC
        INTEGER NWFAKE(3)
C
        INTEGER I, N_TAGS, N_SINGLE_TAGS, N_DOUBLE_TAGS, N_TAGGED_EVENTS
        INTEGER J, IERR, STAT
        INTEGER I_EXP, N_EXP, N_DUMP, I_FAKE
        INTEGER N_TOTAL_WBBB,           N_TOTAL_WCCB, N_TOTAL_WC
        INTEGER N_TOTAL_FAKE
        INTEGER N_TOTAL_WBBB_TAGS,      N_TOTAL_WCCB_TAGS
        INTEGER N_TOTAL_WC_TAGS,        N_TOTAL_FAKE_TAGS
        INTEGER N_TOTAL_WBBB_DTAGS,     N_TOTAL_WCCB_DTAGS
        INTEGER N_TOTAL_WC_DTAGS,       N_TOTAL_FAKE_DTAGS
        INTEGER TALLY(2,4), ITALLY, JTALLY
        INTEGER N_TOTAL_SIG_TAGS, N_TOTAL_SIG_DOUBLE_TAGS
        INTEGER N_TOTAL_SIG_TAGGED_EVENTS
        INTEGER N_SVX_SLT_DIL_OBSERVED(3)
        INTEGER N_SVX_SLT_DIL_DOUBLES_OBSERVED(3)
        INTEGER N_TAGS_OBSERVED, N_DOUBLE_TAGS_OBSERVED
        INTEGER N_TAGGED_EVENTS_OBSERVED
        LOGICAL DUMP
C
        REAL    FREQ_TAGS, FREQ_DOUBLE_TAGS, FREQ_TAGGED_EVENTS
C
        LOGICAL LCHANN(3)
C--
C  External functions
C--
        REAL GAUPRB
        REAL RCDF
C
        WRITE(LUPRT,'('' Which channels [T/F]: SVX/SLT/DIL...'')')
        READ(LUTERM,*) LCHANN(1), LCHANN(2), LCHANN(3)
C
C  Get the real rates and uncertainties for the physical backgrounds.
C  Also the efficiencies
        WRITE(LUPRT,'('' Enter mean Wbbar, Wccbar, Wc event rates:'')')
        READ(LUTERM,*) MUWBBB, MUWCCB, MUWC
        WRITE(LUPRT,'('' Enter sigma Wbbar, Wccbar, Wc event rates'')')
        READ(LUTERM,*) SUWBBB, SUWCCB, SUWC
C
        IF (LCHANN(1)) THEN
          WRITE(LUPRT, '('' SVX Wbbbar tag efficiency + error (%)'')')
          READ(LUTERM,*) EFWBBB(1), SFWBBB(1)
          WRITE(LUPRT, '('' SVX Wccbar tag efficiency + error (%)'')')
          READ(LUTERM,*) EFWCCB(1), SFWCCB(1)
          WRITE(LUPRT, '('' SVX Wc tag efficiency + error (%)'')')
          READ(LUTERM,*) EFWC(1), SFWC(1)
        ENDIF
        IF (LCHANN(2)) THEN
          WRITE(LUPRT, '('' SLT Wbbbar tag efficiency + error (%)'')')
          READ(LUTERM,*) EFWBBB(2), SFWBBB(2)
          WRITE(LUPRT, '('' SLT Wccbar tag efficiency + error (%)'')')
          READ(LUTERM,*) EFWCCB(2), SFWCCB(2)
          WRITE(LUPRT, '('' SLT Wc tag efficiency + error (%)'')')
          READ(LUTERM,*) EFWC(2), SFWC(2)
        ENDIF
C
C  Now get the fake rates (in terms of the percentage of jets tagged
C  as fake)
        WRITE(LUPRT,'('' Enter mean fake rates  for SVX/SLT/DIL '')')
        READ(LUTERM,*) MUFAKE(1), MUFAKE(2), MUFAKE(3)
        WRITE(LUPRT,'('' Enter sigma fake rates for SVX/SLT/DIL '')')
        READ(LUTERM,*) SUFAKE(1), SUFAKE(2), SUFAKE(3)
C
C  How many experiments
C--
        WRITE(LUPRT, '('' How many experiments & how many dumped?'')')
        READ(LUTERM, *) N_EXP, N_DUMP
C--
C  Now read in the observations
C--
        WRITE(LUPRT, '('' Enter Observed SVX, SLT and DIL tags observed:
     &'')')
        READ(LUTERM, *) N_SVX_SLT_DIL_OBSERVED
        WRITE(LUPRT, '('' Enter Observed SVX-SLT SVX-DIL and SLT-DIL '',
     >                  ''double tags.'')')
        READ(LUTERM, *) N_SVX_SLT_DIL_DOUBLES_OBSERVED
C
        N_TAGS_OBSERVED                 = 0
        N_DOUBLE_TAGS_OBSERVED          = 0
        N_TAGGED_EVENTS_OBSERVED        = 0
        IF (LCHANN(1)) THEN
          N_TAGS_OBSERVED = N_TAGS_OBSERVED + N_SVX_SLT_DIL_OBSERVED(1)
          IF (LCHANN(2)) THEN
            N_DOUBLE_TAGS_OBSERVED = N_DOUBLE_TAGS_OBSERVED + 
     >                  N_SVX_SLT_DIL_DOUBLES_OBSERVED(1)
          ENDIF
          IF (LCHANN(3)) THEN
            N_DOUBLE_TAGS_OBSERVED = N_DOUBLE_TAGS_OBSERVED + 
     >                  N_SVX_SLT_DIL_DOUBLES_OBSERVED(2)
          ENDIF
        ENDIF
        IF (LCHANN(2)) THEN
          N_TAGS_OBSERVED = N_TAGS_OBSERVED + N_SVX_SLT_DIL_OBSERVED(2)
          IF (LCHANN(3)) THEN
            N_DOUBLE_TAGS_OBSERVED = N_DOUBLE_TAGS_OBSERVED + 
     >                  N_SVX_SLT_DIL_DOUBLES_OBSERVED(3)
          ENDIF
        ENDIF
        IF (LCHANN(3)) THEN
          N_TAGS_OBSERVED = N_TAGS_OBSERVED + N_SVX_SLT_DIL_OBSERVED(3)
        ENDIF
C
        N_TAGGED_EVENTS_OBSERVED = N_TAGS_OBSERVED - N_DOUBLE_TAGS_OBSER
     &VED
C--
C  Report what we are about to do...
C--
        WRITE(LUPRT, '(//,'' Number of observed tags, double tags and'',
     >          '' tagged events = '',3I6)')
     >          N_TAGS_OBSERVED, N_DOUBLE_TAGS_OBSERVED, 
     >          N_TAGGED_EVENTS_OBSERVED
C--
C  Now start the MC experiments
C--
        N_TOTAL_WBBB            = 0
        N_TOTAL_WCCB            = 0
        N_TOTAL_WC              = 0
        N_TOTAL_FAKE            = 0
        N_TOTAL_WBBB_TAGS               = 0
        N_TOTAL_WCCB_TAGS               = 0
        N_TOTAL_WC_TAGS                 = 0
        N_TOTAL_FAKE_TAGS               = 0
        N_TOTAL_WBBB_DTAGS              = 0
        N_TOTAL_WCCB_DTAGS              = 0
        N_TOTAL_WC_DTAGS                = 0
        N_TOTAL_FAKE_DTAGS              = 0
C
        N_TOTAL_SIG_TAGS                = 0
        N_TOTAL_SIG_DOUBLE_TAGS         = 0
        N_TOTAL_SIG_TAGGED_EVENTS       = 0
C--
C  Set the dump flag on
C--
        DUMP = .TRUE.
C
        DO I_EXP=1, N_EXP
C
          DUMP = (I_EXP.LE.N_DUMP)
C--
C  Now first generate the number of produced background events.
C--
          N_SINGLE_TAGS   = 0
          N_DOUBLE_TAGS   = 0
          N_TAGGED_EVENTS = 0
C  Number of produced Wbbbar
          NWBBB = 0             
C  Number of produced Wccbar
          NWCCB = 0             
C  Number of produced Wc 
          NWC   = 0             
C  Number of fake SVX's
          NWFAKE(1) = 0         
C  Number of fake SLT's
          NWFAKE(2) = 0         
C  Number of fake dileptons
          NWFAKE(3) = 0         
C
          IF (LCHANN(1).OR.LCHANN(2)) THEN
C--
C  Now the physical backgrounds, Wbb, Wcc, and Wc
C--
            IF (SUWBBB.GT.0.0) THEN
              SMEARED_BKGD = -1.0
              DO WHILE(SMEARED_BKGD.LE.0.0) 
                SMEARED_BKGD = GAUPRB(MUWBBB, SUWBBB)
              ENDDO
            ELSE
              SMEARED_BKGD = MUWBBB
            ENDIF
            CALL POISSN(SMEARED_BKGD, NWBBB, IERR) 
C
            IF (SUWCCB.GT.0.0) THEN
              SMEARED_BKGD = -1.0
              DO WHILE(SMEARED_BKGD.LE.0.0) 
                SMEARED_BKGD = GAUPRB(MUWCCB, SUWCCB)
              ENDDO
            ELSE
              SMEARED_BKGD = MUWCCB
            ENDIF
            CALL POISSN(SMEARED_BKGD, NWCCB, IERR)      
C
            IF (SUWC.GT.0.0) THEN
              SMEARED_BKGD = -1.0
              DO WHILE(SMEARED_BKGD.LE.0.0) 
                SMEARED_BKGD = GAUPRB(MUWC, SUWC)
              ENDDO
            ELSE
              SMEARED_BKGD = MUWC  
            ENDIF
            CALL POISSN(SMEARED_BKGD, NWC, IERR)      
C--
            IF (DUMP) THEN
              WRITE(LUPRT, '(/,'' Experiment '',I7,'' N Wbb, Wcc, Wc = '
     &',
     >                  3I4)') I_EXP, NWBBB, NWCCB, NWC
            ENDIF
C--
C  We now have the number of produced physical backgrounds.  Find out
C  which ones get tagged and double tagged.
C--
C  Keep track of how many Wbb, Wcc, Wc and fakes are tagged this experiment
C--
            DO ITALLY = 1,4
C Number of single tags
              TALLY(1,ITALLY) = 0       
C Number of double tags
              TALLY(2,ITALLY) = 0       
            ENDDO
C
            DO I=1, NWBBB
              STAT = 0
              IF (LCHANN(1)) THEN   
                RANNXT = 100.0*RCDF(1)
                IF (RANNXT.LE. EFWBBB(1)) STAT = 1
              ENDIF
              IF (LCHANN(2)) THEN
                RANNXT = 100.0*RCDF(1)
                IF (RANNXT.LE. EFWBBB(2)) STAT = STAT + 2
              ENDIF
C
              IF (STAT.EQ.1 .OR. STAT.EQ.2) THEN
                N_SINGLE_TAGS = N_SINGLE_TAGS + 1
                TALLY(1,1) = TALLY(1,1) + 1
              ELSEIF (STAT.EQ.3) THEN
                N_DOUBLE_TAGS = N_DOUBLE_TAGS + 1
                TALLY(2,1) = TALLY(2,1) + 1
              ENDIF
            ENDDO
C           
            DO I=1, NWCCB
              STAT = 0
              IF (LCHANN(1)) THEN   
                RANNXT = 100.0*RCDF(1)
                IF (RANNXT.LE. EFWCCB(1)) STAT = 1
              ENDIF
              IF (LCHANN(2)) THEN
                RANNXT = 100.0*RCDF(1)
                IF (RANNXT.LE. EFWCCB(2)) STAT = STAT + 2
              ENDIF
C
              IF (STAT.EQ.1 .OR. STAT.EQ.2) THEN
                N_SINGLE_TAGS = N_SINGLE_TAGS + 1
                TALLY(1,2) = TALLY(1,2) + 1
              ELSEIF (STAT.EQ.3) THEN
                N_DOUBLE_TAGS = N_DOUBLE_TAGS + 1
                TALLY(2,2) = TALLY(2,2) + 1
              ENDIF
            ENDDO
C           
            DO I=1, NWC
              STAT = 0
              IF (LCHANN(1)) THEN   
                RANNXT = 100.0*RCDF(1)
                IF (RANNXT.LE. EFWC(1)) STAT = 1
              ENDIF
              IF (LCHANN(2)) THEN
                RANNXT = 100.0*RCDF(1)
                IF (RANNXT.LE. EFWC(2)) STAT = STAT + 2
              ENDIF
C
              IF (STAT.EQ.1 .OR. STAT.EQ.2) THEN
                N_SINGLE_TAGS = N_SINGLE_TAGS + 1
                TALLY(1,3) = TALLY(1,3) + 1
              ELSEIF (STAT.EQ.3) THEN
                N_DOUBLE_TAGS = N_DOUBLE_TAGS + 1
                TALLY(2,3) = TALLY(2,3) + 1
              ENDIF
            ENDDO
C  End of special treatment for SVX/SLT
          ENDIF         
C--
C  Now determine the number of fakes from mistags.  These are uncorrelated
C  with the physical backgrounds.
C--       
C Single tags from fakes
          TALLY(1,4) = 0        
C Double tags from fakes
          TALLY(2,4) = 0        
          IF (LCHANN(1)) THEN
            IF (SUFAKE(1).GT.0.0) THEN
              SMEARED_BKGD = -1.0
              DO WHILE(SMEARED_BKGD.LE.0.0) 
                SMEARED_BKGD = GAUPRB(MUFAKE(1), SUFAKE(1))
              ENDDO
            ELSE
              SMEARED_BKGD = MUFAKE(1)
            ENDIF
            CALL POISSN(SMEARED_BKGD, NWFAKE(1), IERR) 
          ENDIF
C
          IF (LCHANN(2)) THEN
            IF (SUFAKE(2).GT.0.0) THEN
              SMEARED_BKGD = -1.0
              DO WHILE(SMEARED_BKGD.LE.0.0) 
                SMEARED_BKGD = GAUPRB(MUFAKE(2), SUFAKE(2))
              ENDDO
            ELSE
              SMEARED_BKGD = MUFAKE(2)
            ENDIF
            CALL POISSN(SMEARED_BKGD, NWFAKE(2), IERR) 
          ENDIF
C--
C  Now assume that the fake probabilities are uncorrelated...
C--
          DO I_FAKE = 1, NWFAKE(1)+NWFAKE(2)
            RANNXT = RCDF(1)
            IF (LCHANN(1) .AND. MUFAKE(1)/52.0 .GT. RANNXT) THEN
              N_DOUBLE_TAGS = N_DOUBLE_TAGS + 1
              TALLY(2,4) = TALLY(2,4) + 1
            ELSEIF (LCHANN(2) .AND.   
     >                          RANNXT .GT. 1.0-MUFAKE(2)/52.0) THEN
              N_DOUBLE_TAGS = N_DOUBLE_TAGS + 1
              TALLY(2,4) = TALLY(2,4) + 1
            ELSE
              N_SINGLE_TAGS = N_SINGLE_TAGS + 1
              TALLY(1,4) = TALLY(1,4) + 1
            ENDIF
          ENDDO
C
          IF (LCHANN(3)) THEN
            IF (SUFAKE(3).GT.0.0) THEN
              SMEARED_BKGD = -1.0
              DO WHILE(SMEARED_BKGD.LE.0.0) 
                SMEARED_BKGD = GAUPRB(MUFAKE(3), SUFAKE(3))
              ENDDO
            ELSE
              SMEARED_BKGD = MUFAKE(3)
            ENDIF
            CALL POISSN(SMEARED_BKGD, NWFAKE(3), IERR) 
C--
C  Now overlaps...
            DO I_FAKE = 1, NWFAKE(3)
              RANNXT = RCDF(1)
              IF (LCHANN(1) .AND. MUFAKE(1)/52.0 .GT. RANNXT) THEN
                N_DOUBLE_TAGS = N_DOUBLE_TAGS + 1
                TALLY(2,4) = TALLY(2,4) + 1
              ELSEIF (LCHANN(2) .AND. 
     >                          RANNXT .GT. 1.0 - MUFAKE(2)/52.0) THEN
                N_DOUBLE_TAGS = N_DOUBLE_TAGS + 1
                TALLY(2,4) = TALLY(2,4) + 1
              ELSE
                N_SINGLE_TAGS = N_SINGLE_TAGS + 1
                TALLY(1,4) = TALLY(1,4) + 1
              ENDIF
            ENDDO
          ENDIF
C--
C  Now accumulate all the statistics and report if requested.
C--
          N_TOTAL_WBBB = N_TOTAL_WBBB + NWBBB
          N_TOTAL_WCCB = N_TOTAL_WCCB + NWCCB
          N_TOTAL_WC   = N_TOTAL_WC   + NWC
          N_TOTAL_WBBB_TAGS = N_TOTAL_WBBB_TAGS + TALLY(1,1)
          N_TOTAL_WCCB_TAGS = N_TOTAL_WCCB_TAGS + TALLY(1,2)
          N_TOTAL_WC_TAGS   = N_TOTAL_WC_TAGS   + TALLY(1,3)
          N_TOTAL_FAKE_TAGS = N_TOTAL_FAKE_TAGS + TALLY(1,4)
          N_TOTAL_WBBB_DTAGS = N_TOTAL_WBBB_DTAGS + TALLY(2,1)
          N_TOTAL_WCCB_DTAGS = N_TOTAL_WCCB_DTAGS + TALLY(2,2)
          N_TOTAL_WC_DTAGS   = N_TOTAL_WC_DTAGS   + TALLY(2,3)
          N_TOTAL_FAKE_DTAGS = N_TOTAL_FAKE_DTAGS + TALLY(2,4)
C--
C  Test to see how often the most interesting combinations arise...
C--
          N_TAGGED_EVENTS = N_SINGLE_TAGS + N_DOUBLE_TAGS
          N_TAGS = N_SINGLE_TAGS + 2*N_DOUBLE_TAGS
C
          IF (N_TAGS.GE.N_TAGS_OBSERVED) THEN
            N_TOTAL_SIG_TAGS = N_TOTAL_SIG_TAGS + 1
          ENDIF
          IF (N_TAGGED_EVENTS.GE.N_TAGGED_EVENTS_OBSERVED) THEN
            N_TOTAL_SIG_TAGGED_EVENTS = N_TOTAL_SIG_TAGGED_EVENTS + 1
          ENDIF
          IF (N_DOUBLE_TAGS .GE.
     >                   N_DOUBLE_TAGS_OBSERVED) THEN
            N_TOTAL_SIG_DOUBLE_TAGS = 
     >                  N_TOTAL_SIG_DOUBLE_TAGS + 1
          ENDIF
C
          IF (DUMP) THEN
            WRITE(LUPRT, '('' Single tags = '',4I4,''  Double tags = '',
     >                  4I4)') 
     >                          (TALLY(1,ITALLY), ITALLY=1,4),
     >                          (TALLY(2,JTALLY), JTALLY=1,4)
          ENDIF
        ENDDO
C--
C  Now report a summary and a significance
C--
        WRITE(LUPRT, '(/, '' Total Wbb, Wcc, Wc, fakes = '',4I10,
     >                 /, '' Total tags                = '',4I10,
     >                 /, '' Total double tags         = '',4I10)')
     >  N_TOTAL_WBBB, N_TOTAL_WCCB, N_TOTAL_WC, N_TOTAL_FAKE,
     >  N_TOTAL_WBBB_TAGS, N_TOTAL_WCCB_TAGS, N_TOTAL_WC_TAGS, 
     >  N_TOTAL_FAKE_TAGS,
     >  N_TOTAL_WBBB_DTAGS, N_TOTAL_WCCB_DTAGS, N_TOTAL_WC_DTAGS,
     >  N_TOTAL_FAKE_DTAGS
        WRITE(LUPRT, '(/, '' Frequency tags            = '',4G10.3,
     >                 /, '' Frequency double tags     = '',4G10.3)')
     >      FLOAT(N_TOTAL_WBBB_TAGS)/FLOAT(MAX(1, N_TOTAL_WBBB)),
     >      FLOAT(N_TOTAL_WCCB_TAGS)/FLOAT(MAX(1, N_TOTAL_WCCB)),
     >      FLOAT(N_TOTAL_WC_TAGS  )/FLOAT(MAX(1, N_TOTAL_WC)),
     >      FLOAT(N_TOTAL_FAKE_TAGS)/FLOAT(MAX(1, N_EXP)),
     >      FLOAT(N_TOTAL_WBBB_DTAGS)/FLOAT(MAX(1, N_TOTAL_WBBB)),
     >      FLOAT(N_TOTAL_WCCB_DTAGS)/FLOAT(MAX(1, N_TOTAL_WCCB)),
     >      FLOAT(N_TOTAL_WC_DTAGS  )/FLOAT(MAX(1, N_TOTAL_WC)),
     >      FLOAT(N_TOTAL_FAKE_DTAGS)/FLOAT(MAX(1, N_EXP))
C
        WRITE(LUPRT, '(/, '' Average tags/event = '',G10.3)')
     >  FLOAT(N_TOTAL_WBBB_TAGS+N_TOTAL_WCCB_TAGS+
     >        N_TOTAL_WC_TAGS + N_TOTAL_FAKE_TAGS+
     >     2*(N_TOTAL_WBBB_DTAGS+N_TOTAL_WCCB_DTAGS+
     >        N_TOTAL_WC_DTAGS + N_TOTAL_FAKE_DTAGS))/FLOAT(N_EXP)
        WRITE(LUPRT, '(/, '' Physics tags/event = '',G10.3)')
     >  FLOAT(N_TOTAL_WBBB_TAGS+N_TOTAL_WCCB_TAGS+
     >        N_TOTAL_WC_TAGS +
     >     2*(N_TOTAL_WBBB_DTAGS+N_TOTAL_WCCB_DTAGS+
     >        N_TOTAL_WC_DTAGS))/FLOAT(N_EXP)
        WRITE(LUPRT, '(/, '' Fake    tags/event = '',G10.3)')
     >  FLOAT(N_TOTAL_FAKE_TAGS+2*N_TOTAL_FAKE_DTAGS)/FLOAT(N_EXP)
C
        WRITE(LUPRT, '(/,
     >          '' Frequency of .GE. '', I4, '' tags         = '',G10.3,
     >      /,  '' Frequency of .GE. '', I4, '' double-tags  = '',G10.3,
     >      /,  '' Frequency of .GE. '', I4, '' tagged evnts = '',G10.3)
     &')
     >  N_TAGS_OBSERVED, 
     >  FLOAT(N_TOTAL_SIG_TAGS)/FLOAT(MAX(1, N_EXP)),
     >  N_DOUBLE_TAGS_OBSERVED, 
     >  FLOAT(N_TOTAL_SIG_DOUBLE_TAGS)/FLOAT(MAX(1, N_EXP)),
     >  N_TAGGED_EVENTS_OBSERVED, 
     >  FLOAT(N_TOTAL_SIG_TAGGED_EVENTS)/FLOAT(MAX(1, N_EXP))
C
        RETURN
        END
C--
        REAL FUNCTION MC_POISSON_EXP(N_EXPERIMENTS, NEVNTS, 
     >        MU_BKGD, SIGMA_BKGD, EVENTS_PRODUCED,
     >        DUMP)
C--
C  This function calculates the probability that we would see
C  NEVNTS given that we expect a background process with mean rate
C  MU_BKGD and systematic uncertainty SIGMA_BKGD and a signal process
C  with mean rate EVENTS_PRODUCED.  The programme simulates N_EXPERIMENTS
C  experiments.  
C
C  DUMP (Integer*4)     is logical unit number to which resulting prob
C                       distribution is dumped to.  If 0, minimal
C                       information is generated.  If -ve, moderate
C                       amount of info is dumped in unit ABS(DUMP).
C
C--
        IMPLICIT NONE
C--
C  Input variables
C--
        INTEGER         N_EXPERIMENTS
        INTEGER         NEVNTS
        INTEGER         DUMP
        REAL            MU_BKGD, SIGMA_BKGD
        REAL            EVENTS_PRODUCED
C--
C  Local variables
C--
        INTEGER         MAXBIN
C Size of histogram of dist.
        PARAMETER      (MAXBIN = 1000)    
        REAL            THROWN_BKGD
        REAL            PROBABILITY
        REAL            SMEARED_BKGD
        INTEGER         I,      J
        INTEGER         I_EXPERIMENT
        INTEGER         FREQUENCY
        INTEGER         HIST(0:MAXBIN)
        INTEGER         LUPRT
        INTEGER         EVENTS_THROWN
        INTEGER         N_WARNINGS
        INTEGER         IERR
        LOGICAL         LDUMP,  LDEBUG
C--
C  External references.
C--
C CERNLIB Gaussian generator
        REAL            GAUPRB  
C--
C  Initial the Monte Carlo integration.
C--
        MC_POISSON_EXP = 0.0
        LUPRT = ABS(DUMP)
        LDUMP = (DUMP.NE.0)
        LDEBUG= (DUMP.LT.0)
        FREQUENCY = 0
        DO J=0, MAXBIN
          HIST(J) = 0
        ENDDO
        N_WARNINGS = 0
C
        IF (LDUMP) WRITE(LUPRT,'(/,'' MC_POISSON_EXP: Events '',F8.1,
     >                  '' will be produced in '',I8,'' experiments'')')
     >                  EVENTS_PRODUCED, N_EXPERIMENTS
C--
C  Now do a bunch of experiments, each time having a new acceptance
C  for each channel, and a new number of observed events for each
C  channel.
C--
        DO I_EXPERIMENT = 1, N_EXPERIMENTS
C--
C  First smear the expected background rate.
C--
          IF (SIGMA_BKGD.GT.0.000001) THEN
            SMEARED_BKGD = GAUPRB(MU_BKGD, SIGMA_BKGD)
            IF (SMEARED_BKGD.LE.0.0) THEN
              N_WARNINGS = N_WARNINGS + 1
              IF (N_WARNINGS.LT.5) THEN
                WRITE(LUPRT, '('' *** WARNING *** Smeared BKGD <= 0 *** 
     &'')')
              ELSEIF (N_WARNINGS.EQ.5) THEN
                WRITE(LUPRT, '('' *** Further warnings suppressed   *** 
     &'')')
              ENDIF
              SMEARED_BKGD = 0.0
            ENDIF
          ELSE
            SMEARED_BKGD = MU_BKGD
          ENDIF
C--
C  Calculate expected mean number of events then.
C--
          CALL POISSN(SMEARED_BKGD + EVENTS_PRODUCED, 
     >                  EVENTS_THROWN, IERR)
          IF (EVENTS_THROWN.LE.NEVNTS) THEN
              FREQUENCY = FREQUENCY + 1
          ENDIF
          J = MIN(EVENTS_THROWN, MAXBIN)
          HIST(J) = HIST(J) + 1
C--
C  Dump a little info...
C--
          IF (LDEBUG .AND. I_EXPERIMENT.LT.5) THEN
            WRITE(LUPRT, '(1X,''Smeared BKG '',F6.0,''  Events thrown ''
     &, 
     >          I6, ''  Cum freq '',I8)')
     >         SMEARED_BKGD, EVENTS_THROWN, FREQUENCY 
          ENDIF
C
        ENDDO
C--
        PROBABILITY = FLOAT(FREQUENCY)/
     >                        FLOAT(N_EXPERIMENTS)
        MC_POISSON_EXP = PROBABILITY
C
        IF (LDUMP) THEN
C--
C  Dump the individual frequency distributions.
C--
          IF (.NOT.LDEBUG) THEN
            WRITE(LUPRT,'(/,'' Histogram of event distributions'',
     >                    /,'' --------------------------------'')')
            WRITE(LUPRT,'(1X,'' Events  '',''Chan '')')
C
            DO J=1, MAXBIN
              WRITE(LUPRT,'(1X,I6,2X,10(F8.4,2X))')
     >        J, FLOAT(HIST(J))/FLOAT(N_EXPERIMENTS)
            ENDDO
          ENDIF
        ENDIF
C
        RETURN
        END
