HERWIG/src/sasgam.f

   1
   2 CDECK  ID>, STRUCTM.
   3
   4 *CMZ :S        E26/04/91  11.11.54  by  Bryan Webber
   5
   6 *-- Author :    Bryan Webber
   7
   8 C-----------------------------------------------------------------------
   9
  10 C      SUBROUTINE STRUCTM(X,QSCA,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU)
  11
  12 C-----------------------------------------------------------------------
  13
  14 C     DUMMY SUBROUTINE: DELETE IF YOU USE PDFLIB CERN-LIBRARY
  15
  16 C     PACKAGE FOR NUCLEON STRUCTURE FUNCTIONS
  17
  18 C-----------------------------------------------------------------------
  19
  20 C      DOUBLE PRECISION X,QSCA,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU
  21
  22 C      WRITE (6,10)
  23
  24 C  10  FORMAT(/10X,'STRUCTM CALLED BUT NOT LINKED')
  25
  26 C      STOP
  27
  28 C      END
  29
  30 C-----------------------------------------------------------------------
  31
  32 C...SaSgam version 2 - parton distributions of the photon
  33
  34 C...by Gerhard A. Schuler and Torbjorn Sjostrand
  35
  36 C...For further information see Z. Phys. C68 (1995) 607
  37
  38 C...and CERN-TH/96-04 and LU TP 96-2.
  39
  40 C...Program last changed on 18 January 1996.
  41
  42 C
  43
  44 C!!!Note that one further call parameter - IP2 - has been added
  45
  46 C!!!to the SASGAM argument list compared with version 1.
  47
  48 C
  49
  50 C...The user should only need to call the SASGAM routine,
  51
  52 C...which in turn calls the auxiliary routines SASVMD, SASANO,
  53
  54 C...SASBEH and SASDIR. The package is self-contained.
  55
  56 C
  57
  58 C...One particular aspect of these parametrizations is that F2 for
  59
  60 C...the photon is not obtained just as the charge-squared-weighted
  61
  62 C...sum of quark distributions, but differ in the treatment of
  63
  64 C...heavy flavours (in F2 the DIS relation W2 = Q2*(1-x)/x restricts
  65
  66 C...the kinematics range of heavy-flavour production, but the same
  67
  68 C...kinematics is not relevant e.g. for jet production) and, for the
  69
  70 C...'MSbar' fits, in the addition of a Cgamma term related to the
  71
  72 C...separation of direct processes. Schematically:
  73
  74 C...PDF = VMD (rho, omega, phi) + anomalous (d, u, s, c, b).
  75
  76 C...F2  = VMD (rho, omega, phi) + anomalous (d, u, s) +
  77
  78 C...      Bethe-Heitler (c, b) (+ Cgamma (d, u, s)).
  79
  80 C...The J/psi and Upsilon states have not been included in the VMD sum,
  81
  82 C...but low c and b masses in the other components should compensate
  83
  84 C...for this in a duality sense.
  85
  86 C
  87
  88 C...The calling sequence is the following:
  89
  90 C     CALL SASGAM(ISET,X,Q2,P2,IP2,F2GM,XPDFGM)
  91
  92 C...with the following declaration statement:
  93
  94 C     DIMENSION XPDFGM(-6:6)
  95
  96 C...and, optionally, further information in:
  97
  98 C     COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
  99
 100 C    &XPDIR(-6:6)
 101
 102 C     COMMON/SASVAL/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),VXPDGM(-6:6)
 103
 104 C...Input:  ISET = 1 : SaS set 1D ('DIS',   Q0 = 0.6 GeV)
 105
 106 C                = 2 : SaS set 1M ('MSbar', Q0 = 0.6 GeV)
 107
 108 C                = 3 : SaS set 2D ('DIS',   Q0 =  2  GeV)
 109
 110 C                = 4 : SaS set 2M ('MSbar', Q0 =  2  GeV)
 111
 112 C           X : x value.
 113
 114 C           Q2 : Q2 value.
 115
 116 C           P2 : P2 value; should be = 0. for an on-shell photon.
 117
 118 C           IP2 : scheme used to evaluate off-shell anomalous component.
 119
 120 C               = 0 : recommended default, see = 7.
 121
 122 C               = 1 : dipole dampening by integration; very time-consuming.
 123
 124 C               = 2 : P_0^2 = max( Q_0^2, P^2 )
 125
 126 C               = 3 : P'_0^2 = Q_0^2 + P^2.
 127
 128 C               = 4 : P_{eff} that preserves momentum sum.
 129
 130 C               = 5 : P_{int} that preserves momentum and average
 131
 132 C                     evolution range.
 133
 134 C               = 6 : P_{eff}, matched to P_0 in P2 -> Q2 limit.
 135
 136 C               = 7 : P_{eff}, matched to P_0 in P2 -> Q2 limit.
 137
 138 C...Output: F2GM : F2 value of the photon (including factors of alpha_em).
 139
 140 C           XPFDGM :  x times parton distribution functions of the photon,
 141
 142 C               with elements 0 = g, 1 = d, 2 = u, 3 = s, 4 = c, 5 = b,
 143
 144 C               6 = t (always empty!), - for antiquarks (result is same).
 145
 146 C...The breakdown by component is stored in the commonblock SASCOM,
 147
 148 C               with elements as above.
 149
 150 C           XPVMD : rho, omega, phi VMD part only of output.
 151
 152 C           XPANL : d, u, s anomalous part only of output.
 153
 154 C           XPANH : c, b anomalous part only of output.
 155
 156 C           XPBEH : c, b Bethe-Heitler part only of output.
 157
 158 C           XPDIR : Cgamma (direct contribution) part only of output.
 159
 160 C...The above arrays do not distinguish valence and sea contributions,
 161
 162 C...although this information is available internally. The additional
 163
 164 C...commonblock SASVAL provides the valence part only of the above
 165
 166 C...distributions. Array names VXPVMD, VXPANL and VXPANH correspond
 167
 168 C...to XPVMD, XPANL and XPANH, while XPBEH and XPDIR are valence only
 169
 170 C...and therefore not given doubly. VXPDGM gives the sum of valence
 171
 172 C...parts, and so matches XPDFGM. The difference, i.e. XPVMD-VXPVMD
 173
 174 C...and so on, gives the sea part only.
 175
 176 C
 177
 178       SUBROUTINE SASGAM(ISET,X,Q2,P2,IP2,F2GM,XPDFGM)
 179
 180 C...Purpose: to construct the F2 and parton distributions of the photon
 181
 182 C...by summing homogeneous (VMD) and inhomogeneous (anomalous) terms.
 183
 184 C...For F2, c and b are included by the Bethe-Heitler formula;
 185
 186 C...in the 'MSbar' scheme additionally a Cgamma term is added.
 187
 188       DIMENSION XPDFGM(-6:6)
 189
 190       COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
 191
 192      &XPDIR(-6:6)
 193
 194       COMMON/SASVAL/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),VXPDGM(-6:6)
 195
 196       SAVE /SASCOM/,/SASVAL/
 197
 198 C
 199
 200 C...Temporary array.
 201
 202       DIMENSION XPGA(-6:6), VXPGA(-6:6)
 203
 204 C...Charm and bottom masses (low to compensate for J/psi etc.).
 205
 206       DATA PMC/1.3/, PMB/4.6/
 207
 208 C...alpha_em and alpha_em/(2*pi).
 209
 210       DATA AEM/0.007297/, AEM2PI/0.0011614/
 211
 212 C...Lambda value for 4 flavours.
 213
 214       DATA ALAM/0.20/
 215
 216 C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
 217
 218       DATA FRACU/0.8/
 219
 220 C...VMD couplings f_V**2/(4*pi).
 221
 222       DATA FRHO/2.20/, FOMEGA/23.6/, FPHI/18.4/
 223
 224 C...Masses for rho (=omega) and phi.
 225
 226       DATA PMRHO/0.770/, PMPHI/1.020/
 227
 228 C...Number of points in integration for IP2=1.
 229
 230       DATA NSTEP/100/
 231
 232 C
 233
 234 C...Reset output.
 235
 236       F2GM=0.
 237
 238       DO 100 KFL=-6,6
 239
 240       XPDFGM(KFL)=0.
 241
 242       XPVMD(KFL)=0.
 243
 244       XPANL(KFL)=0.
 245
 246       XPANH(KFL)=0.
 247
 248       XPBEH(KFL)=0.
 249
 250       XPDIR(KFL)=0.
 251
 252       VXPVMD(KFL)=0.
 253
 254       VXPANL(KFL)=0.
 255
 256       VXPANH(KFL)=0.
 257
 258       VXPDGM(KFL)=0.
 259
 260   100 CONTINUE
 261
 262 C
 263
 264 C...Check that input sensible.
 265
 266       IF(ISET.LE.0.OR.ISET.GE.5) THEN
 267
 268         WRITE(*,*) ' FATAL ERROR: SaSgam called for unknown set'
 269
 270         WRITE(*,*) ' ISET = ',ISET
 271
 272         STOP
 273
 274       ENDIF
 275
 276       IF(X.LE.0..OR.X.GT.1.) THEN
 277
 278         WRITE(*,*) ' FATAL ERROR: SaSgam called for unphysical x'
 279
 280         WRITE(*,*) ' X = ',X
 281
 282         STOP
 283
 284       ENDIF
 285
 286 C
 287
 288 C...Set Q0 cut-off parameter as function of set used.
 289
 290       IF(ISET.LE.2) THEN
 291
 292         Q0=0.6
 293
 294       ELSE
 295
 296         Q0=2.
 297
 298       ENDIF
 299
 300       Q02=Q0**2
 301
 302 C
 303
 304 C...Scale choice for off-shell photon; common factors.
 305
 306       Q2A=Q2
 307
 308       FACNOR=1.
 309
 310       IF(IP2.EQ.1) THEN
 311
 312         P2MX=P2+Q02
 313
 314         Q2A=Q2+P2*Q02/MAX(Q02,Q2)
 315
 316         FACNOR=LOG(Q2/Q02)/NSTEP
 317
 318       ELSEIF(IP2.EQ.2) THEN
 319
 320         P2MX=MAX(P2,Q02)
 321
 322       ELSEIF(IP2.EQ.3) THEN
 323
 324         P2MX=P2+Q02
 325
 326         Q2A=Q2+P2*Q02/MAX(Q02,Q2)
 327
 328       ELSEIF(IP2.EQ.4) THEN
 329
 330         P2MX=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 331
 332      &  ((Q2+P2)*(Q02+P2)))
 333
 334       ELSEIF(IP2.EQ.5) THEN
 335
 336         P2MXA=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 337
 338      &  ((Q2+P2)*(Q02+P2)))
 339
 340         P2MX=Q0*SQRT(P2MXA)
 341
 342         FACNOR=LOG(Q2/P2MXA)/LOG(Q2/P2MX)
 343
 344       ELSEIF(IP2.EQ.6) THEN
 345
 346         P2MX=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 347
 348      &  ((Q2+P2)*(Q02+P2)))
 349
 350         P2MX=MAX(0.,1.-P2/Q2)*P2MX+MIN(1.,P2/Q2)*MAX(P2,Q02)
 351
 352       ELSE
 353
 354         P2MXA=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
 355
 356      &  ((Q2+P2)*(Q02+P2)))
 357
 358         P2MX=Q0*SQRT(P2MXA)
 359
 360         P2MXB=P2MX
 361
 362         P2MX=MAX(0.,1.-P2/Q2)*P2MX+MIN(1.,P2/Q2)*MAX(P2,Q02)
 363
 364         P2MXB=MAX(0.,1.-P2/Q2)*P2MXB+MIN(1.,P2/Q2)*P2MXA
 365
 366         FACNOR=LOG(Q2/P2MXA)/LOG(Q2/P2MXB)
 367
 368       ENDIF
 369
 370 C
 371
 372 C...Call VMD parametrization for d quark and use to give rho, omega,
 373
 374 C...phi. Note dipole dampening for off-shell photon.
 375
 376       CALL SASVMD(ISET,1,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 377
 378       XFVAL=VXPGA(1)
 379
 380       XPGA(1)=XPGA(2)
 381
 382       XPGA(-1)=XPGA(-2)
 383
 384       FACUD=AEM*(1./FRHO+1./FOMEGA)*(PMRHO**2/(PMRHO**2+P2))**2
 385
 386       FACS=AEM*(1./FPHI)*(PMPHI**2/(PMPHI**2+P2))**2
 387
 388       DO 110 KFL=-5,5
 389
 390       XPVMD(KFL)=(FACUD+FACS)*XPGA(KFL)
 391
 392   110 CONTINUE
 393
 394       XPVMD(1)=XPVMD(1)+(1.-FRACU)*FACUD*XFVAL
 395
 396       XPVMD(2)=XPVMD(2)+FRACU*FACUD*XFVAL
 397
 398       XPVMD(3)=XPVMD(3)+FACS*XFVAL
 399
 400       XPVMD(-1)=XPVMD(-1)+(1.-FRACU)*FACUD*XFVAL
 401
 402       XPVMD(-2)=XPVMD(-2)+FRACU*FACUD*XFVAL
 403
 404       XPVMD(-3)=XPVMD(-3)+FACS*XFVAL
 405
 406       VXPVMD(1)=(1.-FRACU)*FACUD*XFVAL
 407
 408       VXPVMD(2)=FRACU*FACUD*XFVAL
 409
 410       VXPVMD(3)=FACS*XFVAL
 411
 412       VXPVMD(-1)=(1.-FRACU)*FACUD*XFVAL
 413
 414       VXPVMD(-2)=FRACU*FACUD*XFVAL
 415
 416       VXPVMD(-3)=FACS*XFVAL
 417
 418 C
 419
 420       IF(IP2.NE.1) THEN
 421
 422 C...Anomalous parametrizations for different strategies
 423
 424 C...for off-shell photons; except full integration.
 425
 426 C
 427
 428 C...Call anomalous parametrization for d + u + s.
 429
 430         CALL SASANO(-3,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 431
 432         DO 120 KFL=-5,5
 433
 434         XPANL(KFL)=FACNOR*XPGA(KFL)
 435
 436         VXPANL(KFL)=FACNOR*VXPGA(KFL)
 437
 438   120   CONTINUE
 439
 440 C
 441
 442 C...Call anomalous parametrization for c and b.
 443
 444         CALL SASANO(4,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 445
 446         DO 130 KFL=-5,5
 447
 448         XPANH(KFL)=FACNOR*XPGA(KFL)
 449
 450         VXPANH(KFL)=FACNOR*VXPGA(KFL)
 451
 452   130   CONTINUE
 453
 454         CALL SASANO(5,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
 455
 456         DO 140 KFL=-5,5
 457
 458         XPANH(KFL)=XPANH(KFL)+FACNOR*XPGA(KFL)
 459
 460         VXPANH(KFL)=VXPANH(KFL)+FACNOR*VXPGA(KFL)
 461
 462   140   CONTINUE
 463
 464 C
 465
 466       ELSE
 467
 468 C...Special option: loop over flavours and integrate over k2.
 469
 470         DO 170 KF=1,5
 471
 472         DO 160 ISTEP=1,NSTEP
 473
 474         Q2STEP=Q02*(Q2/Q02)**((ISTEP-0.5)/NSTEP)
 475
 476         IF((KF.EQ.4.AND.Q2STEP.LT.PMC**2).OR.
 477
 478      &  (KF.EQ.5.AND.Q2STEP.LT.PMB**2)) GOTO 160
 479
 480         CALL SASVMD(0,KF,X,Q2,Q2STEP,ALAM,XPGA,VXPGA)
 481
 482         FACQ=AEM2PI*(Q2STEP/(Q2STEP+P2))**2*FACNOR
 483
 484         IF(MOD(KF,2).EQ.0) FACQ=FACQ*(8./9.)
 485
 486         IF(MOD(KF,2).EQ.1) FACQ=FACQ*(2./9.)
 487
 488         DO 150 KFL=-5,5
 489
 490         IF(KF.LE.3) XPANL(KFL)=XPANL(KFL)+FACQ*XPGA(KFL)
 491
 492         IF(KF.GE.4) XPANH(KFL)=XPANH(KFL)+FACQ*XPGA(KFL)
 493
 494         IF(KF.LE.3) VXPANL(KFL)=VXPANL(KFL)+FACQ*VXPGA(KFL)
 495
 496         IF(KF.GE.4) VXPANH(KFL)=VXPANH(KFL)+FACQ*VXPGA(KFL)
 497
 498   150   CONTINUE
 499
 500   160   CONTINUE
 501
 502   170   CONTINUE
 503
 504       ENDIF
 505
 506 C
 507
 508 C...Call Bethe-Heitler term expression for charm and bottom.
 509
 510       CALL SASBEH(4,X,Q2,P2,PMC**2,XPBH)
 511
 512       XPBEH(4)=XPBH
 513
 514       XPBEH(-4)=XPBH
 515
 516       CALL SASBEH(5,X,Q2,P2,PMB**2,XPBH)
 517
 518       XPBEH(5)=XPBH
 519
 520       XPBEH(-5)=XPBH
 521
 522 C
 523
 524 C...For MSbar subtraction call C^gamma term expression for d, u, s.
 525
 526       IF(ISET.EQ.2.OR.ISET.EQ.4) THEN
 527
 528         CALL SASDIR(X,Q2,P2,Q02,XPGA)
 529
 530         DO 180 KFL=-5,5
 531
 532         XPDIR(KFL)=XPGA(KFL)
 533
 534   180   CONTINUE
 535
 536       ENDIF
 537
 538 C
 539
 540 C...Store result in output array.
 541
 542       DO 190 KFL=-5,5
 543
 544       CHSQ=1./9.
 545
 546       IF(IABS(KFL).EQ.2.OR.IABS(KFL).EQ.4) CHSQ=4./9.
 547
 548       XPF2=XPVMD(KFL)+XPANL(KFL)+XPBEH(KFL)+XPDIR(KFL)
 549
 550       IF(KFL.NE.0) F2GM=F2GM+CHSQ*XPF2
 551
 552       XPDFGM(KFL)=XPVMD(KFL)+XPANL(KFL)+XPANH(KFL)
 553
 554       VXPDGM(KFL)=VXPVMD(KFL)+VXPANL(KFL)+VXPANH(KFL)
 555
 556   190 CONTINUE
 557
 558 C
 559
 560       RETURN
 561
 562       END