HERWIG/src/hwhigw.f

   1
   2 CDECK  ID>, HWHIGW.
   3
   4 *CMZ :-        -26/04/91  14.55.44  by  Federico Carminati
   5
   6 *-- Author :    Mike Seymour
   7
   8 C-----------------------------------------------------------------------
   9
  10       SUBROUTINE HWHIGW
  11
  12 C-----------------------------------------------------------------------
  13
  14 C     HIGGS PRODUCTION VIA W BOSON FUSION
  15
  16 C     MEAN EVWGT = HIGGS PRODN C-S * BRANCHING FRACTION IN NB
  17
  18 C-----------------------------------------------------------------------
  19
  20       INCLUDE 'HERWIG61.INC'
  21
  22       DOUBLE PRECISION HWULDO,HWRUNI,HWR,HWUAEM,K1MAX2,K1MIN2,K12,
  23
  24      & K2MAX2,K2MIN2,K22,EMW2,EMW,ROOTS,EMH2,EMH,ROOTS2,P1,PHI1,PHI2,
  25
  26      & COSPHI,COSTH1,SINTH1,COSTH2,SINTH2,P2,WEIGHT,TAU,TAULN,CSFAC,
  27
  28      & PSUM,PROB,Q1(5),Q2(5),H(5),A,B,C,TERM2,BRHIGQ,G1WW,G2WW,G1ZZ(6),
  29
  30      & G2ZZ(6),AWW,AZZ(6),PWW,PZZ(6),EMZ,EMZ2,RSUM,GLUSQ,GRUSQ,GLDSQ,
  31
  32      & GRDSQ,GLESQ,GRESQ,CW,CZ,EMFAC,CV,CA,BR,X2,ETA,P1JAC,FACTR,EH2
  33
  34       INTEGER HWRINT,IDEC,I,ID1,ID2,IHAD
  35
  36       LOGICAL EE,EP
  37
  38       EXTERNAL HWULDO,HWRUNI,HWR,HWUAEM,HWRINT
  39
  40       SAVE EMW2,EMZ2,EE,GLUSQ,GRUSQ,GLDSQ,GRDSQ,GLESQ,GRESQ,G1ZZ,G2ZZ,
  41
  42      & G1WW,G2WW,CW,CZ,PSUM,AWW,PWW,AZZ,PZZ,ROOTS,Q1,Q2,H,FACTR
  43
  44       EQUIVALENCE (EMW,RMASS(198)),(EMZ,RMASS(200))
  45
  46       IHAD=2
  47
  48       IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD)
  49
  50       IF (FSTWGT) THEN
  51
  52         EMW2=EMW**2
  53
  54         EMZ2=EMZ**2
  55
  56         GLUSQ=(VFCH(2,1)+AFCH(2,1))**2
  57
  58         GRUSQ=(VFCH(2,1)-AFCH(2,1))**2
  59
  60         GLDSQ=(VFCH(1,1)+AFCH(1,1))**2
  61
  62         GRDSQ=(VFCH(1,1)-AFCH(1,1))**2
  63
  64         GLESQ=(VFCH(11,1)+AFCH(11,1))**2
  65
  66         GRESQ=(VFCH(11,1)-AFCH(11,1))**2
  67
  68         G1ZZ(1)=GLUSQ*GLUSQ+GRUSQ*GRUSQ
  69
  70         G2ZZ(1)=GLUSQ*GRUSQ+GRUSQ*GLUSQ
  71
  72         G1ZZ(2)=GLUSQ*GLDSQ+GRUSQ*GRDSQ
  73
  74         G2ZZ(2)=GLUSQ*GRDSQ+GRUSQ*GLDSQ
  75
  76         G1ZZ(3)=GLDSQ*GLDSQ+GRDSQ*GRDSQ
  77
  78         G2ZZ(3)=GLDSQ*GRDSQ+GRDSQ*GLDSQ
  79
  80         G1ZZ(4)=GLESQ*GLESQ+GRESQ*GRESQ
  81
  82         G2ZZ(4)=GLESQ*GRESQ+GRESQ*GLESQ
  83
  84         G1ZZ(5)=GLESQ*GLUSQ+GRESQ*GRUSQ
  85
  86         G2ZZ(5)=GLESQ*GRUSQ+GRESQ*GLUSQ
  87
  88         G1ZZ(6)=GLESQ*GLDSQ+GRESQ*GRDSQ
  89
  90         G2ZZ(6)=GLESQ*GRDSQ+GRESQ*GLDSQ
  91
  92         G1WW=0.25
  93
  94         G2WW=0
  95
  96         FACTR=GEV2NB/(128.*PIFAC**3)
  97
  98         EH2=RMASS(201)**2
  99
 100         CW=256*(PIFAC*HWUAEM(EH2)/SWEIN)**3*EMW2
 101
 102         CZ=256.*(PIFAC*HWUAEM(EH2))**3*EMZ2/(SWEIN*(1.-SWEIN))
 103
 104       ENDIF
 105
 106       EE=IPRO.LT.10
 107
 108       EP=IPRO.GE.90
 109
 110       IF (.NOT.GENEV) THEN
 111
 112 C---CHOOSE PARAMETERS
 113
 114         EVWGT=0.
 115
 116         CALL HWHIGM(EMH,EMFAC)
 117
 118         IF (EMH.LE.ZERO .OR. EMH.GE.PHEP(5,3)) RETURN
 119
 120         EMSCA=EMH
 121
 122         IF (EE) THEN
 123
 124           ROOTS=PHEP(5,3)
 125
 126         ELSE
 127
 128           TAU=(EMH/PHEP(5,3))**2
 129
 130           TAULN=LOG(TAU)
 131
 132           ROOTS=PHEP(5,3)*SQRT(EXP(HWRUNI(0,-1D-10,TAULN)))
 133
 134         ENDIF
 135
 136         EMH2=EMH**2
 137
 138         ROOTS2=ROOTS**2
 139
 140 C---CHOOSE P1 ACCORDING TO (1-ETA)*(ETA-X2)/ETA**2
 141
 142 C   WHERE ETA=1-2P1/ROOTS AND X2=EMH**2/S
 143
 144         X2=EMH2/ROOTS2
 145
 146  1      ETA=X2**HWR()
 147
 148         IF (HWR()*(1-EMH/ROOTS)**2*ETA.GT.(1-ETA)*(ETA-X2))GOTO 1
 149
 150         P1JAC=0.5*ROOTS*ETA**2/((1-ETA)*(ETA-X2))
 151
 152      &    *(-LOG(X2)*(1+X2)-2*(1-X2))
 153
 154         P1=0.5*ROOTS*(1-ETA)
 155
 156 C---CHOOSE PHI1,2 UNIFORMLY
 157
 158         PHI1=2*PIFAC*HWR()
 159
 160         PHI2=2*PIFAC*HWR()
 161
 162         COSPHI=COS(PHI2-PHI1)
 163
 164 C---CHOOSE K1^2, ON PROPAGATOR FACTOR
 165
 166         K1MAX2=2*P1*ROOTS
 167
 168         K1MIN2=0
 169
 170         K12=EMW2-(EMW2+K1MAX2)*(EMW2+K1MIN2)/
 171
 172      &           ((K1MAX2-K1MIN2)*HWR()+(EMW2+K1MIN2))
 173
 174 C---CALCULATE COSTH1 FROM K1^2
 175
 176         COSTH1=1+K12/(P1*ROOTS)
 177
 178         SINTH1=SQRT(1-COSTH1**2)
 179
 180 C---CHOOSE K2^2
 181
 182         K2MAX2=ROOTS*(ROOTS2-EMH2-2*ROOTS*P1)*(ROOTS-P1-P1*COSTH1)
 183
 184      &        /((ROOTS-P1)**2-(P1*COSTH1)**2-(P1*SINTH1*COSPHI)**2)
 185
 186         K2MIN2=0
 187
 188         K22=EMW2-(EMW2+K2MAX2)*(EMW2+K2MIN2)/
 189
 190      &           ((K2MAX2-K2MIN2)*HWR()+(EMW2+K2MIN2))
 191
 192 C---CALCULATE A,B,C FACTORS, AND...
 193
 194         A=-2*K22*P1*COSTH1 - ROOTS*(ROOTS2-EMH2-2*ROOTS*P1)
 195
 196         B=-2*K22*P1*SINTH1*COSPHI
 197
 198         C=+2*K22*P1 - 2*ROOTS*K22 - ROOTS*(ROOTS2-EMH2-2*ROOTS*P1)
 199
 200 C---SOLVE A*COSTH2 + B*SINTH2 + C = 0 FOR COSTH2
 201
 202         TERM2=B**2 + A**2 - C**2
 203
 204         IF (TERM2.LT.ZERO) RETURN
 205
 206         TERM2=B*SQRT(TERM2)
 207
 208         IF (A.GE.ZERO) RETURN
 209
 210         COSTH2=(-A*C + TERM2)/(A**2+B**2)
 211
 212         SINTH2=SQRT(1-COSTH2**2)
 213
 214 C---FINALLY, GET P2
 215
 216         IF (COSTH2.EQ.-ONE) RETURN
 217
 218         P2=-K22/(ROOTS*(1+COSTH2))
 219
 220 C---LOAD UP CMF MOMENTA
 221
 222         Q1(1)=P1*SINTH1*COS(PHI1)
 223
 224         Q1(2)=P1*SINTH1*SIN(PHI1)
 225
 226         Q1(3)=P1*COSTH1
 227
 228         Q1(4)=P1
 229
 230         Q1(5)=0
 231
 232         Q2(1)=P2*SINTH2*COS(PHI2)
 233
 234         Q2(2)=P2*SINTH2*SIN(PHI2)
 235
 236         Q2(3)=P2*COSTH2
 237
 238         Q2(4)=P2
 239
 240         Q2(5)=0
 241
 242         H(1)=-Q1(1)-Q2(1)
 243
 244         H(2)=-Q1(2)-Q2(2)
 245
 246         H(3)=-Q1(3)-Q2(3)
 247
 248         H(4)=-Q1(4)-Q2(4)+ROOTS
 249
 250         CALL HWUMAS(H)
 251
 252 C---CALCULATE MATRIX ELEMENTS SQUARED
 253
 254         AWW=ENHANC(10)**2 * CW*(ROOTS2/2*HWULDO(Q1,Q2)*G1WW
 255
 256      &         +ROOTS2/4*P1*P2*(1+COSTH1)*(1-COSTH2)*G2WW)
 257
 258         DO 10 I=1,6
 259
 260           AZZ(I)=ENHANC(11)**2 * CZ*(ROOTS2/2*HWULDO(Q1,Q2)*G1ZZ(I)
 261
 262      &               +ROOTS2/4*P1*P2*(1+COSTH1)*(1-COSTH2)*G2ZZ(I))
 263
 264      &          *((K12-EMW2)/(K12-EMZ2)*(K22-EMW2)/(K22-EMZ2))**2
 265
 266  10     CONTINUE
 267
 268 C---CALCULATE WEIGHT IN INTEGRAL
 269
 270         WEIGHT=FACTR*P2*P1JAC/(ROOTS2**2*HWULDO(H,Q2))
 271
 272      &              *(K1MAX2-K1MIN2)/((K1MAX2+EMW2)*(K1MIN2+EMW2))
 273
 274      &              *(K2MAX2-K2MIN2)/((K2MAX2+EMW2)*(K2MIN2+EMW2))
 275
 276      &              * EMFAC
 277
 278         EMSCA=EMW
 279
 280         XXMIN=(ROOTS/PHEP(5,3))**2
 281
 282         XLMIN=LOG(XXMIN)
 283
 284 C---INCLUDE BRANCHING RATIO OF HIGGS
 285
 286         IDEC=MOD(IPROC,100)
 287
 288         IF (IDEC.GT.0.AND.IDEC.LE.12) WEIGHT=WEIGHT*BRHIG(IDEC)
 289
 290         IF (IDEC.EQ.0) THEN
 291
 292           BRHIGQ=0
 293
 294           DO 20 I=1,6
 295
 296  20         BRHIGQ=BRHIGQ+BRHIG(I)
 297
 298           WEIGHT=WEIGHT*BRHIGQ
 299
 300         ENDIF
 301
 302         IF (IDEC.EQ.10) THEN
 303
 304           CALL HWDBOZ(198,ID1,ID2,CV,CA,BR,1)
 305
 306           CALL HWDBOZ(199,ID1,ID2,CV,CA,BR,1)
 307
 308           WEIGHT=WEIGHT*BR
 309
 310         ELSEIF (IDEC.EQ.11) THEN
 311
 312           CALL HWDBOZ(200,ID1,ID2,CV,CA,BR,1)
 313
 314           CALL HWDBOZ(200,ID1,ID2,CV,CA,BR,1)
 315
 316           WEIGHT=WEIGHT*BR
 317
 318         ENDIF
 319
 320         IF (EE) THEN
 321
 322           CSFAC=WEIGHT
 323
 324           PSUM=AWW+AZZ(4)
 325
 326           EVWGT=CSFAC*PSUM
 327
 328         ELSEIF (EP) THEN
 329
 330           CSFAC=-WEIGHT*TAULN
 331
 332           XX(1)=ONE
 333
 334           XX(2)=XXMIN
 335
 336           CALL HWSFUN(XX(2),EMSCA,IDHW(IHAD),NSTRU,DISF(1,2),2)
 337
 338           IF (IDHW(1).LE.126) THEN
 339
 340             PWW=(DISF(2,2)+DISF(4,2)+DISF(7,2)+DISF( 9,2))*AWW
 341
 342           ELSE
 343
 344             PWW=(DISF(1,2)+DISF(3,2)+DISF(8,2)+DISF(10,2))*AWW
 345
 346           ENDIF
 347
 348           PZZ(5)=(DISF(2,2)+DISF(4,2)+DISF(8,2)+DISF(10,2))*AZZ(5)
 349
 350           PZZ(6)=(DISF(1,2)+DISF(3,2)+DISF(7,2)+DISF( 9,2))*AZZ(6)
 351
 352           PSUM=PWW+PZZ(5)+PZZ(6)
 353
 354           EVWGT=CSFAC*PSUM
 355
 356         ELSE
 357
 358           CSFAC=WEIGHT*TAULN*XLMIN
 359
 360           CALL HWSGEN(.TRUE.)
 361
 362           PWW=((DISF(2,1)+DISF(4, 1)+DISF(7,1)+DISF(9,1))
 363
 364      &        *(DISF(8,2)+DISF(10,2)+DISF(1,2)+DISF(3,2))
 365
 366      &        +(DISF(8,1)+DISF(10,1)+DISF(1,1)+DISF(3,1))
 367
 368      &        *(DISF(2,2)+DISF(4, 2)+DISF(7,2)+DISF(9,2)))
 369
 370      &        *AWW
 371
 372           PZZ(1)=((DISF(2,1)+DISF(4,1)+DISF(8,1)+DISF(10,1))
 373
 374      &           *(DISF(2,2)+DISF(4,2)+DISF(8,2)+DISF(10,2)))
 375
 376      &           *AZZ(1)
 377
 378           PZZ(2)=((DISF(2,1)+DISF(4,1)+DISF(8,1)+DISF(10,1))
 379
 380      &           *(DISF(1,2)+DISF(3,2)+DISF(7,2)+DISF(9, 2))
 381
 382      &           +(DISF(1,1)+DISF(3,1)+DISF(7,1)+DISF(9, 1))
 383
 384      &           *(DISF(2,2)+DISF(4,2)+DISF(8,2)+DISF(10,2)))
 385
 386      &           *AZZ(2)
 387
 388           PZZ(3)=((DISF(1,1)+DISF(3,1)+DISF(7,1)+DISF(9,1))
 389
 390      &           *(DISF(1,2)+DISF(3,2)+DISF(7,2)+DISF(9,2)))
 391
 392      &           *AZZ(3)
 393
 394           PSUM=PWW+PZZ(1)+PZZ(2)+PZZ(3)
 395
 396 C---EVENT WEIGHT IS SUM OVER ALL COMBINATIONS
 397
 398           EVWGT=CSFAC*PSUM
 399
 400         ENDIF
 401
 402       ELSE
 403
 404 C---GENERATE EVENT
 405
 406 C---CHOOSE EVENT TYPE
 407
 408         RSUM=PSUM*HWR()
 409
 410 C---ELECTRON BEAMS?
 411
 412         IF (EE) THEN
 413
 414           IDN(1)=IDHW(1)
 415
 416           IDN(2)=IDHW(2)
 417
 418 C---WW FUSION?
 419
 420           IF (RSUM.LT.AWW) THEN
 421
 422             IDN(3)=IDN(1)+1
 423
 424             IDN(4)=IDN(2)+1
 425
 426 C---ZZ FUSION?
 427
 428           ELSE
 429
 430             IDN(3)=IDN(1)
 431
 432             IDN(4)=IDN(2)
 433
 434           ENDIF
 435
 436 C---LEPTON-HADRON COLISION?
 437
 438         ELSEIF (EP) THEN
 439
 440 C---WW FUSION?
 441
 442           IDN(1)=IDHW(1)
 443
 444           IF (RSUM.LT.PWW) THEN
 445
 446  24         IDN(2)=HWRINT(1,8)
 447
 448             IF (IDN(2).GE.5) IDN(2)=IDN(2)+2
 449
 450             IF (ICHRG(IDN(1))*ICHRG(IDN(2)).GT.0) GOTO 24
 451
 452             PROB=DISF(IDN(2),2)*AWW/PWW
 453
 454             IF (HWR().GT.PROB) GOTO 24
 455
 456             IDN(3)=IDN(1)+1
 457
 458             IF (HWR().GT.SCABI) THEN
 459
 460               IDN(4)= 4*INT((IDN(2)-1)/2)-IDN(2)+3
 461
 462             ELSE
 463
 464               IDN(4)=12*INT((IDN(2)-1)/6)-IDN(2)+5
 465
 466             ENDIF
 467
 468 C---ZZ FUSION FROM U-TYPE QUARK?
 469
 470           ELSEIF (RSUM.LT.PWW+PZZ(5)) THEN
 471
 472  26         IDN(2)=2*HWRINT(1,4)
 473
 474             IF (IDN(2).GE.5) IDN(2)=IDN(2)+2
 475
 476             PROB=DISF(IDN(2),2)*AZZ(5)/PZZ(5)
 477
 478             IF (HWR().GT.PROB) GOTO 26
 479
 480             IDN(3)=IDN(1)
 481
 482             IDN(4)=IDN(2)
 483
 484 C---ZZ FUSION FROM D-TYPE QUARK?
 485
 486           ELSE
 487
 488  28         IDN(2)=2*HWRINT(1,4)-1
 489
 490             IF (IDN(2).GE.5) IDN(2)=IDN(2)+2
 491
 492             PROB=DISF(IDN(2),2)*AZZ(6)/PZZ(6)
 493
 494             IF (HWR().GT.PROB) GOTO 28
 495
 496             IDN(3)=IDN(1)
 497
 498             IDN(4)=IDN(2)
 499
 500           ENDIF
 501
 502 C---HADRON BEAMS?
 503
 504         ELSE
 505
 506 C---WW FUSION?
 507
 508           IF (RSUM.LT.PWW) THEN
 509
 510  31         DO 32 I=1,2
 511
 512               IDN(I)=HWRINT(1,8)
 513
 514               IF (IDN(I).GE.5) IDN(I)=IDN(I)+2
 515
 516  32         CONTINUE
 517
 518             IF (ICHRG(IDN(1))*ICHRG(IDN(2)).GT.0) GOTO 31
 519
 520             PROB=DISF(IDN(1),1)*DISF(IDN(2),2)*AWW/PWW
 521
 522             IF (HWR().GT.PROB) GOTO 31
 523
 524 C---CHOOSE OUTGOING QUARKS
 525
 526             DO 33 I=1,2
 527
 528               IF (HWR().GT.SCABI) THEN
 529
 530                 IDN(I+2)=4*INT((IDN(I)-1)/2)-IDN(I)+3
 531
 532               ELSE
 533
 534                 IDN(I+2)=12*INT((IDN(I)-1)/6)-IDN(I)+5
 535
 536               ENDIF
 537
 538  33         CONTINUE
 539
 540 C---ZZ FUSION FROM U-TYPE QUARKS?
 541
 542           ELSEIF (RSUM.LT.PWW+PZZ(1)) THEN
 543
 544  41         DO 42 I=1,2
 545
 546               IDN(I)=2*HWRINT(1,4)
 547
 548               IF (IDN(I).GE.5) IDN(I)=IDN(I)+2
 549
 550  42         CONTINUE
 551
 552             PROB=DISF(IDN(1),1)*DISF(IDN(2),2)*AZZ(1)/PZZ(1)
 553
 554             IF (HWR().GT.PROB) GOTO 41
 555
 556             IDN(3)=IDN(1)
 557
 558             IDN(4)=IDN(2)
 559
 560 C---ZZ FUSION FROM D-TYPE QUARKS?
 561
 562           ELSEIF (RSUM.LT.PWW+PZZ(1)+PZZ(3)) THEN
 563
 564  51         DO 52 I=1,2
 565
 566               IDN(I)=2*HWRINT(1,4)-1
 567
 568               IF (IDN(I).GE.5) IDN(I)=IDN(I)+2
 569
 570  52         CONTINUE
 571
 572             PROB=DISF(IDN(1),1)*DISF(IDN(2),2)*AZZ(3)/PZZ(3)
 573
 574             IF (HWR().GT.PROB) GOTO 51
 575
 576             IDN(3)=IDN(1)
 577
 578             IDN(4)=IDN(2)
 579
 580 C---ZZ FUSION FROM UD-TYPE PAIRS?
 581
 582           ELSE
 583
 584  61         IF (HWR().GT.HALF) THEN
 585
 586               IDN(1)=2*HWRINT(1,4)-1
 587
 588               IDN(2)=2*HWRINT(1,4)
 589
 590             ELSE
 591
 592               IDN(1)=2*HWRINT(1,4)
 593
 594               IDN(2)=2*HWRINT(1,4)-1
 595
 596             ENDIF
 597
 598             DO 62 I=1,2
 599
 600  62           IF (IDN(I).GE.5) IDN(I)=IDN(I)+2
 601
 602             PROB=DISF(IDN(1),1)*DISF(IDN(2),2)*AZZ(2)/PZZ(2)
 603
 604             IF (HWR().GT.PROB) GOTO 61
 605
 606             IDN(3)=IDN(1)
 607
 608             IDN(4)=IDN(2)
 609
 610           ENDIF
 611
 612         ENDIF
 613
 614 C---NOW BOOST TO LAB, AND SET UP STATUS CODES etc
 615
 616         IDCMF=15
 617
 618 C---INCOMING
 619
 620         IF (.NOT.EE) CALL HWEONE
 621
 622 C---CMF POINTERS
 623
 624         JDAHEP(1,NHEP)=NHEP+1
 625
 626         JDAHEP(2,NHEP)=NHEP+3
 627
 628         JMOHEP(1,NHEP+1)=NHEP
 629
 630         JMOHEP(1,NHEP+2)=NHEP
 631
 632         JMOHEP(1,NHEP+3)=NHEP
 633
 634 C---OUTGOING MOMENTA (GIVE QUARKS MASS NON-COVARIANTLY!)
 635
 636         Q1(5)=RMASS(IDN(1))
 637
 638         Q1(4)=SQRT(Q1(4)**2+Q1(5)**2)
 639
 640         Q2(5)=RMASS(IDN(2))
 641
 642         Q2(4)=SQRT(Q2(4)**2+Q2(5)**2)
 643
 644         H(4)=-Q1(4)-Q2(4)+PHEP(5,NHEP)
 645
 646         CALL HWUMAS(H)
 647
 648         CALL HWULOB(PHEP(1,NHEP),Q1,PHEP(1,NHEP+1))
 649
 650         CALL HWULOB(PHEP(1,NHEP),Q2,PHEP(1,NHEP+2))
 651
 652         CALL HWULOB(PHEP(1,NHEP),H,PHEP(1,NHEP+3))
 653
 654 C---STATUS AND IDs
 655
 656         ISTHEP(NHEP+1)=113
 657
 658         ISTHEP(NHEP+2)=114
 659
 660         ISTHEP(NHEP+3)=114
 661
 662         IDHW(NHEP+1)=IDN(3)
 663
 664         IDHEP(NHEP+1)=IDPDG(IDN(3))
 665
 666         IDHW(NHEP+2)=IDN(4)
 667
 668         IDHEP(NHEP+2)=IDPDG(IDN(4))
 669
 670         IDHW(NHEP+3)=201
 671
 672         IDHEP(NHEP+3)=IDPDG(201)
 673
 674 C---COLOUR LABELS
 675
 676         JMOHEP(2,NHEP+1)=NHEP-2
 677
 678         JMOHEP(2,NHEP+2)=NHEP-1
 679
 680         JMOHEP(2,NHEP-1)=NHEP+2
 681
 682         JMOHEP(2,NHEP-2)=NHEP+1
 683
 684         JMOHEP(2,NHEP+3)=NHEP+3
 685
 686         JDAHEP(2,NHEP+1)=NHEP-2
 687
 688         JDAHEP(2,NHEP+2)=NHEP-1
 689
 690         JDAHEP(2,NHEP-1)=NHEP+2
 691
 692         JDAHEP(2,NHEP-2)=NHEP+1
 693
 694         JDAHEP(2,NHEP+3)=NHEP+3
 695
 696         NHEP=NHEP+3
 697
 698       ENDIF
 699
 700   999 END
 701
 702 CDECK  ID>, HWHIGY.
 703
 704 *CMZ :-        -26/04/91  13.37.37  by  Federico Carminati
 705
 706 *-- Author :    Mike Seymour
 707
 708 C-----------------------------------------------------------------------
 709
 710       FUNCTION HWHIGY(A,B,XP)
 711
 712 C-----------------------------------------------------------------------
 713
 714 C     CALCULATE THE INTEGRAL OF BERENDS AND KLEISS APPENDIX B
 715
 716 C-----------------------------------------------------------------------
 717
 718       IMPLICIT NONE
 719
 720       COMPLEX XQ,Z1,Z2,Z3,Z4,C0,C1,C2,C3,C4,C5,C6,C7,C8,FUN,Z
 721
 722       DOUBLE PRECISION HWHIGY,TWO,A,B,XP,Y
 723
 724       PARAMETER (TWO=2.D0)
 725
 726 C---DECLARE ALL THE STATEMENT-FUNCTION DEFINITIONS
 727
 728       C0(Z,A)=(Z**2-A)**2*((Z**2+A)**2-24*Z*(Z**2+A)+8*Z**2*(A+6))/Z**4
 729
 730       C1(Z,A)=A**4/(3*Z)
 731
 732       C2(Z,A)=-A**3*(24*Z-A)/(2*Z**2)
 733
 734       C3(Z,A)=A**2*(8*Z**2*(A+6)-24*A*Z+A**2)/Z**3
 735
 736       C4(Z,A)=-A**2*(24*Z**3+8*Z**2*(A+6)-24*A*Z+A**2)/Z**4
 737
 738       C5(Z,A)=Z**3-24*Z**2+8*Z*(A+6)+24*A
 739
 740       C6(Z,A)=0.5*Z**2-12*Z+4*(A+6)
 741
 742       C7(Z,A)=Z/3-8
 743
 744       C8(Z,A)=0.25
 745
 746       FUN(Z,Y,A)=C0(Z,A)*LOG(Y-Z)
 747
 748      &          +C1(Z,A)/Y**3
 749
 750      &          +C2(Z,A)/Y**2
 751
 752      &          +C3(Z,A)/Y
 753
 754      &          +C4(Z,A)*LOG(Y)
 755
 756      &          +C5(Z,A)*Y
 757
 758      &          +C6(Z,A)*Y**2
 759
 760      &          +C7(Z,A)*Y**3
 761
 762      &          +C8(Z,A)*Y**4
 763
 764 C---NOW EVALUATE THE INTEGRAL
 765
 766       HWHIGY=0
 767
 768       IF (A.GT.4) RETURN
 769
 770       XQ=CMPLX(XP,B)
 771
 772       Z1=XQ+SQRT(XQ**2-A)
 773
 774       Z2=XQ-SQRT(XQ**2-A)
 775
 776       Z3=FUN(Z1,TWO,A)-FUN(Z1,SQRT(A),A)
 777
 778       Z4=FUN(Z2,TWO,A)-FUN(Z2,SQRT(A),A)
 779
 780       HWHIGY=AIMAG((Z3-Z4)/(Z1-Z2))/(8*B)
 781
 782       END