lysspa.F

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lysspa.F

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

      SUBROUTINE lysspa(IPU1,IPU2)

      IMPLICIT NONE

C...NEW X REDEFINITION
C...GENERATES SPACELIKE PARTON SHOWERS
*
* to avoid variable conflictions, a second keep element is necessary
* with the same common block name (see LEPTO2)
*

      COMMON /leptou/ cut(14),lst(40),parl(30),
     &                xlp,ylp,w2lp,q2lp,ulp
      REAL cut,parl,xlp,ylp,w2lp,q2lp,ulp
      INTEGER lst
      SAVE /leptou/

      INTEGER  nlupdm,nplbuf
      parameter(nlupdm=4000,nplbuf=5)
      common/lujets/n,k(nlupdm,5),p(nlupdm,nplbuf),v(nlupdm,5)
      INTEGER  n,k
      REAL  p,v
      SAVE /lujets/

      common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
      INTEGER  mstu,mstj
      REAL  paru,parj
      SAVE /ludat1/

      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
      INTEGER  kchg
      REAL  pmas,parf,vckm
      SAVE /ludat2/

      COMMON /pypara/ ipy(80),pypar(80),pyvar(80)
      REAL pypar,pyvar
      INTEGER ipy
      SAVE /pypara/ 

*
* to avoid variable conflictions, a second keep element is necessary
* with the same common block name (see LYPRO2)
*
      COMMON /lyproc/ isub,kfl(3,2),x(2),sh,th,uh,q2,xsec(0:40)
      REAL    x,sh,th,uh,q2,xsec
      INTEGER isub,kfl
      SAVE /lyproc/ 

      COMMON /lyint1/ xq(2,-6:6),dsig(-6:6,-6:6,5),fsig(10,10,3)
      REAL    xq,dsig,fsig
      SAVE /lyint1/ 
 


      INTEGER ipu1,ipu2,ifls,is,nq,ilep,ifla,ns,ifl,jt,ihfc,ihfx,
     +i,j,jr,iflb,ihft,ipo,i1,i2,itemp,it,ikin,id1,kn1,kd1,id2,
     +ir,jb,ntry,n145, njr
      REAL xs,zs,q2s,tevs,robo,xfs,xfa,xfb,wtap,wtsf,z,xe0,xa,
     +tmax,qmax,q2e,b0,xbmin,xb,q2b,tevb,qmass,xe,wtapq,wtsum,tevxp,
     +q2ref,wtran,wtz,xbnew,rsoft,zu,q2max,alprat,the
      REAL ulmass,rlu,ulangl,xt,searatio
 
      dimension ifls(4),is(2),xs(2),zs(2),q2s(2),tevs(2),robo(5),
     &xfs(2,-6:6),xfa(-6:6),xfb(-6:6),wtap(-6:6),wtsf(-6:6)
      DOUBLE PRECISION dq2(3),dsh,dshz,dshr,dplcm,dpc(3),dpd(4),dms,
     &dmsma,dpt2,dpb(4),dbe1(4),dbe2(4),dbep,dgabep,dpq(4),dpqs(2),
     &dm2,dq2b,drobo(5),dbez,dtemp
C-GI &DQ23,DPH(4),DM2,DQ2B,DQM2
CJR--begin
      LOGICAL seaquark,split
      REAL xft(-6:6)
      REAL xquark,xgluon,xsea,zsplit,zsoft,zmax
      COMMON /seaqte/ xquark,xgluon,xsea,zsplit,zsoft,zmax,split
      INTEGER lastfl,seafl
      COMMON /flavor/ lastfl,seafl
CJR--end
      DATA ifla,nq/0,0/,z,xe0,xa/3*0./,dshz,dmsma,dpt2,dshr/4*0.d0/

C...COMMON CONSTANTS, SET UP INITIAL VALUES
      ilep=0
      IF(ipu1.EQ.0) ilep=1
      IF(ipu2.EQ.0) ilep=2
      q2e=q2
C-GI  IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.3) Q2E=Q2E/PYPAR(26)
      IF(isub.EQ.27) q2e=pmas(23,1)**2
      IF(isub.EQ.28) q2e=pmas(24,1)**2
      tmax=alog(pypar(26)*pypar(27)*q2e/pypar(21)**2)
      IF(ilep.GE.1) THEN
        sh=p(25,5)**2
        IF(n.GE.27) sh=p(27,5)**2
        CALL lscale(-1,qmax)
        q2e=qmax**2
        q2e=max(pypar(21)**2,min(q2e,(0.95/x(3-ilep)-1.)*q2-sh,
     &  q2/2.+sh))
        tmax=alog(q2e/pypar(21)**2)
      ENDIF
CJR--begin
      IF (mod(lst(8),10).EQ.4 .OR. mod(lst(8),10).EQ.5) THEN
         q2e=pypar(22)
         tmax=alog(q2e/pypar(21)**2)
      ELSEIF(pypar(26)*q2e.LT.max(pypar(22),2.*pypar(21)**2).OR.
     &tmax.LT.0.2) THEN
         RETURN
      ENDIF
CJR--end
      IF(ilep.EQ.0) xe0=2.*pypar(23)/pyvar(1)
      b0=(33.-2.*ipy(8))/6.
      ns=n
      mstu(2)=0
CJR--begin
      ntry=0
  100 n=ns
      ntry=ntry+1
      IF (ntry.GT.100) THEN
         lst(21)=17
         RETURN
      ENDIF
CJR--end
CJR  100 N=NS
CJR--begin
      seaquark=.false.
      split=.false.
CJR--end
      IF(ilep.GE.1) THEN
        nq=ipu2-2
        IF(ilep.EQ.2) nq=ipu1+2
        dpqs(1)=dble(p(nq,3))
        dpqs(2)=dble(p(nq,4))
        xbmin=x(3-ilep)*max(0.5,sh/q2)
        CALL lystfu(ipy(43-ilep),xbmin,q2,xfb)
        DO 110 ifl=-6,6
  110   xq(3-ilep,ifl)=xfb(ifl)
      ENDIF
      DO 120 jt=1,2
      ifls(jt)=kfl(2,jt)
      IF(kfl(2,jt).EQ.21) ifls(jt)=0
      ifls(jt+2)=ifls(jt)
      xs(jt)=x(jt)
      zs(jt)=1.
      IF(ilep.EQ.0) q2s(jt)=pypar(26)*q2e
      tevs(jt)=tmax
      DO 120 ifl=-6,6
  120 xfs(jt,ifl)=xq(jt,ifl)
      IF(ilep.GE.1) THEN
        q2s(ilep)=p(nq,5)**2
        dq2(ilep)=q2s(ilep)
        q2s(3-ilep)=q2e
      ENDIF
      dsh=sh
      ihfc=0
      ihfx=0

C...PICK UP LEG WITH HIGHEST VIRTUALITY
  130 CONTINUE
      IF(n.GT.mstu(4)-10) THEN
        WRITE(6,*) ' LYSSPA: no more memory in LUJETS'
        lst(21)=18
        RETURN
      ENDIF
      DO 133 i=n+1,n+8
      DO 133 j=1,5
      k(i,j)=0
  133 p(i,j)=0.
C     CALL GULIST(21,2)
      n=n+2
      jt=1
      IF((n.GT.ns+2.AND.q2s(2).GT.q2s(1).AND.ilep.EQ.0).OR.ilep.EQ.1)
     &jt=2
      jr=3-jt
      iflb=ifls(jt)
      xb=xs(jt)
      IF(ilep.GE.1.AND.n.EQ.ns+2) xb=xs(jt)*max(sh/q2,0.5)
      DO 140 ifl=-6,6
  140 xfb(ifl)=xfs(jt,ifl)
      q2b=q2s(jt)
      tevb=tevs(jt)
      IF(ipy(14).GE.9.AND.n.GT.ns+4) THEN
        q2b=0.5*(1./zs(jt)+1.)*q2s(jt)+0.5*(1./zs(jt)-1.)*(q2s(3-jt)-
     &  sngl(dsh)+sqrt((sngl(dsh)+q2s(1)+q2s(2))**2+8.*q2s(1)*q2s(2)*
     &  zs(jt)/(1.-zs(jt))))
        tevb=alog(pypar(27)*q2b/pypar(21)**2)
      ENDIF
      IF(ilep.EQ.0) THEN
        dshr=2.*dsqrt(dsh)
        dshz=dsh/dble(zs(jt))
      ELSEIF(ilep.GE.1) THEN
        dshz=dsh
        IF(n.GT.ns+4) dshz=(dsh+dq2(jr)-dq2(jt))/zs(jt)-dq2(jr)+
     &  pypar(22)
        dpd(2)=dshz+dq2(jr)+dble(pypar(22))

        mstj(93)=1
        qmass=ulmass(iabs(iflb))
        IF(iabs(iflb).EQ.0) qmass=ulmass(21)
C...CHECK IF QUARK PAIR CREATION ONLY POSSIBILITY
        IF(dq2(jr).LT.4.*qmass**2) THEN
          dm2=qmass**2
          dpc(1)=dq2(jr)*(dble(pypar(22))+dm2)**2
          dpc(2)=dpd(2)*(dpd(2)-2d0*pypar(22))*(pypar(22)+dm2)
          dpc(3)=pypar(22)*(dpd(2)-2d0*pypar(22))**2
          xe0=1d0-(dpc(2)-dsqrt(dpc(2)**2-4d0*dpc(1)*dpc(3)))/
     &    (2d0*dpc(1))
CJR--begin
          zmax=(dpc(2)-dsqrt(dpc(2)**2-4d0*dpc(1)*dpc(3)))/(2d0*dpc(1))
          xe0=xb*(1./zmax-1.)
CJR--end
        ELSE
          xe0=1d0-(dpd(2)-2d0*dble(pypar(22)))*(dpd(2)-dsqrt(dpd(2)**2-
     &    4d0*dq2(jr)*dble(pypar(22))))/(2d0*dq2(jr)*dble(pypar(22)))
