pypdfl.f

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

 
C*********************************************************************
 
C...PYPDFL
C...Gives proton parton distribution at small x and/or Q^2 according to
C...correct limiting behaviour.
 
      SUBROUTINE pypdfl(KF,X,Q2,XPQ)
 
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)
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pyint1/mint(400),vint(400)
      SAVE /pydat1/,/pydat2/,/pypars/,/pyint1/
C...Local arrays.
      dimension xpq(-25:25),xpa(-25:25),xpb(-25:25),wtsb(-3:3)
      DATA rmr/0.92d0/,rmp/0.38d0/,wtsb/0.5d0,1d0,1d0,5d0,1d0,1d0,0.5d0/
 
C...Send everything but protons/neutrons/VMD pions directly to PYPDFU.
      mint(92)=0
      kfa=iabs(kf)
      iacc=0
      IF((kfa.EQ.2212.OR.kfa.EQ.2112).AND.mstp(57).GE.2) iacc=1
      IF(kfa.EQ.211.AND.mstp(57).GE.3) iacc=1
      IF(kfa.EQ.22.AND.mint(109).EQ.2.AND.mstp(57).GE.3) iacc=1
      IF(iacc.EQ.0) THEN
        CALL pypdfu(kf,x,q2,xpq)
        RETURN
      ENDIF
 
C...Reset. Check x.
      DO 100 kfl=-25,25
        xpq(kfl)=0d0
  100 CONTINUE
      IF(x.LE.0d0.OR.x.GE.1d0) THEN
        WRITE(mstu(11),5000) x
        RETURN
      ENDIF
 
C...Define valence content.
      kfc=kf
      nv1=2
      nv2=1
      IF(kf.EQ.2212) THEN
        kfv1=2
        kfv2=1
      ELSEIF(kf.EQ.-2212) THEN
        kfv1=-2
        kfv2=-1
      ELSEIF(kf.EQ.2112) THEN
        kfv1=1
        kfv2=2
      ELSEIF(kf.EQ.-2112) THEN
        kfv1=-1
        kfv2=-2
      ELSEIF(kf.EQ.211) THEN
        nv1=1
        kfv1=2
        kfv2=-1
      ELSEIF(kf.EQ.-211) THEN
        nv1=1
        kfv1=-2
        kfv2=1
      ELSEIF(mint(105).LE.223) THEN
        kfv1=1
        wtv1=0.2d0
        kfv2=2
        wtv2=0.8d0
      ELSEIF(mint(105).EQ.333) THEN
        kfv1=3
        wtv1=1.0d0
        kfv2=1
        wtv2=0.0d0
      ELSEIF(mint(105).EQ.443) THEN
        kfv1=4
        wtv1=1.0d0
        kfv2=1
        wtv2=0.0d0
      ENDIF
 
C...Do naive evaluation and find min Q^2, boundary Q^2 and x_0.
      mint30=mint(30)
      CALL pypdfu(kfc,x,q2,xpa)
      q2mn=max(3d0,vint(231))
      q2b=2d0+0.052d0**2*exp(3.56d0*sqrt(max(0d0,-log(3d0*x))))
      xmn=exp(-(log((q2mn-2d0)/0.052d0**2)/3.56d0)**2)/3d0
 
C...Large Q2 and large x: naive call is enough.
      IF(q2.GT.q2mn.AND.q2.GT.q2b) THEN
        DO 110 kfl=-25,25
          xpq(kfl)=xpa(kfl)
  110   CONTINUE
        mint(92)=1
 
C...Small Q2 and large x: dampen boundary value.
      ELSEIF(x.GT.xmn) THEN
 
C...Evaluate at boundary and define dampening factors.
        mint(30)=mint30
        CALL pypdfu(kfc,x,q2mn,xpa)
        fv=(q2*(q2mn+rmr)/(q2mn*(q2+rmr)))**(0.55d0*(1d0-x)/(1d0-xmn))
        fs=(q2*(q2mn+rmp)/(q2mn*(q2+rmp)))**1.08d0
 
C...Separate valence and sea parts of parton distribution.
        IF(kfa.NE.22) THEN
          xfv1=xpa(kfv1)-xpa(-kfv1)
          xpa(kfv1)=xpa(-kfv1)
          xfv2=xpa(kfv2)-xpa(-kfv2)
          xpa(kfv2)=xpa(-kfv2)
        ELSE
          xpa(kfv1)=xpa(kfv1)-wtv1*vint(232)
          xpa(-kfv1)=xpa(-kfv1)-wtv1*vint(232)
          xpa(kfv2)=xpa(kfv2)-wtv2*vint(232)
          xpa(-kfv2)=xpa(-kfv2)-wtv2*vint(232)
        ENDIF
 
C...Dampen valence and sea separately. Put back together.
        DO 120 kfl=-25,25
          xpq(kfl)=fs*xpa(kfl)
  120   CONTINUE
        IF(kfa.NE.22) THEN
          xpq(kfv1)=xpq(kfv1)+fv*xfv1
          xpq(kfv2)=xpq(kfv2)+fv*xfv2
        ELSE
          xpq(kfv1)=xpq(kfv1)+fv*wtv1*vint(232)
          xpq(-kfv1)=xpq(-kfv1)+fv*wtv1*vint(232)
          xpq(kfv2)=xpq(kfv2)+fv*wtv2*vint(232)
          xpq(-kfv2)=xpq(-kfv2)+fv*wtv2*vint(232)
        ENDIF
        mint(92)=2
 
C...Large Q2 and small x: interpolate behaviour.
      ELSEIF(q2.GT.q2mn) THEN
 
