C...Sample program illustrating how to interface CIRCE with PYTHIA.

C...All real arithmetic in double precision.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
C...Three Pythia functions return integers, so need declaring.
      INTEGER PYK,PYCHGE,PYCOMP
C...Parameter statement to help give large particle numbers
C...(left- and righthanded SUSY, excited fermions).
      PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)

C...EXTERNAL statement links PYDATA on most machines.
      EXTERNAL PYDATA

C...Commonblocks.
C...The event record.
      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
C...Parameters.
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
C...Particle properties + some flavour parameters.
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
C...Decay information.
      COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
C...Selection of hard scattering subprocesses.
      COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
C...Parameters. 
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)

C...CIRCE header declarations from file circe.h.
      integer ELECTR, POSITR, PHOTON
      parameter (ELECTR =  11)
      parameter (POSITR = -11)
      parameter (PHOTON =  22)
      integer SBAND, TESLA, XBAND
      parameter (SBAND  =  1, TESLA  =  2, XBAND  =  3)
      integer SBNDEE, TESLEE, XBNDEE
      parameter (SBNDEE =  4, TESLEE =  5, XBNDEE =  6)
      integer NACC
      parameter (NACC = 6)

C...A declaration to give CIRCE access to random numbers.
      EXTERNAL RNG

C...Number of events, CM energy. 
      NEV=1000
      ECM=500D0
 
C...Example of processes: single Z0 or W pairs.  
      MSEL=0
      MSUB(1)=1
C      MSUB(25)=1

C...Option to switch off initial&final-state showers and fragmentation; 
C...is not needed for cross sections and is hence convenient for tryout.
C      MSTP(61)=0
C      MSTP(71)=0
C      MSTP(111)=0

C...Book histograms.
C...(Replace with your own favourite package and wanted
C...distributions.) 
      CALL PYBOOK(1,'mass fraction, no radiation',100,0.005D0,1.005D0)
      CALL PYBOOK(2,'mass fraction, brems only  ',100,0.005D0,1.005D0)
      CALL PYBOOK(3,'mass fraction, beam only   ',100,0.005D0,1.005D0)
      CALL PYBOOK(4,'mass fraction, beam+brems  ',100,0.005D0,1.005D0)

C...For demonstration: loop over four cases of beamstrahlung and 
C...bremsstrahlung off and on.
C...(In actual applications one might settle for IBEAM=IBREM=1.)          
      DO 240 IBEAM=0,1
      DO 220 IBREM=0,1

C...In case you get many warnings of maximum violations you should
C...increase the assumed maximum by hand here. (A factor of two
C...seems to work fine in current cases.)
      MSTP(121)=1
      PARP(121)=2D0

C...Initialize CIRCE - see CIRCE manual for details and variations.
C...Nonzero x_min value as safety measure; so long as it is reasonably
C...small it should not matter, but if made large then event rates
C...should be multiplied by a reduction factor to allow for the
C...not simulated part of the beams.
      IF(IBEAM.EQ.1) THEN
        XMIN=0.0001D0
        IACC=TESLA
        IVER=5
        IREV=0
        ICHAT=1
        CALL CIRCES(XMIN,XMIN,ECM,IACC,IVER,IREV,ICHAT)

C...Allow variable energies in PYTHIA when CIRCE is used.
        MSTP(171)=1
      ELSE
        MSTP(171)=0
      ENDIF

C...Switch on or off bremsstrahlung.
      MSTP(11)=IBREM
 
C...Initialize PYTHIA without or with beamstrahlung.
      IF(IBEAM.EQ.0) THEN 
        CALL PYINIT('CMS','e+','e-',ECM)
      ELSE
        DO 100 J=1,5
          P(1,J)=0D0
          P(2,J)=0D0
  100   CONTINUE  
        P(1,3)=0.5D0*ECM
        P(2,3)=-0.5D0*ECM
        CALL PYINIT('USER','e+','e-',ECM)
      ENDIF

C...Start loop over events.
      DO 200 IEV=1,NEV
   
C...Event generation without beamstrahlung is trivial.
      IF(IBEAM.EQ.0) THEN
        CALL PYEVNT
      
C...Else let CIRCE select x1 and x2 and feed kinematics to PYTHIA.
      ELSE
        DO 120 J=1,5
          P(1,J)=0D0
          P(2,J)=0D0
  120   CONTINUE  
  140   CALL GIRCEE(DX1,DX2,RNG)
        P(1,3)=0.5D0*ECM*DX1
        P(2,3)=-0.5D0*ECM*DX2

C...Generate event; here rejection possible with need for new CIRCE call.
        CALL PYEVNT
        IF(MSTI(61).EQ.1) GOTO 140
      ENDIF 

C...Look at first few events. Note that bremsstrahlung photons are listed
C...while beamstrahlung ones are not (them being external to Pythia).
      IF(IEV.LE.2) CALL PYLIST(1)

C...Insert event analysis here! Currently only trivial example.
      ECMRED=PARI(13)
      CALL PYFILL(2*IBEAM+IBREM+1,ECMRED/ECM,1D0)

C...End loop over events.
  200 CONTINUE

C...Print cross sections.
      CALL PYSTAT(1)

C...End loop over four beamstrahlung/bremsstrahlung cases.
  220 CONTINUE
  240 CONTINUE

C...Histogram output.
      CALL PYHIST

      END

C**************************************************************************

C...Interface to PYTHIA random number generator for CIRCE.
      SUBROUTINE RNG(DR)
C...All real arithmetic in double precision.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      DR=PYR(0)
      RETURN
      END