CJR--begin
          zmax=(dpd(2)-2d0*dble(pypar(22)))*(dpd(2)-dsqrt(dpd(2)**2-
     &    4d0*dq2(jr)*dble(pypar(22))))/(2d0*dq2(jr)*dble(pypar(22)))
          xe0=xb*(1./zmax-1.)
CJR--end
        ENDIF
CJR--begin
CJR--          radiated parton energy cut
C        XE0=MAX(XE0,2.*PYPAR(23)/SQRT(W2LP))
CJR--end
      ENDIF
CJR  145 XE=MAX(XE0,XB*(1./(1.-PYPAR(24))-1.))
CJR--begin
      n145=0
145   CONTINUE
      n145=n145+1
      IF (n145.GT.100) THEN
CJR         WRITE(*,*) '145'
         goto 100
      ENDIF
      xe=max(xe0,xb*(1./(1.-pypar(24))-1.))
        zmax=xb/(xb+xe)
CJR--end
      IF(xb+xe.GE.0.999) THEN
        q2b=0.
        goto 210
      ENDIF

C...CALCULATE ALTARELLI-PARISI AND STRUCTURE FUNCTION WEIGHTS
      DO 150 ifl=-6,6
      wtap(ifl)=0.
  150 wtsf(ifl)=0.
      IF(iflb.EQ.0) THEN
        wtapq=16.*(1.-sqrt(xb+xe))/(3.*sqrt(xb))
        DO 160 ifl=-ipy(8),ipy(8)
        IF(ifl.EQ.0) wtap(ifl)=6.*alog((1.-xb)/xe)
  160   IF(ifl.NE.0) wtap(ifl)=wtapq
      ELSE
        wtap(0)=0.5*xb*(1./(xb+xe)-1.)
        wtap(iflb)=8.*alog((1.-xb)*(xb+xe)/xe)/3.
      ENDIF
  170 wtsum=0.
      IF(ihfc.EQ.0) THEN
        DO 180 ifl=-ipy(8),ipy(8)
        wtsf(ifl)=xfb(ifl)/max(1e-10,xfb(iflb))
  180   wtsum=wtsum+wtap(ifl)*wtsf(ifl)
        IF(iabs(iflb).GE.4.AND.wtsum.GT.1e3) THEN
          ihfx=1
          DO 185 ifl=-ipy(8),ipy(8)
  185     wtsf(ifl)=wtsf(ifl)*1e3/wtsum
          wtsum=1e3
        ENDIF
      ENDIF

C...CHOOSE NEW T AND FLAVOUR
CJR  190 IF(IPY(14).LE.6.OR.IPY(14).GE.9) THEN
CJR--begin
      njr=0
190   CONTINUE
      seaquark=.false.
      njr=njr+1
      IF (njr.GT.100) THEN
CJR         WRITE(*,*) '190'
         goto 100
      ENDIF
CJR--end
      IF(ipy(14).LE.6.OR.ipy(14).GE.9) THEN
        tevxp=b0/max(0.0001,wtsum)
      ELSE
        tevxp=b0/max(0.0001,5.*wtsum)
      ENDIF
      tevb=tevb*exp(max(-100.,alog(rlu(0))*tevxp))
      q2ref=pypar(21)**2*exp(tevb)/pypar(27)
      q2b=q2ref/pypar(27)
      dq2b=q2b
      IF(ilep.GE.1) THEN
        dshz=dsh
        IF(n.GT.ns+4) dshz=(dsh+dq2(jr)-dq2(jt))/dble(zs(jt))-dq2(jr)+
     &  dq2b
      ENDIF
CJR--begin -- 
CAE--seaquarks up to LST(12), if the quark density != 0
      IF( lst(35).EQ.1) THEN
         IF( q2b.LT.pypar(22) .AND.
     &       (abs(iflb).LE.lst(12).AND.abs(iflb).GE.1)) THEN
            q2ref=min(pypar(22),q2s(jt))
            q2b=min(pypar(22),q2s(jt))
            IF(ilep.GE.1.AND.n.EQ.ns+2) THEN
               xt=x(jt)*(1.+(dsh-q2b)/dq2(jr))
            ELSE
               xt=xb
            ENDIF
            CALL lystfu(ipy(40+jt),xt,q2ref,xft)
            IF(xft(iflb).EQ.0.0.AND.xft(-iflb).EQ.0.0) THEN
              searatio=1.0
            ELSEIF(xft(abs(iflb)).EQ.0.0) THEN
              searatio=1.0
            ELSE
              searatio=xft(-abs(iflb))/xft(abs(iflb))
            ENDIF
CJR-- (protons only)
            IF (rlu(0).LT.searatio) THEN
               seaquark=.true.
               xquark=xt
            ELSE
               q2b=0.
               seaquark=.false.
            ENDIF
         ENDIF
      ENDIF
CJR--end
      IF(q2b.LT.pypar(22)) THEN
         q2b=0.
CJR--begin
         seaquark=.false.
CJR--end
      ELSE
        wtran=rlu(0)*wtsum
        ifla=-ipy(8)-1
  200   ifla=ifla+1
        wtran=wtran-wtap(ifla)*wtsf(ifla)
        IF(ifla.LT.ipy(8).AND.wtran.GT.0.) goto 200

