lystfu.F

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

C **********************************************************************

      SUBROUTINE lystfu(KF,X,Q2,XPQ)

      IMPLICIT NONE

C...Interface to PYSTFU in PYTHIA 5.7 to get parton density distributions,
C...i.e. momentum weighted probability distributions xq(x,Q2), xg(x,Q2).
C...Also gives intrinsic charm and beauty distributions.
*
* to avoid variable conflictions, a second keep element is necessary
* with the same common block name (see LEPTO2)
*
       include "lepto2.inc"
       include "ludat1.inc"
       include "pypars.inc"
       include "arstrf.inc"

c$$$      COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),
c$$$     &                XLP,YLP,W2LP,Q2LP,ULP
c$$$      REAL CUT,PARL,XLP,YLP,W2LP,Q2LP,ULP
c$$$      INTEGER LST
c$$$      SAVE /LEPTOU/
c$$$
c$$$      COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
c$$$      INTEGER  MSTU,MSTJ
c$$$      REAL  PARU,PARJ
c$$$      SAVE /LUDAT1/
c$$$
c$$$      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) 
c$$$      INTEGER MSTP,MSTI
c$$$      REAL PARP,PARI
c$$$      SAVE /PYPARS/
c$$$
c$$$      COMMON /ARSTRF/ KFSAVE(2),XSAVE(2),XQ2SAV(2),
c$$$     +    XPQSAV(2,-6:6),XDPQSAV(2,-6:6)
c$$$      INTEGER KFSAVE
c$$$      REAL XSAVE,XQ2SAV,XPQSAV,XDPQSAV
c$$$      SAVE /ARSTRF/


      REAL xcorr,xpyst,x,q2,xpq,xdpq
      INTEGER kfl,kf
      dimension xdpq(-6:6),xpq(-6:6),xpyst(-25:25)

C...Reset arrays etc.
      DO 100 kfl=-6,6
      xpq(kfl)=0.0
  100 xpqsav(1,kfl)=0.
      xsave(1)=x
      xq2sav(1)=q2
      kfsave(1)=kf
C...Check x and particle species.
      IF(x.LE.0..OR.x.GE.1.) THEN
C        WRITE(MSTU(11),5000) X
         WRITE(*,5000) x
        RETURN
      ENDIF

      IF(lst(15).EQ.-4.OR.lst(15).EQ.-5) THEN
C...Intrinsic charm/bottom quark distribution in the proton...
        IF(q2.LT.1.) RETURN
C...from Phys. Lett 93B (1980) 451, 
C...Amount of intrinsic charm PARL(12)=BETA^2
        xpq(4)=x**3*1800.*parl(12)*
     &         ((1.-x)/3.*(1.+10.*x+x**2)+2.*x*(1.+x)*log(x))
C...plus first order QCD-correction parametrized with polynomia
        IF(x.LT.0.9) THEN
          xcorr=0.22024e-1*x-0.77833e-1*x**2-0.47292*x**3+
     &          2.104*x**4-2.1698*x**5-0.84891*x**6+1.8882*x**7+
     &          0.8989*x**8-2.1072*x**9+0.76351*x**10
        ELSE
          xcorr=-1.
        ENDIF
C...and a Q2 dependence on that
CJR        XCORR=1.125*XCORR*0.190424*EXP(1.15*LOG(LOG(Q2)))
        IF(q2.GT.1) THEN
           xcorr=1.125*xcorr*0.190424*exp(1.15*log(log(q2)))
        ELSE
           xcorr=1.125*xcorr*0.190424
        ENDIF
C...smooth cut-off of the structure function !
        xpq(4)=max(xpq(4)+xcorr,xpq(4)/q2)
        xpq(-4)=xpq(4)
        IF(lst(15).EQ.-5) THEN
C...  Intrinsic bottom assumed to have the same shape as zeroth
C...  approximation but suppressed by (mc/mb)**2=0.1 approximately
           xpq(5)=xpq(4)*0.1
           xpq(-5)=xpq(5)
           xpq(4)=0.
           xpq(-4)=0.
        ENDIF
      ELSE
C...  Parton densities from PYSTFU in PYTHIA 5.7
*     PEPSI>>
         IF(lst(15).LT.100.or.lst(15).gt.1000) THEN !bs> add gt. 1000
          CALL pystfu(kf,x,q2,xpyst)
          DO 110 kfl=-6,6
 110           xpq(kfl)=xpyst(kfl)
         ELSE 
            CALL parton(x,q2,xpq,xdpq)
         ENDIF
*     PEPSI<<
      ENDIF
         
      DO 120 kfl=-6,6
 120     xpqsav(1,kfl)=xpq(kfl)
C...  Formats for error printouts.
 5000 FORMAT(' Error in LYSTFU: x =',1p,e12.4,' outside physical range')
            
      RETURN
      END