pyggam.f

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

 
C*********************************************************************
 
C...PYGGAM
C...Constructs 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.
C...Contains the SaS sets 1D, 1M, 2D and 2M.
C...Adapted from SaSgam library, authors G.A. Schuler and T. Sjostrand.
 
      SUBROUTINE pyggam(ISET,X,Q2,P2,IP2,F2GM,XPDFGM)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pyint8/xpvmd(-6:6),xpanl(-6:6),xpanh(-6:6),xpbeh(-6:6),
     &xpdir(-6:6)
      common/pyint9/vxpvmd(-6:6),vxpanl(-6:6),vxpanh(-6:6),vxpdgm(-6:6)
      SAVE /pyint8/,/pyint9/
C...Local arrays.
      dimension xpdfgm(-6:6),xpga(-6:6), vxpga(-6:6)
C...Charm and bottom masses (low to compensate for J/psi etc.).
      DATA pmc/1.3d0/, pmb/4.6d0/
C...alpha_em and alpha_em/(2*pi).
      DATA aem/0.007297d0/, aem2pi/0.0011614d0/
C...Lambda value for 4 flavours.
      DATA alam/0.20d0/
C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
      DATA fracu/0.8d0/
C...VMD couplings f_V**2/(4*pi).
      DATA frho/2.20d0/, fomega/23.6d0/, fphi/18.4d0/
C...Masses for rho (=omega) and phi.
      DATA pmrho/0.770d0/, pmphi/1.020d0/
C...Number of points in integration for IP2=1.
      DATA nstep/100/
 
C...Reset output.
      f2gm=0d0
      DO 100 kfl=-6,6
        xpdfgm(kfl)=0d0
        xpvmd(kfl)=0d0
        xpanl(kfl)=0d0
        xpanh(kfl)=0d0
        xpbeh(kfl)=0d0
        xpdir(kfl)=0d0
        vxpvmd(kfl)=0d0
        vxpanl(kfl)=0d0
        vxpanh(kfl)=0d0
        vxpdgm(kfl)=0d0
  100 CONTINUE
 
C...Set Q0 cut-off parameter as function of set used.
      IF(iset.LE.2) THEN
        q0=0.6d0
      ELSE
        q0=2d0
      ENDIF
      q02=q0**2
 
C...Scale choice for off-shell photon; common factors.
      q2a=q2
      facnor=1d0
      IF(ip2.EQ.1) THEN
        p2mx=p2+q02
        q2a=q2+p2*q02/max(q02,q2)
        facnor=log(q2/q02)/nstep
      ELSEIF(ip2.EQ.2) THEN
        p2mx=max(p2,q02)
      ELSEIF(ip2.EQ.3) THEN
        p2mx=p2+q02
        q2a=q2+p2*q02/max(q02,q2)
      ELSEIF(ip2.EQ.4) THEN
        p2mx=q2*(q02+p2)/(q2+p2)*exp(p2*(q2-q02)/
     &  ((q2+p2)*(q02+p2)))
      ELSEIF(ip2.EQ.5) THEN
        p2mxa=q2*(q02+p2)/(q2+p2)*exp(p2*(q2-q02)/
     &  ((q2+p2)*(q02+p2)))
        p2mx=q0*sqrt(p2mxa)
        facnor=log(q2/p2mxa)/log(q2/p2mx)
      ELSEIF(ip2.EQ.6) THEN
        p2mx=q2*(q02+p2)/(q2+p2)*exp(p2*(q2-q02)/
     &  ((q2+p2)*(q02+p2)))
        p2mx=max(0d0,1d0-p2/q2)*p2mx+min(1d0,p2/q2)*max(p2,q02)
      ELSE
        p2mxa=q2*(q02+p2)/(q2+p2)*exp(p2*(q2-q02)/
     &  ((q2+p2)*(q02+p2)))
        p2mx=q0*sqrt(p2mxa)
        p2mxb=p2mx
        p2mx=max(0d0,1d0-p2/q2)*p2mx+min(1d0,p2/q2)*max(p2,q02)
        p2mxb=max(0d0,1d0-p2/q2)*p2mxb+min(1d0,p2/q2)*p2mxa
        IF(abs(q2-q02).GT.1d-6) THEN
          facnor=log(q2/p2mxa)/log(q2/p2mxb)
        ELSEIF(p2.LT.q02) THEN
          facnor=q02**3/(q02+p2)/(q02**2-p2**2/2d0)
        ELSE
          facnor=1d0
        ENDIF
      ENDIF
 
