LHAPDF/lhapdf5.3.1/wrapsasg.f

   1       subroutine SASGevolvep(xin,qin,p2in,ip2in,pdf)
   2       real*8 xin,qin,q2in,p2in,pdf(-6:6),xval(45),qcdl4,qcdl5
   3       include 'parmsetup.inc'
   4       character*16 name(nmxset)
   5       integer nmem(nmxset),ndef(nmxset),mmem
   6       common/NAME/name,nmem,ndef,mmem
   7       integer nset
   8       real*8 upv,dnv,usea,dsea,str,chm,bot,top,glu
   9
  10       save
  11       call getnset(iset)
  12       call getnmem(iset,imem)
  13
  14       iimem = imem
  15       if(iimem.eq.0) iimem=6
  16       q2in = qin*qin
  17       call SFSASxx(iimem,xin,Q2in,p2in,ip2,
  18      +                     upv,dnv,usea,dsea,str,chm,bot,top,glu)
  19
  20       pdf(-6)= top
  21       pdf(6)= top
  22       pdf(-5)= bot
  23       pdf(5 )= bot
  24       pdf(-4)= chm
  25       pdf(4 )= chm
  26       pdf(-3)= str
  27       pdf(3 )= str
  28       pdf(-2)= usea
  29       pdf(2 )= upv
  30       pdf(-1)= dsea
  31       pdf(1 )= dnv
  32       pdf(0 )= glu
  33
  34       return
  35 c
  36 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  37       entry SASGread(nset)
  38 c      print *,'calling SASGread'
  39       read(1,*)nmem(nset),ndef(nset)
  40       return
  41 c
  42 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  43       entry SASGalfa(alfas,qalfa)
  44         call getnset(iset)
  45         call getnmem(iset,imem)
  46         call GetOrderAsM(iset,iord)
  47         call Getlam4M(iset,imem,qcdl4)
  48         call Getlam5M(iset,imem,qcdl5)
  49         call aspdflib(alfas,Qalfa,iord,qcdl5)
  50
  51       return
  52 c
  53 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  54       entry SASGinit(Eorder,Q2fit)
  55       return
  56 c
  57 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  58       entry SASGpdf(mem)
  59 c      print *,'calling SASGpdf',mem
  60 c      imem = mem
  61       call getnset(iset)
  62       call setnmem(iset,mem)
  63       return
  64 c
  65  1000 format(5e13.5)
  66       end
  67 c
  68 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  69 C-----------------------------------------------------------------------
  70       SUBROUTINE SFSASxx(iset,DX,DQ2,DP2,ip2,
  71      +                     DUPV,DDNV,DSEA,DSEAD,DSTR,DCHM,DBOT,DTOP,DGL)
  72 C
  73 C   ********************************************************************
  74 C   *                                                                  *
  75 C   *        Interface to SASset of structure functions                *
  76 C   *                                                                  *
  77 C   *        Author:    H. Plothow-Besch (CERN-PPE)                    *
  78 C   *                                                                  *
  79 C   ********************************************************************
  80 C
  81 C :::::::::::: Structure functions from the SAS group version 2
  82 C :::::::::::: Lambda = 0.200 GeV, Q**2 = 0.36 GeV**2 (DIS)
  83 C
  84       double precision
  85      +          DX,DQ2,DP2,
  86      +          DUPV,DDNV,DSEA,DSEAD,DSTR,DCHM,DBOT,DTOP,DGL
  87       DIMENSION XPDFGM(-6:6)
  88       REAL X, Q, Q2, P2, F2GAM, XPDFGM
  89 c      PARAMETER (ISET=1)
  90 C
  91       X  = DX
  92       Q  = SQRT(DQ2)
  93       Q2 = DQ2
  94       P2 = DP2
  95 C
  96 C     generate the individual structure fcn calls
  97 C
  98       if(iset.le.4) then
  99          iiset=iset
 100          CALL LHASASGAM1(iISET,X,Q2,P2,F2GAM,XPDFGM)
 101       else
 102          iiset = iset-4
 103          CALL LHASASGAM2(iISET,X,Q2,P2,ip2,F2GAM,XPDFGM)
 104       endif
 105
 106       UPV = XPDFGM(2)
 107       DUPV = UPV
 108       DNV = XPDFGM(1)
 109       DDNV = DNV
 110       SEAU = XPDFGM(-2)
 111       DSEA = SEAU
 112       SEAD = XPDFGM(-1)
 113       DSEAD = SEAD
 114       STR = XPDFGM(-3)
 115       DSTR = STR
 116       CHM = XPDFGM(-4)
 117       DCHM = CHM
 118       BOT = XPDFGM(-5)
 119       DBOT = BOT
 120       TOP = 0.
 121 C      IF (DSCAL.GT.TMAS) TOP = XPDFGM(6)
 122       DTOP = TOP
 123       GL = XPDFGM(0)
 124       DGL = GL
 125 C
 126       RETURN
 127       END
 128 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 129 c ------------- SASGAM1 ------------------------
 130 C...SaSgam - parton distributions of the photon
 131 C...by Gerhard A. Schuler and Torbjorn Sjostrand
 132 C...For further information see preprint CERN-TH/95-62 and LU TP 95-6:
 133 C...Low- and high-mass components of the photon distribution functions
 134 C...Program last changed on 21 March 1995.
 135
 136 C...The user should only need to call the SASGAM routine,
 137 C...which in turn calls the auxiliary routines SASVM1, SASAN1,
 138 C...SASBEH and SASDIR. The package is self-contained.
 139
 140 C...One particular aspect of these parametrizations is that F2 for
 141 C...the photon is not obtained just as the charge-squared-weighted
 142 C...sum of quark distributions, but differ in the treatment of
 143 C...heavy flavours (in F2 the DIS relation W2 = Q2*(1-x)/x restricts
 144 C...the kinematics range of heavy-flavour production, but the same
 145 C...kinematics is not relevant e.g. for jet production) and, for the
 146 C...'MSbar' fits, in the addition of a Cgamma term related to the
 147 C...separation of direct processes. Schematically:
 148 C...PDF = VMD (rho, omega, phi) + anomalous (d, u, s, c, b).
 149 C...F2  = VMD (rho, omega, phi) + anomalous (d, u, s) +
 150 C...      Bethe-Heitler (c, b) (+ Cgamma (d, u, s)).
 151 C...The J/psi and Upsilon states have not been included in the VMD sum,
 152 C...but low c and b masses in the other components should compensate
 153 C...for this in a duality sense.
 154
 155 C...The calling sequence is the following:
 156 C     CALL SASGAM1(ISET,X,Q2,P2,F2GM,XPDFGM)
 157 C...with the following declaration statement:
 158 C     DIMENSION XPDFGM(-6:6)
 159 C...and, optionally, further information in:
 160 C     COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
 161 C    &XPDIR(-6:6)
 162 C...Input:  ISET = 1 : SaS set 1D ('DIS',   Q0 = 0.6 GeV)
 163 C                = 2 : SaS set 1M ('MSbar', Q0 = 0.6 GeV)
 164 C                = 3 : SaS set 2D ('DIS',   Q0 =  2  GeV)
 165 C                = 4 : SaS set 2M ('MSbar', Q0 =  2  GeV)
 166 C           X : x value.
 167 C           Q2 : Q2 value.
 168 C           P2 : P2 value; should be = 0. for an on-shell photon.
 169 C...Output: F2GM : F2 value of the photon (including factors of alpha_em).
 170 C           XPFDGM :  x times parton distribution functions of the photon,
 171 C               with elements 0 = g, 1 = d, 2 = u, 3 = s, 4 = c, 5 = b,
 172 C               6 = t (always empty!), - for antiquarks (result is same).
 173 C...The breakdown by component is stored in the commonblock SASCOM,
 174 C               with elements as above.
 175 C           XPVMD : rho, omega, phi VMD part only of output.
 176 C           XPANL : d, u, s anomalous part only of output.
 177 C           XPANH : c, b anomalous part only of output.
 178 C           XPBEH : c, b Bethe-Heitler part only of output.
 179 C           XPDIR : Cgamma (direct contribution) part only of output.
 180
 181       SUBROUTINE LHASASGAM1(ISET,X,Q2,P2,F2GM,XPDFGM)
 182 C...Purpose: to construct the F2 and parton distributions of the photon
 183 C...by summing homogeneous (VMD) and inhomogeneous (anomalous) terms.
 184 C...For F2, c and b are included by the Bethe-Heitler formula;
 185 C...in the 'MSbar' scheme additionally a Cgamma term is added.
 186       DIMENSION XPDFGM(-6:6)
 187       COMMON/LHASASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
 188      &XPDIR(-6:6)
 189       SAVE /LHASASCOM/
 190
 191 C...Temporary array.
 192       DIMENSION XPGA(-6:6)
 193 C...Charm and bottom masses (low to compensate for J/psi etc.).
 194       DATA PMC/1.3/, PMB/4.6/
 195 C...alpha_em and alpha_em/(2*pi).
 196       DATA AEM/0.007297/, AEM2PI/0.0011614/
 197 C...Lambda value for 4 flavours.
 198       DATA ALAM/0.20/
 199 C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
 200       DATA FRACU/0.8/
 201 C...VMD couplings f_V**2/(4*pi).
 202       DATA FRHO/2.20/, FOMEGA/23.6/, FPHI/18.4/
 203 C...Masses for rho (=omega) and phi.
 204       DATA PMRHO/0.770/, PMPHI/1.020/
 205
 206 C...Reset output.
 207       F2GM=0.
 208       DO 100 KFL=-6,6
 209       XPDFGM(KFL)=0.
 210       XPVMD(KFL)=0.
 211       XPANL(KFL)=0.
 212       XPANH(KFL)=0.
 213       XPBEH(KFL)=0.
 214       XPDIR(KFL)=0.
 215   100 CONTINUE
 216
 217 C...Check that input sensible.
 218       IF(ISET.LE.0.OR.ISET.GE.5) THEN
 219         WRITE(*,*) ' FATAL ERROR: SaSgam called for unknown set'
 220         WRITE(*,*) ' ISET = ',ISET
 221         STOP
 222       ENDIF
 223       IF(X.LE.0..OR.X.GT.1.) THEN
 224         WRITE(*,*) ' FATAL ERROR: SaSgam called for unphysical x'
 225         WRITE(*,*) ' X = ',X
 226         STOP
 227       ENDIF
 228
 229 C...Set Q0 cut-off parameter as function of set used.
 230       IF(ISET.LE.2) THEN
 231         Q0=0.6
 232       ELSE
 233         Q0=2.
 234       ENDIF
 235       Q02 = Q0**2
 236
 237 C...Call VMD parametrization for d quark and use to give rho, omega, phi.
 238 C...Note scale choice and dipole dampening for off-shell photon.
 239       P2MX=MAX(P2,Q02)
 240       CALL LHASASVM1(ISET,1,X,Q2,P2MX,ALAM,XPGA)
 241       XFVAL=XPGA(1)-XPGA(2)
 242       XPGA(1)=XPGA(2)
 243       XPGA(-1)=XPGA(-2)
 244       FACUD=AEM*(1./FRHO+1./FOMEGA)*(PMRHO**2/(PMRHO**2+P2))**2
 245       FACS=AEM*(1./FPHI)*(PMPHI**2/(PMPHI**2+P2))**2
 246       DO 110 KFL=-5,5
 247       XPVMD(KFL)=(FACUD+FACS)*XPGA(KFL)
 248   110 CONTINUE
 249       XPVMD(1)=XPVMD(1)+(1.-FRACU)*FACUD*XFVAL
 250       XPVMD(2)=XPVMD(2)+FRACU*FACUD*XFVAL
 251       XPVMD(3)=XPVMD(3)+FACS*XFVAL
 252       XPVMD(-1)=XPVMD(-1)+(1.-FRACU)*FACUD*XFVAL
 253       XPVMD(-2)=XPVMD(-2)+FRACU*FACUD*XFVAL
 254       XPVMD(-3)=XPVMD(-3)+FACS*XFVAL
 255
 256 C...Call anomalous parametrization for d + u + s.
 257       CALL LHASASAN1(-3,X,Q2,P2MX,ALAM,XPGA)
 258       DO 120 KFL=-5,5
 259       XPANL(KFL)=XPGA(KFL)
 260   120 CONTINUE
 261
 262 C...Call anomalous parametrization for c and b.
 263       CALL LHASASAN1(4,X,Q2,P2MX,ALAM,XPGA)
 264       DO 130 KFL=-5,5
 265       XPANH(KFL)=XPGA(KFL)
 266   130 CONTINUE
 267       CALL LHASASAN1(5,X,Q2,P2MX,ALAM,XPGA)
 268       DO 140 KFL=-5,5
 269       XPANH(KFL)=XPANH(KFL)+XPGA(KFL)
 270   140 CONTINUE
 271
 272 C...Call Bethe-Heitler term expression for charm and bottom.
 273       CALL LHASASBEH(4,X,Q2,P2,PMC**2,XPBH)
 274       XPBEH(4)=XPBH
 275       XPBEH(-4)=XPBH
 276       CALL LHASASBEH(5,X,Q2,P2,PMB**2,XPBH)
 277       XPBEH(5)=XPBH
 278       XPBEH(-5)=XPBH
 279
 280 C...For MSbar subtraction call C^gamma term expression for d, u, s.
 281       IF(ISET.EQ.2.OR.ISET.EQ.4) THEN
 282         CALL LHASASDIR(X,Q2,P2,Q02,XPGA)
 283         DO 150 KFL=-5,5
 284         XPDIR(KFL)=XPGA(KFL)
 285   150   CONTINUE
 286       ENDIF
 287
 288 C...Store result in output array.
 289       DO 160 KFL=-5,5
 290       CHSQ=1./9.
 291       IF(IABS(KFL).EQ.2.OR.IABS(KFL).EQ.4) CHSQ=4./9.
 292       XPF2=XPVMD(KFL)+XPANL(KFL)+XPBEH(KFL)+XPDIR(KFL)
 293       IF(KFL.NE.0) F2GM=F2GM+CHSQ*XPF2
 294       XPDFGM(KFL)=XPVMD(KFL)+XPANL(KFL)+XPANH(KFL)
 295   160 CONTINUE
 296
 297       RETURN
 298       END
 299 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 300 c ------------------ SASGAM2 ---------------------------
 301 C...SaSgam version 2 - parton distributions of the photon
 302 C...by Gerhard A. Schuler and Torbjorn Sjostrand
 303 C...For further information see Z. Phys. C68 (1995) 607
 304 C...and CERN-TH/96-04 and LU TP 96-2.
 305 C...Program last changed on 18 January 1996.
 306
 307 C!!!Note that one further call parameter - IP2 - has been added
 308 C!!!to the SASGAM argument list compared with version 1.
 309
 310 C...The user should only need to call the SASGAM routine,
 311 C...which in turn calls the auxiliary routines SASVMD, SASANO,
 312 C...SASBEH and SASDIR. The package is self-contained.
 313
 314 C...One particular aspect of these parametrizations is that F2 for
 315 C...the photon is not obtained just as the charge-squared-weighted
 316 C...sum of quark distributions, but differ in the treatment of
 317 C...heavy flavours (in F2 the DIS relation W2 = Q2*(1-x)/x restricts
 318 C...the kinematics range of heavy-flavour production, but the same
 319 C...kinematics is not relevant e.g. for jet production) and, for the
 320 C...'MSbar' fits, in the addition of a Cgamma term related to the
 321 C...separation of direct processes. Schematically:
 322 C...PDF = VMD (rho, omega, phi) + anomalous (d, u, s, c, b).
 323 C...F2  = VMD (rho, omega, phi) + anomalous (d, u, s) +
 324 C...      Bethe-Heitler (c, b) (+ Cgamma (d, u, s)).
 325 C...The J/psi and Upsilon states have not been included in the VMD sum,
 326 C...but low c and b masses in the other components should compensate
 327 C...for this in a duality sense.
 328
 329 C...The calling sequence is the following:
 330 C     CALL SASGAM2(ISET,X,Q2,P2,IP2,F2GM,XPDFGM)
 331 C...with the following declaration statement:
 332 C     DIMENSION XPDFGM(-6:6)
 333 C...and, optionally, further information in:
 334 C     COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
 335 C    &XPDIR(-6:6)
 336 C     COMMON/SASVAL/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),VXPDGM(-6:6)
 337 C...Input:  ISET = 1 : SaS set 1D ('DIS',   Q0 = 0.6 GeV)
 338 C                = 2 : SaS set 1M ('MSbar', Q0 = 0.6 GeV)
 339 C                = 3 : SaS set 2D ('DIS',   Q0 =  2  GeV)
 340 C                = 4 : SaS set 2M ('MSbar', Q0 =  2  GeV)
 341 C           X : x value.
 342 C           Q2 : Q2 value.
 343 C           P2 : P2 value; should be = 0. for an on-shell photon.
 344 C           IP2 : scheme used to evaluate off-shell anomalous component.
 345 C               = 0 : recommended default, see = 7.
 346 C               = 1 : dipole dampening by integration; very time-consuming.
 347 C               = 2 : P_0^2 = max( Q_0^2, P^2 )
 348 C               = 3 : P'_0^2 = Q_0^2 + P^2.
 349 C               = 4 : P_{eff} that preserves momentum sum.
 350 C               = 5 : P_{int} that preserves momentum and average
 351 C                     evolution range.
 352 C               = 6 : P_{eff}, matched to P_0 in P2 -> Q2 limit.
 353 C               = 7 : P_{eff}, matched to P_0 in P2 -> Q2 limit.
 354 C...Output: F2GM : F2 value of the photon (including factors of alpha_em).
 355 C           XPFDGM :  x times parton distribution functions of the photon,
 356 C               with elements 0 = g, 1 = d, 2 = u, 3 = s, 4 = c, 5 = b,
 357 C               6 = t (always empty!), - for antiquarks (result is same).
 358 C...The breakdown by component is stored in the commonblock SASCOM,
 359 C               with elements as above.
 360 C           XPVMD : rho, omega, phi VMD part only of output.
 361 C           XPANL : d, u, s anomalous part only of output.
 362 C           XPANH : c, b anomalous part only of output.
 363 C           XPBEH : c, b Bethe-Heitler part only of output.
 364 C           XPDIR : Cgamma (direct contribution) part only of output.
 365 C...The above arrays do not distinguish valence and sea contributions,
 366 C...although this information is available internally. The additional
 367 C...commonblock SASVAL provides the valence part only of the above
 368 C...distributions. Array names VXPVMD, VXPANL and VXPANH correspond
 369 C...to XPVMD, XPANL and XPANH, while XPBEH and XPDIR are valence only
 370 C...and therefore not given doubly. VXPDGM gives the sum of valence
 371 C...parts, and so matches XPDFGM. The difference, i.e. XPVMD-VXPVMD
 372 C...and so on, gives the sea part only.
 373
 374       SUBROUTINE LHASASGAM2(ISET,X,Q2,P2,IP2,F2GM,XPDFGM)
 375 C...Purpose: to construct the F2 and parton distributions of the photon
 376 C...by summing homogeneous (VMD) and inhomogeneous (anomalous) terms.
 377 C...For F2, c and b are included by the Bethe-Heitler formula;
 378 C...in the 'MSbar' scheme additionally a Cgamma term is added.
 379       DIMENSION XPDFGM(-6:6)
 380       COMMON/LHASASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
 381      &XPDIR(-6:6)
 382       COMMON/LHASASVAL/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),
 383      &VXPDGM(-6:6)
 384       SAVE /LHASASCOM/,/LHASASVAL/
 385
 386 C...Temporary array.
 387       DIMENSION XPGA(-6:6), VXPGA(-6:6)
 388 C...Charm and bottom masses (low to compensate for J/psi etc.).
 389       DATA PMC/1.3/, PMB/4.6/
 390 C...alpha_em and alpha_em/(2*pi).
 391       DATA AEM/0.007297/, AEM2PI/0.0011614/
 392 C...Lambda value for 4 flavours.
 393       DATA ALAM/0.20/
 394 C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
 395       DATA FRACU/0.8/
 396 C...VMD couplings f_V**2/(4*pi).
 397       DATA FRHO/2.20/, FOMEGA/23.6/, FPHI/18.4/
 398 C...Masses for rho (=omega) and phi.
 399       DATA PMRHO/0.770/, PMPHI/1.020/
 400 C...Number of points in integration for IP2=1.
 401       DATA NSTEP/100/
 402
 403 C...Reset output.
 404       F2GM=0.
 405       DO 100 KFL=-6,6
 406       XPDFGM(KFL)=0.
 407       XPVMD(KFL)=0.
 408       XPANL(KFL)=0.
 409       XPANH(KFL)=0.
 410       XPBEH(KFL)=0.
 411       XPDIR(KFL)=0.
 412       VXPVMD(KFL)=0.
 413       VXPANL(KFL)=0.
 414       VXPANH(KFL)=0.
 415       VXPDGM(KFL)=0.
 416   100 CONTINUE
 417
 418 C...Check that input sensible.
 419       IF(ISET.LE.0.OR.ISET.GE.5) THEN
 420         WRITE(*,*) ' FATAL ERROR: SaSgam called for unknown set'
 421         WRITE(*,*) ' ISET = ',ISET
 422         STOP
 423       ENDIF
 424       IF(X.LE.0..OR.X.GT.1.) THEN
 425         WRITE(*,*) ' FATAL ERROR: SaSgam called for unphysical x'
 426         WRITE(*,*) ' X = ',X
 427         STOP
 428       ENDIF
 429
 430 C...Set Q0 cut-off parameter as function of set used.
 431       IF(ISET.LE.2) THEN
 432         Q0=0.6
 433       ELSE
 434         Q0=2.
 435       ENDIF
 436       Q02=Q0**2
 437
 438 C...Scale choice for off-shell photon; common factors.
 439       Q2A=Q2
 440       FACNOR=1.
 441       IF(IP2.EQ.1) THEN
 442         P2MX=P2+Q02
 443         Q2A=Q2+P2*Q02/MAX(Q02,Q2)
 444         FACNOR=LOG(Q2/Q02)/NSTEP
 445       ELSEIF(IP2.EQ.2) THEN
 446         P2MX=MAX(P2,Q02)
 447       ELSEIF(IP2.EQ.3) THEN
 448         P2MX=P2+Q02
 449         Q2A=Q2+P2*Q02/MAX(Q02,Q2)
 450       ELSEIF(IP2.EQ.4) THEN
 451         P2MX=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 452      &  ((Q2+P2)*(Q02+P2)))
 453       ELSEIF(IP2.EQ.5) THEN
 454         P2MXA=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 455      &  ((Q2+P2)*(Q02+P2)))
 456         P2MX=Q0*SQRT(P2MXA)
 457         FACNOR=LOG(Q2/P2MXA)/LOG(Q2/P2MX)
 458       ELSEIF(IP2.EQ.6) THEN
 459         P2MX=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 460      &  ((Q2+P2)*(Q02+P2)))
 461         P2MX=MAX(0.,1.-P2/Q2)*P2MX+MIN(1.,P2/Q2)*MAX(P2,Q02)
 462       ELSE
 463         P2MXA=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 464      &  ((Q2+P2)*(Q02+P2)))
 465         P2MX=Q0*SQRT(P2MXA)
 466         P2MXB=P2MX
 467         P2MX=MAX(0.,1.-P2/Q2)*P2MX+MIN(1.,P2/Q2)*MAX(P2,Q02)
 468         P2MXB=MAX(0.,1.-P2/Q2)*P2MXB+MIN(1.,P2/Q2)*P2MXA
 469         FACNOR=LOG(Q2/P2MXA)/LOG(Q2/P2MXB)
 470       ENDIF
 471
 472 C...Call VMD parametrization for d quark and use to give rho, omega,
 473 C...phi. Note dipole dampening for off-shell photon.
 474       CALL LHASASVMD(ISET,1,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 475       XFVAL=VXPGA(1)
 476       XPGA(1)=XPGA(2)
 477       XPGA(-1)=XPGA(-2)
 478       FACUD=AEM*(1./FRHO+1./FOMEGA)*(PMRHO**2/(PMRHO**2+P2))**2
 479       FACS=AEM*(1./FPHI)*(PMPHI**2/(PMPHI**2+P2))**2
 480       DO 110 KFL=-5,5
 481       XPVMD(KFL)=(FACUD+FACS)*XPGA(KFL)
 482   110 CONTINUE
 483       XPVMD(1)=XPVMD(1)+(1.-FRACU)*FACUD*XFVAL
 484       XPVMD(2)=XPVMD(2)+FRACU*FACUD*XFVAL
 485       XPVMD(3)=XPVMD(3)+FACS*XFVAL
 486       XPVMD(-1)=XPVMD(-1)+(1.-FRACU)*FACUD*XFVAL
 487       XPVMD(-2)=XPVMD(-2)+FRACU*FACUD*XFVAL
 488       XPVMD(-3)=XPVMD(-3)+FACS*XFVAL
 489       VXPVMD(1)=(1.-FRACU)*FACUD*XFVAL
 490       VXPVMD(2)=FRACU*FACUD*XFVAL
 491       VXPVMD(3)=FACS*XFVAL
 492       VXPVMD(-1)=(1.-FRACU)*FACUD*XFVAL
 493       VXPVMD(-2)=FRACU*FACUD*XFVAL
 494       VXPVMD(-3)=FACS*XFVAL
 495
 496       IF(IP2.NE.1) THEN
 497 C...Anomalous parametrizations for different strategies
 498 C...for off-shell photons; except full integration.
 499
 500 C...Call anomalous parametrization for d + u + s.
 501         CALL LHASASANO(-3,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 502         DO 120 KFL=-5,5
 503         XPANL(KFL)=FACNOR*XPGA(KFL)
 504         VXPANL(KFL)=FACNOR*VXPGA(KFL)
 505   120   CONTINUE
 506
 507 C...Call anomalous parametrization for c and b.
 508         CALL LHASASANO(4,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 509         DO 130 KFL=-5,5
 510         XPANH(KFL)=FACNOR*XPGA(KFL)
 511         VXPANH(KFL)=FACNOR*VXPGA(KFL)
 512   130   CONTINUE
 513         CALL LHASASANO(5,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 514         DO 140 KFL=-5,5
 515         XPANH(KFL)=XPANH(KFL)+FACNOR*XPGA(KFL)
 516         VXPANH(KFL)=VXPANH(KFL)+FACNOR*VXPGA(KFL)
 517   140   CONTINUE
 518
 519       ELSE
 520 C...Special option: loop over flavours and integrate over k2.
 521         DO 170 KF=1,5
 522         DO 160 ISTEP=1,NSTEP
 523         Q2STEP=Q02*(Q2/Q02)**((ISTEP-0.5)/NSTEP)
 524         IF((KF.EQ.4.AND.Q2STEP.LT.PMC**2).OR.
 525      &  (KF.EQ.5.AND.Q2STEP.LT.PMB**2)) GOTO 160
 526         CALL LHASASVMD(0,KF,X,Q2,Q2STEP,ALAM,XPGA,VXPGA)
 527         FACQ=AEM2PI*(Q2STEP/(Q2STEP+P2))**2*FACNOR
 528         IF(MOD(KF,2).EQ.0) FACQ=FACQ*(8./9.)
 529         IF(MOD(KF,2).EQ.1) FACQ=FACQ*(2./9.)
 530         DO 150 KFL=-5,5
 531         IF(KF.LE.3) XPANL(KFL)=XPANL(KFL)+FACQ*XPGA(KFL)
 532         IF(KF.GE.4) XPANH(KFL)=XPANH(KFL)+FACQ*XPGA(KFL)
 533         IF(KF.LE.3) VXPANL(KFL)=VXPANL(KFL)+FACQ*VXPGA(KFL)
 534         IF(KF.GE.4) VXPANH(KFL)=VXPANH(KFL)+FACQ*VXPGA(KFL)
 535   150   CONTINUE
 536   160   CONTINUE
 537   170   CONTINUE
 538       ENDIF
 539
 540 C...Call Bethe-Heitler term expression for charm and bottom.
 541       CALL LHASASBEH(4,X,Q2,P2,PMC**2,XPBH)
 542       XPBEH(4)=XPBH
 543       XPBEH(-4)=XPBH
 544       CALL LHASASBEH(5,X,Q2,P2,PMB**2,XPBH)
 545       XPBEH(5)=XPBH
 546       XPBEH(-5)=XPBH
 547
 548 C...For MSbar subtraction call C^gamma term expression for d, u, s.
 549       IF(ISET.EQ.2.OR.ISET.EQ.4) THEN
 550         CALL LHASASDIR(X,Q2,P2,Q02,XPGA)
 551         DO 180 KFL=-5,5
 552         XPDIR(KFL)=XPGA(KFL)
 553   180   CONTINUE
 554       ENDIF
 555
 556 C...Store result in output array.
 557       DO 190 KFL=-5,5
 558       CHSQ=1./9.
 559       IF(IABS(KFL).EQ.2.OR.IABS(KFL).EQ.4) CHSQ=4./9.
 560       XPF2=XPVMD(KFL)+XPANL(KFL)+XPBEH(KFL)+XPDIR(KFL)
 561       IF(KFL.NE.0) F2GM=F2GM+CHSQ*XPF2
 562       XPDFGM(KFL)=XPVMD(KFL)+XPANL(KFL)+XPANH(KFL)
 563       VXPDGM(KFL)=VXPVMD(KFL)+VXPANL(KFL)+VXPANH(KFL)
 564   190 CONTINUE
 565
 566       RETURN
 567       END
 568
 569
 570 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 571       SUBROUTINE LHASASVMD(ISET,KF,X,Q2,P2,ALAM,XPGA,VXPGA)
 572 C...Purpose: to evaluate the VMD parton distributions of a photon,
 573 C...evolved homogeneously from an initial scale P2 to Q2.
 574 C...Does not include dipole suppression factor.
 575 C...ISET is parton distribution set, see above;
 576 C...additionally ISET=0 is used for the evolution of an anomalous photon
 577 C...which branched at a scale P2 and then evolved homogeneously to Q2.
 578 C...ALAM is the 4-flavour Lambda, which is automatically converted
 579 C...to 3- and 5-flavour equivalents as needed.
 580       DIMENSION XPGA(-6:6), VXPGA(-6:6)
 581       DATA PMC/1.3/, PMB/4.6/, AEM/0.007297/, AEM2PI/0.0011614/
 582
 583 C...Reset output.
 584       DO 100 KFL=-6,6
 585       XPGA(KFL)=0.
 586       VXPGA(KFL)=0.
 587   100 CONTINUE
 588       KFA=IABS(KF)
 589
 590 C...Calculate Lambda; protect against unphysical Q2 and P2 input.
 591       ALAM3=ALAM*(PMC/ALAM)**(2./27.)
 592       ALAM5=ALAM*(ALAM/PMB)**(2./23.)
 593       P2EFF=MAX(P2,1.2*ALAM3**2)
 594       IF(KFA.EQ.4) P2EFF=MAX(P2EFF,PMC**2)
 595       IF(KFA.EQ.5) P2EFF=MAX(P2EFF,PMB**2)
 596       Q2EFF=MAX(Q2,P2EFF)
 597
 598 C...Find number of flavours at lower and upper scale.
 599       NFP=4
 600       IF(P2EFF.LT.PMC**2) NFP=3
 601       IF(P2EFF.GT.PMB**2) NFP=5
 602       NFQ=4
 603       IF(Q2EFF.LT.PMC**2) NFQ=3
 604       IF(Q2EFF.GT.PMB**2) NFQ=5
 605
 606 C...Find s as sum of 3-, 4- and 5-flavour parts.
 607       S=0.
 608       IF(NFP.EQ.3) THEN
 609         Q2DIV=PMC**2
 610         IF(NFQ.EQ.3) Q2DIV=Q2EFF
 611         S=S+(6./27.)*LOG(LOG(Q2DIV/ALAM3**2)/LOG(P2EFF/ALAM3**2))
 612       ENDIF
 613       IF(NFP.LE.4.AND.NFQ.GE.4) THEN
 614         P2DIV=P2EFF
 615         IF(NFP.EQ.3) P2DIV=PMC**2
 616         Q2DIV=Q2EFF
 617         IF(NFQ.EQ.5) Q2DIV=PMB**2
 618         S=S+(6./25.)*LOG(LOG(Q2DIV/ALAM**2)/LOG(P2DIV/ALAM**2))
 619       ENDIF
 620       IF(NFQ.EQ.5) THEN
 621         P2DIV=PMB**2
 622         IF(NFP.EQ.5) P2DIV=P2EFF
 623         S=S+(6./23.)*LOG(LOG(Q2EFF/ALAM5**2)/LOG(P2DIV/ALAM5**2))
 624       ENDIF
 625
 626 C...Calculate frequent combinations of x and s.
 627       X1=1.-X
 628       XL=-LOG(X)
 629       S2=S**2
 630       S3=S**3
 631       S4=S**4
 632
 633 C...Evaluate homogeneous anomalous parton distributions below or
 634 C...above threshold.
 635       IF(ISET.EQ.0) THEN
 636       IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
 637      &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
 638         XVAL = X * 1.5 * (X**2+X1**2)
 639         XGLU = 0.
 640         XSEA = 0.
 641       ELSE
 642         XVAL = (1.5/(1.-0.197*S+4.33*S2)*X**2 + (1.5+2.10*S)/
 643      &  (1.+3.29*S)*X1**2 + 5.23*S/(1.+1.17*S+19.9*S3)*X*X1) *
 644      &  X**(1./(1.+1.5*S)) * (1.-X**2)**(2.667*S)
 645         XGLU = 4.*S/(1.+4.76*S+15.2*S2+29.3*S4) *
 646      &  X**(-2.03*S/(1.+2.44*S)) * (X1*XL)**(1.333*S) *
 647      &  ((4.*X**2+7.*X+4.)*X1/3. - 2.*X*(1.+X)*XL)
 648         XSEA = S2/(1.+4.54*S+8.19*S2+8.05*S3) *
 649      &  X**(-1.54*S/(1.+1.29*S)) * X1**(2.667*S) *
 650      &  ((8.-73.*X+62.*X**2)*X1/9. + (3.-8.*X**2/3.)*X*XL +
 651      &  (2.*X-1.)*X*XL**2)
 652       ENDIF
 653
 654 C...Evaluate set 1D parton distributions below or above threshold.
 655       ELSEIF(ISET.EQ.1) THEN
 656       IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
 657      &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
 658         XVAL = 1.294 * X**0.80 * X1**0.76
 659         XGLU = 1.273 * X**0.40 * X1**1.76
 660         XSEA = 0.100 * X1**3.76
 661       ELSE
 662         XVAL = 1.294/(1.+0.252*S+3.079*S2) * X**(0.80-0.13*S) *
 663      &  X1**(0.76+0.667*S) * XL**(2.*S)
 664         XGLU = 7.90*S/(1.+5.50*S) * EXP(-5.16*S) *
 665      &  X**(-1.90*S/(1.+3.60*S)) * X1**1.30 * XL**(0.50+3.*S) +
 666      &  1.273 * EXP(-10.*S) * X**0.40 * X1**(1.76+3.*S)
 667         XSEA = (0.1-0.397*S2+1.121*S3)/(1.+5.61*S2+5.26*S3) *
 668      &  X**(-7.32*S2/(1.+10.3*S2)) *
 669      &  X1**((3.76+15.*S+12.*S2)/(1.+4.*S))
 670         XSEA0 = 0.100 * X1**3.76
 671       ENDIF
 672
 673 C...Evaluate set 1M parton distributions below or above threshold.
 674       ELSEIF(ISET.EQ.2) THEN
 675       IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
 676      &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
 677         XVAL = 0.8477 * X**0.51 * X1**1.37
 678         XGLU = 3.42 * X**0.255 * X1**2.37
 679         XSEA = 0.
 680       ELSE
 681         XVAL = 0.8477/(1.+1.37*S+2.18*S2+3.73*S3) * X**(0.51+0.21*S)
 682      &  * X1**1.37 * XL**(2.667*S)
 683         XGLU = 24.*S/(1.+9.6*S+0.92*S2+14.34*S3) * EXP(-5.94*S) *
 684      &  X**((-0.013-1.80*S)/(1.+3.14*S)) * X1**(2.37+0.4*S) *
 685      &  XL**(0.32+3.6*S) + 3.42 * EXP(-12.*S) * X**0.255 *
 686      &  X1**(2.37+3.*S)
 687         XSEA = 0.842*S/(1.+21.3*S-33.2*S2+229.*S3) *
 688      &  X**((0.13-2.90*S)/(1.+5.44*S)) * X1**(3.45+0.5*S) *
 689      &  XL**(2.8*S)
 690         XSEA0 = 0.
 691       ENDIF
 692
 693 C...Evaluate set 2D parton distributions below or above threshold.
 694       ELSEIF(ISET.EQ.3) THEN
 695       IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
 696      &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
 697         XVAL = X**0.46 * X1**0.64 + 0.76 * X
 698         XGLU = 1.925 * X1**2
 699         XSEA = 0.242 * X1**4
 700       ELSE
 701         XVAL = (1.+0.186*S)/(1.-0.209*S+1.495*S2) * X**(0.46+0.25*S)
 702      &  * X1**((0.64+0.14*S+5.*S2)/(1.+S)) * XL**(1.9*S) +
 703      &  (0.76+0.4*S) * X * X1**(2.667*S)
 704         XGLU = (1.925+5.55*S+147.*S2)/(1.-3.59*S+3.32*S2) *
 705      &  EXP(-18.67*S) * X**((-5.81*S-5.34*S2)/(1.+29.*S-4.26*S2))
 706      &  * X1**((2.-5.9*S)/(1.+1.7*S)) * XL**(9.3*S/(1.+1.7*S))
 707         XSEA = (0.242-0.252*S+1.19*S2)/(1.-0.607*S+21.95*S2) *
 708      &  X**(-12.1*S2/(1.+2.62*S+16.7*S2)) * X1**4 * XL**S
 709         XSEA0 = 0.242 * X1**4
 710       ENDIF
 711
 712 C...Evaluate set 2M parton distributions below or above threshold.
 713       ELSEIF(ISET.EQ.4) THEN
 714       IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
 715      &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
 716         XVAL = 1.168 * X**0.50 * X1**2.60 + 0.965 * X
 717         XGLU = 1.808 * X1**2
 718         XSEA = 0.209 * X1**4
 719       ELSE
 720         XVAL = (1.168+1.771*S+29.35*S2) * EXP(-5.776*S) *
 721      &  X**((0.5+0.208*S)/(1.-0.794*S+1.516*S2)) *
 722      &  X1**((2.6+7.6*S)/(1.+5.*S)) * XL**(5.15*S/(1.+2.*S)) +
 723      &  (0.965+22.35*S)/(1.+18.4*S) * X * X1**(2.667*S)
 724         XGLU = (1.808+29.9*S)/(1.+26.4*S) * EXP(-5.28*S) *
 725      &  X**((-5.35*S-10.11*S2)/(1.+31.71*S)) *
 726      &  X1**((2.-7.3*S+4.*S2)/(1.+2.5*S)) *
 727      &  XL**(10.9*S/(1.+2.5*S))
 728         XSEA = (0.209+0.644*S2)/(1.+0.319*S+17.6*S2) *
 729      &  X**((-0.373*S-7.71*S2)/(1.+0.815*S+11.0*S2)) *
 730      &  X1**(4.+S) * XL**(0.45*S)
 731         XSEA0 = 0.209 * X1**4
 732       ENDIF
 733       ENDIF
 734
 735 C...Threshold factors for c and b sea.
 736       SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
 737       XCHM=0.
 738       IF(Q2.GT.PMC**2.AND.Q2.GT.1.001*P2EFF) THEN
 739         SCH=MAX(0.,LOG(LOG(PMC**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
 740         IF(ISET.EQ.0) THEN
 741           XCHM=XSEA*(1.-(SCH/SLL)**2)
 742         ELSE
 743           XCHM=MAX(0.,XSEA-XSEA0*X1**(2.667*S))*(1.-SCH/SLL)
 744         ENDIF
 745       ENDIF
 746       XBOT=0.
 747       IF(Q2.GT.PMB**2.AND.Q2.GT.1.001*P2EFF) THEN
 748         SBT=MAX(0.,LOG(LOG(PMB**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
 749         IF(ISET.EQ.0) THEN
 750           XBOT=XSEA*(1.-(SBT/SLL)**2)
 751         ELSE
 752           XBOT=MAX(0.,XSEA-XSEA0*X1**(2.667*S))*(1.-SBT/SLL)
 753         ENDIF
 754       ENDIF
 755
 756 C...Fill parton distributions.
 757       XPGA(0)=XGLU
 758       XPGA(1)=XSEA
 759       XPGA(2)=XSEA
 760       XPGA(3)=XSEA
 761       XPGA(4)=XCHM
 762       XPGA(5)=XBOT
 763       XPGA(KFA)=XPGA(KFA)+XVAL
 764       DO 110 KFL=1,5
 765       XPGA(-KFL)=XPGA(KFL)
 766   110 CONTINUE
 767       VXPGA(KFA)=XVAL
 768       VXPGA(-KFA)=XVAL
 769
 770       RETURN
 771       END
 772
 773 C*********************************************************************
 774
 775       SUBROUTINE LHASASANO(KF,X,Q2,P2,ALAM,XPGA,VXPGA)
 776 C...Purpose: to evaluate the parton distributions of the anomalous
 777 C...photon, inhomogeneously evolved from a scale P2 (where it vanishes)
 778 C...to Q2.
 779 C...KF=0 gives the sum over (up to) 5 flavours,
 780 C...KF<0 limits to flavours up to abs(KF),
 781 C...KF>0 is for flavour KF only.
 782 C...ALAM is the 4-flavour Lambda, which is automatically converted
 783 C...to 3- and 5-flavour equivalents as needed.
 784       DIMENSION XPGA(-6:6), VXPGA(-6:6), ALAMSQ(3:5)
 785       DATA PMC/1.3/, PMB/4.6/, AEM/0.007297/, AEM2PI/0.0011614/
 786
 787 C...Reset output.
 788       DO 100 KFL=-6,6
 789       XPGA(KFL)=0.
 790       VXPGA(KFL)=0.
 791   100 CONTINUE
 792       IF(Q2.LE.P2) RETURN
 793       KFA=IABS(KF)
 794
 795 C...Calculate Lambda; protect against unphysical Q2 and P2 input.
 796       ALAMSQ(3)=(ALAM*(PMC/ALAM)**(2./27.))**2
 797       ALAMSQ(4)=ALAM**2
 798       ALAMSQ(5)=(ALAM*(ALAM/PMB)**(2./23.))**2
 799       P2EFF=MAX(P2,1.2*ALAMSQ(3))
 800       IF(KF.EQ.4) P2EFF=MAX(P2EFF,PMC**2)
 801       IF(KF.EQ.5) P2EFF=MAX(P2EFF,PMB**2)
 802       Q2EFF=MAX(Q2,P2EFF)
 803       XL=-LOG(X)
 804
 805 C...Find number of flavours at lower and upper scale.
 806       NFP=4
 807       IF(P2EFF.LT.PMC**2) NFP=3
 808       IF(P2EFF.GT.PMB**2) NFP=5
 809       NFQ=4
 810       IF(Q2EFF.LT.PMC**2) NFQ=3
 811       IF(Q2EFF.GT.PMB**2) NFQ=5
 812
 813 C...Define range of flavour loop.
 814       IF(KF.EQ.0) THEN
 815         KFLMN=1
 816         KFLMX=5
 817       ELSEIF(KF.LT.0) THEN
 818         KFLMN=1
 819         KFLMX=KFA
 820       ELSE
 821         KFLMN=KFA
 822         KFLMX=KFA
 823       ENDIF
 824
 825 C...Loop over flavours the photon can branch into.
 826       DO 110 KFL=KFLMN,KFLMX
 827
 828 C...Light flavours: calculate t range and (approximate) s range.
 829       IF(KFL.LE.3.AND.(KFL.EQ.1.OR.KFL.EQ.KF)) THEN
 830         TDIFF=LOG(Q2EFF/P2EFF)
 831         S=(6./(33.-2.*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
 832      &  LOG(P2EFF/ALAMSQ(NFQ)))
 833         IF(NFQ.GT.NFP) THEN
 834           Q2DIV=PMB**2
 835           IF(NFQ.EQ.4) Q2DIV=PMC**2
 836           SNFQ=(6./(33.-2.*NFQ))*LOG(LOG(Q2DIV/ALAMSQ(NFQ))/
 837      &    LOG(P2EFF/ALAMSQ(NFQ)))
 838           SNFP=(6./(33.-2.*(NFQ-1)))*LOG(LOG(Q2DIV/ALAMSQ(NFQ-1))/
 839      &    LOG(P2EFF/ALAMSQ(NFQ-1)))
 840           S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNFP-SNFQ)
 841         ENDIF
 842         IF(NFQ.EQ.5.AND.NFP.EQ.3) THEN
 843           Q2DIV=PMC**2
 844           SNF4=(6./(33.-2.*4))*LOG(LOG(Q2DIV/ALAMSQ(4))/
 845      &    LOG(P2EFF/ALAMSQ(4)))
 846           SNF3=(6./(33.-2.*3))*LOG(LOG(Q2DIV/ALAMSQ(3))/
 847      &    LOG(P2EFF/ALAMSQ(3)))
 848           S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNF3-SNF4)
 849         ENDIF
 850
 851 C...u and s quark do not need a separate treatment when d has been done.
 852       ELSEIF(KFL.EQ.2.OR.KFL.EQ.3) THEN
 853
 854 C...Charm: as above, but only include range above c threshold.
 855       ELSEIF(KFL.EQ.4) THEN
 856         IF(Q2.LE.PMC**2) GOTO 110
 857         P2EFF=MAX(P2EFF,PMC**2)
 858         Q2EFF=MAX(Q2EFF,P2EFF)
 859         TDIFF=LOG(Q2EFF/P2EFF)
 860         S=(6./(33.-2.*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
 861      &  LOG(P2EFF/ALAMSQ(NFQ)))
 862         IF(NFQ.EQ.5.AND.NFP.EQ.4) THEN
 863           Q2DIV=PMB**2
 864           SNFQ=(6./(33.-2.*NFQ))*LOG(LOG(Q2DIV/ALAMSQ(NFQ))/
 865      &    LOG(P2EFF/ALAMSQ(NFQ)))
 866           SNFP=(6./(33.-2.*(NFQ-1)))*LOG(LOG(Q2DIV/ALAMSQ(NFQ-1))/
 867      &    LOG(P2EFF/ALAMSQ(NFQ-1)))
 868           S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNFP-SNFQ)
 869         ENDIF
 870
 871 C...Bottom: as above, but only include range above b threshold.
 872       ELSEIF(KFL.EQ.5) THEN
 873         IF(Q2.LE.PMB**2) GOTO 110
 874         P2EFF=MAX(P2EFF,PMB**2)
 875         Q2EFF=MAX(Q2,P2EFF)
 876         TDIFF=LOG(Q2EFF/P2EFF)
 877         S=(6./(33.-2.*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
 878      &  LOG(P2EFF/ALAMSQ(NFQ)))
 879       ENDIF
 880
 881 C...Evaluate flavour-dependent prefactor (charge^2 etc.).
 882       CHSQ=1./9.
 883       IF(KFL.EQ.2.OR.KFL.EQ.4) CHSQ=4./9.
 884       FAC=AEM2PI*2.*CHSQ*TDIFF
 885
 886 C...Evaluate parton distributions (normalized to unit momentum sum).
 887       IF(KFL.EQ.1.OR.KFL.EQ.4.OR.KFL.EQ.5.OR.KFL.EQ.KF) THEN
 888         XVAL= ((1.5+2.49*S+26.9*S**2)/(1.+32.3*S**2)*X**2 +
 889      &  (1.5-0.49*S+7.83*S**2)/(1.+7.68*S**2)*(1.-X)**2 +
 890      &  1.5*S/(1.-3.2*S+7.*S**2)*X*(1.-X)) *
 891      &  X**(1./(1.+0.58*S)) * (1.-X**2)**(2.5*S/(1.+10.*S))
 892         XGLU= 2.*S/(1.+4.*S+7.*S**2) *
 893      &  X**(-1.67*S/(1.+2.*S)) * (1.-X**2)**(1.2*S) *
 894      &  ((4.*X**2+7.*X+4.)*(1.-X)/3. - 2.*X*(1.+X)*XL)
 895         XSEA= 0.333*S**2/(1.+4.90*S+4.69*S**2+21.4*S**3) *
 896      &  X**(-1.18*S/(1.+1.22*S)) * (1.-X)**(1.2*S) *
 897      &  ((8.-73.*X+62.*X**2)*(1.-X)/9. + (3.-8.*X**2/3.)*X*XL +
 898      &  (2.*X-1.)*X*XL**2)
 899
 900 C...Threshold factors for c and b sea.
 901         SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
 902         XCHM=0.
 903         IF(Q2.GT.PMC**2.AND.Q2.GT.1.001*P2EFF) THEN
 904           SCH=MAX(0.,LOG(LOG(PMC**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
 905           XCHM=XSEA*(1.-(SCH/SLL)**3)
 906         ENDIF
 907         XBOT=0.
 908         IF(Q2.GT.PMB**2.AND.Q2.GT.1.001*P2EFF) THEN
 909           SBT=MAX(0.,LOG(LOG(PMB**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
 910           XBOT=XSEA*(1.-(SBT/SLL)**3)
 911         ENDIF
 912       ENDIF
 913
 914 C...Add contribution of each valence flavour.
 915       XPGA(0)=XPGA(0)+FAC*XGLU
 916       XPGA(1)=XPGA(1)+FAC*XSEA
 917       XPGA(2)=XPGA(2)+FAC*XSEA
 918       XPGA(3)=XPGA(3)+FAC*XSEA
 919       XPGA(4)=XPGA(4)+FAC*XCHM
 920       XPGA(5)=XPGA(5)+FAC*XBOT
 921       XPGA(KFL)=XPGA(KFL)+FAC*XVAL
 922       VXPGA(KFL)=VXPGA(KFL)+FAC*XVAL
 923   110 CONTINUE
 924       DO 120 KFL=1,5
 925       XPGA(-KFL)=XPGA(KFL)
 926       VXPGA(-KFL)=VXPGA(KFL)
 927   120 CONTINUE
 928
 929       RETURN
 930       END
 931
 932 C*********************************************************************
 933
 934       SUBROUTINE LHASASBEH(KF,X,Q2,P2,PM2,XPBH)
 935 C...Purpose: to evaluate the Bethe-Heitler cross section for
 936 C...heavy flavour production.
 937       DATA AEM2PI/0.0011614/
 938
 939 C...Reset output.
 940       XPBH=0.
 941       SIGBH=0.
 942
 943 C...Check kinematics limits.
 944       IF(X.GE.Q2/(4.*PM2+Q2+P2)) RETURN
 945       W2=Q2*(1.-X)/X-P2
 946       BETA2=1.-4.*PM2/W2
 947       IF(BETA2.LT.1E-10) RETURN
 948       BETA=SQRT(BETA2)
 949       RMQ=4.*PM2/Q2
 950
 951 C...Simple case: P2 = 0.
 952       IF(P2.LT.1E-4) THEN
 953         IF(BETA.LT.0.99) THEN
 954           XBL=LOG((1.+BETA)/(1.-BETA))
 955         ELSE
 956           XBL=LOG((1.+BETA)**2*W2/(4.*PM2))
 957         ENDIF
 958         SIGBH=BETA*(8.*X*(1.-X)-1.-RMQ*X*(1.-X))+
 959      &  XBL*(X**2+(1.-X)**2+RMQ*X*(1.-3.*X)-0.5*RMQ**2*X**2)
 960
 961 C...Complicated case: P2 > 0, based on approximation of
 962 C...C.T. Hill and G.G. Ross, Nucl. Phys. B148 (1979) 373
 963       ELSE
 964         RPQ=1.-4.*X**2*P2/Q2
 965         IF(RPQ.GT.1E-10) THEN
 966           RPBE=SQRT(RPQ*BETA2)
 967           IF(RPBE.LT.0.99) THEN
 968             XBL=LOG((1.+RPBE)/(1.-RPBE))
 969             XBI=2.*RPBE/(1.-RPBE**2)
 970           ELSE
 971             RPBESN=4.*PM2/W2+(4.*X**2*P2/Q2)*BETA2
 972             XBL=LOG((1.+RPBE)**2/RPBESN)
 973             XBI=2.*RPBE/RPBESN
 974           ENDIF
 975           SIGBH=BETA*(6.*X*(1.-X)-1.)+
 976      &    XBL*(X**2+(1.-X)**2+RMQ*X*(1.-3.*X)-0.5*RMQ**2*X**2)+
 977      &    XBI*(2.*X/Q2)*(PM2*X*(2.-RMQ)-P2*X)
 978         ENDIF
 979       ENDIF
 980
 981 C...Multiply by charge-squared etc. to get parton distribution.
 982       CHSQ=1./9.
 983       IF(IABS(KF).EQ.2.OR.IABS(KF).EQ.4) CHSQ=4./9.
 984       XPBH=3.*CHSQ*AEM2PI*X*SIGBH
 985
 986       RETURN
 987       END
 988
 989 C*********************************************************************
 990
 991        SUBROUTINE LHASASDIR(X,Q2,P2,Q02,XPGA)
 992 C...Purpose: to evaluate the direct contribution, i.e. the C^gamma term,
 993 C...as needed in MSbar parametrizations.
 994       DIMENSION XPGA(-6:6)
 995       DATA PMC/1.3/, PMB/4.6/, AEM2PI/0.0011614/
 996
 997 C...Reset output.
 998       DO 100 KFL=-6,6
 999       XPGA(KFL)=0.
1000   100 CONTINUE
1001
1002 C...Evaluate common x-dependent expression.
1003       XTMP = (X**2+(1.-X)**2) * (-LOG(X)) - 1.
1004       CGAM = 3.*AEM2PI*X * (XTMP*(1.+P2/(P2+Q02)) + 6.*X*(1.-X))
1005
1006 C...d, u, s part by simple charge factor.
1007       XPGA(1)=(1./9.)*CGAM
1008       XPGA(2)=(4./9.)*CGAM
1009       XPGA(3)=(1./9.)*CGAM
1010
1011 C...Also fill for antiquarks.
1012       DO 110 KF=1,5
1013       XPGA(-KF)=XPGA(KF)
1014   110 CONTINUE
1015
1016       RETURN
1017       END
1018 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1019        SUBROUTINE LHASASVM1(ISET,KF,X,Q2,P2,ALAM,XPGA)
1020 C...Purpose: to evaluate the VMD parton distributions of a photon,
1021 C...evolved homogeneously from an initial scale P2 to Q2.
1022 C...Does not include dipole suppression factor.
1023 C...ISET is parton distribution set, see above;
1024 C...additionally ISET=0 is used for the evolution of an anomalous photon
1025 C...which branched at a scale P2 and then evolved homogeneously to Q2.
1026 C...ALAM is the 4-flavour Lambda, which is automatically converted
1027 C...to 3- and 5-flavour equivalents as needed.
1028       DIMENSION XPGA(-6:6)
1029       DATA PMC/1.3/, PMB/4.6/, AEM/0.007297/, AEM2PI/0.0011614/
1030
1031 C...Reset output.
1032       DO 100 KFL=-6,6
1033       XPGA(KFL)=0.
1034   100 CONTINUE
1035       KFA=IABS(KF)
1036
1037 C...Calculate Lambda; protect against unphysical Q2 and P2 input.
1038       ALAM3=ALAM*(PMC/ALAM)**(2./27.)
1039       ALAM5=ALAM*(ALAM/PMB)**(2./23.)
1040       P2EFF=MAX(P2,1.2*ALAM3**2)
1041       IF(KFA.EQ.4) P2EFF=MAX(P2EFF,PMC**2)
1042       IF(KFA.EQ.5) P2EFF=MAX(P2EFF,PMB**2)
1043       Q2EFF=MAX(Q2,P2EFF)
1044
1045 C...Find number of flavours at lower and upper scale.
1046       NFP=4
1047       IF(P2EFF.LT.PMC**2) NFP=3
1048       IF(P2EFF.GT.PMB**2) NFP=5
1049       NFQ=4
1050       IF(Q2EFF.LT.PMC**2) NFQ=3
1051       IF(Q2EFF.GT.PMB**2) NFQ=5
1052
1053 C...Find s as sum of 3-, 4- and 5-flavour parts.
1054       S=0.
1055       IF(NFP.EQ.3) THEN
1056         Q2DIV=PMC**2
1057         IF(NFQ.EQ.3) Q2DIV=Q2EFF
1058         S=S+(6./27.)*LOG(LOG(Q2DIV/ALAM3**2)/LOG(P2EFF/ALAM3**2))
1059       ENDIF
1060       IF(NFP.LE.4.AND.NFQ.GE.4) THEN
1061         P2DIV=P2EFF
1062         IF(NFP.EQ.3) P2DIV=PMC**2
1063         Q2DIV=Q2EFF
1064         IF(NFQ.EQ.5) Q2DIV=PMB**2
1065         S=S+(6./25.)*LOG(LOG(Q2DIV/ALAM**2)/LOG(P2DIV/ALAM**2))
1066       ENDIF
1067       IF(NFQ.EQ.5) THEN
1068         P2DIV=PMB**2
1069         IF(NFP.EQ.5) P2DIV=P2EFF
1070         S=S+(6./23.)*LOG(LOG(Q2EFF/ALAM5**2)/LOG(P2DIV/ALAM5**2))
1071       ENDIF
1072
1073 C...Calculate frequent combinations of x and s.
1074       X1=1.-X
1075       XL=-LOG(X)
1076       S2=S**2
1077       S3=S**3
1078       S4=S**4
1079
1080 C...Evaluate homogeneous anomalous parton distributions below or
1081 C...above threshold.
1082       IF(ISET.EQ.0) THEN
1083       IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
1084      &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
1085         XVAL = X * 1.5 * (X**2+X1**2)
1086         XGLU = 0.
1087         XSEA = 0.
1088       ELSE
1089         XVAL = (1.5/(1.-0.197*S+4.33*S2)*X**2 + (1.5+2.10*S)/
1090      &  (1.+3.29*S)*X1**2 + 5.23*S/(1.+1.17*S+19.9*S3)*X*X1) *
1091      &  X**(1./(1.+1.5*S)) * (1.-X**2)**(2.667*S)
1092         XGLU = 4.*S/(1.+4.76*S+15.2*S2+29.3*S4) *
1093      &  X**(-2.03*S/(1.+2.44*S)) * (X1*XL)**(1.333*S) *
1094      &  ((4.*X**2+7.*X+4.)*X1/3. - 2.*X*(1.+X)*XL)
1095         XSEA = S2/(1.+4.54*S+8.19*S2+8.05*S3) *
1096      &  X**(-1.54*S/(1.+1.29*S)) * X1**(2.667*S) *
1097      &  ((8.-73.*X+62.*X**2)*X1/9. + (3.-8.*X**2/3.)*X*XL +
1098      &  (2.*X-1.)*X*XL**2)
1099       ENDIF
1100
1101 C...Evaluate set 1D parton distributions below or above threshold.
1102       ELSEIF(ISET.EQ.1) THEN
1103       IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
1104      &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
1105         XVAL = 1.294 * X**0.80 * X1**0.76
1106         XGLU = 1.273 * X**0.40 * X1**1.76
1107         XSEA = 0.100 * X1**3.76
1108       ELSE
1109         XVAL = 1.294/(1.+0.252*S+3.079*S2) * X**(0.80-0.13*S) *
1110      &  X1**(0.76+0.667*S) * XL**(2.*S)
1111         XGLU = 7.90*S/(1.+5.50*S) * EXP(-5.16*S) *
1112      &  X**(-1.90*S/(1.+3.60*S)) * X1**1.30 * XL**(0.50+3.*S) +
1113      &  1.273 * EXP(-10.*S) * X**0.40 * X1**(1.76+3.*S)
1114         XSEA = (0.1-0.397*S2+1.121*S3)/(1.+5.61*S2+5.26*S3) *
1115      &  X**(-7.32*S2/(1.+10.3*S2)) *
1116      &  X1**((3.76+15.*S+12.*S2)/(1.+4.*S))
1117         XSEA0 = 0.100 * X1**3.76
1118       ENDIF
1119
1120 C...Evaluate set 1M parton distributions below or above threshold.
1121       ELSEIF(ISET.EQ.2) THEN
1122       IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
1123      &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
1124         XVAL = 0.8477 * X**0.51 * X1**1.37
1125         XGLU = 3.42 * X**0.255 * X1**2.37
1126         XSEA = 0.
1127       ELSE
1128         XVAL = 0.8477/(1.+1.37*S+2.18*S2+3.73*S3) * X**(0.51+0.21*S)
1129      &  * X1**1.37 * XL**(2.667*S)
1130         XGLU = 24.*S/(1.+9.6*S+0.92*S2+14.34*S3) * EXP(-5.94*S) *
1131      &  X**((-0.013-1.80*S)/(1.+3.14*S)) * X1**(2.37+0.4*S) *
1132      &  XL**(0.32+3.6*S) + 3.42 * EXP(-12.*S) * X**0.255 *
1133      &  X1**(2.37+3.*S)
1134         XSEA = 0.842*S/(1.+21.3*S-33.2*S2+229.*S3) *
1135      &  X**((0.13-2.90*S)/(1.+5.44*S)) * X1**(3.45+0.5*S) *
1136      &  XL**(2.8*S)
1137         XSEA0 = 0.
1138       ENDIF
1139
1140 C...Evaluate set 2D parton distributions below or above threshold.
1141       ELSEIF(ISET.EQ.3) THEN
1142       IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
1143      &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
1144         XVAL = X**0.46 * X1**0.64 + 0.76 * X
1145         XGLU = 1.925 * X1**2
1146         XSEA = 0.242 * X1**4
1147       ELSE
1148         XVAL = (1.+0.186*S)/(1.-0.209*S+1.495*S2) * X**(0.46+0.25*S)
1149      &  * X1**((0.64+0.14*S+5.*S2)/(1.+S)) * XL**(1.9*S) +
1150      &  (0.76+0.4*S) * X * X1**(2.667*S)
1151         XGLU = (1.925+5.55*S+147.*S2)/(1.-3.59*S+3.32*S2) *
1152      &  EXP(-18.67*S) * X**((-5.81*S-5.34*S2)/(1.+29.*S-4.26*S2))
1153      &  * X1**((2.-5.9*S)/(1.+1.7*S)) * XL**(9.3*S/(1.+1.7*S))
1154         XSEA = (0.242-0.252*S+1.19*S2)/(1.-0.607*S+21.95*S2) *
1155      &  X**(-12.1*S2/(1.+2.62*S+16.7*S2)) * X1**4 * XL**S
1156         XSEA0 = 0.242 * X1**4
1157       ENDIF
1158
1159 C...Evaluate set 2M parton distributions below or above threshold.
1160       ELSEIF(ISET.EQ.4) THEN
1161       IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
1162      &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
1163         XVAL = 1.168 * X**0.50 * X1**2.60 + 0.965 * X
1164         XGLU = 1.808 * X1**2
1165         XSEA = 0.209 * X1**4
1166       ELSE
1167         XVAL = (1.168+1.771*S+29.35*S2) * EXP(-5.776*S) *
1168      &  X**((0.5+0.208*S)/(1.-0.794*S+1.516*S2)) *
1169      &  X1**((2.6+7.6*S)/(1.+5.*S)) * XL**(5.15*S/(1.+2.*S)) +
1170      &  (0.965+22.35*S)/(1.+18.4*S) * X * X1**(2.667*S)
1171         XGLU = (1.808+29.9*S)/(1.+26.4*S) * EXP(-5.28*S) *
1172      &  X**((-5.35*S-10.11*S2)/(1.+31.71*S)) *
1173      &  X1**((2.-7.3*S+4.*S2)/(1.+2.5*S)) *
1174      &  XL**(10.9*S/(1.+2.5*S))
1175         XSEA = (0.209+0.644*S2)/(1.+0.319*S+17.6*S2) *
1176      &  X**((-0.373*S-7.71*S2)/(1.+0.815*S+11.0*S2)) *
1177      &  X1**(4.+S) * XL**(0.45*S)
1178         XSEA0 = 0.209 * X1**4
1179       ENDIF
1180       ENDIF
1181
1182 C...Threshold factors for c and b sea.
1183       SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
1184       XCHM=0.
1185       IF(Q2.GT.PMC**2.AND.Q2.GT.1.001*P2EFF) THEN
1186         SCH=MAX(0.,LOG(LOG(PMC**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
1187         IF(ISET.EQ.0) THEN
1188           XCHM=XSEA*(1.-(SCH/SLL)**2)
1189         ELSE
1190           XCHM=MAX(0.,XSEA-XSEA0*X1**(2.667*S))*(1.-SCH/SLL)
1191         ENDIF
1192       ENDIF
1193       XBOT=0.
1194       IF(Q2.GT.PMB**2.AND.Q2.GT.1.001*P2EFF) THEN
1195         SBT=MAX(0.,LOG(LOG(PMB**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
1196         IF(ISET.EQ.0) THEN
1197           XBOT=XSEA*(1.-(SBT/SLL)**2)
1198         ELSE
1199           XBOT=MAX(0.,XSEA-XSEA0*X1**(2.667*S))*(1.-SBT/SLL)
1200         ENDIF
1201       ENDIF
1202
1203 C...Fill parton distributions.
1204       XPGA(0)=XGLU
1205       XPGA(1)=XSEA
1206       XPGA(2)=XSEA
1207       XPGA(3)=XSEA
1208       XPGA(4)=XCHM
1209       XPGA(5)=XBOT
1210       XPGA(KFA)=XPGA(KFA)+XVAL
1211       DO 110 KFL=1,5
1212       XPGA(-KFL)=XPGA(KFL)
1213   110 CONTINUE
1214
1215       RETURN
1216       END
1217 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1218       SUBROUTINE LHASASAN1(KF,X,Q2,P2,ALAM,XPGA)
1219 C...Purpose: to evaluate the parton distributions of the anomalous
1220 C...photon, inhomogeneously evolved from a scale P2 (where it vanishes)
1221 C...to Q2.
1222 C...KF=0 gives the sum over (up to) 5 flavours,
1223 C...KF<0 limits to flavours up to abs(KF),
1224 C...KF>0 is for flavour KF only.
1225 C...ALAM is the 4-flavour Lambda, which is automatically converted
1226 C...to 3- and 5-flavour equivalents as needed.
1227       DIMENSION XPGA(-6:6),ALAMSQ(3:5)
1228       DATA PMC/1.3/, PMB/4.6/, AEM/0.007297/, AEM2PI/0.0011614/
1229
1230 C...Reset output.
1231       DO 100 KFL=-6,6
1232       XPGA(KFL)=0.
1233   100 CONTINUE
1234       IF(Q2.LE.P2) RETURN
1235       KFA=IABS(KF)
1236
1237 C...Calculate Lambda; protect against unphysical Q2 and P2 input.
1238       ALAMSQ(3)=(ALAM*(PMC/ALAM)**(2./27.))**2
1239       ALAMSQ(4)=ALAM**2
1240       ALAMSQ(5)=(ALAM*(ALAM/PMB)**(2./23.))**2
1241       P2EFF=MAX(P2,1.2*ALAMSQ(3))
1242       IF(KF.EQ.4) P2EFF=MAX(P2EFF,PMC**2)
1243       IF(KF.EQ.5) P2EFF=MAX(P2EFF,PMB**2)
1244       Q2EFF=MAX(Q2,P2EFF)
1245       XL=-LOG(X)
1246
1247 C...Find number of flavours at lower and upper scale.
1248       NFP=4
1249       IF(P2EFF.LT.PMC**2) NFP=3
1250       IF(P2EFF.GT.PMB**2) NFP=5
1251       NFQ=4
1252       IF(Q2EFF.LT.PMC**2) NFQ=3
1253       IF(Q2EFF.GT.PMB**2) NFQ=5
1254
1255 C...Define range of flavour loop.
1256       IF(KF.EQ.0) THEN
1257         KFLMN=1
1258         KFLMX=5
1259       ELSEIF(KF.LT.0) THEN
1260         KFLMN=1
1261         KFLMX=KFA
1262       ELSE
1263         KFLMN=KFA
1264         KFLMX=KFA
1265       ENDIF
1266
1267 C...Loop over flavours the photon can branch into.
1268       DO 110 KFL=KFLMN,KFLMX
1269
1270 C...Light flavours: calculate t range and (approximate) s range.
1271       IF(KFL.LE.3.AND.(KFL.EQ.1.OR.KFL.EQ.KF)) THEN
1272         TDIFF=LOG(Q2EFF/P2EFF)
1273         S=(6./(33.-2.*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
1274      &  LOG(P2EFF/ALAMSQ(NFQ)))
1275         IF(NFQ.GT.NFP) THEN
1276           Q2DIV=PMB**2
1277           IF(NFQ.EQ.4) Q2DIV=PMC**2
1278           SNFQ=(6./(33.-2.*NFQ))*LOG(LOG(Q2DIV/ALAMSQ(NFQ))/
1279      &    LOG(P2EFF/ALAMSQ(NFQ)))
1280           SNFP=(6./(33.-2.*(NFQ-1)))*LOG(LOG(Q2DIV/ALAMSQ(NFQ-1))/
1281      &    LOG(P2EFF/ALAMSQ(NFQ-1)))
1282           S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNFP-SNFQ)
1283         ENDIF
1284         IF(NFQ.EQ.5.AND.NFP.EQ.3) THEN
1285           Q2DIV=PMC**2
1286           SNF4=(6./(33.-2.*4))*LOG(LOG(Q2DIV/ALAMSQ(4))/
1287      &    LOG(P2EFF/ALAMSQ(4)))
1288           SNF3=(6./(33.-2.*3))*LOG(LOG(Q2DIV/ALAMSQ(3))/
1289      &    LOG(P2EFF/ALAMSQ(3)))
1290           S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNF3-SNF4)
1291         ENDIF
1292
1293 C...u and s quark do not need a separate treatment when d has been done.
1294       ELSEIF(KFL.EQ.2.OR.KFL.EQ.3) THEN
1295
1296 C...Charm: as above, but only include range above c threshold.
1297       ELSEIF(KFL.EQ.4) THEN
1298         IF(Q2.LE.PMC**2) GOTO 110
1299         P2EFF=MAX(P2EFF,PMC**2)
1300         Q2EFF=MAX(Q2EFF,P2EFF)
1301         TDIFF=LOG(Q2EFF/P2EFF)
1302         S=(6./(33.-2.*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
1303      &  LOG(P2EFF/ALAMSQ(NFQ)))
1304         IF(NFQ.EQ.5.AND.NFP.EQ.4) THEN
1305           Q2DIV=PMB**2
1306           SNFQ=(6./(33.-2.*NFQ))*LOG(LOG(Q2DIV/ALAMSQ(NFQ))/
1307      &    LOG(P2EFF/ALAMSQ(NFQ)))
1308           SNFP=(6./(33.-2.*(NFQ-1)))*LOG(LOG(Q2DIV/ALAMSQ(NFQ-1))/
1309      &    LOG(P2EFF/ALAMSQ(NFQ-1)))
1310           S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNFP-SNFQ)
1311         ENDIF
1312
1313 C...Bottom: as above, but only include range above b threshold.
1314       ELSEIF(KFL.EQ.5) THEN
1315         IF(Q2.LE.PMB**2) GOTO 110
1316         P2EFF=MAX(P2EFF,PMB**2)
1317         Q2EFF=MAX(Q2,P2EFF)
1318         TDIFF=LOG(Q2EFF/P2EFF)
1319         S=(6./(33.-2.*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
1320      &  LOG(P2EFF/ALAMSQ(NFQ)))
1321       ENDIF
1322
1323 C...Evaluate flavour-dependent prefactor (charge^2 etc.).
1324       CHSQ=1./9.
1325       IF(KFL.EQ.2.OR.KFL.EQ.4) CHSQ=4./9.
1326       FAC=AEM2PI*2.*CHSQ*TDIFF
1327
1328 C...Evaluate parton distributions (normalized to unit momentum sum).
1329       IF(KFL.EQ.1.OR.KFL.EQ.4.OR.KFL.EQ.5.OR.KFL.EQ.KF) THEN
1330         XVAL= ((1.5+2.49*S+26.9*S**2)/(1.+32.3*S**2)*X**2 +
1331      &  (1.5-0.49*S+7.83*S**2)/(1.+7.68*S**2)*(1.-X)**2 +
1332      &  1.5*S/(1.-3.2*S+7.*S**2)*X*(1.-X)) *
1333      &  X**(1./(1.+0.58*S)) * (1.-X**2)**(2.5*S/(1.+10.*S))
1334         XGLU= 2.*S/(1.+4.*S+7.*S**2) *
1335      &  X**(-1.67*S/(1.+2.*S)) * (1.-X**2)**(1.2*S) *
1336      &  ((4.*X**2+7.*X+4.)*(1.-X)/3. - 2.*X*(1.+X)*XL)
1337         XSEA= 0.333*S**2/(1.+4.90*S+4.69*S**2+21.4*S**3) *
1338      &  X**(-1.18*S/(1.+1.22*S)) * (1.-X)**(1.2*S) *
1339      &  ((8.-73.*X+62.*X**2)*(1.-X)/9. + (3.-8.*X**2/3.)*X*XL +
1340      &  (2.*X-1.)*X*XL**2)
1341
1342 C...Threshold factors for c and b sea.
1343         SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
1344         XCHM=0.
1345         IF(Q2.GT.PMC**2.AND.Q2.GT.1.001*P2EFF) THEN
1346           SCH=MAX(0.,LOG(LOG(PMC**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
1347           XCHM=XSEA*(1.-(SCH/SLL)**3)
1348         ENDIF
1349         XBOT=0.
1350         IF(Q2.GT.PMB**2.AND.Q2.GT.1.001*P2EFF) THEN
1351           SBT=MAX(0.,LOG(LOG(PMB**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
1352           XBOT=XSEA*(1.-(SBT/SLL)**3)
1353         ENDIF
1354       ENDIF
1355
1356 C...Add contribution of each valence flavour.
1357       XPGA(0)=XPGA(0)+FAC*XGLU
1358       XPGA(1)=XPGA(1)+FAC*XSEA
1359       XPGA(2)=XPGA(2)+FAC*XSEA
1360       XPGA(3)=XPGA(3)+FAC*XSEA
1361       XPGA(4)=XPGA(4)+FAC*XCHM
1362       XPGA(5)=XPGA(5)+FAC*XBOT
1363       XPGA(KFL)=XPGA(KFL)+FAC*XVAL
1364   110 CONTINUE
1365       DO 120 KFL=1,5
1366       XPGA(-KFL)=XPGA(KFL)
1367   120 CONTINUE
1368
1369       RETURN
1370       END
1371 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc