pythiaeRHIC/pythia/pyggam.f

0001
0002 C*********************************************************************
0003
0004 C...PYGGAM
0005 C...Constructs the F2 and parton distributions of the photon
0006 C...by summing homogeneous (VMD) and inhomogeneous (anomalous) terms.
0007 C...For F2, c and b are included by the Bethe-Heitler formula;
0008 C...in the 'MSbar' scheme additionally a Cgamma term is added.
0009 C...Contains the SaS sets 1D, 1M, 2D and 2M.
0010 C...Adapted from SaSgam library, authors G.A. Schuler and T. Sjostrand.
0011
0012       SUBROUTINE PYGGAM(ISET,X,Q2,P2,IP2,F2GM,XPDFGM)
0013
0014 C...Double precision and integer declarations.
0015       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0016       IMPLICIT INTEGER(I-N)
0017       INTEGER PYK,PYCHGE,PYCOMP
0018 C...Commonblocks.
0019       COMMON/PYINT8/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
0020      &XPDIR(-6:6)
0021       COMMON/PYINT9/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),VXPDGM(-6:6)
0022       SAVE /PYINT8/,/PYINT9/
0023 C...Local arrays.
0024       DIMENSION XPDFGM(-6:6),XPGA(-6:6), VXPGA(-6:6)
0025 C...Charm and bottom masses (low to compensate for J/psi etc.).
0026       DATA PMC/1.3D0/, PMB/4.6D0/
0027 C...alpha_em and alpha_em/(2*pi).
0028       DATA AEM/0.007297D0/, AEM2PI/0.0011614D0/
0029 C...Lambda value for 4 flavours.
0030       DATA ALAM/0.20D0/
0031 C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
0032       DATA FRACU/0.8D0/
0033 C...VMD couplings f_V**2/(4*pi).
0034       DATA FRHO/2.20D0/, FOMEGA/23.6D0/, FPHI/18.4D0/
0035 C...Masses for rho (=omega) and phi.
0036       DATA PMRHO/0.770D0/, PMPHI/1.020D0/
0037 C...Number of points in integration for IP2=1.
0038       DATA NSTEP/100/
0039
0040 C...Reset output.
0041       F2GM=0D0
0042       DO 100 KFL=-6,6
0043         XPDFGM(KFL)=0D0
0044         XPVMD(KFL)=0D0
0045         XPANL(KFL)=0D0
0046         XPANH(KFL)=0D0
0047         XPBEH(KFL)=0D0
0048         XPDIR(KFL)=0D0
0049         VXPVMD(KFL)=0D0
0050         VXPANL(KFL)=0D0
0051         VXPANH(KFL)=0D0
0052         VXPDGM(KFL)=0D0
0053   100 CONTINUE
0054
0055 C...Set Q0 cut-off parameter as function of set used.
0056       IF(ISET.LE.2) THEN
0057         Q0=0.6D0
0058       ELSE
0059         Q0=2D0
0060       ENDIF
0061       Q02=Q0**2
0062
0063 C...Scale choice for off-shell photon; common factors.
0064       Q2A=Q2
0065       FACNOR=1D0
0066       IF(IP2.EQ.1) THEN
0067         P2MX=P2+Q02
0068         Q2A=Q2+P2*Q02/MAX(Q02,Q2)
0069         FACNOR=LOG(Q2/Q02)/NSTEP
0070       ELSEIF(IP2.EQ.2) THEN
0071         P2MX=MAX(P2,Q02)
0072       ELSEIF(IP2.EQ.3) THEN
0073         P2MX=P2+Q02
0074         Q2A=Q2+P2*Q02/MAX(Q02,Q2)
0075       ELSEIF(IP2.EQ.4) THEN
0076         P2MX=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
0077      &  ((Q2+P2)*(Q02+P2)))
0078       ELSEIF(IP2.EQ.5) THEN
0079         P2MXA=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
0080      &  ((Q2+P2)*(Q02+P2)))
0081         P2MX=Q0*SQRT(P2MXA)
0082         FACNOR=LOG(Q2/P2MXA)/LOG(Q2/P2MX)
0083       ELSEIF(IP2.EQ.6) THEN
0084         P2MX=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
0085      &  ((Q2+P2)*(Q02+P2)))
0086         P2MX=MAX(0D0,1D0-P2/Q2)*P2MX+MIN(1D0,P2/Q2)*MAX(P2,Q02)
0087       ELSE
0088         P2MXA=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
0089      &  ((Q2+P2)*(Q02+P2)))
0090         P2MX=Q0*SQRT(P2MXA)
0091         P2MXB=P2MX
0092         P2MX=MAX(0D0,1D0-P2/Q2)*P2MX+MIN(1D0,P2/Q2)*MAX(P2,Q02)
0093         P2MXB=MAX(0D0,1D0-P2/Q2)*P2MXB+MIN(1D0,P2/Q2)*P2MXA
0094         IF(ABS(Q2-Q02).GT.1D-6) THEN
0095           FACNOR=LOG(Q2/P2MXA)/LOG(Q2/P2MXB)
0096         ELSEIF(P2.LT.Q02) THEN
0097           FACNOR=Q02**3/(Q02+P2)/(Q02**2-P2**2/2D0)
0098         ELSE
0099           FACNOR=1D0
0100         ENDIF
0101       ENDIF
0102
0103 C...Call VMD parametrization for d quark and use to give rho, omega,
0104 C...phi. Note dipole dampening for off-shell photon.
0105       CALL PYGVMD(ISET,1,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
0106       XFVAL=VXPGA(1)
0107       XPGA(1)=XPGA(2)
0108       XPGA(-1)=XPGA(-2)
0109       FACUD=AEM*(1D0/FRHO+1D0/FOMEGA)*(PMRHO**2/(PMRHO**2+P2))**2
0110       FACS=AEM*(1D0/FPHI)*(PMPHI**2/(PMPHI**2+P2))**2
0111       DO 110 KFL=-5,5
0112         XPVMD(KFL)=(FACUD+FACS)*XPGA(KFL)
0113   110 CONTINUE
0114       XPVMD(1)=XPVMD(1)+(1D0-FRACU)*FACUD*XFVAL
0115       XPVMD(2)=XPVMD(2)+FRACU*FACUD*XFVAL
0116       XPVMD(3)=XPVMD(3)+FACS*XFVAL
0117       XPVMD(-1)=XPVMD(-1)+(1D0-FRACU)*FACUD*XFVAL
0118       XPVMD(-2)=XPVMD(-2)+FRACU*FACUD*XFVAL
0119       XPVMD(-3)=XPVMD(-3)+FACS*XFVAL
0120       VXPVMD(1)=(1D0-FRACU)*FACUD*XFVAL
0121       VXPVMD(2)=FRACU*FACUD*XFVAL
0122       VXPVMD(3)=FACS*XFVAL
0123       VXPVMD(-1)=(1D0-FRACU)*FACUD*XFVAL
0124       VXPVMD(-2)=FRACU*FACUD*XFVAL
0125       VXPVMD(-3)=FACS*XFVAL
0126
0127       IF(IP2.NE.1) THEN
0128 C...Anomalous parametrizations for different strategies
0129 C...for off-shell photons; except full integration.
0130
0131 C...Call anomalous parametrization for d + u + s.
0132         CALL PYGANO(-3,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
0133         DO 120 KFL=-5,5
0134           XPANL(KFL)=FACNOR*XPGA(KFL)
0135           VXPANL(KFL)=FACNOR*VXPGA(KFL)
0136   120   CONTINUE
0137
0138 C...Call anomalous parametrization for c and b.
0139         CALL PYGANO(4,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
0140         DO 130 KFL=-5,5
0141           XPANH(KFL)=FACNOR*XPGA(KFL)
0142           VXPANH(KFL)=FACNOR*VXPGA(KFL)
0143   130   CONTINUE
0144         CALL PYGANO(5,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
0145         DO 140 KFL=-5,5
0146           XPANH(KFL)=XPANH(KFL)+FACNOR*XPGA(KFL)
0147           VXPANH(KFL)=VXPANH(KFL)+FACNOR*VXPGA(KFL)
0148   140   CONTINUE
0149
0150       ELSE
0151 C...Special option: loop over flavours and integrate over k2.
0152         DO 170 KF=1,5
0153           DO 160 ISTEP=1,NSTEP
0154             Q2STEP=Q02*(Q2/Q02)**((ISTEP-0.5D0)/NSTEP)
0155             IF((KF.EQ.4.AND.Q2STEP.LT.PMC**2).OR.
0156      &      (KF.EQ.5.AND.Q2STEP.LT.PMB**2)) GOTO 160
0157             CALL PYGVMD(0,KF,X,Q2,Q2STEP,ALAM,XPGA,VXPGA)
0158             FACQ=AEM2PI*(Q2STEP/(Q2STEP+P2))**2*FACNOR
0159             IF(MOD(KF,2).EQ.0) FACQ=FACQ*(8D0/9D0)
0160             IF(MOD(KF,2).EQ.1) FACQ=FACQ*(2D0/9D0)
0161             DO 150 KFL=-5,5
0162               IF(KF.LE.3) XPANL(KFL)=XPANL(KFL)+FACQ*XPGA(KFL)
0163               IF(KF.GE.4) XPANH(KFL)=XPANH(KFL)+FACQ*XPGA(KFL)
0164               IF(KF.LE.3) VXPANL(KFL)=VXPANL(KFL)+FACQ*VXPGA(KFL)
0165               IF(KF.GE.4) VXPANH(KFL)=VXPANH(KFL)+FACQ*VXPGA(KFL)
0166   150       CONTINUE
0167   160     CONTINUE
0168   170   CONTINUE
0169       ENDIF
0170
0171 C...Call Bethe-Heitler term expression for charm and bottom.
0172       CALL PYGBEH(4,X,Q2,P2,PMC**2,XPBH)
0173       XPBEH(4)=XPBH
0174       XPBEH(-4)=XPBH
0175       CALL PYGBEH(5,X,Q2,P2,PMB**2,XPBH)
0176       XPBEH(5)=XPBH
0177       XPBEH(-5)=XPBH
0178
0179 C...For MSbar subtraction call C^gamma term expression for d, u, s.
0180       IF(ISET.EQ.2.OR.ISET.EQ.4) THEN
0181         CALL PYGDIR(X,Q2,P2,Q02,XPGA)
0182         DO 180 KFL=-5,5
0183           XPDIR(KFL)=XPGA(KFL)
0184   180   CONTINUE
0185       ENDIF
0186
0187 C...Store result in output array.
0188       DO 190 KFL=-5,5
0189         CHSQ=1D0/9D0
0190         IF(IABS(KFL).EQ.2.OR.IABS(KFL).EQ.4) CHSQ=4D0/9D0
0191         XPF2=XPVMD(KFL)+XPANL(KFL)+XPBEH(KFL)+XPDIR(KFL)
0192         IF(KFL.NE.0) F2GM=F2GM+CHSQ*XPF2
0193         XPDFGM(KFL)=XPVMD(KFL)+XPANL(KFL)+XPANH(KFL)
0194         VXPDGM(KFL)=VXPVMD(KFL)+VXPANL(KFL)+VXPANH(KFL)
0195   190 CONTINUE
0196
0197       RETURN
0198       END