C...Simple example how to run DIS e-gamma collisions in e+e-
C...with variable energies according to some input distribution.
C...Is only valid for sensibly large Q2, and not too small x.
C...An improved version is available from Leif Lonnblad.

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)
      COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3) 
 
C...Number of events, CM energy of machine.
      NEV=200
      ECM=180.

C...Range of invariant masses to consider for e-gamma system.
      EEGMIN=20.
      PARP(2)=EEGMIN
      
C...Generate combinations of resolved and anomalous components.
      MSTP(14)=1

C...Set pT_min, Q2_min, W2min, x_min and x_max.
      CKIN(3)=2.
      CKIN(35)=10.
      CKIN(39)=400. 
      CKIN(23)=0.01
      CKIN(24)=0.9

C...Allow for variable energies in e-gamma collision.
      MSTP(171)=1 

C...Loop over the various photon sets.
      DO 300 ISET=1,9
      IF(ISET.EQ.1) MSTP(55)=1
      IF(ISET.GE.2) MSTP(55)=ISET+3

C...Give weight of photon spectrum.
      MSTP(173)=1
      CALL EGAEN(ECM,EEGMIN,X1,X2,WTEV,WTMX)
      PARP(174)=WTMX

C...Initialize.
      DO 100 I=1,2
      DO 100 J=1,5
  100 P(I,J)=0.
      P(1,3)=0.5*ECM
      P(2,3)=-0.5*ECM
      CALL PYINIT('user','e-','gamma',ECM)
 
C...Generate events.
      DO 200 IEV=1,NEV
  110 CALL EGAEN(ECM,EEGMIN,X1,X2,WTEV,WTMX)
      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
      PARP(173)=WTEV
      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 J=1,5
      K(N+1,J)=0
      P(N+1,J)=0.
  130 V(N+1,J)=0.
      K(N+1,1)=1
      K(N+1,2)=-11
      P(N+1,3)=-0.5*ECM-P(2,3)
      P(N+1,4)=ABS(P(N+1,3))
      N=N+1

C...List first few events. Insert more analysis later.
c      IF(IEV.LE.2) CALL LULIST(1)

C...End of event loop. Print statistics.
  200 CONTINUE
      CALL PYSTAT(1)
  300 CONTINUE
      
      END

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

      SUBROUTINE EGAEN(ECM,EEGMIN,X1,X2,WTEV,WTMX)

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

C...Set electron fraction=1, generate photon fraction.
  100 X1=1.
      X2=XMIN**RLU(0)

C...Evaluate photon flux weights.
      AEM2PI=1./(137.*2.*3.1416)
      Q2RNG=0.77**2/0.0005**2
      WT2=AEM2PI*LOG(Q2RNG*(1.-X2)/X2**2)*(1.+(1.-X2)**2)      

C...Evaluate weight for x selection procedure.
      WTX=-LOG(XMIN)

C...Give total event weight.
      WTEV=WT2*WTX

C...Give maximum event weight possible
      WTIMX=AEM2PI*LOG(Q2RNG/XMIN**2)*2.
      WTMX=WTIMX*WTX

      RETURN
      END