+++ /dev/null
-
-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),TEVCSV(2),TEVESV(2),
- &XFS(2,-25:25),XFA(-25:25),XFB(-25:25),XFN(-25:25),WTAPC(-25:25),
- &WTAPE(-25:25),WTSF(-25:25),THE2(2),ALAM(2),DQ2(3),DPC(3),DPD(4),
- &DPB(4),ROBO(5),MORE(2),KFBEAM(2),Q2MNCS(2),KCFI(2),NFIS(2),
- &THEFIS(2,2),ISFI(2)
- DATA IS/2*0/
-
-C...Read out basic information; set global Q^2 scale.
- IPUS1=IPU1
- IPUS2=IPU2
- ISUB=MINT(1)
- Q2MX=VINT(56)
- IF(ISET(ISUB).EQ.2) Q2MX=PARP(67)*VINT(56)
-
-C...Initialize QCD evolution and check phase space.
- Q2MNC=PARP(62)**2
- Q2MNCS(1)=Q2MNC
- IF(MSTP(66).EQ.1.AND.MINT(107).EQ.3)
- &Q2MNCS(1)=MAX(Q2MNC,VINT(283))
- Q2MNCS(2)=Q2MNC
- IF(MSTP(66).EQ.1.AND.MINT(108).EQ.3)
- &Q2MNCS(2)=MAX(Q2MNC,VINT(284))
- MCEV=0
- XEC0=2.*PARP(65)/VINT(1)
- ALAMS=PARU(112)
- PARU(112)=PARP(61)
- FQ2C=1.
- TCMX=0.
- IF(MINT(47).GE.2.AND.(MINT(47).NE.5.OR.MSTP(12).GE.1)) THEN
- MCEV=1
- IF(MSTP(64).EQ.1) FQ2C=PARP(63)
- IF(MSTP(64).EQ.2) FQ2C=PARP(64)
- TCMX=LOG(FQ2C*Q2MX/PARP(61)**2)
- IF(Q2MX.LT.MAX(Q2MNC,2.*PARP(61)**2).OR.TCMX.LT.0.2)
- & MCEV=0
- ENDIF
-
-C...Initialize QED evolution and check phase space.
- Q2MNE=PARP(68)**2
- MEEV=0
- XEE=1E-6
- SPME=PMAS(11,1)**2
- TEMX=0.
- FWTE=10.
- IF(MINT(45).EQ.3.OR.MINT(46).EQ.3) THEN
- MEEV=1
- TEMX=LOG(Q2MX/SPME)
- IF(Q2MX.LE.Q2MNE.OR.TEMX.LT.0.2) MEEV=0
- ENDIF
- IF(MCEV.EQ.0.AND.MEEV.EQ.0) RETURN
-
-C...Initial values: flavours, momenta, virtualities.
- NS=N
- 100 N=NS
- DO 120 JT=1,2
- MORE(JT)=1
- KFBEAM(JT)=MINT(10+JT)
- IF(MINT(18+JT).EQ.1)KFBEAM(JT)=22
- KFLS(JT)=MINT(14+JT)
- KFLS(JT+2)=KFLS(JT)
- XS(JT)=VINT(40+JT)
- IF(MINT(18+JT).EQ.1) XS(JT)=VINT(40+JT)/VINT(154+JT)
- ZS(JT)=1.
- Q2S(JT)=Q2MX
- TEVCSV(JT)=TCMX
- ALAM(JT)=PARP(61)
- THE2(JT)=100.
- TEVESV(JT)=TEMX
- DO 110 KFL=-25,25
- XFS(JT,KFL)=XSFX(JT,KFL)
- 110 CONTINUE
- 120 CONTINUE
- DSH=VINT(44)
- IF(ISET(ISUB).GE.3.AND.ISET(ISUB).LE.5) DSH=VINT(26)*VINT(2)
-
-C...Find if interference with final state partons.
- MFIS=0
- IF(MSTP(67).GE.1.AND.MSTP(67).LE.3) MFIS=MSTP(67)
- IF(MFIS.NE.0) THEN
- DO 140 I=1,2
- KCFI(I)=0
- KCA=LUCOMP(IABS(KFLS(I)))
- IF(KCA.NE.0) KCFI(I)=KCHG(KCA,2)*ISIGN(1,KFLS(I))
- NFIS(I)=0
- IF(KCFI(I).NE.0) THEN
- IF(I.EQ.1) IPFS=IPUS1
- IF(I.EQ.2) IPFS=IPUS2
- DO 130 J=1,2
- ICSI=MOD(K(IPFS,3+J),MSTU(5))
- IF(ICSI.GT.0.AND.ICSI.NE.IPUS1.AND.ICSI.NE.IPUS2.AND.
- & (KCFI(I).EQ.(-1)**(J+1).OR.KCFI(I).EQ.2)) THEN
- NFIS(I)=NFIS(I)+1
- THEFIS(I,NFIS(I))=ULANGL(P(ICSI,3),SQRT(P(ICSI,1)**2+
- & P(ICSI,2)**2))
- IF(I.EQ.2) THEFIS(I,NFIS(I))=PARU(1)-THEFIS(I,NFIS(I))
- ENDIF
- 130 CONTINUE
- ENDIF
- 140 CONTINUE
- IF(NFIS(1)+NFIS(2).EQ.0) MFIS=0
- ENDIF
-
-C...Pick up leg with highest virtuality.
- 150 N=N+1
- JT=1
- IF(N.GT.NS+1.AND.Q2S(2).GT.Q2S(1)) JT=2
- IF(MORE(JT).EQ.0) JT=3-JT
- KFLB=KFLS(JT)
- XB=XS(JT)
- DO 160 KFL=-25,25
- XFB(KFL)=XFS(JT,KFL)
- 160 CONTINUE
- DSHR=2D0*SQRT(DSH)
- DSHZ=DSH/DBLE(ZS(JT))
-
-C...Check if allowed to branch.
- MCEV=0
- IF(IABS(KFLB).LE.10.OR.KFLB.EQ.21) THEN
- MCEV=1
- XEC=MAX(XEC0,XB*(1./(1.-PARP(66))-1.))
- IF(XB.GE.1.-2.*XEC) MCEV=0
- ENDIF
- MEEV=0
- IF(MINT(44+JT).EQ.3) THEN
- MEEV=1
- IF(XB.GE.1.-2.*XEE) MEEV=0
- IF((IABS(KFLB).LE.10.OR.KFLB.EQ.21).AND.XB.GE.1.-2.*XEC) MEEV=0
-C***Currently kill QED shower for resolved photoproduction.
- IF(MINT(18+JT).EQ.1) MEEV=0
-C***Currently kill shower for W inside electron.
- IF(IABS(KFLB).EQ.24) THEN
- MCEV=0
- MEEV=0
- ENDIF
- ENDIF
- IF(MCEV.EQ.0.AND.MEEV.EQ.0) THEN
- Q2B=0.
- GOTO 250
- ENDIF
-
-C...Maximum Q2 with or without Q2 ordering. Effective Lambda and n_f.
- Q2B=Q2S(JT)
- TEVCB=TEVCSV(JT)
- TEVEB=TEVESV(JT)
- IF(MSTP(62).LE.1) THEN
- IF(ZS(JT).GT.0.99999) THEN
- Q2B=Q2S(JT)
- ELSE
- 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))))
- ENDIF
- IF(MCEV.EQ.1) TEVCB=LOG(FQ2C*Q2B/ALAM(JT)**2)
- IF(MEEV.EQ.1) TEVEB=LOG(Q2B/SPME)
- ENDIF
- IF(MCEV.EQ.1) THEN
- ALSDUM=ULALPS(FQ2C*Q2B)
- TEVCB=TEVCB+2.*LOG(ALAM(JT)/PARU(117))
- ALAM(JT)=PARU(117)
- B0=(33.-2.*MSTU(118))/6.
- ENDIF
- TEVCBS=TEVCB
- TEVEBS=TEVEB
-
-C...Select side for interference with final state partons.
- IF(MFIS.GE.1.AND.N.LE.NS+2) THEN
- IFI=N-NS
- ISFI(IFI)=0
- IF(IABS(KCFI(IFI)).EQ.1.AND.NFIS(IFI).EQ.1) THEN
- ISFI(IFI)=1
- ELSEIF(KCFI(IFI).EQ.2.AND.NFIS(IFI).EQ.1) THEN
- IF(RLU(0).GT.0.5) ISFI(IFI)=1
- ELSEIF(KCFI(IFI).EQ.2.AND.NFIS(IFI).EQ.2) THEN
- ISFI(IFI)=1
- IF(RLU(0).GT.0.5) ISFI(IFI)=2
- ENDIF
- ENDIF
-
-C...Calculate Altarelli-Parisi weights.
- DO 170 KFL=-25,25
- WTAPC(KFL)=0.
- WTAPE(KFL)=0.
- WTSF(KFL)=0.
- 170 CONTINUE
-C...q -> q, g -> q.
- IF(IABS(KFLB).LE.10) THEN
- WTAPC(KFLB)=(8./3.)*LOG((1.-XEC-XB)*(XB+XEC)/(XEC*(1.-XEC)))
- WTAPC(21)=0.5*(XB/(XB+XEC)-XB/(1.-XEC))
-C...f -> f, gamma -> f.
- ELSEIF(IABS(KFLB).LE.20) THEN
- WTAPF1=LOG((1.-XEE-XB)*(XB+XEE)/(XEE*(1.-XEE)))
- WTAPF2=LOG((1.-XEE-XB)*(1.-XEE)/(XEE*(XB+XEE)))
- WTAPE(KFLB)=2.*(WTAPF1+WTAPF2)
- IF(MSTP(12).GE.1) WTAPE(22)=XB/(XB+XEE)-XB/(1.-XEE)
-C...f -> g, g -> g.
- ELSEIF(KFLB.EQ.21) THEN
- WTAPQ=(16./3.)*(SQRT((1.-XEC)/XB)-SQRT((XB+XEC)/XB))
- DO 180 KFL=1,MSTP(58)
- WTAPC(KFL)=WTAPQ
- WTAPC(-KFL)=WTAPQ
- 180 CONTINUE
- WTAPC(21)=6.*LOG((1.-XEC-XB)/XEC)
-C...f -> gamma, W+, W-.
- ELSEIF(KFLB.EQ.22) THEN
- WTAPF=LOG((1.-XEE-XB)*(1.-XEE)/(XEE*(XB+XEE)))/XB
- WTAPE(11)=WTAPF
- WTAPE(-11)=WTAPF
- ELSEIF(KFLB.EQ.24) THEN
- WTAPE(-11)=1./(4.*PARU(102))*LOG((1.-XEE-XB)*(1.-XEE)/
- & (XEE*(XB+XEE)))/XB
- ELSEIF(KFLB.EQ.-24) THEN
- WTAPE(11)=1./(4.*PARU(102))*LOG((1.-XEE-XB)*(1.-XEE)/
- & (XEE*(XB+XEE)))/XB
- ENDIF
-
-C...Calculate structure function weights and sum.
- NTRY=0
- 190 NTRY=NTRY+1
- IF(NTRY.GT.500) THEN
- MINT(51)=1
- RETURN
- ENDIF
- WTSUMC=0.
- WTSUME=0.
- XFBO=MAX(1E-10,XFB(KFLB))
- DO 200 KFL=-25,25
- WTSF(KFL)=XFB(KFL)/XFBO
- WTSUMC=WTSUMC+WTAPC(KFL)*WTSF(KFL)
- WTSUME=WTSUME+WTAPE(KFL)*WTSF(KFL)
- 200 CONTINUE
- WTSUMC=MAX(0.0001,WTSUMC)
- WTSUME=MAX(0.0001/FWTE,WTSUME)
-
-C...Choose new t: fix alpha_s, alpha_s(Q^2), alpha_s(k_T^2).
- NTRY2=0
- 210 NTRY2=NTRY2+1
- IF(NTRY2.GT.500) THEN
- MINT(51)=1
- RETURN
- ENDIF
- IF(MCEV.EQ.1) THEN
- IF(MSTP(64).LE.0) THEN
- TEVCB=TEVCB+LOG(RLU(0))*PARU(2)/(PARU(111)*WTSUMC)
- ELSEIF(MSTP(64).EQ.1) THEN
- TEVCB=TEVCB*EXP(MAX(-50.,LOG(RLU(0))*B0/WTSUMC))
- ELSE
- TEVCB=TEVCB*EXP(MAX(-50.,LOG(RLU(0))*B0/(5.*WTSUMC)))
- ENDIF
- ENDIF
- IF(MEEV.EQ.1) THEN
- TEVEB=TEVEB*EXP(MAX(-50.,LOG(RLU(0))*PARU(2)/
- & (PARU(101)*FWTE*WTSUME*TEMX)))
- ENDIF
-
-C...Translate t into Q2 scale; choose between QCD and QED evolution.
- 220 IF(MCEV.EQ.1) Q2CB=ALAM(JT)**2*EXP(MAX(-50.,TEVCB))/FQ2C
- IF(MEEV.EQ.1) Q2EB=SPME*EXP(MAX(-50.,TEVEB))
- MCE=0
- IF(MCEV.EQ.0.AND.MEEV.EQ.0) THEN
- ELSEIF(MCEV.EQ.1.AND.MEEV.EQ.0) THEN
- IF(Q2CB.GT.Q2MNCS(JT)) MCE=1
- ELSEIF(MCEV.EQ.0.AND.MEEV.EQ.1) THEN
- IF(Q2EB.GT.Q2MNE) MCE=2
- ELSEIF(Q2MNCS(JT).GT.Q2MNE) THEN
- MCE=1
- IF(Q2EB.GT.Q2CB.OR.Q2CB.LE.Q2MNCS(JT)) MCE=2
- IF(MCE.EQ.2.AND.Q2EB.LE.Q2MNE) MCE=0
- ELSE
- MCE=2
- IF(Q2CB.GT.Q2EB.OR.Q2EB.LE.Q2MNE) MCE=1
- IF(MCE.EQ.1.AND.Q2CB.LE.Q2MNCS(JT)) MCE=0
- ENDIF
-
-C...Evolution possibly ended. Update t values.
- IF(MCE.EQ.0) THEN
- Q2B=0.
- GOTO 250
- ELSEIF(MCE.EQ.1) THEN
- Q2B=Q2CB
- Q2REF=FQ2C*Q2B
- IF(MEEV.EQ.1) TEVEB=LOG(Q2B/SPME)
- ELSE
- Q2B=Q2EB
- Q2REF=Q2B
- IF(MCEV.EQ.1) TEVCB=LOG(FQ2C*Q2B/ALAM(JT)**2)
- ENDIF
-
-C...Select flavour for branching parton.
- IF(MCE.EQ.1) WTRAN=RLU(0)*WTSUMC
- IF(MCE.EQ.2) WTRAN=RLU(0)*WTSUME
- KFLA=-25
- 230 KFLA=KFLA+1
- IF(MCE.EQ.1) WTRAN=WTRAN-WTAPC(KFLA)*WTSF(KFLA)
- IF(MCE.EQ.2) WTRAN=WTRAN-WTAPE(KFLA)*WTSF(KFLA)
- IF(KFLA.LE.24.AND.WTRAN.GT.0.) GOTO 230
- IF(KFLA.EQ.25) THEN
- Q2B=0.
- GOTO 250
- ENDIF
-
-C...Choose z value and corrective weight.
- WTZ=0.
-C...q -> q + g.
- IF(IABS(KFLA).LE.10.AND.IABS(KFLB).LE.10) THEN
- Z=1.-((1.-XB-XEC)/(1.-XEC))*
- & (XEC*(1.-XEC)/((XB+XEC)*(1.-XB-XEC)))**RLU(0)
- WTZ=0.5*(1.+Z**2)
-C...q -> g + q.
- ELSEIF(IABS(KFLA).LE.10.AND.KFLB.EQ.21) THEN
- Z=XB/(SQRT(XB+XEC)+RLU(0)*(SQRT(1.-XEC)-SQRT(XB+XEC)))**2
- WTZ=0.5*(1.+(1.-Z)**2)*SQRT(Z)
-C...f -> f + gamma.
- ELSEIF(IABS(KFLA).LE.20.AND.IABS(KFLB).LE.20) THEN
- IF(WTAPF1.GT.RLU(0)*(WTAPF1+WTAPF2)) THEN
- Z=1.-((1.-XB-XEE)/(1.-XEE))*
- & (XEE*(1.-XEE)/((XB+XEE)*(1.-XB-XEE)))**RLU(0)
- ELSE
- Z=XB+XB*(XEE/(1.-XEE))*
- & ((1.-XB-XEE)*(1.-XEE)/(XEE*(XB+XEE)))**RLU(0)
- ENDIF
- WTZ=0.5*(1.+Z**2)*(Z-XB)/(1.-XB)
-C...f -> gamma + f.
- ELSEIF(IABS(KFLA).LE.20.AND.KFLB.EQ.22) THEN
- Z=XB+XB*(XEE/(1.-XEE))*
- & ((1.-XB-XEE)*(1.-XEE)/(XEE*(XB+XEE)))**RLU(0)
- WTZ=0.5*(1.+(1.-Z)**2)*XB*(Z-XB)/Z
-C...f -> W+- + f'.
- ELSEIF(IABS(KFLA).LE.20.AND.IABS(KFLB).EQ.24) THEN
- Z=XB+XB*(XEE/(1.-XEE))*
- & ((1.-XB-XEE)*(1.-XEE)/(XEE*(XB+XEE)))**RLU(0)
- WTZ=0.5*(1.+(1.-Z)**2)*(XB*(Z-XB)/Z)*(Q2B/(Q2B+PMAS(24,1)**2))
-C...g -> q + q~.
- ELSEIF(KFLA.EQ.21.AND.IABS(KFLB).LE.10) THEN
- Z=XB/(1.-XEC)+RLU(0)*(XB/(XB+XEC)-XB/(1.-XEC))
- WTZ=1.-2.*Z*(1.-Z)
-C...g -> g + g.
- ELSEIF(KFLA.EQ.21.AND.KFLB.EQ.21) THEN
- Z=1./(1.+((1.-XEC-XB)/XB)*(XEC/(1.-XEC-XB))**RLU(0))
- WTZ=(1.-Z*(1.-Z))**2
-C...gamma -> f + f~.
- ELSEIF(KFLA.EQ.22.AND.IABS(KFLB).LE.20) THEN
- Z=XB/(1.-XEE)+RLU(0)*(XB/(XB+XEE)-XB/(1.-XEE))
- WTZ=1.-2.*Z*(1.-Z)
- ENDIF
- IF(MCE.EQ.2) WTZ=(WTZ/FWTE)*(TEVEB/TEMX)
-
-C...Option with resummation of soft gluon emission as effective z shift.
- IF(MCE.EQ.1) THEN
- IF(MSTP(65).GE.1) THEN
- RSOFT=6.
- IF(KFLB.NE.21) RSOFT=8./3.
- Z=Z*(TEVCB/TEVCSV(JT))**(RSOFT*XEC/((XB+XEC)*B0))
- IF(Z.LE.XB) GOTO 210
- 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.Q2MNCS(JT)) GOTO 210
- ALPRAT=TEVCB/(TEVCB+LOG(1.-Z))
- IF(ALPRAT.LT.5.*RLU(0)) GOTO 210
- IF(ALPRAT.GT.5.) WTZ=WTZ*ALPRAT/5.
- ENDIF
-
-C...Impose angular constraint in first branching from interference
-C...with final state partons.
- IF(MFIS.GE.1.AND.N.LE.NS+2.AND.NTRY2.LT.200) THEN
- THE2D=(4.*Q2B)/(DSH*(1.-Z))
- IF(N.EQ.NS+1.AND.ISFI(1).GE.1) THEN
- IF(THE2D.GT.THEFIS(1,ISFI(1))**2) GOTO 210
- ELSEIF(N.EQ.NS+2.AND.ISFI(2).GE.1) THEN
- IF(THE2D.GT.THEFIS(2,ISFI(2))**2) GOTO 210
- ENDIF
- ENDIF
-
-C...Option with angular ordering requirement.
- IF(MSTP(62).GE.3.AND.NTRY2.LT.200) THEN
- THE2T=(4.*Z**2*Q2B)/(VINT(2)*(1.-Z)*XB**2)
- IF(THE2T.GT.THE2(JT)) GOTO 210
- ENDIF
- ENDIF
-
-C...Weighting with new structure functions.
- MINT(105)=MINT(102+JT)
- MINT(109)=MINT(106+JT)
- IF(MSTP(57).LE.1) THEN
- CALL PYSTFU(KFBEAM(JT),XB,Q2REF,XFN)
- ELSE
- CALL PYSTFL(KFBEAM(JT),XB,Q2REF,XFN)
- ENDIF
- XFBN=XFN(KFLB)
- IF(XFBN.LT.1E-20) THEN
- IF(KFLA.EQ.KFLB) THEN
- TEVCB=TEVCBS
- TEVEB=TEVEBS
- WTAPC(KFLB)=0.
- WTAPE(KFLB)=0.
- GOTO 190
- ELSEIF(MCE.EQ.1.AND.TEVCBS-TEVCB.GT.0.2) THEN
- TEVCB=0.5*(TEVCBS+TEVCB)
- GOTO 220
- ELSEIF(MCE.EQ.2.AND.TEVEBS-TEVEB.GT.0.2) THEN
- TEVEB=0.5*(TEVEBS+TEVEB)
- GOTO 220
- ELSE
- XFBN=1E-10
- XFN(KFLB)=XFBN
- ENDIF
- ENDIF
- DO 240 KFL=-25,25
- XFB(KFL)=XFN(KFL)
- 240 CONTINUE
- XA=XB/Z
- IF(MSTP(57).LE.1) THEN
- CALL PYSTFU(KFBEAM(JT),XA,Q2REF,XFA)
- ELSE
- CALL PYSTFL(KFBEAM(JT),XA,Q2REF,XFA)
- ENDIF
- XFAN=XFA(KFLA)
- IF(XFAN.LT.1E-20) GOTO 190
- WTSFA=WTSF(KFLA)
- IF(WTZ*XFAN/XFBN.LT.RLU(0)*WTSFA) GOTO 190
-
-C...Define two hard scatterers in their CM-frame.
- 250 IF(N.EQ.NS+2) THEN
- DQ2(JT)=Q2B
- DPLCM=SQRT((DSH+DQ2(1)+DQ2(2))**2-4D0*DQ2(1)*DQ2(2))/DSHR
- DO 270 JR=1,2
- I=NS+JR
- IF(JR.EQ.1) IPO=IPUS1
- IF(JR.EQ.2) IPO=IPUS2
- DO 260 J=1,5
- K(I,J)=0
- P(I,J)=0.
- V(I,J)=0.
- 260 CONTINUE
- 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
- K(IPO,5)=MOD(K(IPO,5),MSTU(5))+MSTU(5)*I
- 270 CONTINUE
-
-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.25*Q2MNC.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 280 J=1,5
- K(IT,J)=0
- P(IT,J)=0.
- V(IT,J)=0.
- 280 CONTINUE
- K(IT,1)=3
-C...f -> f + g (gamma).
- IF(IABS(KFLB).LE.20.AND.IABS(KFLS(JT+2)).LE.20) THEN
- K(IT,2)=21
- IF(IABS(KFLB).GE.11) K(IT,2)=22
-C...f -> g (gamma, W+-) + f.
- ELSEIF(IABS(KFLB).LE.20.AND.IABS(KFLS(JT+2)).GT.20) THEN
- K(IT,2)=KFLB
- IF(KFLS(JT+2).EQ.24) THEN
- K(IT,2)=-12
- ELSEIF(KFLS(JT+2).EQ.-24) THEN
- K(IT,2)=12
- ENDIF
-C...g (gamma) -> f + f~, g + g.
- ELSE
- K(IT,2)=-KFLS(JT+2)
- IF(KFLS(JT+2).GT.20) K(IT,2)=KFLS(JT+2)
- ENDIF
- P(IT,5)=ULMASS(K(IT,2))
- IF(SNGL(DMSMA).LE.P(IT,5)**2) GOTO 100
- IF(MSTP(63).GE.1.AND.MCE.EQ.1) THEN
- MSTJ48=MSTJ(48)
- PARJ85=PARJ(85)
- 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
- Q2TIM=DMSMA
- MSTJ(48)=1
- IF(IKIN.EQ.0) DPT2=DMSMA*(DSHZ+DQ2(3))/(DSH+DQ2(JT))
- IF(IKIN.EQ.1) DPT2=DMSMA*(0.5*DPD(1)*DPD(2)+0.5*DPD(3)*
- & DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)))/(4.*DSH*DPC(3)**2)
- PARJ(85)=SQRT(MAX(0.,SNGL(DPT2)))*
- & (1./P(IT,4)+1./P(IS(JT),4))
- ENDIF
- CALL LUSHOW(IT,0,SQRT(Q2TIM))
- MSTJ(48)=MSTJ48
- PARJ(85)=PARJ85
- 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)=SQRT(DPT2+DPB(1)**2+DMS)
- 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 290 J=1,5
- K(N+1,J)=0
- P(N+1,J)=0.
- V(N+1,J)=0.
- 290 CONTINUE
- 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
- IM1=N+1
- IM2=N+1
-C...f -> f + gamma (Z, W).
- IF(IABS(K(IT,2)).GE.22) THEN
- K(IT,1)=1
- ID1=IS(JT)
- ID2=IS(JT)
-C...f -> gamma (Z, W) + f.
- ELSEIF(IABS(K(IS(JT),2)).GE.22) THEN
- ID1=IT
- ID2=IT
-C...gamma -> q + q~, g + g.
- ELSEIF(K(N+1,2).EQ.22) THEN
- ID1=IS(JT)
- ID2=IT
- IM1=ID2
- IM2=ID1
-C...q -> q + g.
- ELSEIF(K(N+1,2).GT.0.AND.K(N+1,2).NE.21.AND.K(IT,2).EQ.21) THEN
- ID1=IT
- ID2=IS(JT)
-C...q -> g + q.
- ELSEIF(K(N+1,2).GT.0.AND.K(N+1,2).NE.21) THEN
- ID1=IS(JT)
- ID2=IT
-C...q~ -> q~ + g.
- ELSEIF(K(N+1,2).LT.0.AND.K(IT,2).EQ.21) THEN
- ID1=IS(JT)
- ID2=IT
-C...q~ -> g + q~.
- ELSEIF(K(N+1,2).LT.0) THEN
- ID1=IT
- ID2=IS(JT)
-C...g -> g + g; g -> q + q~.
- ELSEIF((K(IT,2).EQ.21.AND.RLU(0).GT.0.5).OR.K(IT,2).LT.0) THEN
- ID1=IS(JT)
- ID2=IT
- ELSE
- ID1=IT
- ID2=IS(JT)
- ENDIF
- IF(IM1.EQ.N+1) K(IM1,4)=K(IM1,4)+ID1
- IF(IM2.EQ.N+1) K(IM2,5)=K(IM2,5)+ID2
- K(ID1,4)=K(ID1,4)+MSTU(5)*IM1
- K(ID2,5)=K(ID2,5)+MSTU(5)*IM2
- IF(ID1.NE.ID2) THEN
- K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
- K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
- ENDIF
- N=N+1
-
-C...Boost to new CM-frame.
- DBSVX=DBLE((P(N,1)+P(IS(JR),1))/(P(N,4)+P(IS(JR),4)))
- DBSVZ=DBLE((P(N,3)+P(IS(JR),3))/(P(N,4)+P(IS(JR),4)))
- IF(DBSVX**2+DBSVZ**2.GE.1D0) GOTO 100
- CALL LUDBRB(NS+1,N,0.,0.,-DBSVX,0D0,-DBSVZ)
- 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...Update kinematics variables.
- IS(JT)=N
- DQ2(JT)=Q2B
- IF(MSTP(62).GE.3) THE2(JT)=THE2T
- DSH=DSHZ
-
-C...Save quantities; loop back.
- Q2S(JT)=Q2B
- IF((MCEV.EQ.1.AND.Q2B.GE.0.25*Q2MNC).OR.
- &(MEEV.EQ.1.AND.Q2B.GE.Q2MNE)) THEN
- KFLS(JT+2)=KFLS(JT)
- KFLS(JT)=KFLA
- XS(JT)=XA
- ZS(JT)=Z
- DO 300 KFL=-25,25
- XFS(JT,KFL)=XFA(KFL)
- 300 CONTINUE
- TEVCSV(JT)=TEVCB
- TEVESV(JT)=TEVEB
- ELSE
- MORE(JT)=0
- 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(MORE(1).EQ.1.OR.MORE(2).EQ.1) GOTO 150
-
-C...Boost hard scattering partons to frame of shower initiators.
- DO 310 J=1,3
- ROBO(J+2)=(P(NS+1,J)+P(NS+2,J))/(P(NS+1,4)+P(NS+2,4))
- 310 CONTINUE
- K(N+2,1)=1
- DO 320 J=1,5
- P(N+2,J)=P(NS+1,J)
- 320 CONTINUE
- 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 330 JT=1,2
- MINT(12+JT)=KFLS(JT)
- VINT(140+JT)=XS(JT)
- IF(MINT(18+JT).EQ.1) VINT(140+JT)=VINT(154+JT)*XS(JT)
- 330 CONTINUE
- PARU(112)=ALAMS
-
- RETURN
- END
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff