C
C batcher.f - a simple 1d MC generator to demonstrate the dangers of doing
C             keep/discard MC generation with a fluctuating maxProb cutoff.
C             This is what is happening in the AmpTools weighted MC generator
C             where events are generated in batches, and maxInten is computed
C             max(I[])*1.5 separately for every batch.
C
      subroutine batcher(mu,sigma)
      real mu,sigma
      include 'batcher.inc'
      logical hexist
      external hexist
      mean=mu
      variance=sigma**2
      id=10
      title='merged sample from all batches'
      if (hexist(id)) then
        call hdelet(id)
      endif
      nbins=100
      xlow=-10
      xhigh=10
      call hbook1(id,title,nbins,xlow,xhigh,0)
      print *, 'generate sample with CALL batch(N)'
      end

      subroutine batch(N)
      integer N
      include 'batcher.inc'
      integer i
      real pmax
      if (N.gt.batch_maxsize) then
        print *, 'Requested batch size',N,' is larger than maximum',
     +           batch_maxsize,'.  Must recompile to increase limit.'
        return
      endif
      call ranlux(rvec,2*N)
      pmax=0
      do i=1,N
        rvec(1,i)=xlow+(xhigh-xlow)*rvec(1,i)
        pvec(i)=pdf(rvec(1,i))
        pmax=max(pmax,pvec(i))
      enddo
      pmax=pmax*1.5
      do i=1,N
        if (rvec(2,i)*pmax.lt.pvec(i)) then
          call hfill(id,rvec(1,i),0.,1.)
        endif
      enddo
      end

      subroutine batches(m,N)
      integer m,N
      integer i
      print *, 'generating',m,' batches of size',N
      do i=1,m
        call batch(N)
      enddo
      end

      real function pdf(x)
      real x
      include 'batcher.inc'
      pdf=exp(-0.5*(x-mean)**2/variance)
      end