pygano.f

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

 
C*********************************************************************
 
C...PYGANO
C...Evaluates the parton distributions of the anomalous photon,
C...inhomogeneously evolved from a scale P2 (where it vanishes) 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.
C...Adapted from SaSgam library, authors G.A. Schuler and T. Sjostrand.
 
      SUBROUTINE pygano(KF,X,Q2,P2,ALAM,XPGA,VXPGA)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Local arrays and data.
      dimension xpga(-6:6), vxpga(-6:6), alamsq(3:5)
      DATA pmc/1.3d0/, pmb/4.6d0/, aem/0.007297d0/, aem2pi/0.0011614d0/
 
C...Reset output.
      DO 100 kfl=-6,6
        xpga(kfl)=0d0
        vxpga(kfl)=0d0
  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)**(2d0/27d0))**2
      alamsq(4)=alam**2
      alamsq(5)=(alam*(alam/pmb)**(2d0/23d0))**2
      p2eff=max(p2,1.2d0*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=(6d0/(33d0-2d0*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=(6d0/(33d0-2d0*nfq))*log(log(q2div/alamsq(nfq))/
     &      log(p2eff/alamsq(nfq)))
            snfp=(6d0/(33d0-2d0*(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=(6d0/(33d0-2d0*4))*log(log(q2div/alamsq(4))/
     &      log(p2eff/alamsq(4)))
            snf3=(6d0/(33d0-2d0*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=(6d0/(33d0-2d0*nfq))*log(log(q2eff/alamsq(nfq))/
     &    log(p2eff/alamsq(nfq)))
          IF(nfq.EQ.5.AND.nfp.EQ.4) THEN
            q2div=pmb**2
            snfq=(6d0/(33d0-2d0*nfq))*log(log(q2div/alamsq(nfq))/
     &      log(p2eff/alamsq(nfq)))
            snfp=(6d0/(33d0-2d0*(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=(6d0/(33d0-2d0*nfq))*log(log(q2eff/alamsq(nfq))/
     &    log(p2eff/alamsq(nfq)))
        ENDIF
 
C...Evaluate flavour-dependent prefactor (charge^2 etc.).
        chsq=1d0/9d0
        IF(kfl.EQ.2.OR.kfl.EQ.4) chsq=4d0/9d0
        fac=aem2pi*2d0*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.5d0+2.49d0*s+26.9d0*s**2)/(1d0+32.3d0*s**2)*x**2 +
     &    (1.5d0-0.49d0*s+7.83d0*s**2)/(1d0+7.68d0*s**2)*(1d0-x)**2 +
     &    1.5d0*s/(1d0-3.2d0*s+7d0*s**2)*x*(1d0-x)) *
     &    x**(1d0/(1d0+0.58d0*s)) * (1d0-x**2)**(2.5d0*s/(1d0+10d0*s))
          xglu= 2d0*s/(1d0+4d0*s+7d0*s**2) *
     &    x**(-1.67d0*s/(1d0+2d0*s)) * (1d0-x**2)**(1.2d0*s) *
     &    ((4d0*x**2+7d0*x+4d0)*(1d0-x)/3d0 - 2d0*x*(1d0+x)*xl)
          xsea= 0.333d0*s**2/(1d0+4.90d0*s+4.69d0*s**2+21.4d0*s**3) *
     &    x**(-1.18d0*s/(1d0+1.22d0*s)) * (1d0-x)**(1.2d0*s) *
     &    ((8d0-73d0*x+62d0*x**2)*(1d0-x)/9d0 +
     &    (3d0-8d0*x**2/3d0)*x*xl + (2d0*x-1d0)*x*xl**2)
 
C...Threshold factors for c and b sea.
          sll=log(log(q2eff/alam**2)/log(p2eff/alam**2))
          xchm=0d0
          IF(q2.GT.pmc**2.AND.q2.GT.1.001d0*p2eff) THEN
            sch=max(0d0,log(log(pmc**2/alam**2)/log(p2eff/alam**2)))
            xchm=xsea*(1d0-(sch/sll)**3)
          ENDIF
          xbot=0d0
          IF(q2.GT.pmb**2.AND.q2.GT.1.001d0*p2eff) THEN
            sbt=max(0d0,log(log(pmb**2/alam**2)/log(p2eff/alam**2)))
            xbot=xsea*(1d0-(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
        vxpga(kfl)=vxpga(kfl)+fac*xval
  110 CONTINUE
      DO 120 kfl=1,5
        xpga(-kfl)=xpga(kfl)
        vxpga(-kfl)=vxpga(kfl)
  120 CONTINUE
 
      RETURN
      END