pythiaeRHIC/pythia/pygvmd.f

0001
0002 C*********************************************************************
0003
0004 C...PYGVMD
0005 C...Evaluates the VMD parton distributions of a photon,
0006 C...evolved homogeneously from an initial scale P2 to Q2.
0007 C...Does not include dipole suppression factor.
0008 C...ISET is parton distribution set, see above;
0009 C...additionally ISET=0 is used for the evolution of an anomalous photon
0010 C...which branched at a scale P2 and then evolved homogeneously to Q2.
0011 C...ALAM is the 4-flavour Lambda, which is automatically converted
0012 C...to 3- and 5-flavour equivalents as needed.
0013 C...Adapted from SaSgam library, authors G.A. Schuler and T. Sjostrand.
0014
0015       SUBROUTINE PYGVMD(ISET,KF,X,Q2,P2,ALAM,XPGA,VXPGA)
0016
0017 C...Double precision and integer declarations.
0018       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0019       IMPLICIT INTEGER(I-N)
0020       INTEGER PYK,PYCHGE,PYCOMP
0021 C...Local arrays and data.
0022       DIMENSION XPGA(-6:6), VXPGA(-6:6)
0023       DATA PMC/1.3D0/, PMB/4.6D0/, AEM/0.007297D0/, AEM2PI/0.0011614D0/
0024
0025 C...Reset output.
0026       DO 100 KFL=-6,6
0027         XPGA(KFL)=0D0
0028         VXPGA(KFL)=0D0
0029   100 CONTINUE
0030       KFA=IABS(KF)
0031
0032 C...Calculate Lambda; protect against unphysical Q2 and P2 input.
0033       ALAM3=ALAM*(PMC/ALAM)**(2D0/27D0)
0034       ALAM5=ALAM*(ALAM/PMB)**(2D0/23D0)
0035       P2EFF=MAX(P2,1.2D0*ALAM3**2)
0036       IF(KFA.EQ.4) P2EFF=MAX(P2EFF,PMC**2)
0037       IF(KFA.EQ.5) P2EFF=MAX(P2EFF,PMB**2)
0038       Q2EFF=MAX(Q2,P2EFF)
0039
0040 C...Find number of flavours at lower and upper scale.
0041       NFP=4
0042       IF(P2EFF.LT.PMC**2) NFP=3
0043       IF(P2EFF.GT.PMB**2) NFP=5
0044       NFQ=4
0045       IF(Q2EFF.LT.PMC**2) NFQ=3
0046       IF(Q2EFF.GT.PMB**2) NFQ=5
0047
0048 C...Find s as sum of 3-, 4- and 5-flavour parts.
0049       S=0D0
0050       IF(NFP.EQ.3) THEN
0051         Q2DIV=PMC**2
0052         IF(NFQ.EQ.3) Q2DIV=Q2EFF
0053         S=S+(6D0/27D0)*LOG(LOG(Q2DIV/ALAM3**2)/LOG(P2EFF/ALAM3**2))
0054       ENDIF
0055       IF(NFP.LE.4.AND.NFQ.GE.4) THEN
0056         P2DIV=P2EFF
0057         IF(NFP.EQ.3) P2DIV=PMC**2
0058         Q2DIV=Q2EFF
0059         IF(NFQ.EQ.5) Q2DIV=PMB**2
0060         S=S+(6D0/25D0)*LOG(LOG(Q2DIV/ALAM**2)/LOG(P2DIV/ALAM**2))
0061       ENDIF
0062       IF(NFQ.EQ.5) THEN
0063         P2DIV=PMB**2
0064         IF(NFP.EQ.5) P2DIV=P2EFF
0065         S=S+(6D0/23D0)*LOG(LOG(Q2EFF/ALAM5**2)/LOG(P2DIV/ALAM5**2))
0066       ENDIF
0067
0068 C...Calculate frequent combinations of x and s.
0069       X1=1D0-X
0070       XL=-LOG(X)
0071       S2=S**2
0072       S3=S**3
0073       S4=S**4
0074
0075 C...Evaluate homogeneous anomalous parton distributions below or
0076 C...above threshold.
0077       IF(ISET.EQ.0) THEN
0078         IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
0079      &  (KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
0080           XVAL = X * 1.5D0 * (X**2+X1**2)
0081           XGLU = 0D0
0082           XSEA = 0D0
0083         ELSE
0084           XVAL = (1.5D0/(1D0-0.197D0*S+4.33D0*S2)*X**2 +
0085      &    (1.5D0+2.10D0*S)/(1D0+3.29D0*S)*X1**2 +
0086      &    5.23D0*S/(1D0+1.17D0*S+19.9D0*S3)*X*X1) *
0087      &    X**(1D0/(1D0+1.5D0*S)) * (1D0-X**2)**(2.667D0*S)
0088           XGLU = 4D0*S/(1D0+4.76D0*S+15.2D0*S2+29.3D0*S4) *
0089      &    X**(-2.03D0*S/(1D0+2.44D0*S)) * (X1*XL)**(1.333D0*S) *
0090      &    ((4D0*X**2+7D0*X+4D0)*X1/3D0 - 2D0*X*(1D0+X)*XL)
0091           XSEA = S2/(1D0+4.54D0*S+8.19D0*S2+8.05D0*S3) *
0092      &    X**(-1.54D0*S/(1D0+1.29D0*S)) * X1**(2.667D0*S) *
0093      &    ((8D0-73D0*X+62D0*X**2)*X1/9D0 + (3D0-8D0*X**2/3D0)*X*XL +
0094      &    (2D0*X-1D0)*X*XL**2)
0095         ENDIF
0096
0097 C...Evaluate set 1D parton distributions below or above threshold.
0098       ELSEIF(ISET.EQ.1) THEN
0099         IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
0100      &  (KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
0101           XVAL = 1.294D0 * X**0.80D0 * X1**0.76D0
0102           XGLU = 1.273D0 * X**0.40D0 * X1**1.76D0
0103           XSEA = 0.100D0 * X1**3.76D0
0104         ELSE
0105           XVAL = 1.294D0/(1D0+0.252D0*S+3.079D0*S2) *
0106      &    X**(0.80D0-0.13D0*S) * X1**(0.76D0+0.667D0*S) * XL**(2D0*S)
0107           XGLU = 7.90D0*S/(1D0+5.50D0*S) * EXP(-5.16D0*S) *
0108      &    X**(-1.90D0*S/(1D0+3.60D0*S)) * X1**1.30D0 *
0109      &    XL**(0.50D0+3D0*S) + 1.273D0 * EXP(-10D0*S) *
0110      &    X**0.40D0 * X1**(1.76D0+3D0*S)
0111           XSEA = (0.1D0-0.397D0*S2+1.121D0*S3)/
0112      &    (1D0+5.61D0*S2+5.26D0*S3) * X**(-7.32D0*S2/(1D0+10.3D0*S2)) *
0113      &    X1**((3.76D0+15D0*S+12D0*S2)/(1D0+4D0*S))
0114           XSEA0 = 0.100D0 * X1**3.76D0
0115         ENDIF
0116
0117 C...Evaluate set 1M parton distributions below or above threshold.
0118       ELSEIF(ISET.EQ.2) THEN
0119         IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
0120      &  (KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
0121           XVAL = 0.8477D0 * X**0.51D0 * X1**1.37D0
0122           XGLU = 3.42D0 * X**0.255D0 * X1**2.37D0
0123           XSEA = 0D0
0124         ELSE
0125           XVAL = 0.8477D0/(1D0+1.37D0*S+2.18D0*S2+3.73D0*S3) *
0126      &    X**(0.51D0+0.21D0*S) * X1**1.37D0 * XL**(2.667D0*S)
0127           XGLU = 24D0*S/(1D0+9.6D0*S+0.92D0*S2+14.34D0*S3) *
0128      &    EXP(-5.94D0*S) * X**((-0.013D0-1.80D0*S)/(1D0+3.14D0*S)) *
0129      &    X1**(2.37D0+0.4D0*S) * XL**(0.32D0+3.6D0*S) + 3.42D0 *
0130      &    EXP(-12D0*S) * X**0.255D0 * X1**(2.37D0+3D0*S)
0131           XSEA = 0.842D0*S/(1D0+21.3D0*S-33.2D0*S2+229D0*S3) *
0132      &    X**((0.13D0-2.90D0*S)/(1D0+5.44D0*S)) * X1**(3.45D0+0.5D0*S) *
0133      &    XL**(2.8D0*S)
0134           XSEA0 = 0D0
0135         ENDIF
0136
0137 C...Evaluate set 2D parton distributions below or above threshold.
0138       ELSEIF(ISET.EQ.3) THEN
0139         IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
0140      &  (KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
0141           XVAL = X**0.46D0 * X1**0.64D0 + 0.76D0 * X
0142           XGLU = 1.925D0 * X1**2
0143           XSEA = 0.242D0 * X1**4
0144         ELSE
0145           XVAL = (1D0+0.186D0*S)/(1D0-0.209D0*S+1.495D0*S2) *
0146      &    X**(0.46D0+0.25D0*S) *
0147      &    X1**((0.64D0+0.14D0*S+5D0*S2)/(1D0+S)) * XL**(1.9D0*S) +
0148      &    (0.76D0+0.4D0*S) * X * X1**(2.667D0*S)
0149           XGLU = (1.925D0+5.55D0*S+147D0*S2)/(1D0-3.59D0*S+3.32D0*S2) *
0150      &    EXP(-18.67D0*S) *
0151      &    X**((-5.81D0*S-5.34D0*S2)/(1D0+29D0*S-4.26D0*S2))
0152      &    * X1**((2D0-5.9D0*S)/(1D0+1.7D0*S)) *
0153      &    XL**(9.3D0*S/(1D0+1.7D0*S))
0154           XSEA = (0.242D0-0.252D0*S+1.19D0*S2)/
0155      &    (1D0-0.607D0*S+21.95D0*S2) *
0156      &    X**(-12.1D0*S2/(1D0+2.62D0*S+16.7D0*S2)) * X1**4 * XL**S
0157           XSEA0 = 0.242D0 * X1**4
0158         ENDIF
0159
0160 C...Evaluate set 2M parton distributions below or above threshold.
0161       ELSEIF(ISET.EQ.4) THEN
0162         IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
0163      &  (KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
0164           XVAL = 1.168D0 * X**0.50D0 * X1**2.60D0 + 0.965D0 * X
0165           XGLU = 1.808D0 * X1**2
0166           XSEA = 0.209D0 * X1**4
0167         ELSE
0168           XVAL = (1.168D0+1.771D0*S+29.35D0*S2) * EXP(-5.776D0*S) *
0169      &    X**((0.5D0+0.208D0*S)/(1D0-0.794D0*S+1.516D0*S2)) *
0170      &    X1**((2.6D0+7.6D0*S)/(1D0+5D0*S)) *
0171      &    XL**(5.15D0*S/(1D0+2D0*S)) +
0172      &    (0.965D0+22.35D0*S)/(1D0+18.4D0*S) * X * X1**(2.667D0*S)
0173           XGLU = (1.808D0+29.9D0*S)/(1D0+26.4D0*S) * EXP(-5.28D0*S) *
0174      &    X**((-5.35D0*S-10.11D0*S2)/(1D0+31.71D0*S)) *
0175      &    X1**((2D0-7.3D0*S+4D0*S2)/(1D0+2.5D0*S)) *
0176      &    XL**(10.9D0*S/(1D0+2.5D0*S))
0177           XSEA = (0.209D0+0.644D0*S2)/(1D0+0.319D0*S+17.6D0*S2) *
0178      &    X**((-0.373D0*S-7.71D0*S2)/(1D0+0.815D0*S+11.0D0*S2)) *
0179      &    X1**(4D0+S) * XL**(0.45D0*S)
0180           XSEA0 = 0.209D0 * X1**4
0181         ENDIF
0182       ENDIF
0183
0184 C...Threshold factors for c and b sea.
0185       SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
0186       XCHM=0D0
0187       IF(Q2.GT.PMC**2.AND.Q2.GT.1.001D0*P2EFF) THEN
0188         SCH=MAX(0D0,LOG(LOG(PMC**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
0189         IF(ISET.EQ.0) THEN
0190           XCHM=XSEA*(1D0-(SCH/SLL)**2)
0191         ELSE
0192           XCHM=MAX(0D0,XSEA-XSEA0*X1**(2.667D0*S))*(1D0-SCH/SLL)
0193         ENDIF
0194       ENDIF
0195       XBOT=0D0
0196       IF(Q2.GT.PMB**2.AND.Q2.GT.1.001D0*P2EFF) THEN
0197         SBT=MAX(0D0,LOG(LOG(PMB**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
0198         IF(ISET.EQ.0) THEN
0199           XBOT=XSEA*(1D0-(SBT/SLL)**2)
0200         ELSE
0201           XBOT=MAX(0D0,XSEA-XSEA0*X1**(2.667D0*S))*(1D0-SBT/SLL)
0202         ENDIF
0203       ENDIF
0204
0205 C...Fill parton distributions.
0206       XPGA(0)=XGLU
0207       XPGA(1)=XSEA
0208       XPGA(2)=XSEA
0209       XPGA(3)=XSEA
0210       XPGA(4)=XCHM
0211       XPGA(5)=XBOT
0212       XPGA(KFA)=XPGA(KFA)+XVAL
0213       DO 110 KFL=1,5
0214         XPGA(-KFL)=XPGA(KFL)
0215   110 CONTINUE
0216       VXPGA(KFA)=XVAL
0217       VXPGA(-KFA)=XVAL
0218
0219       RETURN
0220       END