C DEC/CMS REPLACEMENT HISTORY, Element POISSON.CDF
C *4    12-DEC-1991 23:16:49 GRINDLEY "Split off functions into separate files."
C *3    13-JUL-1991 17:53:03 PEKKA "Change include references"
C *2    22-MAY-1991 04:34:19 PEKKA "Add $$Alias's."
C *1    12-FEB-1990 14:22:15 PEKKA "Convolutes Poisson and Gaussian uncertaintie
C      s to calculate limits."                                                  
C DEC/CMS REPLACEMENT HISTORY, Element POISSON.CDF
        PROGRAM POISSON
C--
C  This program calculates the combined 90 and 95% confidence level
C  upper limits for multiple channels, assuming:
C
C    i)  for each mutually exclusive channel, an efficiency for t-tbar
C        is known.
C   ii)  A known number of events is found for each channel.
C  iii)  There are correlated Gaussian uncertainties (ie uncertainties
C        common to each selection like integrated luminosity).
C   iv)  Each channel has additional uncertainties that are unique to
C        that channel.
C--
C  The program prompts the user for the necessary information.
C--
C       Pekka K. Sinervo
C       University of Pennsylvania and CDF
C
C       November 6, 1989
C--
C
        IMPLICIT NONE
C
        INTEGER         MXCHAN, NCHAN
        PARAMETER      (MXCHAN = 10)
C
        INTEGER         NEVNTS(MXCHAN)
        REAL            CORSYS         
        REAL            UNCSYS           (MXCHAN)
        REAL            ACCEPTANCE(MXCHAN)
        INTEGER         DUMP
        INTEGER         LUPRT, LUTERM
        PARAMETER      (LUPRT = 6)
        PARAMETER      (LUTERM= 5)
        INTEGER         I
        INTEGER         NCL 
        PARAMETER      (NCL  = 2)
        REAL            CLPROB (NCL )
C
        INTEGER         ICL 
        REAL            EVEPRO         ,        CLFUZZ 
        REAL            PREVPR                  
        REAL            TOTACC          
        REAL            TEMP
        REAL            DELTCL  
        REAL            NEWPRO  ,       PREPRO       
        REAL            STEP
        INTEGER         NTRIES ,        MAXTRI   
        INTEGER         NEXPER       ,  OLNEXP           
        REAL            MASLIM    
        LOGICAL         PERMUTE
        LOGICAL         USETOT   
C--
C  External functions referenced
C--
        REAL            CALCCL
        REAL            TOPLIMIT
C--
C  Initialize...
C--
C Confidence levels to calculate...
        CLPROB (1) = 0.05       
        CLPROB (2) = 0.10
C--
C  Get the information from the user.
C--
        WRITE(LUPRT, '('' Enter # of channels being considered: '')')
        READ(LUTERM, *) NCHAN
        IF (NCHAN.LT.1 .OR. NCHAN.GT.MXCHAN) THEN
          WRITE(LUPRT, '('' Illegal number of channels--'',
     >                   ''have to be between 1 and '',I3)') MXCHAN
          STOP 4
        ENDIF
C--
C  Get the correlated systematic uncertainty on the acceptance.
C--
        WRITE(LUPRT, '('' Correlated systematic error (%)?'')')
        READ(LUTERM, *) CORSYS         
        CORSYS          = CORSYS          * 0.01
C--
C  Do we consider permutations of the order of measurements
C  in calculating the upper limit (PERMUTE=.TRUE.)?
C--
        WRITE(LUPRT, '('' Do we consider permutations on the '',
     >     ''order of measurements? (T/F)'')')
        READ(LUTERM, *) PERMUTE
C--
C  Do we only use the overall rate to set the limits.
C--
        WRITE(LUPRT, '('' Do we use the total rate '',
     >     ''to set the limit? (T/F)'')')
        READ(LUTERM, *) USETOT   
C--
C  Now get the information specific to each channel.
C--
        WRITE(LUPRT, '('' # of events for each channel?'')')
        READ(LUTERM, *) (NEVNTS(I), I=1,NCHAN)
        WRITE(LUPRT, '('' Uncorrelated fractional uncertainty '',
     >         ''on the acceptance for each channel (%)?'')')
        READ(LUTERM, *) (UNCSYS           (I),I=1,NCHAN)
        WRITE(LUPRT, '('' Acceptance for '',
     >                  ''each channel (%)?'')')
        READ(LUTERM, *) (ACCEPTANCE(I),I=1,NCHAN)
C--
C  Scale acceptance from  % to fractions.
C--
        DO I=1,NCHAN
          UNCSYS           (I) = UNCSYS           (I) * 0.01
          ACCEPTANCE(I) = ACCEPTANCE(I) * 0.01
        ENDDO
C--
C  Now we go the data--dump it back to the user...
C--
        WRITE(LUPRT, '( //,'' Calculate Poisson CL for '',I3,
     >                                             '' channels'',
     >                  //,'' Correlated systematic error = '',F6.3,
     >                /,('' Channel '',I3,'' # Events '', I4,
     >                          ''  Uncorrelated Systematic '',F6.3,
     >                          ''  Acceptance '',F8.5))')
     >          NCHAN, CORSYS         ,
     >           (I, NEVNTS(I), UNCSYS           (I), ACCEPTANCE(I),
     >                                          I=1, NCHAN)
        IF (PERMUTE) THEN
           WRITE(LUPRT,'(/,'' This probability calculation will'',
     >       '' allow permutations of the measurements.'')')
        ELSE
           WRITE(LUPRT,'(/,'' This probability calculation will'',
     >       '' not permute the order of the measurements.'')')
        ENDIF
C--
C  Now start the MC search for the appropriate confidence levels.
C--
        DO ICL  = 1, NCL 
C--
C  We make an initial guess at the rate to get the right scale
C  for the EVENTS_PRODUCED.
C--
          EVEPRO          = 0.0
          TOTACC           = 0.0
          DO I=1, NCHAN
            EVEPRO          = EVEPRO          +
     >            NEVNTS(I)/ACCEPTANCE(I)
            TOTACC           = TOTACC           + ACCEPTANCE(I)
          ENDDO
C--
C  If we have no observed events, use a rough guess at the overall
C  CL limit.
C--
          IF (EVEPRO         .LE.0.0) THEN
            EVEPRO          = 3.0/TOTACC          
          ENDIF
C--
C  Now verify that with the user.
C--
          WRITE(LUPRT,'('' Initial guess at limit is '',G10.1,
     >            '':  Enter your guess:'')')
     >       EVEPRO         
          READ(LUTERM, *) EVEPRO         
C--
C  Initialize other bookkeeping.
C--
C Prior guess...
          PREVPR                   = 0.0  
C Calculated probability
          NEWPRO   = 0.0                  
C ...and previous step.
          PREPRO        = 0.0             
          DELTCL   = 9999999.
          NTRIES   = 0
C Avoid run-away search....
          MAXTRI    = 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
          CLFUZZ  = 0.0005        
          DO WHILE ((ABS(DELTCL  ).GT.CLFUZZ  .OR.
     >                         (NTRIES .LT.3) .OR.
     >                         (STEP.GT.5.0)) .AND.
     >                   NTRIES .LT.MAXTRI   )
C Store prob from last step
            PREPRO        = NEWPRO      
C and calculate new one.
C
C--
C  Set the number of experiments to do based on how close we are
C  getting...
C--
            OLNEXP            = NEXPER       
            IF (NTRIES .LT.2) THEN
              NEXPER        = 20000
            ELSEIF(DELTCL  .LT.0.001) THEN
              NEXPER        = 500000
            ELSEIF(DELTCL  .LT.0.01) THEN
              NEXPER        = 200000
            ELSEIF(DELTCL  .LT.0.02) THEN
              NEXPER        = 100000
            ELSE
              NEXPER        = 50000
            ENDIF
            NEXPER        = MAX(NEXPER       , OLNEXP           )
C
            NEWPRO   = CALCCL(
     >          NEXPER       , NCHAN, NEVNTS(1), ACCEPTANCE,
     >          UNCSYS           , CORSYS         , EVEPRO         ,
     >          DUMP, PERMUTE, USETOT   )
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 (NTRIES .EQ.0) THEN
              IF (NEWPRO  .LT.CLPROB (ICL )) THEN
                STEP =  -EVEPRO          * 0.3
              ELSE
                STEP =   EVEPRO          * 0.3
              ENDIF
            ELSE
              IF (PREPRO       .EQ.NEWPRO  ) THEN
                WRITE (LUPRT, '(/,'' ALERT...not converging'',
     >              ''...enter next step value: '')')
                READ(LUTERM,*) STEP
              ELSE
                STEP =
     >           ((PREVPR                   - EVEPRO         ) /
     >            (PREPRO                   - NEWPRO  ) ) *
     >            (CLPROB (ICL ) - NEWPRO  )
              ENDIF
            ENDIF
C--
C  Keep the step reasonable.
C--
            IF (ABS(STEP).GT.EVEPRO         *0.2) THEN
              STEP = SIGN(EVEPRO         *0.2, STEP)
            ENDIF
            IF (ABS(STEP).LT.1.0) THEN
              STEP = SIGN(1., STEP)
            ENDIF
            TEMP = EVEPRO          + STEP
C--
C  Update book-keeping.
C--
            PREVPR                   = EVEPRO         
            EVEPRO                   = TEMP
            PREPRO                   = NEWPRO  
C--
C  See how close we are to final answer...
C--
            DELTCL   = NEWPRO   - CLPROB (ICL )
            NTRIES   = NTRIES  + 1
          ENDDO
C--
C  Report results...
C--
          WRITE(LUPRT,'('' +'',65(''-''),''+'')')
          WRITE(LUPRT, 100)
     >          NCHAN, CORSYS         ,
     >           (I, NEVNTS(I), UNCSYS           (I), ACCEPTANCE(I),
     >                                          I=1, NCHAN)
 100      FORMAT(' | Calculated Poisson CL for ',I3,
     >                                 ' channels',26X,'|',
     >         /,' |',65X,'|',
     >         /,' | Correlated systematic error = ',F6.3,27X,' |',
     >         /,' |',65X,'|',
     >        /,(' | Chan',I3,' Evnts ',I3,
     >                          '  Uncorrelated Syst ',F5.3,
     >                          '  Accep ',F8.5,'      |'))
          IF (USETOT   ) THEN
            WRITE(LUPRT,'('' |'',65X,''|'',/,
     >          '' | This calculation was based on the total rat'',
     >          ''e of events          |'')')
          ELSEIF (PERMUTE) THEN
            WRITE(LUPRT,'('' |'',65X,''|'',/,
     >          '' | This calculation permuted the order  of the'',
     >          '' observations        |'')')
          ELSE
            WRITE(LUPRT,'('' |'',65X,''|'',/,
     >          '' | This calculation did not permute the  order'',
     >          '' of observations     |'')')
          ENDIF
          WRITE(LUPRT,'('' |'',65X,''|'',/,
     >                '' | Rate is '',F8.2,'' for '',
     >              F5.2,''% CL   (Ntries ='',I4,'')'',17X,''|'')')
     >         EVEPRO         , 1.0 - CLPROB (ICL ), NTRIES 
          WRITE(LUPRT,'('' |'',65('' ''),''|'')')
C--
C  Calc the top mass limit from this.  This is using the lower limit
C  of the Altarelli calculation.
C--
          WRITE(LUPRT,'('' +'',65(''-''),''+'')')
          MASLIM     = TOPLIMIT(EVEPRO         )
          WRITE(LUPRT, '('' | Top mass limit is '',F8.1,38X,''|'')')
     >                  MASLIM    
          WRITE(LUPRT,'('' +'',65(''-''),''+'')')
C--
C  Dump the resulting probability distribution.
C--
          DUMP = 6
          NEXPER        = 200000
          NEWPRO   = CALCCL(
     >          NEXPER       , NCHAN, NEVNTS(1), ACCEPTANCE,
     >          UNCSYS           , CORSYS         , EVEPRO         ,
     >          DUMP, PERMUTE, USETOT   )
        ENDDO
C--
C  OK, we are done...
C--
        STOP
        END