CJR--begin
        IF (seaquark) THEN
           seafl=-iflb
           ifla=0
CT           XE=XB*(1./(1.-0.001)-1.)
        ENDIF
CJR--end

C...CHOOSE Z VALUE AND CORRECTIVE WEIGHT
        IF(iflb.EQ.0.AND.ifla.EQ.0) THEN
          z=1./(1.+((1.-xb)/xb)*(xe/(1.-xb))**rlu(0))
          wtz=(1.-z*(1.-z))**2
        ELSEIF(iflb.EQ.0) THEN
          z=xb/(1.-rlu(0)*(1.-sqrt(xb+xe)))**2
          wtz=0.5*(1.+(1.-z)**2)*sqrt(z)
        ELSEIF(ifla.EQ.0) THEN
          z=xb*(1.+rlu(0)*(1./(xb+xe)-1.))
          wtz=1.-2.*z*(1.-z)
        ELSE
          z=1.-(1.-xb)*(xe/((xb+xe)*(1.-xb)))**rlu(0)
          wtz=0.5*(1.+z**2)
        ENDIF
CJR--begin
C        IF (SEAQUARK) THEN
C           XSEA=LEXSEA(0.15*XT,Q2B)
C           XE=MIN(XE,XSEA)
C           Z=XT/(XSEA+XT)
C        ENDIF
CJR--end
C...REWEIGHT FIRST LEG BECAUSE OF MODIFIED XB OR CHECK PHASE SPACE
        IF(ilep.GE.1.AND.n.EQ.ns+2) THEN
          xbnew=x(jt)*(1.+(dsh-q2b)/dq2(jr))
          IF(xbnew.GT.min(z,0.999)) goto 190
          xb=xbnew
        ENDIF
C...SUM UP SOFT GLUON EMISSION AS EFFECTIVE Z SHIFT
CJR--       should this realy always be done ??
        IF(ipy(15).GE.1) THEN
          rsoft=6.
          IF(iflb.NE.0) rsoft=8./3.
          z=z*(tevb/tevs(jt))**(rsoft*xe/((xb+xe)*b0))
          IF(z.LE.xb) goto 190
CJR--begin
          zsoft=(tevb/tevs(jt))**(rsoft*xe/((xb+xe)*b0))
          zmax=xb/(xb+xe)
CJR--end
        ENDIF
C...CHECK IF HEAVY FLAVOUR BELOW THRESHOLD
        ihft=0
CIC...Skip for intrinsic charm/bottom simulation, charm quark should
CIC...not come from gluon but is non-perturbative part of proton.
        IF(lst(15).EQ.-4.OR.lst(15).EQ.-5) goto 205
        mstj(93)=1
        IF(ilep.GE.1.AND.iabs(iflb).GE.4.AND.(xfb(iflb).LT.1e-10.OR.
     &    q2b.LT.5.*ulmass(iabs(iflb))**2)) THEN
          ihft=1
          ifla=0
        ENDIF
  205   CONTINUE

C...FOR LEPTOPRODUCTION, CHECK Z AGAINST NEW LIMIT
        IF(ilep.GE.1) THEN
          dpd(2)=dshz+dq2(jr)+dq2b
          mstj(93)=1
          dm2=ulmass(iabs(ifla-iflb))**2
          IF(iabs(ifla-iflb).EQ.0) dm2=ulmass(21)**2
          dpc(1)=dq2(jr)*(dq2b+dm2)**2
          dpc(2)=dpd(2)*(dpd(2)-2d0*dq2b)*(dq2b+dm2)
          dpc(3)=dq2b*(dpd(2)-2d0*dq2b)**2
          zu=(dpc(2)-dsqrt(dpc(2)**2-4d0*dpc(1)*dpc(3)))/(2d0*dpc(1))
          IF(z.GE.zu) goto 190
        ENDIF

C...OPTION WITH EVOLUTION IN KT2=(1-Z)Q2:
        IF(ipy(14).GE.5.AND.ipy(14).LE.6.AND.n.LE.ns+4) THEN
C...CHECK THAT (Q2)LAST BRANCHING < (Q2)HARD
          IF(q2b/(1.-z).GT.pypar(26)*q2) goto 190
        ELSEIF(ipy(14).GE.3.AND.ipy(14).LE.6.AND.n.GE.ns+6) THEN
C...CHECK THAT Z,Q2 COMBINATION IS KINEMATICALLY ALLOWED
          q2max=0.5*(1./zs(jt)+1.)*dq2(jt)+0.5*(1./zs(jt)-1.)*
     &    (dq2(3-jt)-dsh+sqrt((dsh+dq2(1)+dq2(2))**2+8.*dq2(1)*dq2(2)*
     &    zs(jt)/(1.-zs(jt))))
          IF(q2b/(1.-z).GE.q2max) goto 190

        ELSEIF(ipy(14).EQ.7.OR.ipy(14).EQ.8) THEN
