*
* $Id: pysspa.F,v 1.1.1.1 1996/01/11 14:05:27 mclareni Exp $
*
* $Log: pysspa.F,v $
* Revision 1.1.1.1  1996/01/11 14:05:27  mclareni
* Fritiof
*
*
C*********************************************************************
 
      SUBROUTINE PYSSPA(IPU1,IPU2)
 
C...Generates spacelike parton showers.
      IMPLICIT DOUBLE PRECISION(D)
      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/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
      COMMON/PYINT1/MINT(400),VINT(400)
      COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
      COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
      SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
      SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/
      DIMENSION KFLS(4),IS(2),XS(2),ZS(2),Q2S(2),TEVS(2),ROBO(5),
     &XFS(2,-25:25),XFA(-25:25),XFB(-25:25),XFN(-25:25),WTAP(-25:25),
     &WTSF(-25:25),THE2(2),ALAM(2),DQ2(3),DPC(3),DPD(4),DPB(4)
 
C...Calculate maximum virtuality and check that evolution possible.
      IPUS1=IPU1
      IPUS2=IPU2
      ISUB=MINT(1)
      Q2E=VINT(52)
      IF(ISET(ISUB).EQ.1) THEN
        Q2E=Q2E/PARP(67)
      ELSEIF(ISET(ISUB).GE.3.AND.ISET(ISUB).LE.5) THEN
        Q2E=PMAS(23,1)**2
        IF(ISUB.EQ.8.OR.ISUB.EQ.76.OR.ISUB.EQ.77.OR.ISUB.EQ.124.OR.
     &  ISUB.EQ.174.OR.ISUB.EQ.179) Q2E=PMAS(24,1)**2
      ENDIF
      TMAX=LOG(PARP(67)*PARP(63)*Q2E/PARP(61)**2)
      IF(PARP(67)*Q2E.LT.MAX(PARP(62)**2,2.*PARP(61)**2).OR.
     &TMAX.LT.0.2) RETURN
 
C...Common constants and initial values. Save normal Lambda value.
      XE0=2.*PARP(65)/VINT(1)
      ALAMS=PARU(111)
      PARU(111)=PARP(61)
      NS=N
  100 N=NS
      DO 110 JT=1,2
      KFLS(JT)=MINT(14+JT)
      KFLS(JT+2)=KFLS(JT)
      XS(JT)=VINT(40+JT)
      ZS(JT)=1.
      Q2S(JT)=PARP(67)*Q2E
      TEVS(JT)=TMAX
      ALAM(JT)=PARP(61)
      THE2(JT)=100.
      DO 110 KFL=-25,25
  110 XFS(JT,KFL)=XSFX(JT,KFL)
      DSH=VINT(44)
      IF(ISET(ISUB).GE.3.AND.ISET(ISUB).LE.5) DSH=VINT(26)*VINT(2)
 
C...Pick up leg with highest virtuality.
  120 N=N+1
      JT=1
      IF(N.GT.NS+1.AND.Q2S(2).GT.Q2S(1)) JT=2
      KFLB=KFLS(JT)
      XB=XS(JT)
      DO 130 KFL=-25,25
  130 XFB(KFL)=XFS(JT,KFL)
      DSHR=2D0*SQRT(DSH)
      DSHZ=DSH/DBLE(ZS(JT))
      XE=MAX(XE0,XB*(1./(1.-PARP(66))-1.))
      IF(MINT(40+JT).EQ.1.OR.XB+XE.GE.0.999) THEN
        Q2B=0.
        GOTO 220
      ENDIF
 
C...Maximum Q2 without or with Q2 ordering. Effective Lambda and n_f.
      IF(MSTP(62).LE.1) 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=LOG(PARP(63)*Q2B/ALAM(JT)**2)
      ELSE
        Q2B=Q2S(JT)
        TEVB=TEVS(JT)
      ENDIF
      ALSDUM=ULALPS(PARP(63)*Q2B)
      TEVB=TEVB+2.*LOG(ALAM(JT)/PARU(117))
      TEVBSV=TEVB
      ALAM(JT)=PARU(117)
      B0=(33.-2.*MSTU(118))/6.
 
C...Calculate Altarelli-Parisi and structure function weights.
      DO 140 KFL=-25,25
      WTAP(KFL)=0.
  140 WTSF(KFL)=0.
      IF(KFLB.EQ.21) THEN
        WTAPQ=16.*(1.-SQRT(XB+XE))/(3.*SQRT(XB))
        DO 150 KFL=-MSTP(54),MSTP(54)
        IF(KFL.EQ.0) WTAP(KFL)=6.*LOG((1.-XB)/XE)
  150   IF(KFL.NE.0) WTAP(KFL)=WTAPQ
      ELSE
        WTAP(0)=0.5*XB*(1./(XB+XE)-1.)
        WTAP(KFLB)=8.*LOG((1.-XB)*(XB+XE)/XE)/3.
      ENDIF
  160 WTSUM=0.
      IF(KFLB.NE.21) XFBO=XFB(KFLB)
      IF(KFLB.EQ.21) XFBO=XFB(0)
      DO 170 KFL=-MSTP(54),MSTP(54)
      WTSF(KFL)=XFB(KFL)/XFBO
  170 WTSUM=WTSUM+WTAP(KFL)*WTSF(KFL)
      WTSUM=MAX(0.0001,WTSUM)
 