C...Evaluate at extremes and define coefficients for interpolation.
        mint(30)=mint30
        CALL pypdfu(kfc,xmn,q2mn,xpa)
        vi232a=vint(232)
        mint(30)=mint30
        CALL pypdfu(kfc,x,q2b,xpb)
        vi232b=vint(232)
        fla=log(q2b/q2)/log(q2b/q2mn)
        fva=(x/xmn)**0.45d0*fla
        fsa=(x/xmn)**(-0.08d0)*fla
        fb=1d0-fla
 
C...Separate valence and sea parts of parton distribution.
        IF(kfa.NE.22) THEN
          xfva1=xpa(kfv1)-xpa(-kfv1)
          xpa(kfv1)=xpa(-kfv1)
          xfva2=xpa(kfv2)-xpa(-kfv2)
          xpa(kfv2)=xpa(-kfv2)
          xfvb1=xpb(kfv1)-xpb(-kfv1)
          xpb(kfv1)=xpb(-kfv1)
          xfvb2=xpb(kfv2)-xpb(-kfv2)
          xpb(kfv2)=xpb(-kfv2)
        ELSE
          xpa(kfv1)=xpa(kfv1)-wtv1*vi232a
          xpa(-kfv1)=xpa(-kfv1)-wtv1*vi232a
          xpa(kfv2)=xpa(kfv2)-wtv2*vi232a
          xpa(-kfv2)=xpa(-kfv2)-wtv2*vi232a
          xpb(kfv1)=xpb(kfv1)-wtv1*vi232b
          xpb(-kfv1)=xpb(-kfv1)-wtv1*vi232b
          xpb(kfv2)=xpb(kfv2)-wtv2*vi232b
          xpb(-kfv2)=xpb(-kfv2)-wtv2*vi232b
        ENDIF
 
C...Interpolate for valence and sea. Put back together.
        DO 130 kfl=-25,25
          xpq(kfl)=fsa*xpa(kfl)+fb*xpb(kfl)
  130   CONTINUE
        IF(kfa.NE.22) THEN
          xpq(kfv1)=xpq(kfv1)+(fva*xfva1+fb*xfvb1)
          xpq(kfv2)=xpq(kfv2)+(fva*xfva2+fb*xfvb2)
        ELSE
          xpq(kfv1)=xpq(kfv1)+wtv1*(fva*vi232a+fb*vi232b)
          xpq(-kfv1)=xpq(-kfv1)+wtv1*(fva*vi232a+fb*vi232b)
          xpq(kfv2)=xpq(kfv2)+wtv2*(fva*vi232a+fb*vi232b)
          xpq(-kfv2)=xpq(-kfv2)+wtv2*(fva*vi232a+fb*vi232b)
        ENDIF
        mint(92)=3
 
C...Small Q2 and small x: dampen boundary value and add term.
      ELSE
 
C...Evaluate at boundary and define dampening factors.
        mint(30)=mint30
        CALL pypdfu(kfc,xmn,q2mn,xpa)
        fb=(xmn-x)*(q2mn-q2)/(xmn*q2mn)
        fa=1d0-fb
        fvc=(x/xmn)**0.45d0*(q2/(q2+rmr))**0.55d0
        fva=fvc*fa*((q2mn+rmr)/q2mn)**0.55d0
        fvb=fvc*fb*1.10d0*xmn**0.45d0*0.11d0
        fsc=(x/xmn)**(-0.08d0)*(q2/(q2+rmp))**1.08d0
        fsa=fsc*fa*((q2mn+rmp)/q2mn)**1.08d0
        fsb=fsc*fb*0.21d0*xmn**(-0.08d0)*0.21d0
 
C...Separate valence and sea parts of parton distribution.
        IF(kfa.NE.22) THEN
          xfv1=xpa(kfv1)-xpa(-kfv1)
          xpa(kfv1)=xpa(-kfv1)
          xfv2=xpa(kfv2)-xpa(-kfv2)
          xpa(kfv2)=xpa(-kfv2)
        ELSE
          xpa(kfv1)=xpa(kfv1)-wtv1*vint(232)
          xpa(-kfv1)=xpa(-kfv1)-wtv1*vint(232)
          xpa(kfv2)=xpa(kfv2)-wtv2*vint(232)
          xpa(-kfv2)=xpa(-kfv2)-wtv2*vint(232)
        ENDIF
 
C...Dampen valence and sea separately. Add constant terms.
C...Put back together.
        DO 140 kfl=-25,25
          xpq(kfl)=fsa*xpa(kfl)
  140   CONTINUE
        IF(kfa.NE.22) THEN
          DO 150 kfl=-3,3
            xpq(kfl)=xpq(kfl)+fsb*wtsb(kfl)
  150     CONTINUE
          xpq(kfv1)=xpq(kfv1)+(fva*xfv1+fvb*nv1)
          xpq(kfv2)=xpq(kfv2)+(fva*xfv2+fvb*nv2)
        ELSE
          DO 160 kfl=-3,3
            xpq(kfl)=xpq(kfl)+vint(281)*fsb*wtsb(kfl)
  160     CONTINUE
          xpq(kfv1)=xpq(kfv1)+wtv1*(fva*vint(232)+fvb*vint(281))
          xpq(-kfv1)=xpq(-kfv1)+wtv1*(fva*vint(232)+fvb*vint(281))
          xpq(kfv2)=xpq(kfv2)+wtv2*(fva*vint(232)+fvb*vint(281))
          xpq(-kfv2)=xpq(-kfv2)+wtv2*(fva*vint(232)+fvb*vint(281))
        ENDIF
        xpq(21)=xpq(0)
        mint(92)=4
      ENDIF
 
C...Format for error printout.
 5000 FORMAT(' Error: x value outside physical range; x =',1p,d12.3)
 
      RETURN
      END