C...OPTION WITH ALPHAS((1-Z)Q2): DEMAND KT2 > CUTOFF, REWEIGHT
          IF((1.-z)*q2b.LT.pypar(22)) goto 190
          alprat=tevb/(tevb+alog(1.-z))
          IF(alprat.LT.5.*rlu(0)) goto 190
          IF(alprat.GT.5.) wtz=wtz*alprat/5.
        ENDIF

C...WEIGHTING WITH NEW STRUCTURE FUNCTIONS
        CALL lystfu(ipy(40+jt),xb,q2ref,xfb)
        xa=xb/z
        CALL lystfu(ipy(40+jt),xa,q2ref,xfa)
        IF(ihft.EQ.1.OR.ihfx.EQ.1) THEN
           IF(xfa(ifla).LT.1e-10) ihfc=1
           goto 210
        ELSEIF(xfb(iflb).LT.1e-20) THEN
          goto 100
        ENDIF
        IF(wtz*xfa(ifla)/xfb(iflb).LT.rlu(0)*wtsf(ifla)) THEN
          IF(ilep.GE.1.AND.n.EQ.ns+2) goto 145
          goto 170
        ENDIF
CJR--begin
        IF (seaquark) THEN
           split=.true.
           xgluon=xa
           xsea=xa-xb
           zsplit=z
           seaquark=.false.
        ENDIF
CJR--end
      ENDIF

210   CONTINUE
      IF(n.EQ.ns+4-2*min(1,ilep)) THEN
C...DEFINE TWO HARD SCATTERERS IN THEIR CM-FRAME
        dq2(jt)=q2b
        IF(ipy(14).GE.3.AND.ipy(14).LE.6) dq2(jt)=q2b/(1.-z)
        IF(ilep.EQ.0) THEN
          dplcm=dsqrt((dsh+dq2(1)+dq2(2))**2-4.*dq2(1)*dq2(2))/dshr
          DO 220 jr=1,2
          i=ns+2*jr-1
          ipo=19+2*jr
          k(i,1)=14
          k(i,2)=ifls(jr+2)
          IF(ifls(jr+2).EQ.0) k(i,2)=21
          k(i,3)=0
          k(i,4)=ipo
          k(i,5)=ipo
          p(i,1)=0.
          p(i,2)=0.
          p(i,3)=dplcm*(-1)**(jr+1)
          p(i,4)=(dsh+dq2(3-jr)-dq2(jr))/dshr
          p(i,5)=-sqrt(sngl(dq2(jr)))
          k(i+1,1)=-1
          k(i+1,2)=k(ipo+1,2)
          k(i+1,3)=i
          k(i+1,4)=0
          k(i+1,5)=0
          p(i+1,1)=0.
          p(i+1,2)=0.
          p(i+1,3)=ipo
          p(i+1,4)=ipo
          p(i+1,5)=0.
          p(ipo+1,1)=i
          p(ipo+1,2)=i
          k(ipo,4)=mod(k(ipo,4),mstu(5))+i*mstu(5)
          k(ipo,5)=mod(k(ipo,5),mstu(5))+i*mstu(5)
  220     CONTINUE
        ELSE
C..LEPTOPRODUCTION EVENTS: BOSON AND HADRON REST FRAME
          i1=ns+2*ilep-1
          i2=ns-2*ilep+5
          DO 225 itemp=ns+1,ns+4
          DO 225 j=1,5
          k(itemp,j)=0
  225     p(itemp,j)=0.
          DO 230 j=1,5
  230     p(i1,j)=p(nq,j)
          k(ns+1,1)=11
          k(ns+3,1)=14
          IF(ilep.EQ.2) THEN
            k(ns+1,1)=14
            k(ns+3,1)=11
          ENDIF
          k(ns+2,1)=-1
          k(ns+4,1)=-1
          k(ns+1,3)=0
          k(ns+2,3)=ns+1
          k(ns+3,3)=0
          k(ns+4,3)=ns+3
          k(i1,2)=kfl(2,ilep)
          k(i2,2)=kfl(2,3-ilep)
          dpd(1)=dsh+dq2(1)+dq2(2)
          dpd(3)=(3-2*ilep)*dsqrt(dpd(1)**2-4d0*dq2(1)*dq2(2))
          p(i2,3)=(dpqs(2)*dpd(3)-dpqs(1)*dpd(1))/
     &    (2d0*dq2(jr))
          p(i2,4)=(dpqs(1)*dpd(3)-dpqs(2)*dpd(1))/
     &    (2d0*dq2(jr))
          p(i2,5)=-sqrt(sngl(dq2(3-ilep)))
          p(i2+1,3)=max(ipu1,ipu2)
          p(i2+1,4)=max(ipu1,ipu2)
          k(i2,4)=k(i2,4)-mod(k(i2,4),mstu(5))+max(ipu1,ipu2)
          k(i2,5)=k(i2,5)-mod(k(i2,5),mstu(5))+max(ipu1,ipu2)
          p(26-2*ilep,1)=i2
          p(26-2*ilep,2)=i2
          k(25-2*ilep,4)=mod(k(25-2*ilep,4),mstu(5))+i2*mstu(5)
          k(25-2*ilep,5)=mod(k(25-2*ilep,5),mstu(5))+i2*mstu(5)
          n=n+2
        ENDIF

      ELSEIF(n.GT.ns+4) THEN
