subroutine MOMENT_TEST (id1, id2, rmax, chi2, pval, ierr)

c  Computes 1st rmax alrebraic moments of two histograms and performs chi2
c  test of the hypothesis that they are from the same parent distribution.  
c  Must be linked with HBOOK and CERNLIB.

c  Glen Cowan
c  Royal Holloway, University of London
c  7 April, 1999 

      implicit NONE

c  constants

      integer        max_bins
      parameter      (max_bins = 2000)
      integer        max_moment
      parameter      (max_moment = 100)

c  arguments

      integer        id1, id2             ! input histogram ids
      integer        ierr                 ! output: 0 if OK
      integer        rmax                 ! input maximum moment

      real           chi2 (*)             ! output array of chi2 values
      real           pval (*)             ! output array of P-values

c  local variables

      character*80   title

      integer        i
      integer        j
      integer        loc
      integer        nwt
      integer        nx, nx1, nx2
      integer        ny1, ny2

      real           dx
      real           hist1 (max_bins), hist2 (max_bins)
      real           HSTATI
      real           HSUM                ! HBOOK routine
      real           m1 (2*max_moment)
      real           m2 (2*max_moment)
      real           m_ave
      real           ntot1, ntot2
      real           PROB                ! CERNLIB G100, 1-cdf of chi2
      real           scale, sig1, sig2
      real           var_m1 (max_moment)
      real           var_m2 (max_moment)
      real           x
      real           xmi, xmi1, xmi2
      real           xma, xma1, xma2
      real           ymi1, ymi2
      real           yma1, yma2

c  begin

      ierr = 0
      call HGIVE (id1, title, nx1, xmi1, xma1, 
     &            ny1, ymi1, yma1, nwt, loc)
      call HGIVE (id2, title, nx2, xmi2, xma2, 
     &            ny2, ymi2, yma2, nwt, loc) 
      ntot1 = HSUM(id1)
      ntot2 = HSUM(id2)

c  get rms to use as scale factor

      sig1 = HSTATI (id1, 2, ' ', 0)
      sig2 = HSTATI (id2, 2, ' ', 0)
      scale = 0.5*(sig1 + sig2)

c  check binning of histograms and unpack if OK

      if ( nx1 .gt. max_bins .or. nx2 .gt. max_bins ) then
        write (*, *) 'MOMENT_TEST>> too many bins, histograms = ', 
     &                 id1, id2
        write (*, *) title(1:80)
        write (*, *) 'bins = ', nx1, nx2
        ierr = 1
      elseif ( nx1 .ne. nx2 ) then
        write (*, *) 'MOMENT_TEST>> binning not equal:'
        write (*, *) title(1:80)
        write (*, *) 'id1, nx1 = ', id1, nx1
        write (*, *) 'id2 nx2 = ', id2, nx2
        ierr = 2
      elseif ( xmi1 .ne. xmi2 .or. xma1 .ne. xma2 ) then
        write (*, *) 'MOMENT_TEST>> limits not equal:'
        write (*, *) title(1:80)
        write (*, *) 'id1, xmi1, xma1 = ', id1, xmi1, xma1
        write (*, *) 'id2, xmi2, xma2 = ', id2, xmi2, xma2
        ierr = 3
      elseif ( ntot1 .eq. 0 .or. ntot2 .eq. 0 ) then
        write (*, *) 'MOMENT_TEST>> histogram empty:'
        write (*, *) title(1:80)
        write (*, *) 'id1, ntot1 = ', id1, ntot1
        write (*, *) 'id2, ntot2 = ', id2, ntot2
        ierr = 4
      elseif ( scale .le. 0. ) then
        write (*, *) 'MOMENT_TEST>> scale factor (ave rms) zero:'
        write (*, *) title(1:80)
        write (*, *) 'id1, sig1 = ', id1, sig1
        write (*, *) 'id2, sig2 = ', id2, sig2
        ierr = 5
      elseif ( rmax .gt. max_moment ) then
        write (*, *) 'MOMENT_TEST>> too many moments requested: ', rmax
        write (*, *) 'requested number set to ', max_moment
        rmax = max_moment
      endif

c  divide x limits by rms so as to avoid overflow problems

      if ( ierr .eq. 0 ) then
        xmi = xmi1 / scale 
        xma = xma1 / scale
        nx = nx1
      else
        do i = 1, rmax
          chi2(i) = 99999.
          pval(i) = 0.
        end do
        return        ! -------------------------------------> exit
      endif

      call HUNPAK (id1, hist1, ' ', 0)
      call HUNPAK (id2, hist2, ' ', 0)

c  compute sample moments

      dx = (xma - xmi) / FLOAT(nx)
      do i = 1, 2*rmax
        m1(i) = 0.
        m2(i) = 0.
        do j = 1, nx
          x = xmi + dx * (FLOAT(j) - 0.5)            ! middle of bin
          m1(i) = m1(i) + hist1(j) * x**i
          m2(i) = m2(i) + hist2(j) * x**i
        end do
      end do

      do i = 1, 2*rmax
        m1(i) = m1(i) / ntot1
        m2(i) = m2(i) / ntot2 
      end do

      do i = 1, rmax
        var_m1(i) = (m1(2*i) - m1(i)**2) / ntot1
        var_m2(i) = (m2(2*i) - m2(i)**2) / ntot2
      end do

c      do i = 1, rmax
c        write (*, *) 'moment ', i
c        write (*, *) 'sample 1: ', m1(i), ' +-', SQRT(var_m1(i)) 
c        write (*, *) 'sample 2: ', m2(i), ' +-', SQRT(var_m2(i))
c      end do

c  compute least squares average, chi2 and P-values

      do i = 1, rmax
        if ( var_m1(i) .gt. 0. .and. var_m2(i) .gt. 0. ) then
          m_ave = ( m1(i)/var_m1(i) +  m2(i)/var_m2(i) ) /
     &           ( 1./var_m1(i) + 1./var_m2(i) )
          chi2(i) = (m1(i) - m_ave)**2 / var_m1(i)  +
     &              (m2(i) - m_ave)**2 / var_m2(i)
          pval(i) = PROB(chi2(i), 1)
        else
          write (*, *) 'MOMENT_TEST>> unable to compute moment', i 
          write (*, *) 'for histograms ', id1, id2
          write (*, *) 'var_m1, var_m2 = ', var_m1, var_m2
          chi2(i) = 99999.
          pval(i) = 0.
        endif
      end do

      return
      END