PDF/spdf/sasgam2.F

   1 #include "pdf/pilot.h"
   2 C...SaSgam version 2 - parton distributions of the photon
   3 C...by Gerhard A. Schuler and Torbjorn Sjostrand
   4 C...For further information see Z. Phys. C68 (1995) 607
   5 C...and CERN-TH/96-04 and LU TP 96-2.
   6 C...Program last changed on 18 January 1996.
   7
   8 C!!!Note that one further call parameter - IP2 - has been added
   9 C!!!to the SASGAM argument list compared with version 1.
  10
  11 C...The user should only need to call the SASGAM routine,
  12 C...which in turn calls the auxiliary routines SASVMD, SASANO,
  13 C...SASBEH and SASDIR. The package is self-contained.
  14
  15 C...One particular aspect of these parametrizations is that F2 for
  16 C...the photon is not obtained just as the charge-squared-weighted
  17 C...sum of quark distributions, but differ in the treatment of
  18 C...heavy flavours (in F2 the DIS relation W2 = Q2*(1-x)/x restricts
  19 C...the kinematics range of heavy-flavour production, but the same
  20 C...kinematics is not relevant e.g. for jet production) and, for the
  21 C...'MSbar' fits, in the addition of a Cgamma term related to the
  22 C...separation of direct processes. Schematically:
  23 C...PDF = VMD (rho, omega, phi) + anomalous (d, u, s, c, b).
  24 C...F2  = VMD (rho, omega, phi) + anomalous (d, u, s) +
  25 C...      Bethe-Heitler (c, b) (+ Cgamma (d, u, s)).
  26 C...The J/psi and Upsilon states have not been included in the VMD sum,
  27 C...but low c and b masses in the other components should compensate
  28 C...for this in a duality sense.
  29
  30 C...The calling sequence is the following:
  31 C     CALL SASGAM2(ISET,X,Q2,P2,IP2,F2GM,XPDFGM)
  32 C...with the following declaration statement:
  33 C     DIMENSION XPDFGM(-6:6)
  34 C...and, optionally, further information in:
  35 C     COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
  36 C    &XPDIR(-6:6)
  37 C     COMMON/SASVAL/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),VXPDGM(-6:6)
  38 C...Input:  ISET = 1 : SaS set 1D ('DIS',   Q0 = 0.6 GeV)
  39 C                = 2 : SaS set 1M ('MSbar', Q0 = 0.6 GeV)
  40 C                = 3 : SaS set 2D ('DIS',   Q0 =  2  GeV)
  41 C                = 4 : SaS set 2M ('MSbar', Q0 =  2  GeV)
  42 C           X : x value.
  43 C           Q2 : Q2 value.
  44 C           P2 : P2 value; should be = 0. for an on-shell photon.
  45 C           IP2 : scheme used to evaluate off-shell anomalous component.
  46 C               = 0 : recommended default, see = 7.
  47 C               = 1 : dipole dampening by integration; very time-consuming.
  48 C               = 2 : P_0^2 = max( Q_0^2, P^2 )
  49 C               = 3 : P'_0^2 = Q_0^2 + P^2.
  50 C               = 4 : P_{eff} that preserves momentum sum.
  51 C               = 5 : P_{int} that preserves momentum and average
  52 C                     evolution range.
  53 C               = 6 : P_{eff}, matched to P_0 in P2 -> Q2 limit.
  54 C               = 7 : P_{eff}, matched to P_0 in P2 -> Q2 limit.
  55 C...Output: F2GM : F2 value of the photon (including factors of alpha_em).
  56 C           XPFDGM :  x times parton distribution functions of the photon,
  57 C               with elements 0 = g, 1 = d, 2 = u, 3 = s, 4 = c, 5 = b,
  58 C               6 = t (always empty!), - for antiquarks (result is same).
  59 C...The breakdown by component is stored in the commonblock SASCOM,
  60 C               with elements as above.
  61 C           XPVMD : rho, omega, phi VMD part only of output.
  62 C           XPANL : d, u, s anomalous part only of output.
  63 C           XPANH : c, b anomalous part only of output.
  64 C           XPBEH : c, b Bethe-Heitler part only of output.
  65 C           XPDIR : Cgamma (direct contribution) part only of output.
  66 C...The above arrays do not distinguish valence and sea contributions,
  67 C...although this information is available internally. The additional
  68 C...commonblock SASVAL provides the valence part only of the above
  69 C...distributions. Array names VXPVMD, VXPANL and VXPANH correspond
  70 C...to XPVMD, XPANL and XPANH, while XPBEH and XPDIR are valence only
  71 C...and therefore not given doubly. VXPDGM gives the sum of valence
  72 C...parts, and so matches XPDFGM. The difference, i.e. XPVMD-VXPVMD
  73 C...and so on, gives the sea part only.
  74
  75       SUBROUTINE SASGAM2(ISET,X,Q2,P2,IP2,F2GM,XPDFGM)
  76 C...Purpose: to construct the F2 and parton distributions of the photon
  77 C...by summing homogeneous (VMD) and inhomogeneous (anomalous) terms.
  78 C...For F2, c and b are included by the Bethe-Heitler formula;
  79 C...in the 'MSbar' scheme additionally a Cgamma term is added.
  80       DIMENSION XPDFGM(-6:6)
  81       COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
  82      &XPDIR(-6:6)
  83       COMMON/SASVAL/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),VXPDGM(-6:6)
  84       SAVE /SASCOM/,/SASVAL/
  85
  86 C...Temporary array.
  87       DIMENSION XPGA(-6:6), VXPGA(-6:6)
  88 C...Charm and bottom masses (low to compensate for J/psi etc.).
  89       DATA PMC/1.3/, PMB/4.6/
  90 C...alpha_em and alpha_em/(2*pi).
  91       DATA AEM/0.007297/, AEM2PI/0.0011614/
  92 C...Lambda value for 4 flavours.
  93       DATA ALAM/0.20/
  94 C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
  95       DATA FRACU/0.8/
  96 C...VMD couplings f_V**2/(4*pi).
  97       DATA FRHO/2.20/, FOMEGA/23.6/, FPHI/18.4/
  98 C...Masses for rho (=omega) and phi.
  99       DATA PMRHO/0.770/, PMPHI/1.020/
 100 C...Number of points in integration for IP2=1.
 101       DATA NSTEP/100/
 102
 103 C...Reset output.
 104       F2GM=0.
 105       DO 100 KFL=-6,6
 106       XPDFGM(KFL)=0.
 107       XPVMD(KFL)=0.
 108       XPANL(KFL)=0.
 109       XPANH(KFL)=0.
 110       XPBEH(KFL)=0.
 111       XPDIR(KFL)=0.
 112       VXPVMD(KFL)=0.
 113       VXPANL(KFL)=0.
 114       VXPANH(KFL)=0.
 115       VXPDGM(KFL)=0.
 116   100 CONTINUE
 117
 118 C...Check that input sensible.
 119       IF(ISET.LE.0.OR.ISET.GE.5) THEN
 120         WRITE(*,*) ' FATAL ERROR: SaSgam called for unknown set'
 121         WRITE(*,*) ' ISET = ',ISET
 122         STOP
 123       ENDIF
 124       IF(X.LE.0..OR.X.GT.1.) THEN
 125         WRITE(*,*) ' FATAL ERROR: SaSgam called for unphysical x'
 126         WRITE(*,*) ' X = ',X
 127         STOP
 128       ENDIF
 129
 130 C...Set Q0 cut-off parameter as function of set used.
 131       IF(ISET.LE.2) THEN
 132         Q0=0.6
 133       ELSE
 134         Q0=2.
 135       ENDIF
 136       Q02=Q0**2
 137
 138 C...Scale choice for off-shell photon; common factors.
 139       Q2A=Q2
 140       FACNOR=1.
 141       IF(IP2.EQ.1) THEN
 142         P2MX=P2+Q02
 143         Q2A=Q2+P2*Q02/MAX(Q02,Q2)
 144         FACNOR=LOG(Q2/Q02)/NSTEP
 145       ELSEIF(IP2.EQ.2) THEN
 146         P2MX=MAX(P2,Q02)
 147       ELSEIF(IP2.EQ.3) THEN
 148         P2MX=P2+Q02
 149         Q2A=Q2+P2*Q02/MAX(Q02,Q2)
 150       ELSEIF(IP2.EQ.4) THEN
 151         P2MX=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 152      &  ((Q2+P2)*(Q02+P2)))
 153       ELSEIF(IP2.EQ.5) THEN
 154         P2MXA=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 155      &  ((Q2+P2)*(Q02+P2)))
 156         P2MX=Q0*SQRT(P2MXA)
 157         FACNOR=LOG(Q2/P2MXA)/LOG(Q2/P2MX)
 158       ELSEIF(IP2.EQ.6) THEN
 159         P2MX=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 160      &  ((Q2+P2)*(Q02+P2)))
 161         P2MX=MAX(0.,1.-P2/Q2)*P2MX+MIN(1.,P2/Q2)*MAX(P2,Q02)
 162       ELSE
 163         P2MXA=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 164      &  ((Q2+P2)*(Q02+P2)))
 165         P2MX=Q0*SQRT(P2MXA)
 166         P2MXB=P2MX
 167         P2MX=MAX(0.,1.-P2/Q2)*P2MX+MIN(1.,P2/Q2)*MAX(P2,Q02)
 168         P2MXB=MAX(0.,1.-P2/Q2)*P2MXB+MIN(1.,P2/Q2)*P2MXA
 169         FACNOR=LOG(Q2/P2MXA)/LOG(Q2/P2MXB)
 170       ENDIF
 171
 172 C...Call VMD parametrization for d quark and use to give rho, omega,
 173 C...phi. Note dipole dampening for off-shell photon.
 174       CALL SASVMD(ISET,1,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 175       XFVAL=VXPGA(1)
 176       XPGA(1)=XPGA(2)
 177       XPGA(-1)=XPGA(-2)
 178       FACUD=AEM*(1./FRHO+1./FOMEGA)*(PMRHO**2/(PMRHO**2+P2))**2
 179       FACS=AEM*(1./FPHI)*(PMPHI**2/(PMPHI**2+P2))**2
 180       DO 110 KFL=-5,5
 181       XPVMD(KFL)=(FACUD+FACS)*XPGA(KFL)
 182   110 CONTINUE
 183       XPVMD(1)=XPVMD(1)+(1.-FRACU)*FACUD*XFVAL
 184       XPVMD(2)=XPVMD(2)+FRACU*FACUD*XFVAL
 185       XPVMD(3)=XPVMD(3)+FACS*XFVAL
 186       XPVMD(-1)=XPVMD(-1)+(1.-FRACU)*FACUD*XFVAL
 187       XPVMD(-2)=XPVMD(-2)+FRACU*FACUD*XFVAL
 188       XPVMD(-3)=XPVMD(-3)+FACS*XFVAL
 189       VXPVMD(1)=(1.-FRACU)*FACUD*XFVAL
 190       VXPVMD(2)=FRACU*FACUD*XFVAL
 191       VXPVMD(3)=FACS*XFVAL
 192       VXPVMD(-1)=(1.-FRACU)*FACUD*XFVAL
 193       VXPVMD(-2)=FRACU*FACUD*XFVAL
 194       VXPVMD(-3)=FACS*XFVAL
 195
 196       IF(IP2.NE.1) THEN
 197 C...Anomalous parametrizations for different strategies
 198 C...for off-shell photons; except full integration.
 199
 200 C...Call anomalous parametrization for d + u + s.
 201         CALL SASANO(-3,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 202         DO 120 KFL=-5,5
 203         XPANL(KFL)=FACNOR*XPGA(KFL)
 204         VXPANL(KFL)=FACNOR*VXPGA(KFL)
 205   120   CONTINUE
 206
 207 C...Call anomalous parametrization for c and b.
 208         CALL SASANO(4,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 209         DO 130 KFL=-5,5
 210         XPANH(KFL)=FACNOR*XPGA(KFL)
 211         VXPANH(KFL)=FACNOR*VXPGA(KFL)
 212   130   CONTINUE
 213         CALL SASANO(5,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 214         DO 140 KFL=-5,5
 215         XPANH(KFL)=XPANH(KFL)+FACNOR*XPGA(KFL)
 216         VXPANH(KFL)=VXPANH(KFL)+FACNOR*VXPGA(KFL)
 217   140   CONTINUE
 218
 219       ELSE
 220 C...Special option: loop over flavours and integrate over k2.
 221         DO 170 KF=1,5
 222         DO 160 ISTEP=1,NSTEP
 223         Q2STEP=Q02*(Q2/Q02)**((ISTEP-0.5)/NSTEP)
 224         IF((KF.EQ.4.AND.Q2STEP.LT.PMC**2).OR.
 225      &  (KF.EQ.5.AND.Q2STEP.LT.PMB**2)) GOTO 160
 226         CALL SASVMD(0,KF,X,Q2,Q2STEP,ALAM,XPGA,VXPGA)
 227         FACQ=AEM2PI*(Q2STEP/(Q2STEP+P2))**2*FACNOR
 228         IF(MOD(KF,2).EQ.0) FACQ=FACQ*(8./9.)
 229         IF(MOD(KF,2).EQ.1) FACQ=FACQ*(2./9.)
 230         DO 150 KFL=-5,5
 231         IF(KF.LE.3) XPANL(KFL)=XPANL(KFL)+FACQ*XPGA(KFL)
 232         IF(KF.GE.4) XPANH(KFL)=XPANH(KFL)+FACQ*XPGA(KFL)
 233         IF(KF.LE.3) VXPANL(KFL)=VXPANL(KFL)+FACQ*VXPGA(KFL)
 234         IF(KF.GE.4) VXPANH(KFL)=VXPANH(KFL)+FACQ*VXPGA(KFL)
 235   150   CONTINUE
 236   160   CONTINUE
 237   170   CONTINUE
 238       ENDIF
 239
 240 C...Call Bethe-Heitler term expression for charm and bottom.
 241       CALL SASBEH(4,X,Q2,P2,PMC**2,XPBH)
 242       XPBEH(4)=XPBH
 243       XPBEH(-4)=XPBH
 244       CALL SASBEH(5,X,Q2,P2,PMB**2,XPBH)
 245       XPBEH(5)=XPBH
 246       XPBEH(-5)=XPBH
 247
 248 C...For MSbar subtraction call C^gamma term expression for d, u, s.
 249       IF(ISET.EQ.2.OR.ISET.EQ.4) THEN
 250         CALL SASDIR(X,Q2,P2,Q02,XPGA)
 251         DO 180 KFL=-5,5
 252         XPDIR(KFL)=XPGA(KFL)
 253   180   CONTINUE
 254       ENDIF
 255
 256 C...Store result in output array.
 257       DO 190 KFL=-5,5
 258       CHSQ=1./9.
 259       IF(IABS(KFL).EQ.2.OR.IABS(KFL).EQ.4) CHSQ=4./9.
 260       XPF2=XPVMD(KFL)+XPANL(KFL)+XPBEH(KFL)+XPDIR(KFL)
 261       IF(KFL.NE.0) F2GM=F2GM+CHSQ*XPF2
 262       XPDFGM(KFL)=XPVMD(KFL)+XPANL(KFL)+XPANH(KFL)
 263       VXPDGM(KFL)=VXPVMD(KFL)+VXPANL(KFL)+VXPANH(KFL)
 264   190 CONTINUE
 265
 266       RETURN
 267       END