hijing/src/pyhisspa.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE PYHISSPA(IPU1,IPU2)
0005
0006 C...Generates spacelike parton showers.
0007       IMPLICIT DOUBLE PRECISION(D)
0008       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0009       SAVE /LUJETS/
0010       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0011       SAVE /LUDAT1/
0012       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0013       SAVE /LUDAT2/
0014       COMMON/PYHISUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
0015       SAVE /PYHISUBS/
0016       COMMON/PYHIPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0017       SAVE /PYHIPARS/
0018       COMMON/PYHIINT1/MINT(400),VINT(400)
0019       SAVE /PYHIINT1/
0020       COMMON/PYHIINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
0021       SAVE /PYHIINT2/
0022       COMMON/PYHIINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
0023       SAVE /PYHIINT3/
0024       DIMENSION KFLS(4),IS(2),XS(2),ZS(2),Q2S(2),TEVS(2),ROBO(5),
0025      &XFS(2,-6:6),XFA(-6:6),XFB(-6:6),XFN(-6:6),WTAP(-6:6),WTSF(-6:6),
0026      &THE2(2),ALAM(2),DQ2(3),DPC(3),DPD(4),DPB(4)
0027
0028 C...Calculate maximum virtuality and check that evolution possible.
0029       IPUS1=IPU1
0030       IPUS2=IPU2
0031       ISUB=MINT(1)
0032       Q2E=VINT(52)
0033       IF(ISET(ISUB).EQ.1) THEN
0034         Q2E=Q2E/PARP(67)
0035       ELSEIF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) THEN
0036         Q2E=PMAS(23,1)**2
0037         IF(ISUB.EQ.8.OR.ISUB.EQ.76.OR.ISUB.EQ.77) Q2E=PMAS(24,1)**2
0038       ENDIF
0039       TMAX=LOG(PARP(67)*PARP(63)*Q2E/PARP(61)**2)
0040       IF(PARP(67)*Q2E.LT.MAX(PARP(62)**2,2.*PARP(61)**2).OR.
0041      &TMAX.LT.0.2) RETURN
0042
0043 C...Common constants and initial values. Save normal Lambda value.
0044       XE0=2.*PARP(65)/VINT(1)
0045       ALAMS=PARU(111)
0046       PARU(111)=PARP(61)
0047       NS=N
0048   100 N=NS
0049       DO 110 JT=1,2
0050       KFLS(JT)=MINT(14+JT)
0051       KFLS(JT+2)=KFLS(JT)
0052       XS(JT)=VINT(40+JT)
0053       ZS(JT)=1.
0054       Q2S(JT)=PARP(67)*Q2E
0055       TEVS(JT)=TMAX
0056       ALAM(JT)=PARP(61)
0057       THE2(JT)=100.
0058       DO 110 KFL=-6,6
0059   110 XFS(JT,KFL)=XSFX(JT,KFL)
0060       DSH=VINT(44)
0061       IF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) DSH=VINT(26)*VINT(2)
0062
0063 C...Pick up leg with highest virtuality.
0064   120 N=N+1
0065       JT=1
0066       IF(N.GT.NS+1.AND.Q2S(2).GT.Q2S(1)) JT=2
0067       KFLB=KFLS(JT)
0068       XB=XS(JT)
0069       DO 130 KFL=-6,6
0070   130 XFB(KFL)=XFS(JT,KFL)
0071       DSHR=2D0*SQRT(DSH)
0072       DSHZ=DSH/DBLE(ZS(JT))
0073       XE=MAX(XE0,XB*(1./(1.-PARP(66))-1.))
0074       IF(XB+XE.GE.0.999) THEN
0075         Q2B=0.
0076         GOTO 220
0077       ENDIF
0078
0079 C...Maximum Q2 without or with Q2 ordering. Effective Lambda and n_f.
0080       IF(MSTP(62).LE.1) THEN
0081         Q2B=0.5*(1./ZS(JT)+1.)*Q2S(JT)+0.5*(1./ZS(JT)-1.)*(Q2S(3-JT)-
0082      &  SNGL(DSH)+SQRT((SNGL(DSH)+Q2S(1)+Q2S(2))**2+8.*Q2S(1)*Q2S(2)*
0083      &  ZS(JT)/(1.-ZS(JT))))
0084         TEVB=LOG(PARP(63)*Q2B/ALAM(JT)**2)
0085       ELSE
0086         Q2B=Q2S(JT)
0087         TEVB=TEVS(JT)
0088       ENDIF
0089       ALSDUM=ULALPS(PARP(63)*Q2B)
0090       TEVB=TEVB+2.*LOG(ALAM(JT)/PARU(117))
0091       TEVBSV=TEVB
0092       ALAM(JT)=PARU(117)
0093       B0=(33.-2.*MSTU(118))/6.
0094
0095 C...Calculate Altarelli-Parisi and structure function weights.
0096       DO 140 KFL=-6,6
0097       WTAP(KFL)=0.
0098   140 WTSF(KFL)=0.
0099       IF(KFLB.EQ.21) THEN
0100         WTAPQ=16.*(1.-SQRT(XB+XE))/(3.*SQRT(XB))
0101         DO 150 KFL=-MSTP(54),MSTP(54)
0102         IF(KFL.EQ.0) WTAP(KFL)=6.*LOG((1.-XB)/XE)
0103   150   IF(KFL.NE.0) WTAP(KFL)=WTAPQ
0104       ELSE
0105         WTAP(0)=0.5*XB*(1./(XB+XE)-1.)
0106         WTAP(KFLB)=8.*LOG((1.-XB)*(XB+XE)/XE)/3.
0107       ENDIF
0108   160 WTSUM=0.
0109       IF(KFLB.NE.21) XFBO=XFB(KFLB)
0110       IF(KFLB.EQ.21) XFBO=XFB(0)
0111 C***************************************************************
0112 C**********ERROR HAS OCCURED HERE
0113       IF(XFBO.EQ.0.0) THEN
0114                 WRITE(MSTU(11),1000)
0115                 WRITE(MSTU(11),1001) KFLB,XFB(KFLB)
0116                 XFBO=0.00001
0117       ENDIF
0118 C****************************************************************
0119       DO 170 KFL=-MSTP(54),MSTP(54)
0120       WTSF(KFL)=XFB(KFL)/XFBO
0121   170 WTSUM=WTSUM+WTAP(KFL)*WTSF(KFL)
0122       WTSUM=MAX(0.0001,WTSUM)
0123
0124 C...Choose new t: fix alpha_s, alpha_s(Q2), alpha_s(k_T2).
0125   180 IF(MSTP(64).LE.0) THEN
0126         TEVB=TEVB+LOG(RLU(0))*PARU(2)/(PARU(111)*WTSUM)
0127       ELSEIF(MSTP(64).EQ.1) THEN
0128         TEVB=TEVB*EXP(MAX(-100.,LOG(RLU(0))*B0/WTSUM))
0129       ELSE
0130         TEVB=TEVB*EXP(MAX(-100.,LOG(RLU(0))*B0/(5.*WTSUM)))
0131       ENDIF
0132   190 Q2REF=ALAM(JT)**2*EXP(TEVB)
0133       Q2B=Q2REF/PARP(63)
0134
0135 C...Evolution ended or select flavour for branching parton.
0136       IF(Q2B.LT.PARP(62)**2) THEN
0137         Q2B=0.
0138       ELSE
0139         WTRAN=RLU(0)*WTSUM
0140         KFLA=-MSTP(54)-1
0141   200   KFLA=KFLA+1
0142         WTRAN=WTRAN-WTAP(KFLA)*WTSF(KFLA)
0143         IF(KFLA.LT.MSTP(54).AND.WTRAN.GT.0.) GOTO 200
0144         IF(KFLA.EQ.0) KFLA=21
0145
0146 C...Choose z value and corrective weight.
0147         IF(KFLB.EQ.21.AND.KFLA.EQ.21) THEN
0148           Z=1./(1.+((1.-XB)/XB)*(XE/(1.-XB))**RLU(0))
0149           WTZ=(1.-Z*(1.-Z))**2
0150         ELSEIF(KFLB.EQ.21) THEN
0151           Z=XB/(1.-RLU(0)*(1.-SQRT(XB+XE)))**2
0152           WTZ=0.5*(1.+(1.-Z)**2)*SQRT(Z)
0153         ELSEIF(KFLA.EQ.21) THEN
0154           Z=XB*(1.+RLU(0)*(1./(XB+XE)-1.))
0155           WTZ=1.-2.*Z*(1.-Z)
0156         ELSE
0157           Z=1.-(1.-XB)*(XE/((XB+XE)*(1.-XB)))**RLU(0)
0158           WTZ=0.5*(1.+Z**2)
0159         ENDIF
0160
0161 C...Option with resummation of soft gluon emission as effective z shift.
0162         IF(MSTP(65).GE.1) THEN
0163           RSOFT=6.
0164           IF(KFLB.NE.21) RSOFT=8./3.
0165           Z=Z*(TEVB/TEVS(JT))**(RSOFT*XE/((XB+XE)*B0))
0166           IF(Z.LE.XB) GOTO 180
0167         ENDIF
0168
0169 C...Option with alpha_s(k_T2)Q2): demand k_T2 > cutoff, reweight.
0170         IF(MSTP(64).GE.2) THEN
0171           IF((1.-Z)*Q2B.LT.PARP(62)**2) GOTO 180
0172           ALPRAT=TEVB/(TEVB+LOG(1.-Z))
0173           IF(ALPRAT.LT.5.*RLU(0)) GOTO 180
0174           IF(ALPRAT.GT.5.) WTZ=WTZ*ALPRAT/5.
0175         ENDIF
0176
0177 C...Option with angular ordering requirement.
0178         IF(MSTP(62).GE.3) THEN
0179           THE2T=(4.*Z**2*Q2B)/(VINT(2)*(1.-Z)*XB**2)
0180           IF(THE2T.GT.THE2(JT)) GOTO 180
0181         ENDIF
0182
0183 C...Weighting with new structure functions.
0184         CALL PYHISTFU(MINT(10+JT),XB,Q2REF,XFN,JT)
0185         IF(KFLB.NE.21) XFBN=XFN(KFLB)
0186         IF(KFLB.EQ.21) XFBN=XFN(0)
0187         IF(XFBN.LT.1E-20) THEN
0188           IF(KFLA.EQ.KFLB) THEN
0189             TEVB=TEVBSV
0190             WTAP(KFLB)=0.
0191             GOTO 160
0192           ELSEIF(TEVBSV-TEVB.GT.0.2) THEN
0193             TEVB=0.5*(TEVBSV+TEVB)
0194             GOTO 190
0195           ELSE
0196             XFBN=1E-10
0197           ENDIF
0198         ENDIF
0199         DO 210 KFL=-MSTP(54),MSTP(54)
0200   210   XFB(KFL)=XFN(KFL)
0201         XA=XB/Z
0202         CALL PYHISTFU(MINT(10+JT),XA,Q2REF,XFA,JT)
0203         IF(KFLA.NE.21) XFAN=XFA(KFLA)
0204         IF(KFLA.EQ.21) XFAN=XFA(0)
0205         IF(XFAN.LT.1E-20) GOTO 160
0206         IF(KFLA.NE.21) WTSFA=WTSF(KFLA)
0207         IF(KFLA.EQ.21) WTSFA=WTSF(0)
0208         IF(WTZ*XFAN/XFBN.LT.RLU(0)*WTSFA) GOTO 160
0209       ENDIF
0210
0211 C...Define two hard scatterers in their CM-frame.
0212   220 IF(N.EQ.NS+2) THEN
0213         DQ2(JT)=Q2B
0214         DPLCM=SQRT((DSH+DQ2(1)+DQ2(2))**2-4D0*DQ2(1)*DQ2(2))/DSHR
0215         DO 240 JR=1,2
0216         I=NS+JR
0217         IF(JR.EQ.1) IPO=IPUS1
0218         IF(JR.EQ.2) IPO=IPUS2
0219         DO 230 J=1,5
0220         K(I,J)=0
0221         P(I,J)=0.
0222   230   V(I,J)=0.
0223         K(I,1)=14
0224         K(I,2)=KFLS(JR+2)
0225         K(I,4)=IPO
0226         K(I,5)=IPO
0227         P(I,3)=DPLCM*(-1)**(JR+1)
0228         P(I,4)=(DSH+DQ2(3-JR)-DQ2(JR))/DSHR
0229         P(I,5)=-SQRT(SNGL(DQ2(JR)))
0230         K(IPO,1)=14
0231         K(IPO,3)=I
0232         K(IPO,4)=MOD(K(IPO,4),MSTU(5))+MSTU(5)*I
0233   240   K(IPO,5)=MOD(K(IPO,5),MSTU(5))+MSTU(5)*I
0234
0235 C...Find maximum allowed mass of timelike parton.
0236       ELSEIF(N.GT.NS+2) THEN
0237         JR=3-JT
0238         DQ2(3)=Q2B
0239         DPC(1)=P(IS(1),4)
0240         DPC(2)=P(IS(2),4)
0241         DPC(3)=0.5*(ABS(P(IS(1),3))+ABS(P(IS(2),3)))
0242         DPD(1)=DSH+DQ2(JR)+DQ2(JT)
0243         DPD(2)=DSHZ+DQ2(JR)+DQ2(3)
0244         DPD(3)=SQRT(DPD(1)**2-4D0*DQ2(JR)*DQ2(JT))
0245         DPD(4)=SQRT(DPD(2)**2-4D0*DQ2(JR)*DQ2(3))
0246         IKIN=0
0247         IF(Q2S(JR).GE.(0.5*PARP(62))**2.AND.DPD(1)-DPD(3).GE.
0248      &  1D-10*DPD(1)) IKIN=1
0249         IF(IKIN.EQ.0) DMSMA=(DQ2(JT)/DBLE(ZS(JT))-DQ2(3))*(DSH/
0250      &  (DSH+DQ2(JT))-DSH/(DSHZ+DQ2(3)))
0251         IF(IKIN.EQ.1) DMSMA=(DPD(1)*DPD(2)-DPD(3)*DPD(4))/(2.*
0252      &  DQ2(JR))-DQ2(JT)-DQ2(3)
0253
0254 C...Generate timelike parton shower (if required).
0255         IT=N
0256         DO 250 J=1,5
0257         K(IT,J)=0
0258         P(IT,J)=0.
0259   250   V(IT,J)=0.
0260         K(IT,1)=3
0261         K(IT,2)=21
0262         IF(KFLB.EQ.21.AND.KFLS(JT+2).NE.21) K(IT,2)=-KFLS(JT+2)
0263         IF(KFLB.NE.21.AND.KFLS(JT+2).EQ.21) K(IT,2)=KFLB
0264         P(IT,5)=ULMASS(K(IT,2))
0265         IF(SNGL(DMSMA).LE.P(IT,5)**2) GOTO 100
0266         IF(MSTP(63).GE.1) THEN
0267           P(IT,4)=(DSHZ-DSH-P(IT,5)**2)/DSHR
0268           P(IT,3)=SQRT(P(IT,4)**2-P(IT,5)**2)
0269           IF(MSTP(63).EQ.1) THEN
0270             Q2TIM=DMSMA
0271           ELSEIF(MSTP(63).EQ.2) THEN
0272             Q2TIM=MIN(SNGL(DMSMA),PARP(71)*Q2S(JT))
0273           ELSE
0274 C'''Here remains to introduce angular ordering in first branching.
0275             Q2TIM=DMSMA
0276           ENDIF
0277           CALL LUSHOW(IT,0,SQRT(Q2TIM))
0278           IF(N.GE.IT+1) P(IT,5)=P(IT+1,5)
0279         ENDIF
0280
0281 C...Reconstruct kinematics of branching: timelike parton shower.
0282         DMS=P(IT,5)**2
0283         IF(IKIN.EQ.0) DPT2=(DMSMA-DMS)*(DSHZ+DQ2(3))/(DSH+DQ2(JT))
0284         IF(IKIN.EQ.1) DPT2=(DMSMA-DMS)*(0.5*DPD(1)*DPD(2)+0.5*DPD(3)*
0285      &  DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)+DMS))/(4.*DSH*DPC(3)**2)
0286         IF(DPT2.LT.0.) GOTO 100
0287         DPB(1)=(0.5*DPD(2)-DPC(JR)*(DSHZ+DQ2(JR)-DQ2(JT)-DMS)/
0288      &  DSHR)/DPC(3)-DPC(3)
0289         P(IT,1)=SQRT(SNGL(DPT2))
0290         P(IT,3)=DPB(1)*(-1)**(JT+1)
0291         P(IT,4)=(DSHZ-DSH-DMS)/DSHR
0292         IF(N.GE.IT+1) THEN
0293           DPB(1)=SQRT(DPB(1)**2+DPT2)
0294           DPB(2)=SQRT(DPB(1)**2+DMS)
0295           DPB(3)=P(IT+1,3)
0296           DPB(4)=SQRT(DPB(3)**2+DMS)
0297           DBEZ=(DPB(4)*DPB(1)-DPB(3)*DPB(2))/(DPB(4)*DPB(2)-DPB(3)*
0298      &    DPB(1))
0299           CALL LUDBRB(IT+1,N,0.,0.,0D0,0D0,DBEZ)
0300           THE=ULANGL(P(IT,3),P(IT,1))
0301           CALL LUDBRB(IT+1,N,THE,0.,0D0,0D0,0D0)
0302         ENDIF
0303
0304 C...Reconstruct kinematics of branching: spacelike parton.
0305         DO 260 J=1,5
0306         K(N+1,J)=0
0307         P(N+1,J)=0.
0308   260   V(N+1,J)=0.
0309         K(N+1,1)=14
0310         K(N+1,2)=KFLB
0311         P(N+1,1)=P(IT,1)
0312         P(N+1,3)=P(IT,3)+P(IS(JT),3)
0313         P(N+1,4)=P(IT,4)+P(IS(JT),4)
0314         P(N+1,5)=-SQRT(SNGL(DQ2(3)))
0315
0316 C...Define colour flow of branching.
0317         K(IS(JT),3)=N+1
0318         K(IT,3)=N+1
0319         ID1=IT
0320         IF((K(N+1,2).GT.0.AND.K(N+1,2).NE.21.AND.K(ID1,2).GT.0.AND.
0321      &  K(ID1,2).NE.21).OR.(K(N+1,2).LT.0.AND.K(ID1,2).EQ.21).OR.
0322      &  (K(N+1,2).EQ.21.AND.K(ID1,2).EQ.21.AND.RLU(0).GT.0.5).OR.
0323      &  (K(N+1,2).EQ.21.AND.K(ID1,2).LT.0)) ID1=IS(JT)
0324         ID2=IT+IS(JT)-ID1
0325         K(N+1,4)=K(N+1,4)+ID1
0326         K(N+1,5)=K(N+1,5)+ID2
0327         K(ID1,4)=K(ID1,4)+MSTU(5)*(N+1)
0328         K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
0329         K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
0330         K(ID2,5)=K(ID2,5)+MSTU(5)*(N+1)
0331         N=N+1
0332
0333 C...Boost to new CM-frame.
0334         CALL LUDBRB(NS+1,N,0.,0.,-DBLE((P(N,1)+P(IS(JR),1))/(P(N,4)+
0335      &  P(IS(JR),4))),0D0,-DBLE((P(N,3)+P(IS(JR),3))/(P(N,4)+
0336      &  P(IS(JR),4))))
0337         IR=N+(JT-1)*(IS(1)-N)
0338         CALL LUDBRB(NS+1,N,-ULANGL(P(IR,3),P(IR,1)),PARU(2)*RLU(0),
0339      &  0D0,0D0,0D0)
0340       ENDIF
0341
0342 C...Save quantities, loop back.
0343       IS(JT)=N
0344       Q2S(JT)=Q2B
0345       DQ2(JT)=Q2B
0346       IF(MSTP(62).GE.3) THE2(JT)=THE2T
0347       DSH=DSHZ
0348       IF(Q2B.GE.(0.5*PARP(62))**2) THEN
0349         KFLS(JT+2)=KFLS(JT)
0350         KFLS(JT)=KFLA
0351         XS(JT)=XA
0352         ZS(JT)=Z
0353         DO 270 KFL=-6,6
0354   270   XFS(JT,KFL)=XFA(KFL)
0355         TEVS(JT)=TEVB
0356       ELSE
0357         IF(JT.EQ.1) IPU1=N
0358         IF(JT.EQ.2) IPU2=N
0359       ENDIF
0360       IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
0361         CALL LUERRM(11,'(PYHISSPA:) no more memory left in LUJETS')
0362         IF(MSTU(21).GE.1) N=NS
0363         IF(MSTU(21).GE.1) RETURN
0364       ENDIF
0365       IF(MAX(Q2S(1),Q2S(2)).GE.(0.5*PARP(62))**2.OR.N.LE.NS+1) GOTO 120
0366
0367 C...Boost hard scattering partons to frame of shower initiators.
0368       DO 280 J=1,3
0369   280 ROBO(J+2)=(P(NS+1,J)+P(NS+2,J))/(P(NS+1,4)+P(NS+2,4))
0370       DO 290 J=1,5
0371   290 P(N+2,J)=P(NS+1,J)
0372       ROBOT=ROBO(3)**2+ROBO(4)**2+ROBO(5)**2
0373       IF(ROBOT.GE.0.999999) THEN
0374         ROBOT=1.00001*SQRT(ROBOT)
0375         ROBO(3)=ROBO(3)/ROBOT
0376         ROBO(4)=ROBO(4)/ROBOT
0377         ROBO(5)=ROBO(5)/ROBOT
0378       ENDIF
0379       CALL LUDBRB(N+2,N+2,0.,0.,-DBLE(ROBO(3)),-DBLE(ROBO(4)),
0380      &-DBLE(ROBO(5)))
0381       ROBO(2)=ULANGL(P(N+2,1),P(N+2,2))
0382       ROBO(1)=ULANGL(P(N+2,3),SQRT(P(N+2,1)**2+P(N+2,2)**2))
0383       CALL LUDBRB(MINT(83)+5,NS,ROBO(1),ROBO(2),DBLE(ROBO(3)),
0384      &DBLE(ROBO(4)),DBLE(ROBO(5)))
0385
0386 C...Store user information. Reset Lambda value.
0387       K(IPU1,3)=MINT(83)+3
0388       K(IPU2,3)=MINT(83)+4
0389       DO 300 JT=1,2
0390       MINT(12+JT)=KFLS(JT)
0391   300 VINT(140+JT)=XS(JT)
0392       PARU(111)=ALAMS
0393  1000 FORMAT(5X,'structure function has a zero point here')
0394  1001 FORMAT(5X,'xf(x,i=',I5,')=',F10.5)
0395
0396       RETURN
0397       END