        REAL FUNCTION CALCCL(N_EXPERIMENTS, NCHAN, NEVNTS, ACCEPTANCE,
     >        UNCORRELATED_SYST, CORRELATED_SYST, EVENTS_PRODUCED,
     >        DUMP, PERMUTE, USE_TOTAL)
C--
C  This function calculates the probability that we would see
C  EVNTS(i) events for NCHAN channels, given the acceptance
C  ACCEPTANCE(i), uncorrelated systematic uncertainty of each
C  acceptance UNCORRELATED_SYST(i), the correlated systematic
C  uncertainty CORRELATED_SYST common to all channels, and the
C  supposed total rate of events produced, EVENTS_PRODUCED.
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  PERMUTE (Logical)    Allow permutations of the order of measurements
C                       if .TRUE..
C                       Otherwise, consider each measurement unique.
C
C  USE_TOTAL (Logical)  Use total rate to set limit instead of
C                       individual rates when .TRUE..
C--
        IMPLICIT NONE
C--
C  Input variables
C--
        INTEGER         N_EXPERIMENTS
        INTEGER         NCHAN
        INTEGER         NEVNTS(*)
        INTEGER         DUMP
        REAL            CORRELATED_SYST
        REAL            UNCORRELATED_SYST(*)
        REAL            ACCEPTANCE(*)
        REAL            EVENTS_PRODUCED
        LOGICAL         PERMUTE
        LOGICAL         USE_TOTAL
C--
C  Local variables
C--
C
        INTEGER         MXCHAN
        PARAMETER      (MXCHAN = 10)
        INTEGER         MAXBIN
C Size of histogram of dist.
        PARAMETER      (MAXBIN = 20)    
C Thrown acceptance for channel
        REAL            EPSILON(MXCHAN) 
C How often we throw <=
        INTEGER         FREQUENCY(MXCHAN)       
C # events observed
C
        REAL            PROB_FOR_CHANNEL(MXCHAN)
        INTEGER         N_PASS_PERMUTE
        REAL            PROB_WITH_PERMUTE
        REAL            PROB_WITHOUT_PERMUTE
        INTEGER         I,      J
        INTEGER         I_EXPERIMENT
        REAL            OVERALL_UNCERTAINTY
        INTEGER         EVENTS_THROWN(MXCHAN)
        REAL            MEAN_EVENTS_FOR_CHANNEL(MXCHAN)
        INTEGER         N_OBSERVED(0:MAXBIN,MXCHAN)
        INTEGER         LUPRT
        LOGICAL         LDUMP,  LDEBUG
        LOGICAL         FAILED
        LOGICAL         LSUSPICIOUS
        INTEGER         ITEMP
        INTEGER         IERR
        REAL            RTEMP
CBackground counts
        real            bkgnd(mxchan),sbkgnd(mxchan)    
CUncertainty on bkgnd
        real            bkunc(mxchan)           
C1=Low Lim, 2=Up Lim
        integer         round                   
C
        REAL            TOTAL_RATE_OBSERVED
        REAL            TOTAL_RATE_THROWN
        INTEGER         N_PASS_TOTAL_RATE
        REAL            PROB_USING_TOTAL_RATE
        REAL            BEVENTS_PRODUCED
        REAL            REVENTS_THROWN(MXCHAN)
C--
C  External references.
C--
C CERNLIB Gaussian generator
        REAL            GAUPRB  
C
        common /baksub/ BEVENTS_PRODUCED,bkunc,round,bkgnd
C--
C  Initial the Monte Carlo integration.
C--
        LSUSPICIOUS = .FALSE.
        CALCCL = 0.0
        LUPRT = ABS(DUMP)
        LDUMP = (DUMP.NE.0)
        LDEBUG= (DUMP.LT.0)
        DO I=1, NCHAN
          PROB_FOR_CHANNEL(I) = 0.0
          FREQUENCY(I) = 0
          DO J=0, MAXBIN
            N_OBSERVED(J,I) = 0
          ENDDO
        ENDDO
C
        IF (PERMUTE) N_PASS_PERMUTE = 0
        N_PASS_TOTAL_RATE = 0
C
        IF (LDUMP) WRITE(LUPRT,'(/,'' CALCCL: 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
           do i=1,nchan
                sbkgnd(i)=bkgnd(i)
           end do
C--
C  First get smearing by overall systematic uncertainty.
C--
          IF (CORRELATED_SYST.GT.0.000001) THEN
            OVERALL_UNCERTAINTY = GAUPRB(1.0, CORRELATED_SYST)
          ELSE
            OVERALL_UNCERTAINTY = 1.0
          ENDIF
C--
C  Fold in additional Gaussian systematic smearings to individual
C  channels.
C--
          DO I = 1, NCHAN
            IF (UNCORRELATED_SYST(I).GT.0.000001) THEN
              EPSILON(I) = OVERALL_UNCERTAINTY *
     >                          GAUPRB(1.0, UNCORRELATED_SYST(I))
            ELSE
              EPSILON(I) = OVERALL_UNCERTAINTY
            ENDIF
C
C  Now smear backgrounds 
C
C
            sbkgnd(i)=sbkgnd(i)*gauprb(1.0,bkunc(i))
C
C--
C  Calculate expected mean number of events then.
C--
            EPSILON(I) = EPSILON(I) * ACCEPTANCE(I)
            MEAN_EVENTS_FOR_CHANNEL(I) = EPSILON(I) * EVENTS_PRODUCED
     >      +SBKGND(I)
C--
C  Throw the observed number for this channel, and bin it.
C--
            CALL POISSN(MEAN_EVENTS_FOR_CHANNEL(I), EVENTS_THROWN(I),
     >                                  IERR)
C--
C  Now see if this is less than or equal to the number observed
C  for each channel.  Also, if we are allowing permutations, keep
C  track of the number that satisfy some permutation of the order
C  of the measurements.
C--
           if(round.eq.1) then
            IF (EVENTS_THROWN(I).LE.NEVNTS(I)) THEN
              FREQUENCY(I) = FREQUENCY(I) + 1
            ENDIF
           elseif(round.eq.2) then
            IF (EVENTS_THROWN(I).GE.NEVNTS(I)) THEN
              FREQUENCY(I) = FREQUENCY(I) + 1
            ENDIF
           endif
            J = MIN(EVENTS_THROWN(I), MAXBIN)
            N_OBSERVED(J,I) = N_OBSERVED(J,I) + 1
          ENDDO
C--
C  And see if it is less than or equal to the total rate observed.
C--
          IF (USE_TOTAL) THEN
            TOTAL_RATE_OBSERVED = 0.0
            TOTAL_RATE_THROWN   = 0.0
            DO I=1, NCHAN
              TOTAL_RATE_OBSERVED =  TOTAL_RATE_OBSERVED +
     >                                          FLOAT(NEVNTS(I))
              TOTAL_RATE_THROWN   =  TOTAL_RATE_THROWN   +
     >                                          FLOAT(EVENTS_THROWN(I))
            ENDDO
            if(round.eq.1) then
            IF (TOTAL_RATE_THROWN .LE. TOTAL_RATE_OBSERVED) THEN
              N_PASS_TOTAL_RATE = N_PASS_TOTAL_RATE + 1
            ENDIF
            elseif(round.eq.2) then
            IF (TOTAL_RATE_THROWN .GE. TOTAL_RATE_OBSERVED) THEN
              N_PASS_TOTAL_RATE = N_PASS_TOTAL_RATE + 1
            ENDIF
            endif
            IF (LDEBUG .AND. I_EXPERIMENT.LT.5) THEN
              WRITE(LUPRT, '('' Total Observed and Thrown = '',
     >                     2G10.2)')
     >              TOTAL_RATE_OBSERVED, TOTAL_RATE_THROWN
            ENDIF
          ENDIF
C--
C  To check against permutations, we sort the number of thrown
C  events AND the number observed using a crude bubble-sort.  For
C  each sorted channel, we see if we would fall within the interval
C  we are calculating the probability for.
C--
          IF (PERMUTE) THEN
            FAILED = .FALSE.
            I = 1
            DO WHILE(I.LE. NCHAN .AND. .NOT.FAILED)
              DO J=I+1, NCHAN
                IF (EVENTS_THROWN(I).LT.EVENTS_THROWN(J)) THEN
C  Sort events
                  ITEMP = EVENTS_THROWN(I)             
C  thrown in this
                  EVENTS_THROWN(I) = EVENTS_THROWN(J)  
C  experiment
                  EVENTS_THROWN(J) = ITEMP             
                ENDIF
C  Sort # events
                IF (NEVNTS(I).LT.NEVNTS(J)) THEN       
C  used to set CL.
                  ITEMP = NEVNTS(I)                    
                  NEVNTS(I) = NEVNTS(J)
                  NEVNTS(J) = ITEMP
                ENDIF
              ENDDO
              if(round.eq.1) then
C  Does it pass
              IF (EVENTS_THROWN(I).GT.NEVNTS(I)) THEN  
C  this cut.
                FAILED = .TRUE.                        
              ENDIF
              I = I+1
              elseif(round.eq.2) then
C  Does it pass
              IF (EVENTS_THROWN(I).LT.NEVNTS(I)) THEN  
C  this cut.
                FAILED = .TRUE.                        
              ENDIF
              I = I+1
              endif
            ENDDO
C
C  Add up exper.
            IF (.NOT.FAILED) THEN                      
C  that pass.
              N_PASS_PERMUTE = N_PASS_PERMUTE + 1      
            ENDIF
C
            IF (LDEBUG .AND. I_EXPERIMENT.LT.5) THEN
              WRITE(LUPRT, '('' FAILED permute = '',L2)') FAILED
            ENDIF
C--
C  Look for events that pass the TOTAL_RATE selection but not
C  when individual rates are being calculated.
C--
            IF (LSUSPICIOUS .AND.
     >          TOTAL_RATE_THROWN .LE. TOTAL_RATE_OBSERVED .AND.
     >              FAILED) THEN
              WRITE(LUPRT, '('' ***Smaller total'',
     >         '' rate but larger individual rates '')')
              WRITE(LUPRT, '(1X,F5.3,4(4('' /'',F7.4,1X,F7.2,I4),/))')
     >           OVERALL_UNCERTAINTY, (EPSILON(J),
     >            MEAN_EVENTS_FOR_CHANNEL(J),
     >            EVENTS_THROWN(J), J=1, NCHAN)
              WRITE(LUPRT, '('' Total Observed and Thrown = '',
     >                     2G10.2)')
     >              TOTAL_RATE_OBSERVED, TOTAL_RATE_THROWN
            ENDIF
C
          ENDIF
C--
C  Dump a little info...
C--
          IF (LDEBUG .AND. I_EXPERIMENT.LT.5) THEN
            WRITE(LUPRT, '(1X,F5.3,4(4('' /'',F7.4,1X,F7.2,I4),/))')
     >         OVERALL_UNCERTAINTY, (EPSILON(J),
     >          MEAN_EVENTS_FOR_CHANNEL(J),
     >          EVENTS_THROWN(J), J=1, NCHAN)
          ENDIF
C
        ENDDO
C--
C  Now we can calculate the individual probabilities for each channel.
C--
        PROB_WITHOUT_PERMUTE = 1.0
        DO I=1, NCHAN
          PROB_FOR_CHANNEL(I) = FLOAT(FREQUENCY(I))/FLOAT(N_EXPERIMENTS)
          PROB_WITHOUT_PERMUTE =
     >             PROB_WITHOUT_PERMUTE * PROB_FOR_CHANNEL(I)
        ENDDO
C--
C  If we are allowing permutations, total probability comes from
C  the probability allowing for permutations.
C--
        IF (USE_TOTAL) THEN
          PROB_USING_TOTAL_RATE = FLOAT(N_PASS_TOTAL_RATE)/
     >                      FLOAT(N_EXPERIMENTS)
          CALCCL = PROB_USING_TOTAL_RATE
        ELSEIF (PERMUTE) THEN
          PROB_WITH_PERMUTE = FLOAT(N_PASS_PERMUTE)/
     >                        FLOAT(N_EXPERIMENTS)
          CALCCL = PROB_WITH_PERMUTE
        ELSE
           CALCCL = PROB_WITHOUT_PERMUTE
        ENDIF
C
        IF (LDUMP) THEN
          WRITE(LUPRT, '('' For '',F8.1,'' events produced, '',
     >                       ''probs are ''/,6X,10(1X,G8.3))')
     >          EVENTS_PRODUCED, (PROB_FOR_CHANNEL(I), I=1, NCHAN)
          WRITE(LUPRT, '(31X,''freqs are '',/,6X,10I9)')
     >          (FREQUENCY(I), I=1, NCHAN)
          IF (PERMUTE) THEN
            PROB_WITH_PERMUTE = FLOAT(N_PASS_PERMUTE)/
     >                          FLOAT(N_EXPERIMENTS)
            WRITE(LUPRT, '(/,2X,I6,'' Events pass allowing for '',
     >          ''permutations; Prob = '',G9.4)')
     >                  N_PASS_PERMUTE, PROB_WITH_PERMUTE
          ENDIF
          WRITE(LUPRT, '(/,8X,'' Probability without permu'',
     >          ''tations is =         '',G9.4)')
     >                  PROB_WITHOUT_PERMUTE
          WRITE(LUPRT, '(/,2X,I6,'' Events pass total rate te'',
     >          ''st; Prob =           '',G9.4)')
     >                  N_PASS_TOTAL_RATE, PROB_USING_TOTAL_RATE
C--
C  Dump the individual frequency distributions.
C--
          IF (.NOT.LDEBUG) THEN
            WRITE(LUPRT,'(/,'' Histogram of event distributions'',
     >                    /,'' --------------------------------'')')
            WRITE(LUPRT,'(1X,'' Events  '',10(''Chan '',I2,3X))')
     >                          (I,I=1,NCHAN)
            DO J=1, MAXBIN
              WRITE(LUPRT,'(1X,I6,2X,10(F8.4,2X))')
     >        J, (FLOAT(N_OBSERVED(J,I))/FLOAT(N_EXPERIMENTS),
     >                                          I=1,NCHAN)
            ENDDO
          ENDIF
        ENDIF
C
C
        BEVENTS_PRODUCED=EVENTS_PRODUCED
C       do i=1,nchan
C         bevents_produced=bevents_produced-bkgnd(i)/acceptance(i)
C       end do
C
        RETURN
        END
