pythiaeRHIC/pythia/pyx3jt.f

0001
0002 C*********************************************************************
0003
0004 C...PYX3JT
0005 C...Selects the kinematical variables of three-jet events.
0006
0007       SUBROUTINE PYX3JT(NJET,CUT,KFL,ECM,X1,X2)
0008
0009 C...Double precision and integer declarations.
0010       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0011       IMPLICIT INTEGER(I-N)
0012       INTEGER PYK,PYCHGE,PYCOMP
0013 C...Commonblocks.
0014       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0015       SAVE /PYDAT1/
0016 C...Local array.
0017       DIMENSION ZHUP(5,12)
0018
0019 C...Coefficients of Zhu second order parametrization.
0020       DATA ((ZHUP(IC1,IC2),IC2=1,12),IC1=1,5)/
0021      &18.29D0,  89.56D0,  4.541D0,  -52.09D0, -109.8D0,  24.90D0,
0022      &11.63D0,  3.683D0,  17.50D0,0.002440D0, -1.362D0,-0.3537D0,
0023      &11.42D0,  6.299D0, -22.55D0,  -8.915D0,  59.25D0, -5.855D0,
0024      &-32.85D0, -1.054D0, -16.90D0,0.006489D0,-0.8156D0,0.01095D0,
0025      &7.847D0, -3.964D0, -35.83D0,   1.178D0,  29.39D0, 0.2806D0,
0026      &47.82D0, -12.36D0, -56.72D0, 0.04054D0,-0.4365D0, 0.6062D0,
0027      &5.441D0, -56.89D0, -50.27D0,   15.13D0,  114.3D0, -18.19D0,
0028      &97.05D0, -1.890D0, -139.9D0, 0.08153D0,-0.4984D0, 0.9439D0,
0029      &-17.65D0,  51.44D0, -58.32D0,   70.95D0, -255.7D0, -78.99D0,
0030      &476.9D0,  29.65D0, -239.3D0,  0.4745D0, -1.174D0,  6.081D0/
0031
0032 C...Dilogarithm of x for x<0.5 (x>0.5 obtained by analytic trick).
0033       DILOG(X)=X+X**2/4D0+X**3/9D0+X**4/16D0+X**5/25D0+X**6/36D0+
0034      &X**7/49D0
0035
0036 C...Event type. Mass effect factors and other common constants.
0037       MSTJ(120)=2
0038       MSTJ(121)=0
0039       PMQ=PYMASS(KFL)
0040       QME=(2D0*PMQ/ECM)**2
0041       IF(MSTJ(109).NE.1) THEN
0042         CUTL=LOG(CUT)
0043         CUTD=LOG(1D0/CUT-2D0)
0044         IF(MSTJ(109).EQ.0) THEN
0045           CF=4D0/3D0
0046           CN=3D0
0047           TR=2D0
0048           WTMX=MIN(20D0,37D0-6D0*CUTD)
0049           IF(MSTJ(110).EQ.2) WTMX=2D0*(7.5D0+80D0*CUT)
0050         ELSE
0051           CF=1D0
0052           CN=0D0
0053           TR=12D0
0054           WTMX=0D0
0055         ENDIF
0056
0057 C...Alpha_strong and effects of optimized Q^2 scale. Maximum weight.
0058         ALS2PI=PARU(118)/PARU(2)
0059         WTOPT=0D0
0060         IF(MSTJ(111).EQ.1) WTOPT=(33D0-2D0*MSTU(112))/6D0*
0061      &  LOG(PARJ(169))*ALS2PI
0062         WTMAX=MAX(0D0,1D0+WTOPT+ALS2PI*WTMX)
0063
0064 C...Choose three-jet events in allowed region.
0065   100   NJET=3
0066   110   Y13L=CUTL+CUTD*PYR(0)
0067         Y23L=CUTL+CUTD*PYR(0)
0068         Y13=EXP(Y13L)
0069         Y23=EXP(Y23L)
0070         Y12=1D0-Y13-Y23
0071         IF(Y12.LE.CUT) GOTO 110
0072         IF(Y13**2+Y23**2+2D0*Y12.LE.2D0*PYR(0)) GOTO 110
0073
0074 C...Second order corrections.
0075         IF(MSTJ(101).EQ.2.AND.MSTJ(110).LE.1) THEN
0076           Y12L=LOG(Y12)
0077           Y13M=LOG(1D0-Y13)
0078           Y23M=LOG(1D0-Y23)
0079           Y12M=LOG(1D0-Y12)
0080           IF(Y13.LE.0.5D0) Y13I=DILOG(Y13)
0081           IF(Y13.GE.0.5D0) Y13I=1.644934D0-Y13L*Y13M-DILOG(1D0-Y13)
0082           IF(Y23.LE.0.5D0) Y23I=DILOG(Y23)
0083           IF(Y23.GE.0.5D0) Y23I=1.644934D0-Y23L*Y23M-DILOG(1D0-Y23)
0084           IF(Y12.LE.0.5D0) Y12I=DILOG(Y12)
0085           IF(Y12.GE.0.5D0) Y12I=1.644934D0-Y12L*Y12M-DILOG(1D0-Y12)
0086           WT1=(Y13**2+Y23**2+2D0*Y12)/(Y13*Y23)
0087           WT2=CF*(-2D0*(CUTL-Y12L)**2-3D0*CUTL-1D0+3.289868D0+
0088      &    2D0*(2D0*CUTL-Y12L)*CUT/Y12)+
0089      &    CN*((CUTL-Y12L)**2-(CUTL-Y13L)**2-(CUTL-Y23L)**2-
0090      &    11D0*CUTL/6D0+67D0/18D0+1.644934D0-(2D0*CUTL-Y12L)*CUT/Y12+
0091      &    (2D0*CUTL-Y13L)*CUT/Y13+(2D0*CUTL-Y23L)*CUT/Y23)+
0092      &    TR*(2D0*CUTL/3D0-10D0/9D0)+
0093      &    CF*(Y12/(Y12+Y13)+Y12/(Y12+Y23)+(Y12+Y23)/Y13+(Y12+Y13)/Y23+
0094      &    Y13L*(4D0*Y12**2+2D0*Y12*Y13+4D0*Y12*Y23+Y13*Y23)/
0095      &    (Y12+Y23)**2+Y23L*(4D0*Y12**2+2D0*Y12*Y23+4D0*Y12*Y13+
0096      &    Y13*Y23)/(Y12+Y13)**2)/WT1+
0097      &    CN*(Y13L*Y13/(Y12+Y23)+Y23L*Y23/(Y12+Y13))/WT1+(CN-2D0*CF)*
0098      &    ((Y12**2+(Y12+Y13)**2)*(Y12L*Y23L-Y12L*Y12M-Y23L*
0099      &    Y23M+1.644934D0-Y12I-Y23I)/(Y13*Y23)+(Y12**2+(Y12+Y23)**2)*
0100      &    (Y12L*Y13L-Y12L*Y12M-Y13L*Y13M+1.644934D0-Y12I-Y13I)/
0101      &    (Y13*Y23)+(Y13**2+Y23**2)/(Y13*Y23*(Y13+Y23))-
0102      &    2D0*Y12L*Y12**2/(Y13+Y23)**2-4D0*Y12L*Y12/(Y13+Y23))/WT1-
0103      &    CN*(Y13L*Y23L-Y13L*Y13M-Y23L*Y23M+1.644934D0-Y13I-Y23I)
0104           IF(1D0+WTOPT+ALS2PI*WT2.LE.0D0) MSTJ(121)=1
0105           IF(1D0+WTOPT+ALS2PI*WT2.LE.WTMAX*PYR(0)) GOTO 110
0106           PARJ(156)=(WTOPT+ALS2PI*WT2)/(1D0+WTOPT+ALS2PI*WT2)
0107
0108         ELSEIF(MSTJ(101).EQ.2.AND.MSTJ(110).EQ.2) THEN
0109 C...Second order corrections; Zhu parametrization of ERT.
0110           ZX=(Y23-Y13)**2
0111           ZY=1D0-Y12
0112           IZA=0
0113           DO 120 IY=1,5
0114             IF(ABS(CUT-0.01D0*IY).LT.0.0001D0) IZA=IY
0115   120     CONTINUE
0116           IF(IZA.NE.0) THEN
0117             IZ=IZA
0118             WT2=ZHUP(IZ,1)+ZHUP(IZ,2)*ZX+ZHUP(IZ,3)*ZX**2+(ZHUP(IZ,4)+
0119      &      ZHUP(IZ,5)*ZX)*ZY+(ZHUP(IZ,6)+ZHUP(IZ,7)*ZX)*ZY**2+
0120      &      (ZHUP(IZ,8)+ZHUP(IZ,9)*ZX)*ZY**3+ZHUP(IZ,10)/(ZX-ZY**2)+
0121      &      ZHUP(IZ,11)/(1D0-ZY)+ZHUP(IZ,12)/ZY
0122           ELSE
0123             IZ=100D0*CUT
0124             WTL=ZHUP(IZ,1)+ZHUP(IZ,2)*ZX+ZHUP(IZ,3)*ZX**2+(ZHUP(IZ,4)+
0125      &      ZHUP(IZ,5)*ZX)*ZY+(ZHUP(IZ,6)+ZHUP(IZ,7)*ZX)*ZY**2+
0126      &      (ZHUP(IZ,8)+ZHUP(IZ,9)*ZX)*ZY**3+ZHUP(IZ,10)/(ZX-ZY**2)+
0127      &      ZHUP(IZ,11)/(1D0-ZY)+ZHUP(IZ,12)/ZY
0128             IZ=IZ+1
0129             WTU=ZHUP(IZ,1)+ZHUP(IZ,2)*ZX+ZHUP(IZ,3)*ZX**2+(ZHUP(IZ,4)+
0130      &      ZHUP(IZ,5)*ZX)*ZY+(ZHUP(IZ,6)+ZHUP(IZ,7)*ZX)*ZY**2+
0131      &      (ZHUP(IZ,8)+ZHUP(IZ,9)*ZX)*ZY**3+ZHUP(IZ,10)/(ZX-ZY**2)+
0132      &      ZHUP(IZ,11)/(1D0-ZY)+ZHUP(IZ,12)/ZY
0133             WT2=WTL+(WTU-WTL)*(100D0*CUT+1D0-IZ)
0134           ENDIF
0135           IF(1D0+WTOPT+2D0*ALS2PI*WT2.LE.0D0) MSTJ(121)=1
0136           IF(1D0+WTOPT+2D0*ALS2PI*WT2.LE.WTMAX*PYR(0)) GOTO 110
0137           PARJ(156)=(WTOPT+2D0*ALS2PI*WT2)/(1D0+WTOPT+2D0*ALS2PI*WT2)
0138         ENDIF
0139
0140 C...Impose mass cuts (gives two jets). For fixed jet number new try.
0141         X1=1D0-Y23
0142         X2=1D0-Y13
0143         X3=1D0-Y12
0144         IF(4D0*Y23*Y13*Y12/X3**2.LE.QME) NJET=2
0145         IF(MOD(MSTJ(103),4).GE.2.AND.IABS(MSTJ(101)).LE.1.AND.QME*X3+
0146      &  0.5D0*QME**2+(0.5D0*QME+0.25D0*QME**2)*((1D0-X2)/(1D0-X1)+
0147      &  (1D0-X1)/(1D0-X2)).GT.(X1**2+X2**2)*PYR(0)) NJET=2
0148         IF(MSTJ(101).EQ.-1.AND.NJET.EQ.2) GOTO 100
0149
0150 C...Scalar gluon model (first order only, no mass effects).
0151       ELSE
0152   130   NJET=3
0153   140   X3=SQRT(4D0*CUT**2+PYR(0)*((1D0-CUT)**2-4D0*CUT**2))
0154         IF(LOG((X3-CUT)/CUT).LE.PYR(0)*LOG((1D0-2D0*CUT)/CUT)) GOTO 140
0155         YD=SIGN(2D0*CUT*((X3-CUT)/CUT)**PYR(0)-X3,PYR(0)-0.5D0)
0156         X1=1D0-0.5D0*(X3+YD)
0157         X2=1D0-0.5D0*(X3-YD)
0158         IF(4D0*(1D0-X1)*(1D0-X2)*(1D0-X3)/X3**2.LE.QME) NJET=2
0159         IF(MSTJ(102).GE.2) THEN
0160           IF(X3**2-2D0*(1D0+X3)*(1D0-X1)*(1D0-X2)*PARJ(171).LT.
0161      &    X3**2*PYR(0)) NJET=2
0162         ENDIF
0163         IF(MSTJ(101).EQ.-1.AND.NJET.EQ.2) GOTO 130
0164       ENDIF
0165
0166       RETURN
0167       END