C...Choose new t: fix alpha_s, alpha_s(Q^2), alpha_s(k_T^2).
  180 IF(MSTP(64).LE.0) THEN
        TEVB=TEVB+LOG(RLU(0))*PARU(2)/(PARU(111)*WTSUM)
      ELSEIF(MSTP(64).EQ.1) THEN
        TEVB=TEVB*EXP(MAX(-70.,LOG(RLU(0))*B0/WTSUM))
      ELSE
        TEVB=TEVB*EXP(MAX(-70.,LOG(RLU(0))*B0/(5.*WTSUM)))
      ENDIF
  190 Q2REF=ALAM(JT)**2*EXP(TEVB)
      Q2B=Q2REF/PARP(63)
 
C...Evolution ended or select flavour for branching parton.
      IF(Q2B.LT.PARP(62)**2) THEN
        Q2B=0.
      ELSE
        WTRAN=RLU(0)*WTSUM
        KFLA=-MSTP(54)-1
  200   KFLA=KFLA+1
        WTRAN=WTRAN-WTAP(KFLA)*WTSF(KFLA)
        IF(KFLA.LT.MSTP(54).AND.WTRAN.GT.0.) GOTO 200
        IF(KFLA.EQ.0) KFLA=21
 
C...Choose z value and corrective weight.
        IF(KFLB.EQ.21.AND.KFLA.EQ.21) THEN
          Z=1./(1.+((1.-XB)/XB)*(XE/(1.-XB))**RLU(0))
          WTZ=(1.-Z*(1.-Z))**2
        ELSEIF(KFLB.EQ.21) THEN
          Z=XB/(1.-RLU(0)*(1.-SQRT(XB+XE)))**2
          WTZ=0.5*(1.+(1.-Z)**2)*SQRT(Z)
        ELSEIF(KFLA.EQ.21) 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
 
C...Option with resummation of soft gluon emission as effective z shift.
        IF(MSTP(65).GE.1) THEN
          RSOFT=6.
          IF(KFLB.NE.21) RSOFT=8./3.
          Z=Z*(TEVB/TEVS(JT))**(RSOFT*XE/((XB+XE)*B0))
          IF(Z.LE.XB) GOTO 180
        ENDIF
 
C...Option with alpha_s(k_T^2): demand k_T^2 > cutoff, reweight.
        IF(MSTP(64).GE.2) THEN
          IF((1.-Z)*Q2B.LT.PARP(62)**2) GOTO 180
          ALPRAT=TEVB/(TEVB+LOG(1.-Z))
          IF(ALPRAT.LT.5.*RLU(0)) GOTO 180
          IF(ALPRAT.GT.5.) WTZ=WTZ*ALPRAT/5.
        ENDIF
 
C...Option with angular ordering requirement.
        IF(MSTP(62).GE.3) THEN
          THE2T=(4.*Z**2*Q2B)/(VINT(2)*(1.-Z)*XB**2)
          IF(THE2T.GT.THE2(JT)) GOTO 180
        ENDIF
 
C...Weighting with new structure functions.
        CALL PYSTFU(MINT(10+JT),XB,Q2REF,XFN)
        IF(KFLB.NE.21) XFBN=XFN(KFLB)
        IF(KFLB.EQ.21) XFBN=XFN(0)
        IF(XFBN.LT.1E-20) THEN
          IF(KFLA.EQ.KFLB) THEN
            TEVB=TEVBSV
            WTAP(KFLB)=0.
            GOTO 160
          ELSEIF(TEVBSV-TEVB.GT.0.2) THEN
            TEVB=0.5*(TEVBSV+TEVB)
            GOTO 190
          ELSE
            XFBN=1E-10
            IF(KFLB.NE.21) XFN(KFLB)=XFBN
            IF(KFLB.EQ.21) XFN(0)=XFBN
          ENDIF
        ENDIF
        DO 210 KFL=-MSTP(54),MSTP(54)
  210   XFB(KFL)=XFN(KFL)
        XA=XB/Z
        CALL PYSTFU(MINT(10+JT),XA,Q2REF,XFA)
        IF(KFLA.NE.21) XFAN=XFA(KFLA)
        IF(KFLA.EQ.21) XFAN=XFA(0)
        IF(XFAN.LT.1E-20) GOTO 160
        IF(KFLA.NE.21) WTSFA=WTSF(KFLA)
        IF(KFLA.EQ.21) WTSFA=WTSF(0)
        IF(WTZ*XFAN/XFBN.LT.RLU(0)*WTSFA) GOTO 160
      ENDIF
 
C...Define two hard scatterers in their CM-frame.
  220 IF(N.EQ.NS+2) THEN
        DQ2(JT)=Q2B
        DPLCM=SQRT((DSH+DQ2(1)+DQ2(2))**2-4D0*DQ2(1)*DQ2(2))/DSHR
        DO 240 JR=1,2
        I=NS+JR
        IF(JR.EQ.1) IPO=IPUS1
        IF(JR.EQ.2) IPO=IPUS2
        DO 230 J=1,5
        K(I,J)=0
        P(I,J)=0.
  230   V(I,J)=0.
        K(I,1)=14
        K(I,2)=KFLS(JR+2)
        K(I,4)=IPO
        K(I,5)=IPO
        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(IPO,1)=14
        K(IPO,3)=I
        K(IPO,4)=MOD(K(IPO,4),MSTU(5))+MSTU(5)*I
  240   K(IPO,5)=MOD(K(IPO,5),MSTU(5))+MSTU(5)*I
 
C...Find maximum allowed mass of timelike parton.
      ELSEIF(N.GT.NS+2) THEN
        JR=3-JT
        DQ2(3)=Q2B
        DPC(1)=P(IS(1),4)
        DPC(2)=P(IS(2),4)
        DPC(3)=0.5*(ABS(P(IS(1),3))+ABS(P(IS(2),3)))
        DPD(1)=DSH+DQ2(JR)+DQ2(JT)
        DPD(2)=DSHZ+DQ2(JR)+DQ2(3)
        DPD(3)=SQRT(DPD(1)**2-4D0*DQ2(JR)*DQ2(JT))
        DPD(4)=SQRT(DPD(2)**2-4D0*DQ2(JR)*DQ2(3))
        IKIN=0
        IF(Q2S(JR).GE.(0.5*PARP(62))**2.AND.DPD(1)-DPD(3).GE.
     &  1D-10*DPD(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
        DO 250 J=1,5
        K(IT,J)=0
        P(IT,J)=0.
  250   V(IT,J)=0.
        K(IT,1)=3
        K(IT,2)=21
        IF(KFLB.EQ.21.AND.KFLS(JT+2).NE.21) K(IT,2)=-KFLS(JT+2)
        IF(KFLB.NE.21.AND.KFLS(JT+2).EQ.21) K(IT,2)=KFLB
        P(IT,5)=ULMASS(K(IT,2))
        IF(SNGL(DMSMA).LE.P(IT,5)**2) GOTO 100
        IF(MSTP(63).GE.1) THEN
          P(IT,4)=(DSHZ-DSH-P(IT,5)**2)/DSHR
          P(IT,3)=SQRT(P(IT,4)**2-P(IT,5)**2)
          IF(MSTP(63).EQ.1) THEN
            Q2TIM=DMSMA
          ELSEIF(MSTP(63).EQ.2) THEN
            Q2TIM=MIN(SNGL(DMSMA),PARP(71)*Q2S(JT))
          ELSE
C'''Here remains to introduce angular ordering in first branching.
            Q2TIM=DMSMA
          ENDIF
          CALL LUSHOW(IT,0,SQRT(Q2TIM))
          IF(N.GE.IT+1) P(IT,5)=P(IT+1,5)
        ENDIF
 
C...Reconstruct kinematics of branching: timelike parton shower.
        DMS=P(IT,5)**2
        IF(IKIN.EQ.0) DPT2=(DMSMA-DMS)*(DSHZ+DQ2(3))/(DSH+DQ2(JT))
        IF(IKIN.EQ.1) 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(DPT2.LT.0.) GOTO 100
        DPB(1)=(0.5*DPD(2)-DPC(JR)*(DSHZ+DQ2(JR)-DQ2(JT)-DMS)/
     &  DSHR)/DPC(3)-DPC(3)
        P(IT,1)=SQRT(SNGL(DPT2))
        P(IT,3)=DPB(1)*(-1)**(JT+1)
        P(IT,4)=(DSHZ-DSH-DMS)/DSHR
        IF(N.GE.IT+1) THEN
          DPB(1)=SQRT(DPB(1)**2+DPT2)
          DPB(2)=SQRT(DPB(1)**2+DMS)
          DPB(3)=P(IT+1,3)
          DPB(4)=SQRT(DPB(3)**2+DMS)
          DBEZ=(DPB(4)*DPB(1)-DPB(3)*DPB(2))/(DPB(4)*DPB(2)-DPB(3)*
     &    DPB(1))
          CALL LUDBRB(IT+1,N,0.,0.,0D0,0D0,DBEZ)
          THE=ULANGL(P(IT,3),P(IT,1))
          CALL LUDBRB(IT+1,N,THE,0.,0D0,0D0,0D0)
        ENDIF
 
C...Reconstruct kinematics of branching: spacelike parton.
        DO 260 J=1,5
        K(N+1,J)=0
        P(N+1,J)=0.
  260   V(N+1,J)=0.
        K(N+1,1)=14
        K(N+1,2)=KFLB
        P(N+1,1)=P(IT,1)
        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)))
 
C...Define colour flow of branching.
        K(IS(JT),3)=N+1
        K(IT,3)=N+1
        ID1=IT
        IF((K(N+1,2).GT.0.AND.K(N+1,2).NE.21.AND.K(ID1,2).GT.0.AND.
     &  K(ID1,2).NE.21).OR.(K(N+1,2).LT.0.AND.K(ID1,2).EQ.21).OR.
     &  (K(N+1,2).EQ.21.AND.K(ID1,2).EQ.21.AND.RLU(0).GT.0.5).OR.
     &  (K(N+1,2).EQ.21.AND.K(ID1,2).LT.0)) ID1=IS(JT)
        ID2=IT+IS(JT)-ID1
        K(N+1,4)=K(N+1,4)+ID1
        K(N+1,5)=K(N+1,5)+ID2
        K(ID1,4)=K(ID1,4)+MSTU(5)*(N+1)
        K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
        K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
        K(ID2,5)=K(ID2,5)+MSTU(5)*(N+1)
        N=N+1
 
C...Boost to new CM-frame.
        CALL LUDBRB(NS+1,N,0.,0.,-DBLE((P(N,1)+P(IS(JR),1))/(P(N,4)+
     &  P(IS(JR),4))),0D0,-DBLE((P(N,3)+P(IS(JR),3))/(P(N,4)+
     &  P(IS(JR),4))))
        IR=N+(JT-1)*(IS(1)-N)
        CALL LUDBRB(NS+1,N,-ULANGL(P(IR,3),P(IR,1)),PARU(2)*RLU(0),
     &  0D0,0D0,0D0)
      ENDIF
 
C...Save quantities, loop back.
      IS(JT)=N
      Q2S(JT)=Q2B
      DQ2(JT)=Q2B
      IF(MSTP(62).GE.3) THE2(JT)=THE2T
      DSH=DSHZ
      IF(Q2B.GE.(0.5*PARP(62))**2) THEN
        KFLS(JT+2)=KFLS(JT)
        KFLS(JT)=KFLA
        XS(JT)=XA
        ZS(JT)=Z
        DO 270 KFL=-25,25
  270   XFS(JT,KFL)=XFA(KFL)
        TEVS(JT)=TEVB
      ELSE
        IF(JT.EQ.1) IPU1=N
        IF(JT.EQ.2) IPU2=N
      ENDIF
      IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
        CALL LUERRM(11,'(PYSSPA:) no more memory left in LUJETS')
        IF(MSTU(21).GE.1) N=NS
        IF(MSTU(21).GE.1) RETURN
      ENDIF
      IF(MAX(Q2S(1),Q2S(2)).GE.(0.5*PARP(62))**2.OR.N.LE.NS+1) GOTO 120
 
C...Boost hard scattering partons to frame of shower initiators.
      DO 280 J=1,3
  280 ROBO(J+2)=(P(NS+1,J)+P(NS+2,J))/(P(NS+1,4)+P(NS+2,4))
      K(N+2,1)=1
      DO 290 J=1,5
  290 P(N+2,J)=P(NS+1,J)
      ROBOT=ROBO(3)**2+ROBO(4)**2+ROBO(5)**2
      IF(ROBOT.GE.0.999999) THEN
        ROBOT=1.00001*SQRT(ROBOT)
        ROBO(3)=ROBO(3)/ROBOT
        ROBO(4)=ROBO(4)/ROBOT
        ROBO(5)=ROBO(5)/ROBOT
      ENDIF
      CALL LUDBRB(N+2,N+2,0.,0.,-DBLE(ROBO(3)),-DBLE(ROBO(4)),
     &-DBLE(ROBO(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))
      CALL LUDBRB(MINT(83)+5,NS,ROBO(1),ROBO(2),DBLE(ROBO(3)),
     &DBLE(ROBO(4)),DBLE(ROBO(5)))
 
C...Store user information. Reset Lambda value.
      K(IPU1,3)=MINT(83)+3
      K(IPU2,3)=MINT(83)+4
      DO 300 JT=1,2
      MINT(12+JT)=KFLS(JT)
  300 VINT(140+JT)=XS(JT)
      PARU(111)=ALAMS
 
      RETURN
      END
