LHAPDF 5.2.2 source code.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / 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