pyggam.F

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

      SUBROUTINE pyggam(ISET,X,Q2,P2,F2GM,XPDFGM)
C...Purpose: to construct the F2 and parton distributions of the photon
C...by summing homogeneous (VMD) and inhomogeneous (anomalous) terms.
C...For F2, c and b are included by the Bethe-Heitler formula;
C...in the 'MSbar' scheme additionally a Cgamma term is added.
      dimension xpdfgm(-6:6)
      common/pyint8/xpvmd(-6:6),xpanl(-6:6),xpanh(-6:6),xpbeh(-6:6),
     &xpdir(-6:6)
      SAVE /pyint8/
 
C...Temporary array.
      dimension xpga(-6:6)
C...Charm and bottom masses (low to compensate for J/psi etc.).
      DATA pmc/1.3/, pmb/4.6/
C...alpha_em and alpha_em/(2*pi).
      DATA aem/0.007297/, aem2pi/0.0011614/
C...Lambda value for 4 flavours.
      DATA alam/0.20/
C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
      DATA fracu/0.8/
C...VMD couplings f_V**2/(4*pi).
      DATA frho/2.20/, fomega/23.6/, fphi/18.4/
C...Masses for rho (=omega) and phi.
      DATA pmrho/0.770/, pmphi/1.020/
 
C...Reset output.
      f2gm=0.
      DO 100 kfl=-6,6
      xpdfgm(kfl)=0.
      xpvmd(kfl)=0.
      xpanl(kfl)=0.
      xpanh(kfl)=0.
      xpbeh(kfl)=0.
      xpdir(kfl)=0.
  100 CONTINUE
 
C...Set k0 cut-off parameter as function of set used.
      IF(iset.LE.2) THEN
        ak0=0.6
      ELSE
        ak0=2.
      ENDIF
 
C...Call VMD parametrization for d quark and use to give rho+omega+ phi.
C...Note scale choice and dipole dampening for off-shell photon.
      p2mx=max(p2,ak0**2)
      CALL pygvmd(iset,1,x,q2,p2mx,alam,xpga)
      xfval=xpga(1)-xpga(2)
      xpga(1)=xpga(2)
      xpga(-1)=xpga(-2)
      facud=aem*(1./frho+1./fomega)*(pmrho**2/(pmrho**2+p2))**2
      facs=aem*(1./fphi)*(pmphi**2/(pmphi**2+p2))**2
      DO 110 kfl=-5,5
      xpvmd(kfl)=(facud+facs)*xpga(kfl)
  110 CONTINUE
      xpvmd(1)=xpvmd(1)+(1.-fracu)*facud*xfval
      xpvmd(2)=xpvmd(2)+fracu*facud*xfval
      xpvmd(3)=xpvmd(3)+facs*xfval
      xpvmd(-1)=xpvmd(-1)+(1.-fracu)*facud*xfval
      xpvmd(-2)=xpvmd(-2)+fracu*facud*xfval
      xpvmd(-3)=xpvmd(-3)+facs*xfval
 
C...Call anomalous parametrization for d + u + s.
      CALL pygano(-3,x,q2,p2mx,alam,xpga)
      DO 120 kfl=-5,5
      xpanl(kfl)=xpga(kfl)
  120 CONTINUE
 
C...Call anomalous parametrization for c and b.
      CALL pygano(4,x,q2,p2mx,alam,xpga)
      DO 130 kfl=-5,5
      xpanh(kfl)=xpga(kfl)
  130 CONTINUE
      CALL pygano(5,x,q2,p2mx,alam,xpga)
      DO 140 kfl=-5,5
      xpanh(kfl)=xpanh(kfl)+xpga(kfl)
  140 CONTINUE
 
C...Call Bethe-Heitler term expression for charm and bottom.
      CALL pygbeh(4,x,q2,p2,pmc**2,xpbh)
      xpbeh(4)=xpbh
      xpbeh(-4)=xpbh
      CALL pygbeh(5,x,q2,p2,pmb**2,xpbh)
      xpbeh(5)=xpbh
      xpbeh(-5)=xpbh
 
C...For MSbar subtraction call C^gamma term expression for d, u, s.
      IF(iset.EQ.2.OR.iset.EQ.4) THEN
        CALL pygdir(x,q2,p2,ak0,xpga)
        DO 150 kfl=-5,5
        xpdir(kfl)=xpga(kfl)
  150   CONTINUE
      ENDIF
 
C...Store result in output array.
      DO 160 kfl=-5,5
      chsq=1./9.
      IF(iabs(kfl).EQ.2.OR.iabs(kfl).EQ.4) chsq=4./9.
      xpf2=xpvmd(kfl)+xpanl(kfl)+xpbeh(kfl)+xpdir(kfl)
      IF(kfl.NE.0) f2gm=f2gm+chsq*xpf2
      xpdfgm(kfl)=xpvmd(kfl)+xpanl(kfl)+xpanh(kfl)
  160 CONTINUE
 
      RETURN
      END