pypdel.f

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

 
C*********************************************************************
 
C...PYPDEL
C...Gives electron (or muon, or tau) parton distribution.
 
      SUBROUTINE pypdel(KFA,X,Q2,XPEL)
 
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 xpel(-25:25),xpga(-6:6),sxp(0:6)
 
C...Interface to PDFLIB.
      common/w50513/xmin,xmax,q2min,q2max
      SAVE /w50513/
      DOUBLE PRECISION xx,qq,upv,dnv,usea,dsea,str,chm,bot,top,glu,
     &value(20),xmin,xmax,q2min,q2max
      CHARACTER*20 parm(20)
      DATA value/20*0d0/,parm/20*' '/
 
C...Some common constants.
      DO 100 kfl=-25,25
        xpel(kfl)=0d0
  100 CONTINUE
      aem=paru(101)
      pme=pmas(11,1)
      IF(kfa.EQ.13) pme=pmas(13,1)
      IF(kfa.EQ.15) pme=pmas(15,1)
      xl=log(max(1d-10,x))
      x1l=log(max(1d-10,1d0-x))
      hle=log(max(3d0,q2/pme**2))
      hbe2=(aem/paru(1))*(hle-1d0)
 
C...Electron inside electron, see R. Kleiss et al., in Z physics at
C...LEP 1, CERN 89-08, p. 34
      IF(mstp(59).LE.1) THEN
        hde=1d0+(aem/paru(1))*(1.5d0*hle+1.289868d0)+(aem/paru(1))**2*
     &  (-2.164868d0*hle**2+9.840808d0*hle-10.130464d0)
        hee=hbe2*(1d0-x)**(hbe2-1d0)*sqrt(max(0d0,hde))-
     &  0.5d0*hbe2*(1d0+x)+hbe2**2/8d0*((1d0+x)*(-4d0*x1l+3d0*xl)-
     &  4d0*xl/(1d0-x)-5d0-x)
      ELSE
        hee=hbe2*(1d0-x)**(hbe2-1d0)*exp(0.172784d0*hbe2)/
     &  pygamm(1d0+hbe2)-0.5d0*hbe2*(1d0+x)+hbe2**2/8d0*((1d0+x)*
     &  (-4d0*x1l+3d0*xl)-4d0*xl/(1d0-x)-5d0-x)
      ENDIF
C...Zero distribution for very large x and rescale it for intermediate.
      IF(x.GT.1d0-1d-10) THEN
        hee=0d0
      ELSEIF(x.GT.1d0-1d-7) THEN
        hee=hee*1000d0**hbe2/(1000d0**hbe2-1d0)
      ENDIF
      xpel(kfa)=x*hee
 
C...Photon and (transverse) W- inside electron.
      aemp=pyalem(pme*sqrt(max(0d0,q2)))/paru(2)
      IF(mstp(13).LE.1) THEN
        hlg=hle
      ELSE
        hlg=log(max(1d0,(parp(13)/pme**2)*(1d0-x)/x**2))
      ENDIF
      xpel(22)=aemp*hlg*(1d0+(1d0-x)**2)
      hlw=log(1d0+q2/pmas(24,1)**2)/(4d0*paru(102))
      xpel(-24)=aemp*hlw*(1d0+(1d0-x)**2)
 
C...Electron or positron inside photon inside electron.
      IF(kfa.EQ.11.AND.mstp(12).EQ.1) THEN
        xfsea=0.5d0*(aemp*(hle-1d0))**2*(4d0/3d0+x-x**2-4d0*x**3/3d0+
     &  2d0*x*(1d0+x)*xl)
        xpel(11)=xpel(11)+xfsea
        xpel(-11)=xfsea
 
C...Initialize PDFLIB photon parton distributions.
        IF(mstp(56).EQ.2) THEN
          parm(1)='NPTYPE'
          value(1)=3
          parm(2)='NGROUP'
          value(2)=mstp(55)/1000
          parm(3)='NSET'
          value(3)=mod(mstp(55),1000)
          IF(mint(93).NE.3000000+mstp(55)) THEN
C             CALL PDFSET(PARM,VALUE)
            mint(93)=3000000+mstp(55)
          ENDIF
        ENDIF
 
