pyxjet.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyxjet.f

 
C*********************************************************************
 
C...PYXJET
C...Selects number of jets in matrix element approach.
 
      SUBROUTINE pyxjet(ECM,NJET,CUT)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      SAVE /pydat1/
C...Local array and data.
      dimension zhut(5)
      DATA zhut/3.0922d0, 6.2291d0, 7.4782d0, 7.8440d0, 8.2560d0/
 
C...Trivial result for two-jets only, including parton shower.
      IF(mstj(101).EQ.0.OR.mstj(101).EQ.5) THEN
        cut=0d0
 
C...QCD and Abelian vector gluon theory: Q^2 for jet rate and R.
      ELSEIF(mstj(109).EQ.0.OR.mstj(109).EQ.2) THEN
        cf=4d0/3d0
        IF(mstj(109).EQ.2) cf=1d0
        IF(mstj(111).EQ.0) THEN
          q2=ecm**2
          q2r=ecm**2
        ELSEIF(mstu(111).EQ.0) THEN
          parj(169)=min(1d0,parj(129))
          q2=parj(169)*ecm**2
          parj(168)=min(1d0,max(parj(128),exp(-12d0*paru(1)/
     &    ((33d0-2d0*mstu(112))*paru(111)))))
          q2r=parj(168)*ecm**2
        ELSE
          parj(169)=min(1d0,max(parj(129),(2d0*paru(112)/ecm)**2))
          q2=parj(169)*ecm**2
          parj(168)=min(1d0,max(parj(128),paru(112)/ecm,
     &    (2d0*paru(112)/ecm)**2))
          q2r=parj(168)*ecm**2
        ENDIF
 
C...alpha_strong for R and R itself.
        alspi=(3d0/4d0)*cf*pyalps(q2r)/paru(1)
        IF(iabs(mstj(101)).EQ.1) THEN
          rqcd=1d0+alspi
        ELSEIF(mstj(109).EQ.0) THEN
          rqcd=1d0+alspi+(1.986d0-0.115d0*mstu(118))*alspi**2
          IF(mstj(111).EQ.1) rqcd=max(1d0,rqcd+
     &    (33d0-2d0*mstu(112))/12d0*log(parj(168))*alspi**2)
        ELSE
          rqcd=1d0+alspi-(3d0/32d0+0.519d0*mstu(118))*(4d0*alspi/3d0)**2
        ENDIF
 
C...alpha_strong for jet rate. Initial value for y cut.
        alspi=(3d0/4d0)*cf*pyalps(q2)/paru(1)
        cut=max(0.001d0,parj(125),(parj(126)/ecm)**2)
        IF(iabs(mstj(101)).LE.1.OR.(mstj(109).EQ.0.AND.mstj(111).EQ.0))
     &  cut=max(cut,exp(-sqrt(0.75d0/alspi))/2d0)
        IF(mstj(110).EQ.2) cut=max(0.01d0,min(0.05d0,cut))
 
C...Parametrization of first order three-jet cross-section.
  100   IF(mstj(101).EQ.0.OR.cut.GE.0.25d0) THEN
          parj(152)=0d0
        ELSE
          parj(152)=(2d0*alspi/3d0)*((3d0-6d0*cut+2d0*log(cut))*
     &    log(cut/(1d0-2d0*cut))+(2.5d0+1.5d0*cut-6.571d0)*
     &    (1d0-3d0*cut)+5.833d0*(1d0-3d0*cut)**2-3.894d0*
     &    (1d0-3d0*cut)**3+1.342d0*(1d0-3d0*cut)**4)/rqcd
          IF(mstj(109).EQ.2.AND.(mstj(101).EQ.2.OR.mstj(101).LE.-2))
     &    parj(152)=0d0
        ENDIF
 
C...Parametrization of second order three-jet cross-section.
        IF(iabs(mstj(101)).LE.1.OR.mstj(101).EQ.3.OR.mstj(109).EQ.2.OR.
     &  cut.GE.0.25d0) THEN
          parj(153)=0d0
        ELSEIF(mstj(110).LE.1) THEN
          ct=log(1d0/cut-2d0)
          parj(153)=alspi**2*ct**2*(2.419d0+0.5989d0*ct+0.6782d0*ct**2-
     &    0.2661d0*ct**3+0.01159d0*ct**4)/rqcd
 
C...Interpolation in second/first order ratio for Zhu parametrization.
        ELSEIF(mstj(110).EQ.2) THEN
          iza=0
          DO 110 iy=1,5
            IF(abs(cut-0.01d0*iy).LT.0.0001d0) iza=iy
  110     CONTINUE
          IF(iza.NE.0) THEN
            zhurat=zhut(iza)
          ELSE
            iz=100d0*cut
            zhurat=zhut(iz)+(100d0*cut-iz)*(zhut(iz+1)-zhut(iz))
          ENDIF
          parj(153)=alspi*parj(152)*zhurat
        ENDIF
 
