HERWIG/src/hwhbsg.f

   1
   2 CDECK  ID>, HWHBSG.
   3
   4 *CMZ :-        -03/07/95  19.02.12  by  Giovanni Abbiendi
   5
   6 *-- Author :    Giovanni Abbiendi & Luca Stanco
   7
   8 C----------------------------------------------------------------------
   9
  10       SUBROUTINE HWHBSG
  11
  12 C----------------------------------------------------------------------
  13
  14 C     Returns differential cross section DSIGMA in (Y,Q2,ETA,Z,PHI)
  15
  16 C     Scale for structure functions and alpha_s selected by BGSHAT
  17
  18 C----------------------------------------------------------------------
  19
  20       INCLUDE 'HERWIG61.INC'
  21
  22       DOUBLE PRECISION HWUALF,HWUAEM,LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,
  23
  24      & ME,MP,ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT,
  25
  26      & SFUN(13),ALPHA,LDSIG,DLQ(7),SG,XG,MF1,MF2,MSUM,MDIF,MPRO,FFUN,
  27
  28      & GFUN,H43,H41,H11,H12,H14,H16,H21,H22,G11,G12,G1A,G1B,G21,G22,G3,
  29
  30      & GC,A11,A12,A44,ALPHAS,PDENS,AFACT,BFACT,CFACT,DFACT,GAMMA,S,T,U,
  31
  32      & MREMIN,POL,CCOL,ETA
  33
  34       INTEGER IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,IPROO,IHAD,ILEPT,IQ,IS
  35
  36       LOGICAL CHARGD,INCLUD(18),INSIDE(18)
  37
  38       EXTERNAL HWUALF,HWUAEM
  39
  40       COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF,
  41
  42      & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,
  43
  44      & IPROO,CHARGD,INCLUD,INSIDE
  45
  46 C
  47
  48       IHAD=2
  49
  50       IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD)
  51
  52 C---set masses
  53
  54       IF (CHARGD) THEN
  55
  56         MREMIN=MP
  57
  58         IF (LEP.EQ.ONE) THEN
  59
  60           MF1=RMASS(IFLAVD)
  61
  62           MF2=RMASS(IFLAVU)
  63
  64         ELSE
  65
  66           MF1=RMASS(IFLAVU)
  67
  68           MF2=RMASS(IFLAVD)
  69
  70         ENDIF
  71
  72       ELSE
  73
  74         IS=IFL
  75
  76         IF (IFL.EQ.164) IS=IQK
  77
  78         MREMIN = MREMIF(IS)
  79
  80         MF1 = MFIN1(IS)
  81
  82         MF2 = MFIN2(IS)
  83
  84       ENDIF
  85
  86 C---choose subprocess scale
  87
  88       IF (BGSHAT) THEN
  89
  90         EMSCA = RSHAT
  91
  92       ELSE
  93
  94         S=SHAT+Q2
  95
  96         IF (IFL.GE.7.AND.IFL.LE.18) S=SHAT+Q2-MF1**2
  97
  98         T=-S*Z
  99
 100         U=-S-T
 101
 102         IF (IFL.GE.7.AND.IFL.LE.18) U=-S-T-2*MF1**2
 103
 104         EMSCA = SQRT(TWO*S*T*U/(S**2+T**2+U**2))
 105
 106         IF (IFL.EQ.164) EMSCA=SQRT(-U)
 107
 108       ENDIF
 109
 110       ALPHAS = HWUALF(1,EMSCA)
 111
 112       IF (ALPHAS.GE.ONE.OR.ALPHAS.LE.ZERO) CALL HWWARN('HWHBSG',51,*888)
 113
 114 C---structure functions
 115
 116       ETA = (SHAT+Q2)/SMA/Y
 117
 118       IF (ETA.GT.ONE) ETA=ONE
 119
 120       CALL HWSFUN (ETA,EMSCA,IDHW(IHAD),NSTRU,SFUN,2)
 121
 122       XG = Q2/(SHAT + Q2)
 123
 124       SG = ETA*SMA
 125
 126       IF (SG.LE.(RSHAT+ML)**2.OR.SG.GE.(RS-MREMIN)**2) GOTO 888
 127
 128 C
 129
 130       IF (IFL.EQ.164) GOTO 200
 131
 132 C
 133
 134 C---Electroweak couplings
 135
 136       ALPHA=HWUAEM(-Q2)
 137
 138       IF (CHARGD) THEN
 139
 140         POL = PPOLN(3) - EPOLN(3)
 141
 142         DLQ(1)=.0625*VCKM(IFLAVU/2,(IFLAVD+1)/2)/SWEIN**2 *
 143
 144      +         Q2**2/((Q2+RMASS(198)**2)**2+(RMASS(198)*GAMW)**2) *
 145
 146      +         (ONE + POL)
 147
 148         DLQ(2)=ZERO
 149
 150         DLQ(3)=DLQ(1)
 151
 152       ELSE
 153
 154         IQ=MOD(IFL-1,6)+1
 155
 156         ILEPT=MOD(IDHW(1)-121,6)+11
 157
 158         CALL HWUCFF(ILEPT,IQ,-Q2,DLQ(1))
 159
 160       ENDIF
 161
 162 C
 163
 164       IF (IFL.LE.6) THEN
 165
 166 C---For Boson-Gluon Fusion
 167
 168         PDENS = SFUN(13)/ETA
 169
 170         CCOL = HALF
 171
 172         MSUM = (MF1**2 + MF2**2) / (Y*SG)
 173
 174         MDIF = (MF1**2 - MF2**2) / (Y*SG)
 175
 176         MPRO = MF1*MF2 / (Y*SG)
 177
 178 C
 179
 180         FFUN = (1.D0-XG)*Z*(1.D0-Z) + (MDIF*(2.D0*Z-1.D0)-MSUM)/2.D0
 181
 182         GFUN = (1.D0-XG)*(1.D0-Z) + XG*Z + MDIF
 183
 184         IF ( FFUN .LT. ZERO ) FFUN = ZERO
 185
 186         H43 = (8.D0*(2.D0*Z**2*XG-Z**2-2.D0*Z*XG+2.D0*Z*MDIF+Z-MDIF
 187
 188      &         -MSUM)) / (Z*(1.D0-Z))**2
 189
 190 C
 191
 192         H41 = (8.D0*(Z**2-Z*XG+Z*MDIF-MDIF-MSUM)) / (Z**2*(1.D0-Z))
 193
 194 C
 195
 196         H11 = (4.D0*(2.D0*Z**4-4.D0*Z**3+2.D0*Z**2*MSUM*XG
 197
 198      &         -2.D0*Z**2*MSUM+2.D0*Z**2*XG**2-2.D0*Z**2*XG+3.D0*Z**2
 199
 200      &         +2.D0*Z*MDIF*MSUM+2.D0*Z*MDIF*XG-2.D0*Z*MSUM*XG
 201
 202      &         +2.D0*Z*MSUM-2.D0*Z*XG**2+2.D0*Z*XG-Z-MDIF*MSUM-MDIF*XG
 203
 204      &         -MSUM**2-MSUM*XG)) / (Z*(1.D0-Z))**2
 205
 206 C
 207
 208         H12 = (16.D0*(-Z*MDIF+Z*XG+MDIF+MSUM))/(Z**2*(1.D0-Z))
 209
 210 C
 211
 212         H14 = (16.D0*(-2.D0*Z**2*XG-2.D0*Z*MDIF+2.D0*Z*XG+MDIF+MSUM))
 213
 214      &        / (Z*(1.D0-Z))**2
 215
 216 C
 217
 218         H16 = (32.D0*(Z*MDIF-Z*XG-MDIF-MSUM)) / (Z**2*(1.D0-Z))
 219
 220 C
 221
 222         H21 = (8.D0*MPRO*(-2.D0*Z**2*XG+2.D0*Z**2-2.D0*Z*MDIF+2.D0*Z*XG
 223
 224      +         -2.D0*Z+MDIF+MSUM)) / (Z*(1.D0-Z))**2
 225
 226 C
 227
 228         H22 = (-32.D0*MPRO) / (Z*(1.D0-Z))
 229
 230 C
 231
 232         G11 = -2.D0*H11 + FFUN*H14
 233
 234         G12 = 2.D0*XG*FFUN*H14 + H12 + GFUN * ( H16+GFUN*H14 )
 235
 236         G1A = SQRT( XG*FFUN ) * ( H16 + 2.D0*GFUN*H14 )
 237
 238         G1B = FFUN*H14
 239
 240         G21 = -2.D0*H21
 241
 242         G22 = H22
 243
 244         G3  = H41 - GFUN*H43
 245
 246         GC  = SQRT( XG*FFUN ) * (-2.D0*XG*H43 )
 247
 248       ELSE
 249
 250 C---for QCD Compton, massless matrix element
 251
 252         PDENS = SFUN(IFL-6)/ETA
 253
 254         CCOL = CFFAC
 255
 256         FFUN = XG*(ONE-XG)*Z*(ONE-Z)
 257
 258         GFUN = (ONE-XG)*(ONE-Z)+XG*Z
 259
 260         G11 = 8.D0*((Z**2+XG**2)/(ONE-XG)/(ONE-Z)+TWO*(XG*Z+ONE))
 261
 262         G12 = 64.D0*XG**2*Z+TWO*XG*G11
 263
 264         G1A = 32.D0*XG*GFUN*SQRT(FFUN)/((ONE-XG)*(ONE-Z))
 265
 266         G1B = 16.D0*XG*Z
 267
 268         G3  = -16.D0*(ONE-XG)*(ONE-Z)+G11
 269
 270         GC  = -16.D0*XG*SQRT(FFUN)*(ONE-Z-XG)/((ONE-XG)*(ONE-Z))
 271
 272         G21 = ZERO
 273
 274         G22 = ZERO
 275
 276       ENDIF
 277
 278 C
 279
 280       A11 = XG * Y**2 * G11  +  (1.D0-Y) * G12
 281
 282      &      - (2.D0-Y) * SQRT( 1.D0-Y ) * G1A  *  COS( PHI )
 283
 284      &      + 2.D0 * XG * (1.D0-Y) * G1B  *  COS( 2.D0*PHI )
 285
 286 C
 287
 288       A12 = XG * Y**2 * G21  +  (1.D0-Y) * G22
 289
 290 C
 291
 292       A44 = XG * Y * (2.D0-Y) * G3
 293
 294      &      - 2.D0 * Y * SQRT( 1.D0-Y ) * GC  *  COS( PHI )
 295
 296 C
 297
 298       IF ( Y*Q2**2 .LT. 1D-38 ) THEN
 299
 300 C---prevent numerical uncertainties in DSIGMA computation
 301
 302         DSIGMA = PDENS*ALPHA**2*ALPHAS*GEV2NB*CCOL/(16.D0*PIFAC)
 303
 304      &           *(DLQ(1)*A11 + DLQ(2)*A12 + LEP*DLQ(3)*A44)
 305
 306         IF ( DSIGMA .LE. ZERO ) GOTO 888
 307
 308         LDSIG = LOG (DSIGMA) - LOG (Y) - 2.D0 * LOG (Q2)
 309
 310         DSIGMA = EXP (LDSIG)
 311
 312       ELSE
 313
 314         DSIGMA = PDENS*ALPHA**2*ALPHAS*GEV2NB*CCOL
 315
 316      &         * (DLQ(1)*A11 + DLQ(2)*A12 + LEP*DLQ(3)*A44)
 317
 318      &         / (16.D0*PIFAC*Y*Q2**2)
 319
 320       ENDIF
 321
 322       IF (DSIGMA.LT.ZERO) GOTO 888
 323
 324       RETURN
 325
 326 C
 327
 328   200 CONTINUE
 329
 330 C--- J/psi production
 331
 332       ALPHA = HWUAEM(-Q2)
 333
 334       GAMMA = 4.8D-6
 335
 336       PDENS = SFUN(13)/ETA
 337
 338       AFACT = (8.D0*PIFAC*ALPHAS**2*RMASS(164)**3*GAMMA)/(3.D0*ALPHA)
 339
 340       BFACT = ONE/(Y*SG*Z**2*((Z-ONE)*Y*SG-RMASS(164)**2)**2)
 341
 342       CFACT = (RMASS(164)**2-Z*Y*SG)**2/(Y*SG*(ONE-XG)**2*
 343
 344      &        ((ONE-XG)*Y*SG-RMASS(164)**2)**2*
 345
 346      &        ((Z-ONE)*Y*SG-RMASS(164)**2)**2)
 347
 348       DFACT = ((Z-ONE)*Y*SG)**2/(Y*SG*(ONE-XG)**2*
 349
 350      &          ((ONE-XG)*Y*SG-RMASS(164)**2)**2*(Z*Y*SG)**2)
 351
 352       DSIGMA = GEV2NB*ALPHA/(TWO*PIFAC)*AFACT*(BFACT+CFACT+DFACT)*PDENS
 353
 354       IF (DSIGMA.LT.ZERO ) GOTO 888
 355
 356       RETURN
 357
 358  888  DSIGMA=ZERO
 359
 360       END