pystel.F

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

      SUBROUTINE pystel(X,Q2,XPEL)
 
C...Gives electron structure function.
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
      common/ludat2/kchg(500,3),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 /ludat1/,/ludat2/
      SAVE /pypars/,/pyint1/
      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)=0.
  100 CONTINUE
      aem=paru(101)
      pme=pmas(11,1)
      xl=log(max(1e-10,x))
      x1l=log(max(1e-10,1.-x))
      hle=log(max(3.,q2/pme**2))
      hbe2=(aem/paru(1))*(hle-1.)
 
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=1.+(aem/paru(1))*(1.5*hle+1.289868)+(aem/paru(1))**2*
     &  (-2.164868*hle**2+9.840808*hle-10.130464)
        hee=hbe2*(1.-x)**(hbe2-1.)*sqrt(max(0.,hde))-
     &  0.5*hbe2*(1.+x)+hbe2**2/8.*((1.+x)*(-4.*x1l+3.*xl)-
     &  4.*xl/(1.-x)-5.-x)
      ELSE
        hee=hbe2*(1.-x)**(hbe2-1.)*exp(0.172784*hbe2)/pygamm(1.+hbe2)-
     &  0.5*hbe2*(1.+x)+hbe2**2/8.*((1.+x)*(-4.*x1l+3.*xl)-
     &  4.*xl/(1.-x)-5.-x)        
      ENDIF
      IF(x.GT.0.9999.AND.x.LE.0.999999) THEN
        hee=hee*100.**hbe2/(100.**hbe2-1.)
      ELSEIF(x.GT.0.999999) THEN
        hee=0.
      ENDIF
      xpel(11)=x*hee
 
C...Photon and (transverse) W- inside electron.
      aemp=ulalem(pme*sqrt(max(0.,q2)))/paru(2)
      IF(mstp(13).LE.1) THEN
        hlg=hle
      ELSE
        hlg=log(max(1.,(parp(13)/pme**2)*(1.-x)/x**2))
      ENDIF
      xpel(22)=aemp*hlg*(1.+(1.-x)**2)
      hlw=log(1.+q2/pmas(24,1)**2)/(4.*paru(102))
      xpel(-24)=aemp*hlw*(1.+(1.-x)**2)
 
C...Electron or positron inside photon inside electron.
      IF(mstp(12).EQ.1) THEN
        xfsea=0.5*(aemp*(hle-1.))**2*(4./3.+x-x**2-4.*x**3/3.+
     &  2.*x*(1.+x)*xl)
        xpel(11)=xpel(11)+xfsea
        xpel(-11)=xfsea
 
C...Initialize PDFLIB photon structure functions.
        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
            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)=0.
  110   CONTINUE
        sumxpp=0.
        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.5*sxp(kfl)
  130   CONTINUE
        wtstp=0.5/nstp
        IF(iter.EQ.0) wtstp=0.5
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.5)/nstp
        ENDIF
        xe=min(0.999999,exp(xle))
        xg=min(0.999999,x/xe)
C...Evaluate photon inside electron structure function for convolution.
        xpgp=1.+(1.-xe)**2
        IF(mstp(13).LE.1) THEN
          xpgp=xpgp*hle
        ELSE
          xpgp=xpgp*log(max(1.,(parp(13)/pme**2)*(1.-xe)/xe**2))
        ENDIF
C...Evaluate photon structure functions for convolution.
        IF(mstp(56).EQ.1) THEN
          CALL pystga(xg,q2,xpga)
          DO 140 kfl=0,5
          sxp(kfl)=sxp(kfl)+wtstp*xpgp*xpga(kfl)
  140     CONTINUE
        ELSEIF(mstp(56).EQ.2) THEN
C...Call PDFLIB structure functions.
          xx=xg
          qq=sqrt(max(0.,sngl(q2min),q2))
          IF(mstp(57).EQ.0) qq=sqrt(q2min)
          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)+2.*sxp(1)+2.*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