program tstpyt
C
C      test program for stdhep/pythia interface
C      compare event list
C      write events and cross-section
C
#include "stdlun.inc"

      integer istr,mstat
      integer nevt,i,nprnt,lok
C
      nevt=100
      nprnt=1
      call upyinit(nevt,istr)

C          Fill Stdhep common block 1 (only valid for Pythia)
      call stdflpyxsec(nevt)
C          Write Stdhep begin-run record   
      call stdxwrt(100,istr,lok)
      if(lok.ne.0) write(lnhout,*)
     1        ' Problem writing stdhep begin run record'
c
      do 100 i=1,nevt
        call pyevnt
        call lunhep(1)
        if(i.le.nprnt)then
          call pylist(2)
          call heplst(1)
        endif
        call pytflem(1,i)
        call stdxwrt(1,istr,lok)
C          Fill Stdhep common block 1 (only valid for Pythia)
        call stdflpyxsec(nevt)
        call pytflem(2,i)
 100  continue
C
C          Fill Stdhep common block 1
      call stdflpyxsec(nevt)
C          Write end-of-run record  
      call stdxwrt(200,istr,lok)
      if(lok.ne.0) write(lnhout,*) ' Problem writing end run record'

C...   Print end of run statistics.
      mstat=1
      call pystat(mstat)

c...close event file
      call stdxend(istr)
C...print histograms
      call houtpu(lnhout)
      call histdo
      close(unit=lnhout)
      stop
      end

      subroutine upyinit(nevt,istr)
C
C   Pythia initialization
C
#include "pydat1.inc"
#include "stdlun.inc"

C...   PYTHIA initialization variables.
      character*4 frame
      character*5 beam,target
      integer nevt,istr
      double precision win
C
C initialize HEP logical units
      lnhwrt=23
      lnhrd=0
      lnhdcy=0
      lnhout=22
      mstu(11)=lnhout
      open(unit=lnhout,file='stdtstpx.lpt',status='new')
      call stdxwinit('stdtstpx.io','StdHep/Pythia example',
     1               nevt,istr,lok)
C
      CALL PYTBKEM
C
C...   Initialize PYTHIA
      frame  = 'cms'
      beam   = 'pbar'
      target = 'p'
      win    = 1800.0d0
      call pygive('msel=5')
      call pygive('mstj(21)=0')
      call pyinit(frame,beam,target,win)

      return
      end

      subroutine pytflem(iblk,iev)
C... fill histograms
C   STDXSEC  - cross-section
C
#include "stdhep.inc"
#include "stdcm1.inc"

      integer iblk,i,iev
      real bn

      if(iblk.eq.1)then
        call hfill(1,float(nevhep),0.,1.)
        call hfill(2,float(nhep),0.,1.)
        do i=1,nhep
          call hfill(3,float(idhep(i)),0.,1.)
          call hfill(4,float(isthep(i)),0.,1.)
          call hfill(5,sngl(phep(1,i)),0.,1.)
          call hfill(6,sngl(phep(2,i)),0.,1.)
          call hfill(7,sngl(phep(3,i)),0.,1.)
          call hfill(8,sngl(phep(4,i)),0.,1.)
          call hfill(9,sngl(phep(5,i)),0.,1.)
        enddo
      elseif(iblk.eq.2)then
        bn = float(iev)+0.5
        call hfill(10,bn,0.,stdxsec)
      endif
      return
      end
      subroutine pytbkem
C--- book histograms
C
      implicit none
C---   The HBOOK common
      real hmemor
      common/pawc/hmemor(50000)
C---   Setup the HBOOK memory limit.
      call hlimit(50000)
C---   Book histograms.
      call hbook1(1,'PYTHIA: NEVHEP',100,0.,100.,0.)
      call hbook1(2,'PYTHIA: NHEP',100,0.,200.,0.)
      call hbook1(3,'PYTHIA: IDHEP',100,0.,100.,0.)
      call hbook1(4,'PYTHIA: ISTHEP',20,0.,20.,0.)
      call hbook1(5,'PYTHIA: PHEP(1,)',100,-100.,100.,0.)
      call hbook1(6,'PYTHIA: PHEP(2,)',100,-100.,100.,0.)
      call hbook1(7,'PYTHIA: PHEP(3,)',100,-100.,100.,0.)
      call hbook1(8,'PYTHIA: PHEP(4,)',100,-100.,100.,0.)
      call hbook1(9,'PYTHIA: PHEP(5,)',100,-100.,100.,0.)
      call hbook1(10,'STDHEP: cross-section',100,0.,100.,0.)
      return
      end