SUBROUTINE PYUPEV(ISUB,SIGEV)

 
C...  All real arithmetic in double precision.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)

C...  Cross-section information.
      COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)

C...  Information on user-selected processes.
      COMMON/PYUPPR/NUP,KUP(20,7),NFUP,IFUP(10,2),PUP(20,5),Q2UP(0:10)

C...  The user's own transfer information
      common/mycomm/ECM,PTMIN,INIT,hlmass

C...  Local arrays and parameters.
      DIMENSION XP1(-25:25),XP2(-25:25)
      DATA PI/3.141592653589793D0/
      DATA CONV/0.3894D0/
      DATA SI2TW/0.232D0/
      DATA WMASS/80.419D0/
      DIMENSION QMASS(6)
      DATA QMASS/0.002,0.006,0.0475,1.25,4.2,175.0/

C-----------------------------------------------------------------

C...Default return values.
      ISUB=200
      SIGEV=0D0

C...Default flag for event rejection.
      IREJ=0

      PTMIN=20D0
      ECM=14000D0

C...
      X1MIN=2D0*hlmass/ECM
      X2MIN=X1MIN

      UHTMIN=hlmass**2
      UHTMAX=ECM**2      

      IF(INIT.EQ.0) THEN

         X1=X1MIN**PYR(0)
         X2=X2MIN**PYR(0)
         UHAT=-UHTMIN**PYR(0)

      ELSE

         X1=2.2D0*hlmass/ECM
         X2=X1
         UHAT=-hlmass**2

      ENDIF

      SHAT=X1*X2*ECM**2
      THAT=-SHAT-UHAT+2*hlmass**2

      PT2=THAT*UHAT/SHAT

c     phase space volum 
      PHSPV=(-LOG(X1MIN))*(-LOG(X2MIN))*LOG(UHTMAX/UHTMIN)

      Q2=PT2
      AEM=PYALEM(Q2)
      ALPS=PYALPS(Q2) 

c...  start of the cross section formula calculation
c...
c...  4 the cross section and the amplitude I use the following articles:
c...  Phys. Lett. B 156 (1985) 429
c...  Phys. Rev. D 55 (1997) 68
c...  

      IF(SHAT.LT.122500.) THEN 

c...  ! shat < 4 mtop^2 = 122500. 
c...  XPAR(2) = 175.0 is the top quark mass
c...  this is just for further studies
c...  shat < 4 mtop^2 will never take place because X1MIN=2*hlmass/ECM
c...  hlmass > 200 GeV/c^2 

         rlamu=QMASS(1)**2 /SHAT
         rlamd=QMASS(2)**2 /SHAT
         rlams=QMASS(3)**2 /SHAT
         rlamc=QMASS(4)**2 /SHAT
         rlamb=QMASS(5)**2 /SHAT
         rlamt=QMASS(6)**2 /SHAT
         
         amp2=16*(
     |         rlamd*asin(1/(2*sqrt(rlamd)))**2
     |        +rlamc*asin(1/(2*sqrt(rlamc)))**2
     |        +rlamt*asin(1/(2*sqrt(rlamt)))**2
     |        -rlamu*asin(1/(2*sqrt(rlamu)))**2
     |        -rlams*asin(1/(2*sqrt(rlams)))**2
     |        -rlamb*asin(1/(2*sqrt(rlamb)))**2
     |        )**2

      ELSE
c...  ! shat > 4 mtop
      rlamu=QMASS(1)**2 /SHAT
      rlamd=QMASS(2)**2 /SHAT
      rlams=QMASS(3)**2 /SHAT
      rlamc=QMASS(4)**2 /SHAT
      rlamb=QMASS(5)**2 /SHAT
      rlamt=QMASS(6)**2 /SHAT

      rpu=1.0+sqrt(1.0-4.0*rlamu)
      rpd=1.0+sqrt(1.0-4.0*rlamd)
      rps=1.0+sqrt(1.0-4.0*rlams)
      rpc=1.0+sqrt(1.0-4.0*rlamc)
      rpb=1.0+sqrt(1.0-4.0*rlamb)
      rpt=1.0+sqrt(1.0-4.0*rlamt)

      rmu=1.0-sqrt(1.0-4.0*rlamu)
      rmd=1.0-sqrt(1.0-4.0*rlamd)
      rms=1.0-sqrt(1.0-4.0*rlams)
      rmc=1.0-sqrt(1.0-4.0*rlamc)
      rmb=1.0-sqrt(1.0-4.0*rlamb)
      rmt=1.0-sqrt(1.0-4.0*rlamt)

      a1=PI**2*(rlamt+rlamc+rlamd-rlamb-rlams-rlamu)
     |     +rlamu*(log(rpu/rmu))**2
     |     +rlams*(log(rps/rms))**2
     |     +rlamb*(log(rpb/rmb))**2
     |     -rlamd*(log(rpd/rmd))**2
     |     -rlamc*(log(rpc/rmc))**2
     |     -rlamt*(log(rpt/rmt))**2

      a2=2.0*PI*(
     |     +rlamu*log(rpu/rmu)
     |     +rlams*log(rps/rms)   
     |     +rlamb*log(rpb/rmb)
     |     -rlamd*log(rpd/rmd)   
     |     -rlamc*log(rpc/rmc)
     |     -rlamt*log(rpt/rmt)
     |     )

         amp2=a1**2+a2**2

      ENDIF

