pyxtee.f

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

 
C*********************************************************************
 
C...PYXTEE
C...Calculates total cross-section, including initial state
C...radiation effects.
 
      SUBROUTINE pyxtee(KFL,ECM,XTOT)
 
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)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      SAVE /pydat1/,/pydat2/
 
C...Status, (optimized) Q^2 scale, alpha_strong.
      parj(151)=ecm
      mstj(119)=10*mstj(102)+kfl
      IF(mstj(111).EQ.0) THEN
        q2r=ecm**2
      ELSEIF(mstu(111).EQ.0) THEN
        parj(168)=min(1d0,max(parj(128),exp(-12d0*paru(1)/
     &  ((33d0-2d0*mstu(112))*paru(111)))))
        q2r=parj(168)*ecm**2
      ELSE
        parj(168)=min(1d0,max(parj(128),paru(112)/ecm,
     &  (2d0*paru(112)/ecm)**2))
        q2r=parj(168)*ecm**2
      ENDIF
      alspi=pyalps(q2r)/paru(1)
 
C...QCD corrections factor in R.
      IF(mstj(101).EQ.0.OR.mstj(109).EQ.1) THEN
        rqcd=1d0
      ELSEIF(iabs(mstj(101)).EQ.1.AND.mstj(109).EQ.0) 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)
      ELSEIF(iabs(mstj(101)).EQ.1) THEN
        rqcd=1d0+(3d0/4d0)*alspi
      ELSE
        rqcd=1d0+(3d0/4d0)*alspi-(3d0/32d0+0.519d0*mstu(118))*alspi**2
      ENDIF
 
C...Calculate Z0 width if default value not acceptable.
      IF(mstj(102).GE.3) THEN
        rva=3d0*(3d0+(4d0*paru(102)-1d0)**2)+6d0*rqcd*(2d0+
     &  (1d0-8d0*paru(102)/3d0)**2+(4d0*paru(102)/3d0-1d0)**2)
        DO 100 kflc=5,6
          vq=1d0
          IF(mod(mstj(103),2).EQ.1) vq=sqrt(max(0d0,1d0-
     &    (2d0*pymass(kflc)/ ecm)**2))
          IF(kflc.EQ.5) vf=4d0*paru(102)/3d0-1d0
          IF(kflc.EQ.6) vf=1d0-8d0*paru(102)/3d0
          rva=rva+3d0*rqcd*(0.5d0*vq*(3d0-vq**2)*vf**2+vq**3)
  100   CONTINUE
        parj(124)=paru(101)*parj(123)*rva/(48d0*paru(102)*
     &  (1d0-paru(102)))
      ENDIF
 
C...Calculate propagator and related constants for QFD case.
      poll=1d0-parj(131)*parj(132)
      IF(mstj(102).GE.2) THEN
        sff=1d0/(16d0*paru(102)*(1d0-paru(102)))
        sfw=ecm**4/((ecm**2-parj(123)**2)**2+(parj(123)*parj(124))**2)
        sfi=sfw*(1d0-(parj(123)/ecm)**2)
        ve=4d0*paru(102)-1d0
        sf1i=sff*(ve*poll+parj(132)-parj(131))
        sf1w=sff**2*((ve**2+1d0)*poll+2d0*ve*(parj(132)-parj(131)))
        hf1i=sfi*sf1i
        hf1w=sfw*sf1w
      ENDIF
 
C...Loop over different flavours: charge, velocity.
      rtot=0d0
      rqq=0d0
      rqv=0d0
      rva=0d0
      DO 110 kflc=1,max(mstj(104),kfl)
        IF(kfl.GT.0.AND.kflc.NE.kfl) goto 110
        mstj(93)=1
        pmq=pymass(kflc)
        IF(ecm.LT.2d0*pmq+parj(127)) goto 110
        qf=kchg(kflc,1)/3d0
        vq=1d0
        IF(mod(mstj(103),2).EQ.1) vq=sqrt(1d0-(2d0*pmq/ecm)**2)
 
C...Calculate R and sum of charges for QED or QFD case.
        rqq=rqq+3d0*qf**2*poll
        IF(mstj(102).LE.1) THEN
          rtot=rtot+3d0*0.5d0*vq*(3d0-vq**2)*qf**2*poll
        ELSE
          vf=sign(1d0,qf)-4d0*qf*paru(102)
          rqv=rqv-6d0*qf*vf*sf1i
          rva=rva+3d0*(vf**2+1d0)*sf1w
          rtot=rtot+3d0*(0.5d0*vq*(3d0-vq**2)*(qf**2*poll-
     &    2d0*qf*vf*hf1i+vf**2*hf1w)+vq**3*hf1w)
        ENDIF
  110 CONTINUE
      rsum=rqq
      IF(mstj(102).GE.2) rsum=rqq+sfi*rqv+sfw*rva
 
C...Calculate cross-section, including QCD corrections.
      parj(141)=rqq
      parj(142)=rtot
      parj(143)=rtot*rqcd
      parj(144)=parj(143)
      parj(145)=parj(141)*86.8d0/ecm**2
      parj(146)=parj(142)*86.8d0/ecm**2
      parj(147)=parj(143)*86.8d0/ecm**2
      parj(148)=parj(147)
      parj(157)=rsum*rqcd
      parj(158)=0d0
      parj(159)=0d0
      xtot=parj(147)
      IF(mstj(107).LE.0) RETURN
 
C...Virtual cross-section.
      xkl=parj(135)
      xku=min(parj(136),1d0-(2d0*parj(127)/ecm)**2)
      ale=2d0*log(ecm/pymass(11))-1d0
      sigv=ale/3d0+2d0*log(ecm**2/(pymass(13)*pymass(15)))/3d0-4d0/3d0+
     &1.526d0*log(ecm**2/0.932d0)
 
C...Soft and hard radiative cross-section in QED case.
      IF(mstj(102).LE.1) THEN
        sigv=1.5d0*ale-0.5d0+paru(1)**2/3d0+2d0*sigv
        sigs=ale*(2d0*log(xkl)-log(1d0-xkl)-xkl)
        sigh=ale*(2d0*log(xku/xkl)-log((1d0-xku)/(1d0-xkl))-(xku-xkl))
 
C...Soft and hard radiative cross-section in QFD case.
      ELSE
        szm=1d0-(parj(123)/ecm)**2
        szw=parj(123)*parj(124)/ecm**2
        parj(161)=-rqq/rsum
        parj(162)=-(rqq+rqv+rva)/rsum
        parj(163)=(rqv*(1d0-0.5d0*szm-sfi)+rva*(1.5d0-szm-sfw))/rsum
        parj(164)=(rqv*szw**2*(1d0-2d0*sfw)+rva*(2d0*sfi+szw**2-
     &  4d0+3d0*szm-szm**2))/(szw*rsum)
        sigv=1.5d0*ale-0.5d0+paru(1)**2/3d0+((2d0*rqq+sfi*rqv)/
     &  rsum)*sigv+(szw*sfw*rqv/rsum)*paru(1)*20d0/9d0
        sigs=ale*(2d0*log(xkl)+parj(161)*log(1d0-xkl)+parj(162)*xkl+
     &  parj(163)*log(((xkl-szm)**2+szw**2)/(szm**2+szw**2))+
     &  parj(164)*(atan((xkl-szm)/szw)-atan(-szm/szw)))
        sigh=ale*(2d0*log(xku/xkl)+parj(161)*log((1d0-xku)/
     &  (1d0-xkl))+parj(162)*(xku-xkl)+parj(163)*
     &  log(((xku-szm)**2+szw**2)/((xkl-szm)**2+szw**2))+
     &  parj(164)*(atan((xku-szm)/szw)-atan((xkl-szm)/szw)))
      ENDIF
 
C...Total cross-section and fraction of hard photon events.
      parj(160)=sigh/(paru(1)/paru(101)+sigv+sigs+sigh)
      parj(157)=rsum*(1d0+(paru(101)/paru(1))*(sigv+sigs+sigh))*rqcd
      parj(144)=parj(157)
      parj(148)=parj(144)*86.8d0/ecm**2
      xtot=parj(148)
 
      RETURN
      END