C...Quarks and gluons inside photon inside electron:
C...numerical convolution required.
        DO 110 kfl=0,6
          sxp(kfl)=0d0
  110   CONTINUE
        sumxpp=0d0
        iter=-1
  120   iter=iter+1
        sumxp=sumxpp
        nstp=2**(iter-1)
        IF(iter.EQ.0) nstp=2
        DO 130 kfl=0,6
          sxp(kfl)=0.5d0*sxp(kfl)
  130   CONTINUE
        wtstp=0.5d0/nstp
        IF(iter.EQ.0) wtstp=0.5d0
C...Pick grid of x_{gamma} values logarithmically even.
        DO 150 istp=1,nstp
          IF(iter.EQ.0) THEN
            xle=xl*(istp-1)
          ELSE
            xle=xl*(istp-0.5d0)/nstp
          ENDIF
          xe=min(1d0-1d-10,exp(xle))
          xg=min(1d0-1d-10,x/xe)
C...Evaluate photon inside electron parton distribution for convolution.
          xpgp=1d0+(1d0-xe)**2
          IF(mstp(13).LE.1) THEN
            xpgp=xpgp*hle
          ELSE
            xpgp=xpgp*log(max(1d0,(parp(13)/pme**2)*(1d0-xe)/xe**2))
          ENDIF
C...Evaluate photon parton distributions for convolution.
          IF(mstp(56).EQ.1) THEN
            IF(mstp(55).EQ.1) THEN
              CALL pypdga(xg,q2,xpga)
            ELSEIF(mstp(55).GE.5.AND.mstp(55).LE.8) THEN
              q2mx=q2
              p2mx=0.36d0
              IF(mstp(55).GE.7) p2mx=4.0d0
              IF(mstp(57).EQ.0) q2mx=p2mx
              p2=0d0
              IF(vint(120).LT.0d0) p2=vint(120)**2
              CALL pyggam(mstp(55)-4,xg,q2mx,p2,mstp(60),f2gam,xpga)
              vint(231)=p2mx
            ELSEIF(mstp(55).GE.9.AND.mstp(55).LE.12) THEN
              q2mx=q2
              p2mx=0.36d0
              IF(mstp(55).GE.11) p2mx=4.0d0
              IF(mstp(57).EQ.0) q2mx=p2mx
              p2=0d0
              IF(vint(120).LT.0d0) p2=vint(120)**2
              CALL pyggam(mstp(55)-8,xg,q2mx,p2,mstp(60),f2gam,xpga)
              vint(231)=p2mx
            ENDIF
            DO 140 kfl=0,5
              sxp(kfl)=sxp(kfl)+wtstp*xpgp*xpga(kfl)
  140       CONTINUE
          ELSEIF(mstp(56).EQ.2) THEN
C...Call PDFLIB parton distributions.
            xx=xg
            qq=sqrt(max(0d0,q2min,q2))
            IF(mstp(57).EQ.0) qq=sqrt(q2min)
C            CALL STRUCTM(XX,QQ,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU)
            sxp(0)=sxp(0)+wtstp*xpgp*glu
            sxp(1)=sxp(1)+wtstp*xpgp*dnv
            sxp(2)=sxp(2)+wtstp*xpgp*upv
            sxp(3)=sxp(3)+wtstp*xpgp*str
            sxp(4)=sxp(4)+wtstp*xpgp*chm
            sxp(5)=sxp(5)+wtstp*xpgp*bot
            sxp(6)=sxp(6)+wtstp*xpgp*top
          ENDIF
  150   CONTINUE
        sumxpp=sxp(0)+2d0*sxp(1)+2d0*sxp(2)
        IF(iter.LE.2.OR.(iter.LE.7.AND.abs(sumxpp-sumxp).GT.
     &  parp(14)*(sumxpp+sumxp))) goto 120
 
C...Put convolution into output arrays.
        fconv=aemp*(-xl)
        xpel(0)=fconv*sxp(0)
        DO 160 kfl=1,6
          xpel(kfl)=fconv*sxp(kfl)
          xpel(-kfl)=xpel(kfl)
  160   CONTINUE
      ENDIF
 
      RETURN
      END