C...Call VMD parametrization for d quark and use to give rho, omega,
C...phi. Note dipole dampening for off-shell photon.
      CALL pygvmd(iset,1,x,q2a,p2mx,alam,xpga,vxpga)
      xfval=vxpga(1)
      xpga(1)=xpga(2)
      xpga(-1)=xpga(-2)
      facud=aem*(1d0/frho+1d0/fomega)*(pmrho**2/(pmrho**2+p2))**2
      facs=aem*(1d0/fphi)*(pmphi**2/(pmphi**2+p2))**2
      DO 110 kfl=-5,5
        xpvmd(kfl)=(facud+facs)*xpga(kfl)
  110 CONTINUE
      xpvmd(1)=xpvmd(1)+(1d0-fracu)*facud*xfval
      xpvmd(2)=xpvmd(2)+fracu*facud*xfval
      xpvmd(3)=xpvmd(3)+facs*xfval
      xpvmd(-1)=xpvmd(-1)+(1d0-fracu)*facud*xfval
      xpvmd(-2)=xpvmd(-2)+fracu*facud*xfval
      xpvmd(-3)=xpvmd(-3)+facs*xfval
      vxpvmd(1)=(1d0-fracu)*facud*xfval
      vxpvmd(2)=fracu*facud*xfval
      vxpvmd(3)=facs*xfval
      vxpvmd(-1)=(1d0-fracu)*facud*xfval
      vxpvmd(-2)=fracu*facud*xfval
      vxpvmd(-3)=facs*xfval
 
      IF(ip2.NE.1) THEN
C...Anomalous parametrizations for different strategies
C...for off-shell photons; except full integration.
 
C...Call anomalous parametrization for d + u + s.
        CALL pygano(-3,x,q2a,p2mx,alam,xpga,vxpga)
        DO 120 kfl=-5,5
          xpanl(kfl)=facnor*xpga(kfl)
          vxpanl(kfl)=facnor*vxpga(kfl)
  120   CONTINUE
 
C...Call anomalous parametrization for c and b.
        CALL pygano(4,x,q2a,p2mx,alam,xpga,vxpga)
        DO 130 kfl=-5,5
          xpanh(kfl)=facnor*xpga(kfl)
          vxpanh(kfl)=facnor*vxpga(kfl)
  130   CONTINUE
        CALL pygano(5,x,q2a,p2mx,alam,xpga,vxpga)
        DO 140 kfl=-5,5
          xpanh(kfl)=xpanh(kfl)+facnor*xpga(kfl)
          vxpanh(kfl)=vxpanh(kfl)+facnor*vxpga(kfl)
  140   CONTINUE
 
      ELSE
C...Special option: loop over flavours and integrate over k2.
        DO 170 kf=1,5
          DO 160 istep=1,nstep
            q2step=q02*(q2/q02)**((istep-0.5d0)/nstep)
            IF((kf.EQ.4.AND.q2step.LT.pmc**2).OR.
     &      (kf.EQ.5.AND.q2step.LT.pmb**2)) goto 160
            CALL pygvmd(0,kf,x,q2,q2step,alam,xpga,vxpga)
            facq=aem2pi*(q2step/(q2step+p2))**2*facnor
            IF(mod(kf,2).EQ.0) facq=facq*(8d0/9d0)
            IF(mod(kf,2).EQ.1) facq=facq*(2d0/9d0)
            DO 150 kfl=-5,5
              IF(kf.LE.3) xpanl(kfl)=xpanl(kfl)+facq*xpga(kfl)
              IF(kf.GE.4) xpanh(kfl)=xpanh(kfl)+facq*xpga(kfl)
              IF(kf.LE.3) vxpanl(kfl)=vxpanl(kfl)+facq*vxpga(kfl)
              IF(kf.GE.4) vxpanh(kfl)=vxpanh(kfl)+facq*vxpga(kfl)
  150       CONTINUE
  160     CONTINUE
  170   CONTINUE
      ENDIF
 
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,q02,xpga)
        DO 180 kfl=-5,5
          xpdir(kfl)=xpga(kfl)
  180   CONTINUE
      ENDIF
 
C...Store result in output array.
      DO 190 kfl=-5,5
        chsq=1d0/9d0
        IF(iabs(kfl).EQ.2.OR.iabs(kfl).EQ.4) chsq=4d0/9d0
        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)
        vxpdgm(kfl)=vxpvmd(kfl)+vxpanl(kfl)+vxpanh(kfl)
  190 CONTINUE
 
      RETURN
      END