pygano.F

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

      SUBROUTINE pygano(KF,X,Q2,P2,ALAM,XPGA)
C...Purpose: to evaluate the parton distributions of the anomalous
C...photon, inhomogeneously evolved from a scale P2 (where it vanishes)
C...to Q2.
C...KF=0 gives the sum over (up to) 5 flavours,
C...KF<0 limits to flavours up to abs(KF),
C...KF>0 is for flavour KF only.
C...ALAM is the 4-flavour Lambda, which is automatically converted
C...to 3- and 5-flavour equivalents as needed.
      dimension xpga(-6:6),alamsq(3:5)
      DATA pmc/1.3/, pmb/4.6/, aem/0.007297/, aem2pi/0.0011614/
 
C...Reset output.
      DO 100 kfl=-6,6
      xpga(kfl)=0.
  100 CONTINUE
      IF(q2.LE.p2) RETURN
      kfa=iabs(kf)
 
C...Calculate Lambda; protect against unphysical Q2 and P2 input.
      alamsq(3)=(alam*(pmc/alam)**(2./27.))**2
      alamsq(4)=alam**2
      alamsq(5)=(alam*(alam/pmb)**(2./23.))**2
      p2eff=max(p2,1.2*alamsq(3))
      IF(kf.EQ.4) p2eff=max(p2eff,pmc**2)
      IF(kf.EQ.5) p2eff=max(p2eff,pmb**2)
      q2eff=max(q2,p2eff)
      xl=-log(x)
 
C...Find number of flavours at lower and upper scale.
      nfp=4
      IF(p2eff.LT.pmc**2) nfp=3
      IF(p2eff.GT.pmb**2) nfp=5
      nfq=4
      IF(q2eff.LT.pmc**2) nfq=3
      IF(q2eff.GT.pmb**2) nfq=5
 
C...Define range of flavour loop.
      IF(kf.EQ.0) THEN
        kflmn=1
        kflmx=5
      ELSEIF(kf.LT.0) THEN
        kflmn=1
        kflmx=kfa
      ELSE
        kflmn=kfa
        kflmx=kfa
      ENDIF
 
C...Loop over flavours the photon can branch into.
      DO 110 kfl=kflmn,kflmx
 
C...Light flavours: calculate t range and (approximate) s range.
      IF(kfl.LE.3.AND.(kfl.EQ.1.OR.kfl.EQ.kf)) THEN
        tdiff=log(q2eff/p2eff)
        s=(6./(33.-2.*nfq))*log(log(q2eff/alamsq(nfq))/
     &  log(p2eff/alamsq(nfq)))
        IF(nfq.GT.nfp) THEN
          q2div=pmb**2
          IF(nfq.EQ.4) q2div=pmc**2
          snfq=(6./(33.-2.*nfq))*log(log(q2div/alamsq(nfq))/
     &    log(p2eff/alamsq(nfq)))
          snfp=(6./(33.-2.*(nfq-1)))*log(log(q2div/alamsq(nfq-1))/
     &    log(p2eff/alamsq(nfq-1)))
          s=s+(log(q2div/p2eff)/log(q2eff/p2eff))*(snfp-snfq)
        ENDIF
        IF(nfq.EQ.5.AND.nfp.EQ.3) THEN
          q2div=pmc**2
          snf4=(6./(33.-2.*4))*log(log(q2div/alamsq(4))/
     &    log(p2eff/alamsq(4)))
          snf3=(6./(33.-2.*3))*log(log(q2div/alamsq(3))/
     &    log(p2eff/alamsq(3)))
          s=s+(log(q2div/p2eff)/log(q2eff/p2eff))*(snf3-snf4)
        ENDIF
 
C...u and s quark do not need a separate treatment when d has been done.
      ELSEIF(kfl.EQ.2.OR.kfl.EQ.3) THEN
 
