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...Common blocks.
      COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
      COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
      COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
      COMMON/PYSUBS/MSEL,MSUB(200),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=180.

C...Range of invariant masses to consider for gamma-gamma.
      EGGMIN=10.
      EGGMAX=100.
      PARP(2)=0.999*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)=0.
      P(1,3)=0.5*EGGMAX
      P(2,3)=-0.5*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)=0.
      P(1,3)=0.5*ECM*X1
      P(2,3)=-0.5*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)=0.
  130 V(I,J)=0.
      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.5*ECM-P(1,3)
      P(N+2,3)=-0.5*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 LULIST(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)

C...Range of possible photon energy fractions.
      XMIN=(EGGMIN/ECM)**2

C...Generate x1 and x2 independently.
  100 X1=XMIN**RLU(0)
      IF(1.+(1.-X1)**2.LT.2.*RLU(0)) GOTO 100    
  110 X2=XMIN**RLU(0)
      IF(1.+(1.-X2)**2.LT.2.*RLU(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