C...FIND MAXIMUM ALLOWED MASS OF TIMELIKE PARTON
        dq2(3)=q2b
        IF(ipy(14).GE.3.AND.ipy(14).LE.6) dq2(3)=q2b/(1.-z)
        IF(is(1).GE.1.AND.is(1).LE.mstu(4)) THEN
          dpc(1)=p(is(1),4)
          dpc(3)=0.5*(abs(p(is(1),3))+abs(p(is(2),3)))
        ELSE
C...IS(1) not initialized
          dpc(1)=0.
          dpc(3)=0.5*(       0.      +abs(p(is(2),3)))
        ENDIF
        dpc(2)=p(is(2),4)
        dpd(1)=dsh+dq2(jr)+dq2(jt)
        dpd(2)=dshz+dq2(jr)+dq2(3)
        dpd(3)=dsqrt(dpd(1)**2-4.*dq2(jr)*dq2(jt))
        dpd(4)=dsqrt(dpd(2)**2-4.*dq2(jr)*dq2(3))
        ikin=0
        IF((q2s(jr).GE.0.5*pypar(22).AND.dpd(1)-dpd(3).GE.1d-10*dpd(1))
     &  .OR.ilep.GE.1) ikin=1
        IF(ikin.EQ.0) dmsma=(dq2(jt)/dble(zs(jt))-dq2(3))*(dsh/
     &  (dsh+dq2(jt))-dsh/(dshz+dq2(3)))
        IF(ikin.EQ.1) dmsma=(dpd(1)*dpd(2)-dpd(3)*dpd(4))/(2.*
     &  dq2(jr))-dq2(jt)-dq2(3)

C...GENERATE TIMELIKE PARTON SHOWER (IF REQUIRED)
        it=n-1
        k(it,1)=3
        k(it,2)=iflb-ifls(jt+2)
        IF(iflb-ifls(jt+2).EQ.0) k(it,2)=21
        mstj(93)=1
        p(it,5)=ulmass(k(it,2))
        IF(sngl(dmsma).LE.p(it,5)**2) goto 100
        p(it,2)=0.
        DO 240 j=1,5
        k(it+1,j)=0
  240   p(it+1,j)=0.
        k(it+1,1)=-1
        k(it+1,2)=k(is(jt)+1,2)
        k(it+1,3)=it
        IF(mod(ipy(14),2).EQ.0) THEN
          p(it,1)=0.
          IF(ilep.EQ.0) p(it,4)=(dshz-dsh-p(it,5)**2)/dshr
          IF(ilep.GE.1) p(it,4)=0.5*(p(is(jt),3)*dpd(2)+
     &    dpqs(1)*(dq2(jt)+dq2(3)+p(it,5)**2))/(p(is(jt),3)*dpqs(2)-
     &    p(is(jt),4)*dpqs(1))-dpc(jt)
          p(it,3)=sqrt(max(0.,p(it,4)**2-p(it,5)**2))
          CALL lushow(it,0,sqrt(min(sngl(dmsma),pypar(25)*q2)))
          IF(n.GE.it+2) p(it,5)=p(it+2,5)
          IF(n.GT.mstu(4)-10) THEN
            WRITE(6,*) ' LYSSPA: no more memory in LUJETS'
            lst(21)=19
            RETURN
          ENDIF
          DO 243 i=n+1,n+8
          DO 243 j=1,5
          k(i,j)=0
  243     p(i,j)=0.
        ENDIF

C...RECONSTRUCT KINEMATICS OF BRANCHING: TIMELIKE PARTON SHOWER
        dms=p(it,5)**2
        IF(ikin.EQ.0.AND.ilep.EQ.0) dpt2=(dmsma-dms)*(dshz+dq2(3))/
     &  (dsh+dq2(jt))
        IF(ikin.EQ.1.AND.ilep.EQ.0) dpt2=(dmsma-dms)*(0.5*dpd(1)*
     &  dpd(2)+0.5*dpd(3)*dpd(4)-dq2(jr)*(dq2(jt)+dq2(3)+dms))/
     &  (4.*dsh*dpc(3)**2)
        IF(ikin.EQ.1.AND.ilep.GE.1) dpt2=(dmsma-dms)*(0.5*dpd(1)*
     &  dpd(2)+0.5*dpd(3)*dpd(4)-dq2(jr)*(dq2(jt)+dq2(3)+dms))/
     &  dpd(3)**2
        IF(dpt2.LT.0.) goto 100
        k(it,3)=n+1
        p(it,1)=sqrt(sngl(dpt2))
        IF(ilep.EQ.0) THEN
          dpb(1)=(0.5*dpd(2)-dpc(jr)*(dshz+dq2(jr)-dq2(jt)-dms)/
     &    dshr)/dpc(3)-dpc(3)
          p(it,3)=dpb(1)*(-1)**(jt+1)
          p(it,4)=(dshz-dsh-dms)/dshr
        ELSE
          dpc(3)=dq2(jt)+dq2(3)+dms
          dpb(2)=dpqs(2)*dble(p(is(jt),3))-dpqs(1)*dpc(jt)
          dpb(1)=0.5d0*(dpc(jt)*dpd(2)+dpqs(2)*dpc(3))/dpb(2)-
     &    dble(p(is(jt),3))
          p(it,3)=dpb(1)
          p(it,4)=0.5d0*(dble(p(is(jt),3))*dpd(2)+
     &    dpqs(1)*dpc(3))/dpb(2)-dpc(jt)
        ENDIF
        IF(n.GE.it+2) THEN
          mstu(1)=it+2
          dpb(1)=dsqrt(dpb(1)**2+dpt2)
          dpb(2)=dsqrt(dpb(1)**2+dms)
          dpb(3)=p(it+2,3)
          dpb(4)=dsqrt(dpb(3)**2+dms)
          dbez=(dpb(4)*dpb(1)-dpb(3)*dpb(2))/(dpb(4)*dpb(2)-dpb(3)*
     &    dpb(1))
          CALL ludbrb(mstu(1),mstu(2),0.,0.,0.d0,0.d0,dbez)
          the=ulangl(p(it,3),p(it,1))
          CALL ludbrb(mstu(1),mstu(2),the,0.,0.d0,0.d0,0.d0)
        ENDIF