C...Charm: as above, but only include range above c threshold.
      ELSEIF(kfl.EQ.4) THEN
        IF(q2.LE.pmc**2) goto 110
        p2eff=max(p2eff,pmc**2)
        q2eff=max(q2eff,p2eff)
        tdiff=log(q2eff/p2eff)
        s=(6./(33.-2.*nfq))*log(log(q2eff/alamsq(nfq))/
     &  log(p2eff/alamsq(nfq)))
        IF(nfq.EQ.5.AND.nfp.EQ.4) THEN
          q2div=pmb**2
          snfq=(6./(33.-2.*nfq))*log(log(q2div/alamsq(nfq))/
     &    log(p2eff/alamsq(nfq)))
          snfp=(6./(33.-2.*(nfq-1)))*log(log(q2div/alamsq(nfq-1))/
     &    log(p2eff/alamsq(nfq-1)))
          s=s+(log(q2div/p2eff)/log(q2eff/p2eff))*(snfp-snfq)
        ENDIF
 
C...Bottom: as above, but only include range above b threshold.
      ELSEIF(kfl.EQ.5) THEN
        IF(q2.LE.pmb**2) goto 110
        p2eff=max(p2eff,pmb**2)
        q2eff=max(q2,p2eff)
        tdiff=log(q2eff/p2eff)
        s=(6./(33.-2.*nfq))*log(log(q2eff/alamsq(nfq))/
     &  log(p2eff/alamsq(nfq)))
      ENDIF
 
C...Evaluate flavour-dependent prefactor (charge^2 etc.).
      chsq=1./9.
      IF(kfl.EQ.2.OR.kfl.EQ.4) chsq=4./9.
      fac=aem2pi*2.*chsq*tdiff
 
C...Evaluate parton distributions (normalized to unit momentum sum).
      IF(kfl.EQ.1.OR.kfl.EQ.4.OR.kfl.EQ.5.OR.kfl.EQ.kf) THEN
        xval= ((1.5+2.49*s+26.9*s**2)/(1.+32.3*s**2)*x**2 +
     &  (1.5-0.49*s+7.83*s**2)/(1.+7.68*s**2)*(1.-x)**2 +
     &  1.5*s/(1.-3.2*s+7.*s**2)*x*(1.-x)) *
     &  x**(1./(1.+0.58*s)) * (1.-x**2)**(2.5*s/(1.+10.*s))
        xglu= 2.*s/(1.+4.*s+7.*s**2) *
     &  x**(-1.67*s/(1.+2.*s)) * (1.-x**2)**(1.2*s) *
     &  ((4.*x**2+7.*x+4.)*(1.-x)/3. - 2.*x*(1.+x)*xl)
        xsea= 0.333*s**2/(1.+4.90*s+4.69*s**2+21.4*s**3) *
     &  x**(-1.18*s/(1.+1.22*s)) * (1.-x)**(1.2*s) *
     &  ((8.-73.*x+62.*x**2)*(1.-x)/9. + (3.-8.*x**2/3.)*x*xl +
     &  (2.*x-1.)*x*xl**2)
 
C...Threshold factors for c and b sea.
        sll=log(log(q2eff/alam**2)/log(p2eff/alam**2))
        xchm=0.
        IF(q2.GT.pmc**2.AND.q2.GT.1.001*p2eff) THEN
          sch=max(0.,log(log(pmc**2/alam**2)/log(p2eff/alam**2)))
          xchm=xsea*(1.-(sch/sll)**3)
        ENDIF
        xbot=0.
        IF(q2.GT.pmb**2.AND.q2.GT.1.001*p2eff) THEN
          sbt=max(0.,log(log(pmb**2/alam**2)/log(p2eff/alam**2)))
          xbot=xsea*(1.-(sbt/sll)**3)
        ENDIF
      ENDIF
 
C...Add contribution of each valence flavour.
      xpga(0)=xpga(0)+fac*xglu
      xpga(1)=xpga(1)+fac*xsea
      xpga(2)=xpga(2)+fac*xsea
      xpga(3)=xpga(3)+fac*xsea
      xpga(4)=xpga(4)+fac*xchm
      xpga(5)=xpga(5)+fac*xbot
      xpga(kfl)=xpga(kfl)+fac*xval
  110 CONTINUE
      DO 120 kfl=1,5
      xpga(-kfl)=xpga(kfl)
  120 CONTINUE
 
      RETURN
      END