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
8 real*8 upv,dnv,usea,dsea,str,chm,bot,top,glu
12 call getnmem(iset,imem)
15 if(iimem.eq.0) iimem=6
17 call SFSASxx(iimem,xin,Q2in,p2in,ip2,
18 + upv,dnv,usea,dsea,str,chm,bot,top,glu)
36 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
38 c print *,'calling SASGread'
39 read(1,*)nmem(nset),ndef(nset)
42 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
43 entry SASGalfa(alfas,qalfa)
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)
53 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
54 entry SASGinit(Eorder,Q2fit)
57 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
59 c print *,'calling SASGpdf',mem
62 call setnmem(iset,mem)
68 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
69 C-----------------------------------------------------------------------
70 SUBROUTINE SFSASxx(iset,DX,DQ2,DP2,ip2,
71 + DUPV,DDNV,DSEA,DSEAD,DSTR,DCHM,DBOT,DTOP,DGL)
73 C ********************************************************************
75 C * Interface to SASset of structure functions *
77 C * Author: H. Plothow-Besch (CERN-PPE) *
79 C ********************************************************************
81 C :::::::::::: Structure functions from the SAS group version 2
82 C :::::::::::: Lambda = 0.200 GeV, Q**2 = 0.36 GeV**2 (DIS)
86 + DUPV,DDNV,DSEA,DSEAD,DSTR,DCHM,DBOT,DTOP,DGL
87 DIMENSION XPDFGM(-6:6)
88 REAL X, Q, Q2, P2, F2GAM, XPDFGM
96 C generate the individual structure fcn calls
100 CALL LHASASGAM1(iISET,X,Q2,P2,F2GAM,XPDFGM)
103 CALL LHASASGAM2(iISET,X,Q2,P2,ip2,F2GAM,XPDFGM)
121 C IF (DSCAL.GT.TMAS) TOP = XPDFGM(6)
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.
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.
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.
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),
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)
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.
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),
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.
199 C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
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/
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
223 IF(X.LE.0..OR.X.GT.1.) THEN
224 WRITE(*,*) ' FATAL ERROR: SaSgam called for unphysical x'
229 C...Set Q0 cut-off parameter as function of set used.
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.
240 CALL LHASASVM1(ISET,1,X,Q2,P2MX,ALAM,XPGA)
241 XFVAL=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
247 XPVMD(KFL)=(FACUD+FACS)*XPGA(KFL)
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
256 C...Call anomalous parametrization for d + u + s.
257 CALL LHASASAN1(-3,X,Q2,P2MX,ALAM,XPGA)
262 C...Call anomalous parametrization for c and b.
263 CALL LHASASAN1(4,X,Q2,P2MX,ALAM,XPGA)
267 CALL LHASASAN1(5,X,Q2,P2MX,ALAM,XPGA)
269 XPANH(KFL)=XPANH(KFL)+XPGA(KFL)
272 C...Call Bethe-Heitler term expression for charm and bottom.
273 CALL LHASASBEH(4,X,Q2,P2,PMC**2,XPBH)
276 CALL LHASASBEH(5,X,Q2,P2,PMB**2,XPBH)
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)
288 C...Store result in output array.
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)
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.
307 C!!!Note that one further call parameter - IP2 - has been added
308 C!!!to the SASGAM argument list compared with version 1.
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.
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.
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),
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)
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
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.
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),
382 COMMON/LHASASVAL/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),
384 SAVE /LHASASCOM/,/LHASASVAL/
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.
394 C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
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.
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
424 IF(X.LE.0..OR.X.GT.1.) THEN
425 WRITE(*,*) ' FATAL ERROR: SaSgam called for unphysical x'
430 C...Set Q0 cut-off parameter as function of set used.
438 C...Scale choice for off-shell photon; common factors.
443 Q2A=Q2+P2*Q02/MAX(Q02,Q2)
444 FACNOR=LOG(Q2/Q02)/NSTEP
445 ELSEIF(IP2.EQ.2) THEN
447 ELSEIF(IP2.EQ.3) THEN
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)))
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)
463 P2MXA=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
464 & ((Q2+P2)*(Q02+P2)))
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)
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)
478 FACUD=AEM*(1./FRHO+1./FOMEGA)*(PMRHO**2/(PMRHO**2+P2))**2
479 FACS=AEM*(1./FPHI)*(PMPHI**2/(PMPHI**2+P2))**2
481 XPVMD(KFL)=(FACUD+FACS)*XPGA(KFL)
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
492 VXPVMD(-1)=(1.-FRACU)*FACUD*XFVAL
493 VXPVMD(-2)=FRACU*FACUD*XFVAL
494 VXPVMD(-3)=FACS*XFVAL
497 C...Anomalous parametrizations for different strategies
498 C...for off-shell photons; except full integration.
500 C...Call anomalous parametrization for d + u + s.
501 CALL LHASASANO(-3,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
503 XPANL(KFL)=FACNOR*XPGA(KFL)
504 VXPANL(KFL)=FACNOR*VXPGA(KFL)
507 C...Call anomalous parametrization for c and b.
508 CALL LHASASANO(4,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
510 XPANH(KFL)=FACNOR*XPGA(KFL)
511 VXPANH(KFL)=FACNOR*VXPGA(KFL)
513 CALL LHASASANO(5,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
515 XPANH(KFL)=XPANH(KFL)+FACNOR*XPGA(KFL)
516 VXPANH(KFL)=VXPANH(KFL)+FACNOR*VXPGA(KFL)
520 C...Special option: loop over flavours and integrate over k2.
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.)
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)
540 C...Call Bethe-Heitler term expression for charm and bottom.
541 CALL LHASASBEH(4,X,Q2,P2,PMC**2,XPBH)
544 CALL LHASASBEH(5,X,Q2,P2,PMB**2,XPBH)
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)
556 C...Store result in output array.
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)
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/
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)
598 C...Find number of flavours at lower and upper scale.
600 IF(P2EFF.LT.PMC**2) NFP=3
601 IF(P2EFF.GT.PMB**2) NFP=5
603 IF(Q2EFF.LT.PMC**2) NFQ=3
604 IF(Q2EFF.GT.PMB**2) NFQ=5
606 C...Find s as sum of 3-, 4- and 5-flavour parts.
610 IF(NFQ.EQ.3) Q2DIV=Q2EFF
611 S=S+(6./27.)*LOG(LOG(Q2DIV/ALAM3**2)/LOG(P2EFF/ALAM3**2))
613 IF(NFP.LE.4.AND.NFQ.GE.4) THEN
615 IF(NFP.EQ.3) P2DIV=PMC**2
617 IF(NFQ.EQ.5) Q2DIV=PMB**2
618 S=S+(6./25.)*LOG(LOG(Q2DIV/ALAM**2)/LOG(P2DIV/ALAM**2))
622 IF(NFP.EQ.5) P2DIV=P2EFF
623 S=S+(6./23.)*LOG(LOG(Q2EFF/ALAM5**2)/LOG(P2DIV/ALAM5**2))
626 C...Calculate frequent combinations of x and s.
633 C...Evaluate homogeneous anomalous parton distributions below or
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)
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 +
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
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
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
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 *
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) *
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
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
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
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
735 C...Threshold factors for c and b sea.
736 SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
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)))
741 XCHM=XSEA*(1.-(SCH/SLL)**2)
743 XCHM=MAX(0.,XSEA-XSEA0*X1**(2.667*S))*(1.-SCH/SLL)
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)))
750 XBOT=XSEA*(1.-(SBT/SLL)**2)
752 XBOT=MAX(0.,XSEA-XSEA0*X1**(2.667*S))*(1.-SBT/SLL)
756 C...Fill parton distributions.
763 XPGA(KFA)=XPGA(KFA)+XVAL
773 C*********************************************************************
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)
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/
795 C...Calculate Lambda; protect against unphysical Q2 and P2 input.
796 ALAMSQ(3)=(ALAM*(PMC/ALAM)**(2./27.))**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)
805 C...Find number of flavours at lower and upper scale.
807 IF(P2EFF.LT.PMC**2) NFP=3
808 IF(P2EFF.GT.PMB**2) NFP=5
810 IF(Q2EFF.LT.PMC**2) NFQ=3
811 IF(Q2EFF.GT.PMB**2) NFQ=5
813 C...Define range of flavour loop.
825 C...Loop over flavours the photon can branch into.
826 DO 110 KFL=KFLMN,KFLMX
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)))
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)
842 IF(NFQ.EQ.5.AND.NFP.EQ.3) THEN
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)
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
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
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)
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)
876 TDIFF=LOG(Q2EFF/P2EFF)
877 S=(6./(33.-2.*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
878 & LOG(P2EFF/ALAMSQ(NFQ)))
881 C...Evaluate flavour-dependent prefactor (charge^2 etc.).
883 IF(KFL.EQ.2.OR.KFL.EQ.4) CHSQ=4./9.
884 FAC=AEM2PI*2.*CHSQ*TDIFF
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 +
900 C...Threshold factors for c and b sea.
901 SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
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)
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)
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
926 VXPGA(-KFL)=VXPGA(KFL)
932 C*********************************************************************
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/
943 C...Check kinematics limits.
944 IF(X.GE.Q2/(4.*PM2+Q2+P2)) RETURN
947 IF(BETA2.LT.1E-10) RETURN
951 C...Simple case: P2 = 0.
953 IF(BETA.LT.0.99) THEN
954 XBL=LOG((1.+BETA)/(1.-BETA))
956 XBL=LOG((1.+BETA)**2*W2/(4.*PM2))
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)
961 C...Complicated case: P2 > 0, based on approximation of
962 C...C.T. Hill and G.G. Ross, Nucl. Phys. B148 (1979) 373
965 IF(RPQ.GT.1E-10) THEN
967 IF(RPBE.LT.0.99) THEN
968 XBL=LOG((1.+RPBE)/(1.-RPBE))
969 XBI=2.*RPBE/(1.-RPBE**2)
971 RPBESN=4.*PM2/W2+(4.*X**2*P2/Q2)*BETA2
972 XBL=LOG((1.+RPBE)**2/RPBESN)
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)
981 C...Multiply by charge-squared etc. to get parton distribution.
983 IF(IABS(KF).EQ.2.OR.IABS(KF).EQ.4) CHSQ=4./9.
984 XPBH=3.*CHSQ*AEM2PI*X*SIGBH
989 C*********************************************************************
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.
995 DATA PMC/1.3/, PMB/4.6/, AEM2PI/0.0011614/
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))
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
1011 C...Also fill for antiquarks.
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/
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)
1045 C...Find number of flavours at lower and upper scale.
1047 IF(P2EFF.LT.PMC**2) NFP=3
1048 IF(P2EFF.GT.PMB**2) NFP=5
1050 IF(Q2EFF.LT.PMC**2) NFQ=3
1051 IF(Q2EFF.GT.PMB**2) NFQ=5
1053 C...Find s as sum of 3-, 4- and 5-flavour parts.
1057 IF(NFQ.EQ.3) Q2DIV=Q2EFF
1058 S=S+(6./27.)*LOG(LOG(Q2DIV/ALAM3**2)/LOG(P2EFF/ALAM3**2))
1060 IF(NFP.LE.4.AND.NFQ.GE.4) THEN
1062 IF(NFP.EQ.3) P2DIV=PMC**2
1064 IF(NFQ.EQ.5) Q2DIV=PMB**2
1065 S=S+(6./25.)*LOG(LOG(Q2DIV/ALAM**2)/LOG(P2DIV/ALAM**2))
1069 IF(NFP.EQ.5) P2DIV=P2EFF
1070 S=S+(6./23.)*LOG(LOG(Q2EFF/ALAM5**2)/LOG(P2DIV/ALAM5**2))
1073 C...Calculate frequent combinations of x and s.
1080 C...Evaluate homogeneous anomalous parton distributions below or
1081 C...above threshold.
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)
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)
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
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
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
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 *
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) *
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
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
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
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
1182 C...Threshold factors for c and b sea.
1183 SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
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)))
1188 XCHM=XSEA*(1.-(SCH/SLL)**2)
1190 XCHM=MAX(0.,XSEA-XSEA0*X1**(2.667*S))*(1.-SCH/SLL)
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)))
1197 XBOT=XSEA*(1.-(SBT/SLL)**2)
1199 XBOT=MAX(0.,XSEA-XSEA0*X1**(2.667*S))*(1.-SBT/SLL)
1203 C...Fill parton distributions.
1210 XPGA(KFA)=XPGA(KFA)+XVAL
1212 XPGA(-KFL)=XPGA(KFL)
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)
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/
1237 C...Calculate Lambda; protect against unphysical Q2 and P2 input.
1238 ALAMSQ(3)=(ALAM*(PMC/ALAM)**(2./27.))**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)
1247 C...Find number of flavours at lower and upper scale.
1249 IF(P2EFF.LT.PMC**2) NFP=3
1250 IF(P2EFF.GT.PMB**2) NFP=5
1252 IF(Q2EFF.LT.PMC**2) NFQ=3
1253 IF(Q2EFF.GT.PMB**2) NFQ=5
1255 C...Define range of flavour loop.
1259 ELSEIF(KF.LT.0) THEN
1267 C...Loop over flavours the photon can branch into.
1268 DO 110 KFL=KFLMN,KFLMX
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)))
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)
1284 IF(NFQ.EQ.5.AND.NFP.EQ.3) THEN
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)
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
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
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)
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)
1318 TDIFF=LOG(Q2EFF/P2EFF)
1319 S=(6./(33.-2.*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
1320 & LOG(P2EFF/ALAMSQ(NFQ)))
1323 C...Evaluate flavour-dependent prefactor (charge^2 etc.).
1325 IF(KFL.EQ.2.OR.KFL.EQ.4) CHSQ=4./9.
1326 FAC=AEM2PI*2.*CHSQ*TDIFF
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)
1342 C...Threshold factors for c and b sea.
1343 SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
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)
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)
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
1366 XPGA(-KFL)=XPGA(KFL)
1371 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc