program COMPUTE_CHI2_DIST

C  Author:  Glen Cowan
C  Date:    7-OCT-1997

c  Program computes sampling p.d.f. of a chi2 statistic by Monte Carlo.
c  Because of the small numbers of entries in various bins, this is
c  expected to differ from the usual chi2 distribution.
c  The usual chi2 distribution is also computed for comparison.

c  Must be linked with the CERN program library.

      implicit        NONE

c  Needed for HBOOK routines

      integer         hsize
      parameter       (hsize = 100000)

      integer         hmemor (hsize)
      common  /pawc/  hmemor

c  Parameters and local variables

      integer         max_points
      parameter       (max_points = 200)      ! increase if necessary

      character*80    infile
      character*80    outfile
      character*80    line

      integer         i
      integer         icycle
      integer         ierror
      integer         ios
      integer         istat
      integer         j
      integer         lun
      integer         n (max_points)
      integer         n_dof
      integer         num_exper
      integer         num_points

      real            bin_width
      real            chi2
      real            f
      real            GAMMA
      real            hist (max_points)
      real            nu (max_points)
      real            x
      real            x_min (max_points)
      real            x_max (max_points)

c  Initialize HBOOK, open histogram file, book histograms.

      call HLIMIT (hsize)
      outfile = 'chi2.his'
      call HROPEN (20, 'histog', outfile, 'N', 1024, istat)
      if (istat .ne. 0) write (*, *) 'HROPEN error, istat = ', istat
      call HBOOK1 (1, 'test statistic', 100, 0., 100., 0.)
      call HBOOK1 (2, 'chi2', 100, 0., 100., 0.)

c  Read theoretical prediction from file.

      write (*, *) 'enter file name for theoretical prediction'
      read (*, '(a80)') infile
      lun = 30
      open (unit=lun, file = infile,
     1  form = 'formatted', status = 'old')

      i = 0
      ios = 0
      do while ( ios .EQ. 0 )
        read (lun, fmt = '(a80)', iostat = ios) line
        if (ios .eq. 0 .and. .not. line(1:1) .eq. '!') then
          i = i + 1
          read (line, *) x_min(i), x_max(i), nu(i)
        endif
      end do
      close (lun)
      num_points = i

c  Perform Monte Carlo experiments to determine distribution of
c  chi2 statistic.

      write (*, *) 'enter number of MC experiments'
      read (*, *) num_exper

      do i = 1, num_exper

c  Generate data -- each bin is treated as a Poisson variable with mean
c  nu(i), generated with CERN routine RNPSSN (V136).

        do j = 1, num_points
          call RNPSSN (nu(j), n(j), ierror)
        end do

c  Compute chi2 and enter into histogram.

        chi2 = 0.
        do j = 1, num_points
          chi2 = chi2 + (FLOAT(n(j)) - nu(j))**2 / nu(j)
        end do

        call HF1 (1, chi2, 1.)

      end do

c  For comparison, compute the usual chi2 distribution.  This requires
c  a gamma function, taken from CERNLIB routine GAMMA (C305).
c  Multiply by the number of entries and bin width so that the curve
c  has the same area as the MC histogram.

      n_dof = num_points
      bin_width = 1.0
      do i = 1, 100
        x = FLOAT(i)
        f = x**(n_dof/2. - 1.) * EXP (-x/2.) /
     1      ( 2.**(n_dof/2.) * GAMMA(n_dof/2.) )
        hist(i) = f * FLOAT(num_exper) * bin_width
      end do
      call HPAK (2, hist)

c  Store histogram and close

      call HROUT (0, icycle, ' ')
      call HREND ('histog')

      stop
      END