| 4e9e3152 |
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 |
git://git.uio.no / u / mrichter / AliRoot.git / blame