+PATCH,//MY_INCLUDE_W19/PYTHIA
+DECK,pyupev.
*CMZ :          27/05/98  09.13.22  by  Dugan C ONeil
*-- Author :
C*********************************************************************

c...PYUPEV
c...Dummy routine, to be replaced by user. When called from PYTHIA
c...the subprocess number ISUB will be given, and PYUPEV is supposed
c...to generate an event of this type, to be stored in the PYUPPR
c...commonblock. SIGEV gives the differential cross-section associated
c...with the event, i.e. the acceptance probability of the event is
c...taken to be SIGEV/SIGMAX, where SIGMAX was given in the PYUPIN
c...call.

      SUBROUTINE PYUPEV(ISUB,SIGEV)
c
c...Double precision and integer declarations (must be present for
c PYTHIA 6 common blocks to work.
c
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      INTEGER PYK,PYCHGE,PYCOMP
c
c my stuff (local variables)
c
      double precision sigma0, sigev
c
      double precision rapid1,rapid2
      double precision x1,x2,mw
c
      double precision xpp(-25:25)
c
      integer i,j
      double precision raplim
      double precision roots
      double precision pi
      parameter(pi=3.1415927)
c
      double precision comb
      integer icomb,ncomb
      integer gen1,gen2
      parameter(ncomb=6)
      integer p1(ncomb),p2(ncomb)
      data p1/2,2,4,-1,-3,-3/
      data p2/-1,-3,-3,2,2,4/
c
c
C...Commonblocks.
+SEQ,LUDAT2.
c
c include PYTHIA common block to be filled by routine
c
      COMMON/PYUPPR/NUP,KUP(20,7),NFUP,IFUP(10,2),PUP(20,5),Q2UP(0:10)
      SAVE /PYDAT1/,/PYUPPR/
c
c get values of useful parameters (could read these from PYTHIA or datacard)
c
c roots = centre of mass energy
c mw    = W mass
c gf    = fermi coupling
c
      roots = 14000.d0
      mw=80.3
      gf=1.17e-5
c
      if (isub .eq. 195) then
c
c randomly choose flavour of quark
c we have six combinations to choose from to give W+
c   p1            p2
c 1 u             dbar
c 2 u             sbar
c 3 c             sbar
c 4 dbar          u
c 5 sbar          u
c 6 sbar          c
c
c pyr is PYTHIA random number generator
c
        comb = min((pyr(0)*6. + 1.),6.)
        icomb = int(comb)
c
c translate particle code in p1 and p2 into generation number
c for later use in VCKM matrix. gen1=up-type number, gen2=down-type number
c
        if (mod(p1(icomb),2) .eq. 0)then
          gen1 = abs(p1(icomb))/2
        else
          gen2 = (abs(p1(icomb)) + 1)/2
        endif
        if (mod(p2(icomb),2) .eq. 0)then
          gen1 = abs(p2(icomb))/2
        else
          gen2 = (abs(p2(icomb)) + 1)/2
        endif
ccc
ccc        print *,'icomb,gen1,gen2',icomb,gen1,gen2
ccc
c
c choose numbers flat in rapidity from -raplim to raplim
c
        raplim = 5.16
        rapid1 = pyr(0)*2*raplim - raplim
        rapid2 = pyr(0)*2*raplim - raplim
c
c convert rapidity to momentum fraction
c
        x1 = (mw/roots)*exp(rapid1)
c
c the momentum fraction of one is related to the other via
c x1*x2=mw**2/roots**2 in order to produce a W
c
        x2 = (mw**2/(roots**2*x1))
c
c use internal pythia routine to calculate structure functions
c at these two momentum fractions
c 2212 = proton
c
        call pypdfu(2212,x1,mw**2,xpp)
        f1 = xpp(p1(icomb))/x1
c
        call pypdfu(2212,x2,mw**2,xpp)
        f2 = xpp(p2(icomb))/x2
ccc
ccc        print *,'p1,p2,icomb ',p1(icomb),p2(icomb),icomb
ccc
c
c evaluate constant factors in differential cross section calculation
c .3894 is conversion from GeV^(-2) to mb
c
        sigma0=0.3894*(2./3.)*pi*(mw**2/roots**2)*gf/(sqrt(2.))
c
c evaluate for chosen quark combination
c
        sigev = sigma0*(vckm(gen1,gen2)*f1*f2)
ccc
ccc        print *,'vckm,f1,f2 ',vckm(gen1,gen2),f1,f2
ccc        print *,'sigev  ',sigev
ccc
c
c multiply by integral of chosen variable (in this case rapidity from
c -raplim to raplim) and sum over number of possible combinations of incoming
c quarks
c
        sigev = sigev*2*raplim*ncomb
c
ccc        print *,' sigma0,sigev --> ',sigma0,sigev
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                   fill event common block                         c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c number of particles
c
        nup = 3
c
c zero event common except for flavour and status codes
c
        do i = 1, nup
          do j = 3, 7
            kup(i,j)=0
          enddo
          do j = 1, 5
            pup(i,j)=0.d0
          enddo
        enddo
c
c allow all particles to be decayed by pythia
c
c particle 1&2 are incoming partons, particle 3 is the W
c
        kup(1,1) = 1
        kup(2,1) = 1
        kup(3,1) = 1
c
c define flavour codes of particles (choice of quarks is made above)
c
        kup(1,2) = p1(icomb)
        kup(2,2) = p2(icomb)
        kup(3,2) = 24
c
c define position of mother (not needed)
c
        kup(1,3) = 0
        kup(2,3) = 0
        kup(3,3) = 0
c
c which final state particle carries colour....none
c
        kup(1,4) = 0
        kup(2,4) = 0
        kup(3,4) = 0
c
c  which final state particle carries anti-colour....none
c
        kup(1,5) = 0
        kup(2,5) = 0
        kup(3,5) = 0
c
c where does initial state colour go to from each particle
c make quarks carry colour, anti-quarks carry anti-colour
c
c eg. kup(1,6) = 2 means that colour from particle 1 goes to particle 2
c
        if (kup(1,2) .gt. 0) then
          kup(1,6) = 2
          kup(2,6) = 0
          kup(3,6) = 0
c...anticolour
          kup(1,7) = 0
          kup(2,7) = 1
          kup(3,7) = 0
        else
          kup(1,6) = 0
          kup(2,6) = 1
          kup(3,6) = 0
c...anticolour
          kup(1,7) = 2
          kup(2,7) = 0
          kup(3,7) = 0
       endif
c
c define momenta in COM frame
c
c...x
        pup(1,1) = 0.
        pup(2,1) = 0.
        pup(3,1) = 0.
c...y
        pup(1,2) = 0.
        pup(2,2) = 0.
        pup(3,2) = 0.
c...z
        pup(1,3) = x1*roots/2.
        pup(2,3) = -x2*roots/2.
        pup(3,3) = pup(1,3) + pup(2,3)
c...energy
        pup(1,4) = pup(1,3)
        pup(2,4) = -pup(2,3)
        pup(3,4) = pup(1,4) + pup(2,4)
c...mass
        pup(1,5) = 0.
        pup(2,5) = 0.
        pup(3,5) = mw
c
      endif
c
 999  continue
c
      RETURN
      END