C...  Parton distributions (multiplied by x) of the protons
      CALL PYPDFU(2212,X1,Q2,XP1)
      CALL PYPDFU(2212,X2,Q2,XP2)

C...  CONV=PARU(5)=0.3894 provides the conversion from GeV^-2 to mb
      IF(
     |     IREJ.EQ.0
     |     ) THEN
         const=(AEM*ALPS/SI2TW)**2/(512.0*PI)
         SIGEV=
     |        (hlmass/WMASS**2)**2*
     |        SQRT(1-4.0*hlmass**2/SHAT)*
     |        amp2*XP1(21)*XP2(21)*PHSPV*CONV
         SIGEV=SIGEV*const
      ENDIF

c...  I'll return SIGMAX
      IF(INIT.EQ.1)  RETURN 

C...Define number of partons involved and reset to zero
      NUP=4
      DO 130 I=1,NUP
      KUP(I,1)=1
      DO 120 J=2,7
  120 KUP(I,J)=0
      DO 130 J=1,5
  130 PUP(I,J)=0D0
      
      KUP(1,2)=21
      KUP(2,2)=21
      KUP(3,2)=17
      KUP(4,2)=-17

C...  masses for the outgoing particles
      PUP(3,5)=hlmass
      PUP(4,5)=PUP(3,5)

C...  4-momenta of the incoming partons
      PUP(1,4)=0.5D0*X1*ECM
      PUP(1,3)=PUP(1,4)
      PUP(2,4)=0.5D0*X2*ECM
      PUP(2,3)=-PUP(2,4)

C...Energies and absolute momentum of the outgoing partons in 
C...the subsystem frame.
      RTSHAT=SQRT(X1*X2)*ECM
      PUP(3,4)=0.5D0*( RTSHAT**2 + PUP(3,5)**2 - PUP(4,5)**2 )/RTSHAT
      PUP(4,4)=RTSHAT-PUP(3,4)
      PABS=0.5D0*SQRT(
     |     MAX(
     |     0D0,
     |     ( RTSHAT**2 - PUP(3,5)**2 - PUP(4,5)**2 )**2
     |     -4D0*PUP(3,5)**2*PUP(4,5)**2
     |     )
     |     )/RTSHAT

C...Subsystem scattering angle defined neglecting quark mass.
      COSTHE=(THAT-UHAT)/SHAT
      SINTHE=SQRT(MAX(0D0,1D0-COSTHE**2))

C...Azimuthal angle at random.
      PHI=2D0*PI*PYR(0)

C...Momenta of outgoing partons in the subsystem frame.
      PUP(3,1)=PABS*SINTHE*COS(PHI)
      PUP(3,2)=PABS*SINTHE*SIN(PHI)
      PUP(3,3)=PABS*COSTHE
      PUP(4,1)=-PUP(3,1)
      PUP(4,2)=-PUP(3,2)
      PUP(4,3)=-PUP(3,3)

C...Longitudinal boost of outgoing partons to cm frame.
      BETA=(X1-X2)/(X1+X2)
      GAMMA=0.5D0*(X1+X2)/SQRT(X1*X2)
      PUP(3,3)=GAMMA*(PUP(3,3)+BETA*PUP(3,4))
      PUP(3,4)=GAMMA*(PUP(3,4)+BETA*PUP(3,3))
      PUP(4,3)=GAMMA*(PUP(4,3)+BETA*PUP(4,4))
      PUP(4,4)=GAMMA*(PUP(4,4)+BETA*PUP(4,3))

C...  Pick Q2 scales for initial state shower evolution.
      Q2UP(0)=Q2

      END