pythiaeRHIC/pythia/pygano.f

0001
0002 C*********************************************************************
0003
0004 C...PYGANO
0005 C...Evaluates the parton distributions of the anomalous photon,
0006 C...inhomogeneously evolved from a scale P2 (where it vanishes) to Q2.
0007 C...KF=0 gives the sum over (up to) 5 flavours,
0008 C...KF<0 limits to flavours up to abs(KF),
0009 C...KF>0 is for flavour KF only.
0010 C...ALAM is the 4-flavour Lambda, which is automatically converted
0011 C...to 3- and 5-flavour equivalents as needed.
0012 C...Adapted from SaSgam library, authors G.A. Schuler and T. Sjostrand.
0013
0014       SUBROUTINE PYGANO(KF,X,Q2,P2,ALAM,XPGA,VXPGA)
0015
0016 C...Double precision and integer declarations.
0017       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0018       IMPLICIT INTEGER(I-N)
0019       INTEGER PYK,PYCHGE,PYCOMP
0020 C...Local arrays and data.
0021       DIMENSION XPGA(-6:6), VXPGA(-6:6), ALAMSQ(3:5)
0022       DATA PMC/1.3D0/, PMB/4.6D0/, AEM/0.007297D0/, AEM2PI/0.0011614D0/
0023
0024 C...Reset output.
0025       DO 100 KFL=-6,6
0026         XPGA(KFL)=0D0
0027         VXPGA(KFL)=0D0
0028   100 CONTINUE
0029       IF(Q2.LE.P2) RETURN
0030       KFA=IABS(KF)
0031
0032 C...Calculate Lambda; protect against unphysical Q2 and P2 input.
0033       ALAMSQ(3)=(ALAM*(PMC/ALAM)**(2D0/27D0))**2
0034       ALAMSQ(4)=ALAM**2
0035       ALAMSQ(5)=(ALAM*(ALAM/PMB)**(2D0/23D0))**2
0036       P2EFF=MAX(P2,1.2D0*ALAMSQ(3))
0037       IF(KF.EQ.4) P2EFF=MAX(P2EFF,PMC**2)
0038       IF(KF.EQ.5) P2EFF=MAX(P2EFF,PMB**2)
0039       Q2EFF=MAX(Q2,P2EFF)
0040       XL=-LOG(X)
0041
0042 C...Find number of flavours at lower and upper scale.
0043       NFP=4
0044       IF(P2EFF.LT.PMC**2) NFP=3
0045       IF(P2EFF.GT.PMB**2) NFP=5
0046       NFQ=4
0047       IF(Q2EFF.LT.PMC**2) NFQ=3
0048       IF(Q2EFF.GT.PMB**2) NFQ=5
0049
0050 C...Define range of flavour loop.
0051       IF(KF.EQ.0) THEN
0052         KFLMN=1
0053         KFLMX=5
0054       ELSEIF(KF.LT.0) THEN
0055         KFLMN=1
0056         KFLMX=KFA
0057       ELSE
0058         KFLMN=KFA
0059         KFLMX=KFA
0060       ENDIF
0061
0062 C...Loop over flavours the photon can branch into.
0063       DO 110 KFL=KFLMN,KFLMX
0064
0065 C...Light flavours: calculate t range and (approximate) s range.
0066         IF(KFL.LE.3.AND.(KFL.EQ.1.OR.KFL.EQ.KF)) THEN
0067           TDIFF=LOG(Q2EFF/P2EFF)
0068           S=(6D0/(33D0-2D0*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
0069      &    LOG(P2EFF/ALAMSQ(NFQ)))
0070           IF(NFQ.GT.NFP) THEN
0071             Q2DIV=PMB**2
0072             IF(NFQ.EQ.4) Q2DIV=PMC**2
0073             SNFQ=(6D0/(33D0-2D0*NFQ))*LOG(LOG(Q2DIV/ALAMSQ(NFQ))/
0074      &      LOG(P2EFF/ALAMSQ(NFQ)))
0075             SNFP=(6D0/(33D0-2D0*(NFQ-1)))*LOG(LOG(Q2DIV/ALAMSQ(NFQ-1))/
0076      &      LOG(P2EFF/ALAMSQ(NFQ-1)))
0077             S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNFP-SNFQ)
0078           ENDIF
0079           IF(NFQ.EQ.5.AND.NFP.EQ.3) THEN
0080             Q2DIV=PMC**2
0081             SNF4=(6D0/(33D0-2D0*4))*LOG(LOG(Q2DIV/ALAMSQ(4))/
0082      &      LOG(P2EFF/ALAMSQ(4)))
0083             SNF3=(6D0/(33D0-2D0*3))*LOG(LOG(Q2DIV/ALAMSQ(3))/
0084      &      LOG(P2EFF/ALAMSQ(3)))
0085             S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNF3-SNF4)
0086           ENDIF
0087
0088 C...u and s quark do not need a separate treatment when d has been done.
0089         ELSEIF(KFL.EQ.2.OR.KFL.EQ.3) THEN
0090
0091 C...Charm: as above, but only include range above c threshold.
0092         ELSEIF(KFL.EQ.4) THEN
0093           IF(Q2.LE.PMC**2) GOTO 110
0094           P2EFF=MAX(P2EFF,PMC**2)
0095           Q2EFF=MAX(Q2EFF,P2EFF)
0096           TDIFF=LOG(Q2EFF/P2EFF)
0097           S=(6D0/(33D0-2D0*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
0098      &    LOG(P2EFF/ALAMSQ(NFQ)))
0099           IF(NFQ.EQ.5.AND.NFP.EQ.4) THEN
0100             Q2DIV=PMB**2
0101             SNFQ=(6D0/(33D0-2D0*NFQ))*LOG(LOG(Q2DIV/ALAMSQ(NFQ))/
0102      &      LOG(P2EFF/ALAMSQ(NFQ)))
0103             SNFP=(6D0/(33D0-2D0*(NFQ-1)))*LOG(LOG(Q2DIV/ALAMSQ(NFQ-1))/
0104      &      LOG(P2EFF/ALAMSQ(NFQ-1)))
0105             S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNFP-SNFQ)
0106           ENDIF
0107
0108 C...Bottom: as above, but only include range above b threshold.
0109         ELSEIF(KFL.EQ.5) THEN
0110           IF(Q2.LE.PMB**2) GOTO 110
0111           P2EFF=MAX(P2EFF,PMB**2)
0112           Q2EFF=MAX(Q2,P2EFF)
0113           TDIFF=LOG(Q2EFF/P2EFF)
0114           S=(6D0/(33D0-2D0*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
0115      &    LOG(P2EFF/ALAMSQ(NFQ)))
0116         ENDIF
0117
0118 C...Evaluate flavour-dependent prefactor (charge^2 etc.).
0119         CHSQ=1D0/9D0
0120         IF(KFL.EQ.2.OR.KFL.EQ.4) CHSQ=4D0/9D0
0121         FAC=AEM2PI*2D0*CHSQ*TDIFF
0122
0123 C...Evaluate parton distributions (normalized to unit momentum sum).
0124         IF(KFL.EQ.1.OR.KFL.EQ.4.OR.KFL.EQ.5.OR.KFL.EQ.KF) THEN
0125           XVAL= ((1.5D0+2.49D0*S+26.9D0*S**2)/(1D0+32.3D0*S**2)*X**2 +
0126      &    (1.5D0-0.49D0*S+7.83D0*S**2)/(1D0+7.68D0*S**2)*(1D0-X)**2 +
0127      &    1.5D0*S/(1D0-3.2D0*S+7D0*S**2)*X*(1D0-X)) *
0128      &    X**(1D0/(1D0+0.58D0*S)) * (1D0-X**2)**(2.5D0*S/(1D0+10D0*S))
0129           XGLU= 2D0*S/(1D0+4D0*S+7D0*S**2) *
0130      &    X**(-1.67D0*S/(1D0+2D0*S)) * (1D0-X**2)**(1.2D0*S) *
0131      &    ((4D0*X**2+7D0*X+4D0)*(1D0-X)/3D0 - 2D0*X*(1D0+X)*XL)
0132           XSEA= 0.333D0*S**2/(1D0+4.90D0*S+4.69D0*S**2+21.4D0*S**3) *
0133      &    X**(-1.18D0*S/(1D0+1.22D0*S)) * (1D0-X)**(1.2D0*S) *
0134      &    ((8D0-73D0*X+62D0*X**2)*(1D0-X)/9D0 +
0135      &    (3D0-8D0*X**2/3D0)*X*XL + (2D0*X-1D0)*X*XL**2)
0136
0137 C...Threshold factors for c and b sea.
0138           SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
0139           XCHM=0D0
0140           IF(Q2.GT.PMC**2.AND.Q2.GT.1.001D0*P2EFF) THEN
0141             SCH=MAX(0D0,LOG(LOG(PMC**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
0142             XCHM=XSEA*(1D0-(SCH/SLL)**3)
0143           ENDIF
0144           XBOT=0D0
0145           IF(Q2.GT.PMB**2.AND.Q2.GT.1.001D0*P2EFF) THEN
0146             SBT=MAX(0D0,LOG(LOG(PMB**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
0147             XBOT=XSEA*(1D0-(SBT/SLL)**3)
0148           ENDIF
0149         ENDIF
0150
0151 C...Add contribution of each valence flavour.
0152         XPGA(0)=XPGA(0)+FAC*XGLU
0153         XPGA(1)=XPGA(1)+FAC*XSEA
0154         XPGA(2)=XPGA(2)+FAC*XSEA
0155         XPGA(3)=XPGA(3)+FAC*XSEA
0156         XPGA(4)=XPGA(4)+FAC*XCHM
0157         XPGA(5)=XPGA(5)+FAC*XBOT
0158         XPGA(KFL)=XPGA(KFL)+FAC*XVAL
0159         VXPGA(KFL)=VXPGA(KFL)+FAC*XVAL
0160   110 CONTINUE
0161       DO 120 KFL=1,5
0162         XPGA(-KFL)=XPGA(KFL)
0163         VXPGA(-KFL)=VXPGA(KFL)
0164   120 CONTINUE
0165
0166       RETURN
0167       END