C...RECONSTRUCT KINEMATICS OF BRANCHING: SPACELIKE PARTON
        k(n+1,1)=14
        k(n+1,2)=iflb
        IF(iflb.EQ.0) k(n+1,2)=21
        k(n+1,3)=0
CJR--begin 
CJR--       give all radiated partons 5 as mother particle 
        k(n+1,3)=5
CJR--end
        p(n+1,1)=p(it,1)
        p(n+1,2)=0.
        p(n+1,3)=p(it,3)+p(is(jt),3)
        p(n+1,4)=p(it,4)+p(is(jt),4)
        p(n+1,5)=-sqrt(sngl(dq2(3)))
        DO 250 j=1,5
        k(n+2,j)=0
  250   p(n+2,j)=0.
        k(n+2,1)=-1
        k(n+2,2)=k(is(jt)+1,2)
        k(n+2,3)=n+1

C...DEFINE COLOUR FLOW OF BRANCHING
        k(is(jt),1)=14
        k(is(jt),3)=n+1
        id1=it
      kn1=isign(500+iabs(k(n+1,2)),2*k(n+1,2)+1)
      kd1=isign(500+iabs(k(id1,2)),2*k(id1,2)+1)
      IF(k(n+1,2).EQ.21) kn1=500
      IF(k(id1,2).EQ.21) kd1=500
        IF((kn1.GE.501.AND.kd1.GE.501).OR.(kn1.LT.0.AND.
     &  kd1.EQ.500).OR.(kn1.EQ.500.AND.kd1.EQ.500.AND.
     &  rlu(0).GT.0.5).OR.(kn1.EQ.500.AND.kd1.LT.0))
     &  id1=is(jt)
        id2=it+is(jt)-id1
        p(n+2,3)=id1
        p(n+2,4)=id2
        p(id1+1,1)=n+1
        p(id1+1,2)=id2
        p(id2+1,1)=id1
        p(id2+1,2)=n+1
        k(n+1,4)=k(n+1,4)-mod(k(n+1,4),mstu(5))+id1
        k(n+1,5)=k(n+1,5)-mod(k(n+1,5),mstu(5))+id2
        k(id1,4)=mod(k(id1,4),mstu(5))+(n+1)*mstu(5)
        k(id1,5)=mod(k(id1,5),mstu(5))+id2*mstu(5)
        k(id2,4)=mod(k(id2,4),mstu(5))+id1*mstu(5)
        k(id2,5)=mod(k(id2,5),mstu(5))+(n+1)*mstu(5)
        n=n+2
C     CALL GULIST(22,2)

C...BOOST TO NEW CM-FRAME
        mstu(1)=ns+1
        IF(ilep.EQ.0) THEN
          CALL ludbrb(mstu(1),mstu(2),0.,0.,
     &    -dble(p(n-1,1)+p(is(jr),1))/dble(p(n-1,4)+p(is(jr),4)),
     &    0.d0,-dble(p(n-1,3)+p(is(jr),3))/dble(p(n-1,4)+p(is(jr),4)))
          ir=n-1+(jt-1)*(is(1)-n+1)
          CALL ludbrb(mstu(1),mstu(2),
     &    -ulangl(p(ir,3),p(ir,1)),paru(2)*rlu(0),0.d0,0.d0,0.d0)
        ELSE
