pythiaeRHIC/pythia/pyx4jt.f

0001
0002 C*********************************************************************
0003
0004 C...PYX4JT
0005 C...Selects the kinematical variables of four-jet events.
0006
0007       SUBROUTINE PYX4JT(NJET,CUT,KFL,ECM,KFLN,X1,X2,X4,X12,X14)
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 arrays.
0017       DIMENSION WTA(4),WTB(4),WTC(4),WTD(4),WTE(4)
0018
0019 C...Common constants. Colour factors for QCD and Abelian gluon theory.
0020       PMQ=PYMASS(KFL)
0021       QME=(2D0*PMQ/ECM)**2
0022       CT=LOG(1D0/CUT-5D0)
0023       IF(MSTJ(109).EQ.0) THEN
0024         CF=4D0/3D0
0025         CN=3D0
0026         TR=2.5D0
0027       ELSE
0028         CF=1D0
0029         CN=0D0
0030         TR=15D0
0031       ENDIF
0032
0033 C...Choice of process (qqbargg or qqbarqqbar).
0034   100 NJET=4
0035       IT=1
0036       IF(PARJ(155).GT.PYR(0)) IT=2
0037       IF(MSTJ(101).LE.-3) IT=-MSTJ(101)-2
0038       IF(IT.EQ.1) WTMX=0.7D0/CUT**2
0039       IF(IT.EQ.1.AND.MSTJ(109).EQ.2) WTMX=0.6D0/CUT**2
0040       IF(IT.EQ.2) WTMX=0.1125D0*CF*TR/CUT**2
0041       ID=1
0042
0043 C...Sample the five kinematical variables (for qqgg preweighted in y34).
0044   110 Y134=3D0*CUT+(1D0-6D0*CUT)*PYR(0)
0045       Y234=3D0*CUT+(1D0-6D0*CUT)*PYR(0)
0046       IF(IT.EQ.1) Y34=(1D0-5D0*CUT)*EXP(-CT*PYR(0))
0047       IF(IT.EQ.2) Y34=CUT+(1D0-6D0*CUT)*PYR(0)
0048       IF(Y34.LE.Y134+Y234-1D0.OR.Y34.GE.Y134*Y234) GOTO 110
0049       VT=PYR(0)
0050       CP=COS(PARU(1)*PYR(0))
0051       Y14=(Y134-Y34)*VT
0052       Y13=Y134-Y14-Y34
0053       VB=Y34*(1D0-Y134-Y234+Y34)/((Y134-Y34)*(Y234-Y34))
0054       Y24=0.5D0*(Y234-Y34)*(1D0-4D0*SQRT(MAX(0D0,VT*(1D0-VT)*
0055      &VB*(1D0-VB)))*CP-(1D0-2D0*VT)*(1D0-2D0*VB))
0056       Y23=Y234-Y34-Y24
0057       Y12=1D0-Y134-Y23-Y24
0058       IF(MIN(Y12,Y13,Y14,Y23,Y24).LE.CUT) GOTO 110
0059       Y123=Y12+Y13+Y23
0060       Y124=Y12+Y14+Y24
0061
0062 C...Calculate matrix elements for qqgg or qqqq process.
0063       IC=0
0064       WTTOT=0D0
0065   120 IC=IC+1
0066       IF(IT.EQ.1) THEN
0067         WTA(IC)=(Y12*Y34**2-Y13*Y24*Y34+Y14*Y23*Y34+3D0*Y12*Y23*Y34+
0068      &  3D0*Y12*Y14*Y34+4D0*Y12**2*Y34-Y13*Y23*Y24+2D0*Y12*Y23*Y24-
0069      &  Y13*Y14*Y24-2D0*Y12*Y13*Y24+2D0*Y12**2*Y24+Y14*Y23**2+2D0*Y12*
0070      &  Y23**2+Y14**2*Y23+4D0*Y12*Y14*Y23+4D0*Y12**2*Y23+2D0*Y12*Y14**2+
0071      &  2D0*Y12*Y13*Y14+4D0*Y12**2*Y14+2D0*Y12**2*Y13+2D0*Y12**3)/
0072      &  (2D0*Y13*Y134*Y234*Y24)+(Y24*Y34+Y12*Y34+Y13*Y24-
0073      &  Y14*Y23+Y12*Y13)/(Y13*Y134**2)+2D0*Y23*(1D0-Y13)/
0074      &  (Y13*Y134*Y24)+Y34/(2D0*Y13*Y24)
0075         WTB(IC)=(Y12*Y24*Y34+Y12*Y14*Y34-Y13*Y24**2+Y13*Y14*Y24+2D0*Y12*
0076      &  Y14*Y24)/(Y13*Y134*Y23*Y14)+Y12*(1D0+Y34)*Y124/(Y134*Y234*Y14*
0077      &  Y24)-(2D0*Y13*Y24+Y14**2+Y13*Y23+2D0*Y12*Y13)/(Y13*Y134*Y14)+
0078      &  Y12*Y123*Y124/(2D0*Y13*Y14*Y23*Y24)
0079         WTC(IC)=-(5D0*Y12*Y34**2+2D0*Y12*Y24*Y34+2D0*Y12*Y23*Y34+
0080      &  2D0*Y12*Y14*Y34+2D0*Y12*Y13*Y34+4D0*Y12**2*Y34-Y13*Y24**2+
0081      &  Y14*Y23*Y24+Y13*Y23*Y24+Y13*Y14*Y24-Y12*Y14*Y24-Y13**2*Y24-
0082      &  3D0*Y12*Y13*Y24-Y14*Y23**2-Y14**2*Y23+Y13*Y14*Y23-
0083      &  3D0*Y12*Y14*Y23-Y12*Y13*Y23)/(4D0*Y134*Y234*Y34**2)+
0084      &  (3D0*Y12*Y34**2-3D0*Y13*Y24*Y34+3D0*Y12*Y24*Y34+
0085      &  3D0*Y14*Y23*Y34-Y13*Y24**2-Y12*Y23*Y34+6D0*Y12*Y14*Y34+
0086      &  2D0*Y12*Y13*Y34-2D0*Y12**2*Y34+Y14*Y23*Y24-3D0*Y13*Y23*Y24-
0087      &  2D0*Y13*Y14*Y24+4D0*Y12*Y14*Y24+2D0*Y12*Y13*Y24+
0088      &  3D0*Y14*Y23**2+2D0*Y14**2*Y23+2D0*Y14**2*Y12+
0089      &  2D0*Y12**2*Y14+6D0*Y12*Y14*Y23-2D0*Y12*Y13**2-
0090      &  2D0*Y12**2*Y13)/(4D0*Y13*Y134*Y234*Y34)
0091         WTC(IC)=WTC(IC)+(2D0*Y12*Y34**2-2D0*Y13*Y24*Y34+Y12*Y24*Y34+
0092      &  4D0*Y13*Y23*Y34+4D0*Y12*Y14*Y34+2D0*Y12*Y13*Y34+2D0*Y12**2*Y34-
0093      &  Y13*Y24**2+3D0*Y14*Y23*Y24+4D0*Y13*Y23*Y24-2D0*Y13*Y14*Y24+
0094      &  4D0*Y12*Y14*Y24+2D0*Y12*Y13*Y24+2D0*Y14*Y23**2+4D0*Y13*Y23**2+
0095      &  2D0*Y13*Y14*Y23+2D0*Y12*Y14*Y23+4D0*Y12*Y13*Y23+2D0*Y12*Y14**2+
0096      &  4D0*Y12**2*Y13+4D0*Y12*Y13*Y14+2D0*Y12**2*Y14)/
0097      &  (4D0*Y13*Y134*Y24*Y34)-(Y12*Y34**2-2D0*Y14*Y24*Y34-
0098      &  2D0*Y13*Y24*Y34-Y14*Y23*Y34+Y13*Y23*Y34+Y12*Y14*Y34+
0099      &  2D0*Y12*Y13*Y34-2D0*Y14**2*Y24-4D0*Y13*Y14*Y24-
0100      &  4D0*Y13**2*Y24-Y14**2*Y23-Y13**2*Y23+Y12*Y13*Y14-
0101      &  Y12*Y13**2)/(2D0*Y13*Y34*Y134**2)+(Y12*Y34**2-
0102      &  4D0*Y14*Y24*Y34-2D0*Y13*Y24*Y34-2D0*Y14*Y23*Y34-
0103      &  4D0*Y13*Y23*Y34-4D0*Y12*Y14*Y34-4D0*Y12*Y13*Y34-
0104      &  2D0*Y13*Y14*Y24+2D0*Y13**2*Y24+2D0*Y14**2*Y23-
0105      &  2D0*Y13*Y14*Y23-Y12*Y14**2-6D0*Y12*Y13*Y14-
0106      &  Y12*Y13**2)/(4D0*Y34**2*Y134**2)
0107         WTTOT=WTTOT+Y34*CF*(CF*WTA(IC)+(CF-0.5D0*CN)*WTB(IC)+
0108      &  CN*WTC(IC))/8D0
0109       ELSE
0110         WTD(IC)=(Y13*Y23*Y34+Y12*Y23*Y34-Y12**2*Y34+Y13*Y23*Y24+2D0*Y12*
0111      &  Y23*Y24-Y14*Y23**2+Y12*Y13*Y24+Y12*Y14*Y23+Y12*Y13*Y14)/(Y13**2*
0112      &  Y123**2)-(Y12*Y34**2-Y13*Y24*Y34+Y12*Y24*Y34-Y14*Y23*Y34-Y12*
0113      &  Y23*Y34-Y13*Y24**2+Y14*Y23*Y24-Y13*Y23*Y24-Y13**2*Y24+Y14*
0114      &  Y23**2)/(Y13**2*Y123*Y134)+(Y13*Y14*Y12+Y34*Y14*Y12-Y34**2*Y12+
0115      &  Y13*Y14*Y24+2D0*Y34*Y14*Y24-Y23*Y14**2+Y34*Y13*Y24+Y34*Y23*Y14+
0116      &  Y34*Y13*Y23)/(Y13**2*Y134**2)-(Y34*Y12**2-Y13*Y24*Y12+Y34*Y24*
0117      &  Y12-Y23*Y14*Y12-Y34*Y14*Y12-Y13*Y24**2+Y23*Y14*Y24-Y13*Y14*Y24-
0118      &  Y13**2*Y24+Y23*Y14**2)/(Y13**2*Y134*Y123)
0119         WTE(IC)=(Y12*Y34*(Y23-Y24+Y14+Y13)+Y13*Y24**2-Y14*Y23*Y24+Y13*
0120      &  Y23*Y24+Y13*Y14*Y24+Y13**2*Y24-Y14*Y23*(Y14+Y23+Y13))/(Y13*Y23*
0121      &  Y123*Y134)-Y12*(Y12*Y34-Y23*Y24-Y13*Y24-Y14*Y23-Y14*Y13)/(Y13*
0122      &  Y23*Y123**2)-(Y14+Y13)*(Y24+Y23)*Y34/(Y13*Y23*Y134*Y234)+
0123      &  (Y12*Y34*(Y14-Y24+Y23+Y13)+Y13*Y24**2-Y23*Y14*Y24+Y13*Y14*Y24+
0124      &  Y13*Y23*Y24+Y13**2*Y24-Y23*Y14*(Y14+Y23+Y13))/(Y13*Y14*Y134*
0125      &  Y123)-Y34*(Y34*Y12-Y14*Y24-Y13*Y24-Y23*Y14-Y23*Y13)/(Y13*Y14*
0126      &  Y134**2)-(Y23+Y13)*(Y24+Y14)*Y12/(Y13*Y14*Y123*Y124)
0127         WTTOT=WTTOT+CF*(TR*WTD(IC)+(CF-0.5D0*CN)*WTE(IC))/16D0
0128       ENDIF
0129
0130 C...Permutations of momenta in matrix element. Weighting.
0131   130 IF(IC.EQ.1.OR.IC.EQ.3.OR.ID.EQ.2.OR.ID.EQ.3) THEN
0132         YSAV=Y13
0133         Y13=Y14
0134         Y14=YSAV
0135         YSAV=Y23
0136         Y23=Y24
0137         Y24=YSAV
0138         YSAV=Y123
0139         Y123=Y124
0140         Y124=YSAV
0141       ENDIF
0142       IF(IC.EQ.2.OR.IC.EQ.4.OR.ID.EQ.3.OR.ID.EQ.4) THEN
0143         YSAV=Y13
0144         Y13=Y23
0145         Y23=YSAV
0146         YSAV=Y14
0147         Y14=Y24
0148         Y24=YSAV
0149         YSAV=Y134
0150         Y134=Y234
0151         Y234=YSAV
0152       ENDIF
0153       IF(IC.LE.3) GOTO 120
0154       IF(ID.EQ.1.AND.WTTOT.LT.PYR(0)*WTMX) GOTO 110
0155       IC=5
0156
0157 C...qqgg events: string configuration and event type.
0158       IF(IT.EQ.1) THEN
0159         IF(MSTJ(109).EQ.0.AND.ID.EQ.1) THEN
0160           PARJ(156)=Y34*(2D0*(WTA(1)+WTA(2)+WTA(3)+WTA(4))+4D0*(WTC(1)+
0161      &    WTC(2)+WTC(3)+WTC(4)))/(9D0*WTTOT)
0162           IF(WTA(2)+WTA(4)+2D0*(WTC(2)+WTC(4)).GT.PYR(0)*(WTA(1)+WTA(2)+
0163      &    WTA(3)+WTA(4)+2D0*(WTC(1)+WTC(2)+WTC(3)+WTC(4)))) ID=2
0164           IF(ID.EQ.2) GOTO 130
0165         ELSEIF(MSTJ(109).EQ.2.AND.ID.EQ.1) THEN
0166           PARJ(156)=Y34*(WTA(1)+WTA(2)+WTA(3)+WTA(4))/(8D0*WTTOT)
0167           IF(WTA(2)+WTA(4).GT.PYR(0)*(WTA(1)+WTA(2)+WTA(3)+WTA(4))) ID=2
0168           IF(ID.EQ.2) GOTO 130
0169         ENDIF
0170         MSTJ(120)=3
0171         IF(MSTJ(109).EQ.0.AND.0.5D0*Y34*(WTC(1)+WTC(2)+WTC(3)+
0172      &  WTC(4)).GT.PYR(0)*WTTOT) MSTJ(120)=4
0173         KFLN=21
0174
0175 C...Mass cuts. Kinematical variables out.
0176         IF(Y12.LE.CUT+QME) NJET=2
0177         IF(NJET.EQ.2) GOTO 150
0178         Q12=0.5D0*(1D0-SQRT(1D0-QME/Y12))
0179         X1=1D0-(1D0-Q12)*Y234-Q12*Y134
0180         X4=1D0-(1D0-Q12)*Y134-Q12*Y234
0181         X2=1D0-Y124
0182         X12=(1D0-Q12)*Y13+Q12*Y23
0183         X14=Y12-0.5D0*QME
0184         IF(Y134*Y234/((1D0-X1)*(1D0-X4)).LE.PYR(0)) NJET=2
0185
0186 C...qqbarqqbar events: string configuration, choose new flavour.
0187       ELSE
0188         IF(ID.EQ.1) THEN
0189           WTR=PYR(0)*(WTD(1)+WTD(2)+WTD(3)+WTD(4))
0190           IF(WTR.LT.WTD(2)+WTD(3)+WTD(4)) ID=2
0191           IF(WTR.LT.WTD(3)+WTD(4)) ID=3
0192           IF(WTR.LT.WTD(4)) ID=4
0193           IF(ID.GE.2) GOTO 130
0194         ENDIF
0195         MSTJ(120)=5
0196         PARJ(156)=CF*TR*(WTD(1)+WTD(2)+WTD(3)+WTD(4))/(16D0*WTTOT)
0197   140   KFLN=1+INT(5D0*PYR(0))
0198         IF(KFLN.NE.KFL.AND.0.2D0*PARJ(156).LE.PYR(0)) GOTO 140
0199         IF(KFLN.EQ.KFL.AND.1D0-0.8D0*PARJ(156).LE.PYR(0)) GOTO 140
0200         IF(KFLN.GT.MSTJ(104)) NJET=2
0201         PMQN=PYMASS(KFLN)
0202         QMEN=(2D0*PMQN/ECM)**2
0203
0204 C...Mass cuts. Kinematical variables out.
0205         IF(Y24.LE.CUT+QME.OR.Y13.LE.1.1D0*QMEN) NJET=2
0206         IF(NJET.EQ.2) GOTO 150
0207         Q24=0.5D0*(1D0-SQRT(1D0-QME/Y24))
0208         Q13=0.5D0*(1D0-SQRT(1D0-QMEN/Y13))
0209         X1=1D0-(1D0-Q24)*Y123-Q24*Y134
0210         X4=1D0-(1D0-Q24)*Y134-Q24*Y123
0211         X2=1D0-(1D0-Q13)*Y234-Q13*Y124
0212         X12=(1D0-Q24)*((1D0-Q13)*Y14+Q13*Y34)+Q24*((1D0-Q13)*Y12+
0213      &  Q13*Y23)
0214         X14=Y24-0.5D0*QME
0215         X34=(1D0-Q24)*((1D0-Q13)*Y23+Q13*Y12)+Q24*((1D0-Q13)*Y34+
0216      &  Q13*Y14)
0217         IF(PMQ**2+PMQN**2+MIN(X12,X34)*ECM**2.LE.
0218      &  (PARJ(127)+PMQ+PMQN)**2) NJET=2
0219         IF(Y123*Y134/((1D0-X1)*(1D0-X4)).LE.PYR(0)) NJET=2
0220       ENDIF
0221   150 IF(MSTJ(101).LE.-2.AND.NJET.EQ.2) GOTO 100
0222
0223       RETURN
0224       END