pythiaeRHIC/pythia/pyxjet.f

0001
0002 C*********************************************************************
0003
0004 C...PYXJET
0005 C...Selects number of jets in matrix element approach.
0006
0007       SUBROUTINE PYXJET(ECM,NJET,CUT)
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 and data.
0017       DIMENSION ZHUT(5)
0018       DATA ZHUT/3.0922D0, 6.2291D0, 7.4782D0, 7.8440D0, 8.2560D0/
0019
0020 C...Trivial result for two-jets only, including parton shower.
0021       IF(MSTJ(101).EQ.0.OR.MSTJ(101).EQ.5) THEN
0022         CUT=0D0
0023
0024 C...QCD and Abelian vector gluon theory: Q^2 for jet rate and R.
0025       ELSEIF(MSTJ(109).EQ.0.OR.MSTJ(109).EQ.2) THEN
0026         CF=4D0/3D0
0027         IF(MSTJ(109).EQ.2) CF=1D0
0028         IF(MSTJ(111).EQ.0) THEN
0029           Q2=ECM**2
0030           Q2R=ECM**2
0031         ELSEIF(MSTU(111).EQ.0) THEN
0032           PARJ(169)=MIN(1D0,PARJ(129))
0033           Q2=PARJ(169)*ECM**2
0034           PARJ(168)=MIN(1D0,MAX(PARJ(128),EXP(-12D0*PARU(1)/
0035      &    ((33D0-2D0*MSTU(112))*PARU(111)))))
0036           Q2R=PARJ(168)*ECM**2
0037         ELSE
0038           PARJ(169)=MIN(1D0,MAX(PARJ(129),(2D0*PARU(112)/ECM)**2))
0039           Q2=PARJ(169)*ECM**2
0040           PARJ(168)=MIN(1D0,MAX(PARJ(128),PARU(112)/ECM,
0041      &    (2D0*PARU(112)/ECM)**2))
0042           Q2R=PARJ(168)*ECM**2
0043         ENDIF
0044
0045 C...alpha_strong for R and R itself.
0046         ALSPI=(3D0/4D0)*CF*PYALPS(Q2R)/PARU(1)
0047         IF(IABS(MSTJ(101)).EQ.1) THEN
0048           RQCD=1D0+ALSPI
0049         ELSEIF(MSTJ(109).EQ.0) THEN
0050           RQCD=1D0+ALSPI+(1.986D0-0.115D0*MSTU(118))*ALSPI**2
0051           IF(MSTJ(111).EQ.1) RQCD=MAX(1D0,RQCD+
0052      &    (33D0-2D0*MSTU(112))/12D0*LOG(PARJ(168))*ALSPI**2)
0053         ELSE
0054           RQCD=1D0+ALSPI-(3D0/32D0+0.519D0*MSTU(118))*(4D0*ALSPI/3D0)**2
0055         ENDIF
0056
0057 C...alpha_strong for jet rate. Initial value for y cut.
0058         ALSPI=(3D0/4D0)*CF*PYALPS(Q2)/PARU(1)
0059         CUT=MAX(0.001D0,PARJ(125),(PARJ(126)/ECM)**2)
0060         IF(IABS(MSTJ(101)).LE.1.OR.(MSTJ(109).EQ.0.AND.MSTJ(111).EQ.0))
0061      &  CUT=MAX(CUT,EXP(-SQRT(0.75D0/ALSPI))/2D0)
0062         IF(MSTJ(110).EQ.2) CUT=MAX(0.01D0,MIN(0.05D0,CUT))
0063
0064 C...Parametrization of first order three-jet cross-section.
0065   100   IF(MSTJ(101).EQ.0.OR.CUT.GE.0.25D0) THEN
0066           PARJ(152)=0D0
0067         ELSE
0068           PARJ(152)=(2D0*ALSPI/3D0)*((3D0-6D0*CUT+2D0*LOG(CUT))*
0069      &    LOG(CUT/(1D0-2D0*CUT))+(2.5D0+1.5D0*CUT-6.571D0)*
0070      &    (1D0-3D0*CUT)+5.833D0*(1D0-3D0*CUT)**2-3.894D0*
0071      &    (1D0-3D0*CUT)**3+1.342D0*(1D0-3D0*CUT)**4)/RQCD
0072           IF(MSTJ(109).EQ.2.AND.(MSTJ(101).EQ.2.OR.MSTJ(101).LE.-2))
0073      &    PARJ(152)=0D0
0074         ENDIF
0075
0076 C...Parametrization of second order three-jet cross-section.
0077         IF(IABS(MSTJ(101)).LE.1.OR.MSTJ(101).EQ.3.OR.MSTJ(109).EQ.2.OR.
0078      &  CUT.GE.0.25D0) THEN
0079           PARJ(153)=0D0
0080         ELSEIF(MSTJ(110).LE.1) THEN
0081           CT=LOG(1D0/CUT-2D0)
0082           PARJ(153)=ALSPI**2*CT**2*(2.419D0+0.5989D0*CT+0.6782D0*CT**2-
0083      &    0.2661D0*CT**3+0.01159D0*CT**4)/RQCD
0084
0085 C...Interpolation in second/first order ratio for Zhu parametrization.
0086         ELSEIF(MSTJ(110).EQ.2) THEN
0087           IZA=0
0088           DO 110 IY=1,5
0089             IF(ABS(CUT-0.01D0*IY).LT.0.0001D0) IZA=IY
0090   110     CONTINUE
0091           IF(IZA.NE.0) THEN
0092             ZHURAT=ZHUT(IZA)
0093           ELSE
0094             IZ=100D0*CUT
0095             ZHURAT=ZHUT(IZ)+(100D0*CUT-IZ)*(ZHUT(IZ+1)-ZHUT(IZ))
0096           ENDIF
0097           PARJ(153)=ALSPI*PARJ(152)*ZHURAT
0098         ENDIF
0099
0100 C...Shift in second order three-jet cross-section with optimized Q^2.
0101         IF(MSTJ(111).EQ.1.AND.IABS(MSTJ(101)).GE.2.AND.MSTJ(101).NE.3
0102      &  .AND.CUT.LT.0.25D0) PARJ(153)=PARJ(153)+
0103      &  (33D0-2D0*MSTU(112))/12D0*LOG(PARJ(169))*ALSPI*PARJ(152)
0104
0105 C...Parametrization of second order four-jet cross-section.
0106         IF(IABS(MSTJ(101)).LE.1.OR.CUT.GE.0.125D0) THEN
0107           PARJ(154)=0D0
0108         ELSE
0109           CT=LOG(1D0/CUT-5D0)
0110           IF(CUT.LE.0.018D0) THEN
0111             XQQGG=6.349D0-4.330D0*CT+0.8304D0*CT**2
0112             IF(MSTJ(109).EQ.2) XQQGG=(4D0/3D0)**2*(3.035D0-2.091D0*CT+
0113      &      0.4059D0*CT**2)
0114             XQQQQ=1.25D0*(-0.1080D0+0.01486D0*CT+0.009364D0*CT**2)
0115             IF(MSTJ(109).EQ.2) XQQQQ=8D0*XQQQQ
0116           ELSE
0117             XQQGG=-0.09773D0+0.2959D0*CT-0.2764D0*CT**2+0.08832D0*CT**3
0118             IF(MSTJ(109).EQ.2) XQQGG=(4D0/3D0)**2*(-0.04079D0+
0119      &      0.1340D0*CT-0.1326D0*CT**2+0.04365D0*CT**3)
0120             XQQQQ=1.25D0*(0.003661D0-0.004888D0*CT-0.001081D0*CT**2+
0121      &      0.002093D0*CT**3)
0122             IF(MSTJ(109).EQ.2) XQQQQ=8D0*XQQQQ
0123           ENDIF
0124           PARJ(154)=ALSPI**2*CT**2*(XQQGG+XQQQQ)/RQCD
0125           PARJ(155)=XQQQQ/(XQQGG+XQQQQ)
0126         ENDIF
0127
0128 C...If negative three-jet rate, change y' optimization parameter.
0129         IF(MSTJ(111).EQ.1.AND.PARJ(152)+PARJ(153).LT.0D0.AND.
0130      &  PARJ(169).LT.0.99D0) THEN
0131           PARJ(169)=MIN(1D0,1.2D0*PARJ(169))
0132           Q2=PARJ(169)*ECM**2
0133           ALSPI=(3D0/4D0)*CF*PYALPS(Q2)/PARU(1)
0134           GOTO 100
0135         ENDIF
0136
0137 C...If too high cross-section, use harder cuts, or fail.
0138         IF(PARJ(152)+PARJ(153)+PARJ(154).GE.1) THEN
0139           IF(MSTJ(110).EQ.2.AND.CUT.GT.0.0499D0.AND.MSTJ(111).EQ.1.AND.
0140      &    PARJ(169).LT.0.99D0) THEN
0141             PARJ(169)=MIN(1D0,1.2D0*PARJ(169))
0142             Q2=PARJ(169)*ECM**2
0143             ALSPI=(3D0/4D0)*CF*PYALPS(Q2)/PARU(1)
0144             GOTO 100
0145           ELSEIF(MSTJ(110).EQ.2.AND.CUT.GT.0.0499D0) THEN
0146             CALL PYERRM(26,
0147      &      '(PYXJET:) no allowed y cut value for Zhu parametrization')
0148           ENDIF
0149           CUT=0.26D0*(4D0*CUT)**(PARJ(152)+PARJ(153)+
0150      &    PARJ(154))**(-1D0/3D0)
0151           IF(MSTJ(110).EQ.2) CUT=MAX(0.01D0,MIN(0.05D0,CUT))
0152           GOTO 100
0153         ENDIF
0154
0155 C...Scalar gluon (first order only).
0156       ELSE
0157         ALSPI=PYALPS(ECM**2)/PARU(1)
0158         CUT=MAX(0.001D0,PARJ(125),(PARJ(126)/ECM)**2,EXP(-3D0/ALSPI))
0159         PARJ(152)=0D0
0160         IF(CUT.LT.0.25D0) PARJ(152)=(ALSPI/3D0)*((1D0-2D0*CUT)*
0161      &  LOG((1D0-2D0*CUT)/CUT)+0.5D0*(9D0*CUT**2-1D0))
0162         PARJ(153)=0D0
0163         PARJ(154)=0D0
0164       ENDIF
0165
0166 C...Select number of jets.
0167       PARJ(150)=CUT
0168       IF(MSTJ(101).EQ.0.OR.MSTJ(101).EQ.5) THEN
0169         NJET=2
0170       ELSEIF(MSTJ(101).LE.0) THEN
0171         NJET=MIN(4,2-MSTJ(101))
0172       ELSE
0173         RNJ=PYR(0)
0174         NJET=2
0175         IF(PARJ(152)+PARJ(153)+PARJ(154).GT.RNJ) NJET=3
0176         IF(PARJ(154).GT.RNJ) NJET=4
0177       ENDIF
0178
0179       RETURN
0180       END