C...REORIENTATE EVENT WITHOUT CHANGING THE BOSON FOUR MOMENTUM
          DO 260 j=1,4
  260     dpq(j)=p(nq,j)
          dbe1(4)=dpq(4)+dble(p(n-1,4))
          DO 270 j=1,3,2
  270     dbe1(j)=-(dpq(j)+dble(p(n-1,j)))/dbe1(4)
          dtemp=1.d0-dbe1(1)**2-dbe1(3)**2
          IF(dtemp.LE.0.d0) THEN
            lst(21)=20
            IF(lst(3).GE.1) WRITE(6,*) ' Warning from LYSSPA: sqrt of',
     &      dtemp,' New event generated.'
            RETURN
          ENDIF
          dbe1(4)=1.d0/dsqrt(dtemp)
          dbep=dbe1(1)*dpq(1)+dbe1(3)*dpq(3)
          dgabep=dbe1(4)*(dbe1(4)*dbep/(1d0+dbe1(4))+dpq(4))
          DO 280 j=1,3,2
  280     dpq(j)=dpq(j)+dgabep*dbe1(j)
          dpq(4)=dbe1(4)*(dpq(4)+dbep)
          dpc(1)=dsqrt(dpq(1)**2+dpq(3)**2)
          dbe2(4)=-(dpq(4)*dpc(1)-dpqs(2)*dsqrt(dpqs(2)**2+dpc(1)**2-
     &    dpq(4)**2))/(dpc(1)**2+dpqs(2)**2)
          the=ulangl(sngl(dpq(3)),sngl(dpq(1)))
          dbe2(1)=dbe2(4)*dsin(dble(the))
          dbe2(3)=dbe2(4)*dcos(dble(the))
          dbe2(4)=1d0/(1d0-dbe2(1)**2-dbe2(3)*2)

C...CONSTRUCT THE COMBINED BOOST
          dpb(1)=dbe1(4)**2*dbe2(4)/(1d0+dbe1(4))
          dpb(2)=dbe1(1)*dbe2(1)+dbe1(3)*dbe2(3)
          dpb(3)=dbe1(4)*dbe2(4)*(1d0+dpb(2))
          DO 290 jb=1,3,2
  290     drobo(jb+2)=(dbe1(4)*dbe2(4)*dbe1(jb)+dbe2(4)*dbe2(jb)+
     &    dpb(1)*dbe1(jb)*dpb(2))/dpb(3)
          CALL ludbrb(mstu(1),mstu(2),0.,0.,drobo(3),0.d0,drobo(5))
          IF(ilep.EQ.1) the=ulangl(p(ns+1,3),p(ns+1,1))
          IF(ilep.EQ.2) the=paru(1)+ulangl(p(ns+3,3),p(ns+3,1))
          CALL ludbrb(mstu(1),mstu(2),-the,paru(2)*rlu(0),0d0,0d0,0d0)
        ENDIF
        mstu(1)=0
      ENDIF

C...SAVE QUANTITIES, LOOP BACK
      is(jt)=n-1
      IF(ilep.EQ.2.AND.n.EQ.ns+4) is(jt)=n-3
      q2s(jt)=q2b
      dq2(jt)=q2b
      IF(ipy(14).GE.3.AND.ipy(14).LE.6) dq2(jt)=q2b/(1.-z)
      dsh=dshz
      IF(q2b.GE.0.5*pypar(22)) THEN
        ifls(jt+2)=ifls(jt)
        ifls(jt)=ifla
        xs(jt)=xa
        zs(jt)=z
        DO 300 ifl=-6,6
  300   xfs(jt,ifl)=xfa(ifl)
        tevs(jt)=tevb
      ELSE
        IF(jt.EQ.1) ipu1=n-1
        IF(jt.EQ.2) ipu2=n-1
      ENDIF
      IF(max(iabs(1-ilep)*q2s(1),min(1,2-ilep)*q2s(2)).GE.0.5*pypar(22)
     &.OR.n.LE.ns+2) goto 130

      IF(ilep.EQ.0) THEN
C...BOOST HARD SCATTERING PARTONS TO FRAME OF SHOWER INITIATORS
        DO 310 j=1,3
  310   drobo(j+2)=(p(ns+1,j)+p(ns+3,j))/(p(ns+1,4)+p(ns+3,4))
        DO 320 j=1,5
  320   p(n+2,j)=p(ns+1,j)
        mstu(1)=n+2
        mstu(2)=n+2
        CALL ludbrb(n+2,n+2,0.,0.,-drobo(3),-drobo(4),-drobo(5))
        robo(2)=ulangl(p(n+2,1),p(n+2,2))
        robo(1)=ulangl(p(n+2,3),sqrt(p(n+2,1)**2+p(n+2,2)**2))
        mstu(1)=4
        mstu(2)=ns
        CALL ludbrb(4,ns,robo(1),robo(2),drobo(3),drobo(4),drobo(5))
        mstu(1)=0
        mstu(2)=0
      ENDIF

C...STORE USER INFORMATION
      k(21,1)=14
      IF(ilep.NE.0) k(21,1)=11
      k(23,1)=14
      k(21,3)=ns+1
      k(23,3)=ns+3
      DO 330 jt=1,2
      kfl(1,jt)=ifls(jt)
      IF(ifls(jt).EQ.0) kfl(1,jt)=21
  330 pyvar(30+jt)=xs(jt)

      DO 340 i=ns+1,n
      DO 340 j=1,5
  340 v(i,j)=0.
  
CJR--begin
      lastfl=ifla
CJR--end

      RETURN
      END