C...Simple example how to run gamma-gamma collisions in e+e-
C...with variable energies according to some input distribution.
C...Could be generalized to allow not quite back-to-back photons,
C...but remember the physics of the gamma-gamma collisions is for
C...real incoming photons.
C...Note: program has been converted automatically from Pythia 5.7.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      INTEGER PYK,PYCHGE,PYCOMP
      PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)

C...EXTERNAL statement links PYDATA on most machines.
      EXTERNAL PYDATA
 
C...Common blocks.
      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
      COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
      COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
 
C...Number of events, CM energy of machine.
      NEV=1000
      ECM=180D0
 
C...Range of invariant masses to consider for gamma-gamma.
      EGGMIN=10D0
      EGGMAX=100D0
      PARP(2)=0.999D0*EGGMIN
 
C...Generate combinations of resolved, anomalous and direct components.
      MSEL=2
      MSTP(14)=10
 
C...Allow for variable energies in gamma-gamma collision.
      MSTP(171)=1
 
C...Initialize.
      DO 100 I=1,2
      DO 100 J=1,5
  100 P(I,J)=0D0
      P(1,3)=0.5D0*EGGMAX
      P(2,3)=-0.5D0*EGGMAX
      CALL PYINIT('user','gamma','gamma',ECM)
      NGEN=0
 
C...Generate events.
      DO 200 IEV=1,NEV
  110 CALL GAGAEN(ECM,EGGMIN,EGGMAX,X1,X2)
      NGEN=NGEN+1
      DO 120 I=1,2
      DO 120 J=1,5
  120 P(I,J)=0D0
      P(1,3)=0.5D0*ECM*X1
      P(2,3)=-0.5D0*ECM*X2
      CALL PYEVNT
C...Note that rejection here is IMPORTANT ingredient.
      IF(MSTI(61).EQ.1) GOTO 110
 
C...Here one can add scattered electron/positron by hand.
      DO 130 I=N+1,N+2
      DO 130 J=1,5
      K(I,J)=0
      P(I,J)=0D0
  130 V(I,J)=0D0
      K(N+1,1)=1
      K(N+2,1)=1
      K(N+1,2)=11
      K(N+2,2)=-11
      P(N+1,3)=0.5D0*ECM-P(1,3)
      P(N+2,3)=-0.5D0*ECM-P(2,3)
      P(N+1,4)=ABS(P(N+1,3))
      P(N+2,4)=ABS(P(N+2,3))
      N=N+2
 
C...List first few events. Insert more analysis later.
      IF(IEV.LE.2) CALL PYLIST(1)
 
C...End of event loop. Print statistics.
  200 CONTINUE
      CALL PYSTAT(1)
      WRITE(6,*) 'Number of calls to GAGAEN:',NGEN
 
      END
 
C********************************************************************
 
      SUBROUTINE GAGAEN(ECM,EGGMIN,EGGMAX,X1,X2)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      INTEGER PYK,PYCHGE,PYCOMP
      PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KEXCIT=4000000)
 
C...Range of possible photon energy fractions.
      XMIN=(EGGMIN/ECM)**2
 
C...Generate x1 and x2 independently.
  100 X1=XMIN**PYR(0)
      IF(1D0+(1D0-X1)**2.LT.2D0*PYR(0)) GOTO 100
  110 X2=XMIN**PYR(0)
      IF(1D0+(1D0-X2)**2.LT.2D0*PYR(0)) GOTO 110
 
C...Check that product is in allowed range.
      IF(X1*X2*ECM**2.LT.EGGMIN**2) GOTO 100
      IF(X1*X2*ECM**2.GT.EGGMAX**2) GOTO 100
 
      RETURN
      END