C...Shift in second order three-jet cross-section with optimized Q^2.
        IF(mstj(111).EQ.1.AND.iabs(mstj(101)).GE.2.AND.mstj(101).NE.3
     &  .AND.cut.LT.0.25d0) parj(153)=parj(153)+
     &  (33d0-2d0*mstu(112))/12d0*log(parj(169))*alspi*parj(152)
 
C...Parametrization of second order four-jet cross-section.
        IF(iabs(mstj(101)).LE.1.OR.cut.GE.0.125d0) THEN
          parj(154)=0d0
        ELSE
          ct=log(1d0/cut-5d0)
          IF(cut.LE.0.018d0) THEN
            xqqgg=6.349d0-4.330d0*ct+0.8304d0*ct**2
            IF(mstj(109).EQ.2) xqqgg=(4d0/3d0)**2*(3.035d0-2.091d0*ct+
     &      0.4059d0*ct**2)
            xqqqq=1.25d0*(-0.1080d0+0.01486d0*ct+0.009364d0*ct**2)
            IF(mstj(109).EQ.2) xqqqq=8d0*xqqqq
          ELSE
            xqqgg=-0.09773d0+0.2959d0*ct-0.2764d0*ct**2+0.08832d0*ct**3
            IF(mstj(109).EQ.2) xqqgg=(4d0/3d0)**2*(-0.04079d0+
     &      0.1340d0*ct-0.1326d0*ct**2+0.04365d0*ct**3)
            xqqqq=1.25d0*(0.003661d0-0.004888d0*ct-0.001081d0*ct**2+
     &      0.002093d0*ct**3)
            IF(mstj(109).EQ.2) xqqqq=8d0*xqqqq
          ENDIF
          parj(154)=alspi**2*ct**2*(xqqgg+xqqqq)/rqcd
          parj(155)=xqqqq/(xqqgg+xqqqq)
        ENDIF
 
C...If negative three-jet rate, change y' optimization parameter.
        IF(mstj(111).EQ.1.AND.parj(152)+parj(153).LT.0d0.AND.
     &  parj(169).LT.0.99d0) THEN
          parj(169)=min(1d0,1.2d0*parj(169))
          q2=parj(169)*ecm**2
          alspi=(3d0/4d0)*cf*pyalps(q2)/paru(1)
          goto 100
        ENDIF
 
C...If too high cross-section, use harder cuts, or fail.
        IF(parj(152)+parj(153)+parj(154).GE.1) THEN
          IF(mstj(110).EQ.2.AND.cut.GT.0.0499d0.AND.mstj(111).EQ.1.AND.
     &    parj(169).LT.0.99d0) THEN
            parj(169)=min(1d0,1.2d0*parj(169))
            q2=parj(169)*ecm**2
            alspi=(3d0/4d0)*cf*pyalps(q2)/paru(1)
            goto 100
          ELSEIF(mstj(110).EQ.2.AND.cut.GT.0.0499d0) THEN
            CALL pyerrm(26,
     &      '(PYXJET:) no allowed y cut value for Zhu parametrization')
          ENDIF
          cut=0.26d0*(4d0*cut)**(parj(152)+parj(153)+
     &    parj(154))**(-1d0/3d0)
          IF(mstj(110).EQ.2) cut=max(0.01d0,min(0.05d0,cut))
          goto 100
        ENDIF
 
C...Scalar gluon (first order only).
      ELSE
        alspi=pyalps(ecm**2)/paru(1)
        cut=max(0.001d0,parj(125),(parj(126)/ecm)**2,exp(-3d0/alspi))
        parj(152)=0d0
        IF(cut.LT.0.25d0) parj(152)=(alspi/3d0)*((1d0-2d0*cut)*
     &  log((1d0-2d0*cut)/cut)+0.5d0*(9d0*cut**2-1d0))
        parj(153)=0d0
        parj(154)=0d0
      ENDIF
 
C...Select number of jets.
      parj(150)=cut
      IF(mstj(101).EQ.0.OR.mstj(101).EQ.5) THEN
        njet=2
      ELSEIF(mstj(101).LE.0) THEN
        njet=min(4,2-mstj(101))
      ELSE
        rnj=pyr(0)
        njet=2
        IF(parj(152)+parj(153)+parj(154).GT.rnj) njet=3
        IF(parj(154).GT.rnj) njet=4
      ENDIF
 
      RETURN
      END