HERWIG/src/hwegam.f

   1
   2 CDECK  ID>, HWEGAM.
   3
   4 *CMZ :-        -26/04/91  11.11.55  by  Bryan Webber
   5
   6 *-- Author :    Bryan Webber & Luca Stanco
   7
   8 C-----------------------------------------------------------------------
   9
  10       SUBROUTINE HWEGAM(IHEP,ZMI,ZMA,WWA)
  11
  12 C-----------------------------------------------------------------------
  13
  14 C     GENERATES A PHOTON IN WEIZSACKER-WILLIAMS (WWA=.TRUE.) OR
  15
  16 C     ELSE EQUIVALENT PHOTON APPROX FROM INCOMING E+, E-, MU+ OR MU-
  17
  18 C-----------------------------------------------------------------------
  19
  20       INCLUDE 'HERWIG61.INC'
  21
  22       DOUBLE PRECISION HWR,HWRUNI,EGMIN,ZMIN,ZMAX,ZGAM,SS,ZMI,ZMA,
  23
  24      & PPL,PMI,QT2,Q2,QQMIN,QQMAX,S0,A,RPM(2)
  25
  26       INTEGER IHEP,IHADIS,HQ,I
  27
  28       LOGICAL WWA
  29
  30       EXTERNAL HWR,HWRUNI
  31
  32       DATA EGMIN/5.D0/
  33
  34       IF (IERROR.NE.0)  RETURN
  35
  36       IF (IHEP.LT.1.OR.IHEP.GT.2) CALL HWWARN('HWEGAM',500,*999)
  37
  38       SS=PHEP(5,3)
  39
  40       IF (IHEP.EQ.1) THEN
  41
  42         IHADIS=2
  43
  44       ELSE
  45
  46         IHADIS=1
  47
  48         IF (JDAHEP(1,IHADIS).NE.0) IHADIS=JDAHEP(1,IHADIS)
  49
  50       ENDIF
  51
  52 C---DEFINE LIMITS FOR GAMMA MOMENTUM FRACTION
  53
  54       IF (ZMI.LE.ZERO .OR. ZMA.GT.ONE) THEN
  55
  56         IF (IPRO.EQ.13.OR.IPRO.EQ.14) THEN
  57
  58           S0 = EMMIN**2
  59
  60         ELSEIF(IPRO.EQ.15.OR.IPRO.EQ.18.OR.IPRO.EQ.22.OR.IPRO.EQ.24.OR.
  61
  62      &         IPRO.EQ.50.OR.IPRO.EQ.53.OR.IPRO.EQ.55)THEN
  63
  64           S0 = 4.D0*PTMIN**2
  65
  66         ELSEIF (IPRO.EQ.17.OR.IPRO.EQ.51) THEN
  67
  68           HQ = MOD(IPROC,100)
  69
  70           S0 = 4.D0*(PTMIN**2+RMASS(HQ)**2)
  71
  72         ELSEIF (IPRO.EQ.16.OR.IPRO.EQ.19.OR.IPRO.EQ.95) THEN
  73
  74           S0 = MAX(2*RMASS(1),RMASS(201)-GAMMAX*GAMH)**2
  75
  76         ELSEIF (IPRO.EQ.23) THEN
  77
  78           S0 = MAX(2*RMASS(1),RMASS(201)-GAMMAX*GAMH)**2
  79
  80           S0 = (PTMIN+SQRT(PTMIN**2+S0))**2
  81
  82         ELSEIF (IPRO.EQ.20) THEN
  83
  84           S0 = RMASS(201)**2
  85
  86         ELSEIF (IPRO.EQ.21) THEN
  87
  88           S0 = (PTMIN+SQRT(PTMIN**2+RMASS(198)**2))**2
  89
  90 C--PR MOD 7/7/99
  91
  92         ELSEIF(IPRO.EQ.30) THEN
  93
  94           S0 = 4.0D0*(PTMIN**2+RMMNSS**2)
  95
  96         ELSEIF(IPRO.EQ.40.OR.IPRO.EQ.41) THEN
  97
  98           HQ = IPROC-100*IPRO
  99
 100           RPM(1) = RMMNSS
 101
 102           RPM(2) = ZERO
 103
 104           IF(HQ.GE.10.AND.HQ.LT.20) THEN
 105
 106             RPM(1) = ABS(RMASS(450))
 107
 108             IF(HQ.GT.10) RPM(1) = ABS(RMASS(449+MOD(HQ,10)))
 109
 110           ELSEIF(HQ.GE.20.AND.HQ.LT.30) THEN
 111
 112             RPM(1) = ABS(RMASS(454))
 113
 114             IF(HQ.GT.20) RPM(1) = ABS(RMASS(453+MOD(HQ,20)))
 115
 116           ELSEIF(HQ.EQ.30) THEN
 117
 118             RPM(1) = RMASS(449)
 119
 120           ELSEIF(HQ.EQ.40) THEN
 121
 122             IF(IPRO.EQ.40) THEN
 123
 124               RPM(1) = RMASS(425)
 125
 126               DO I=1,5
 127
 128                 RPM(1) = MIN(RPM(1),RMASS(425+I))
 129
 130               ENDDO
 131
 132             ELSE
 133
 134               RPM(1) = MIN(RMASS(405),RMASS(406))
 135
 136             ENDIF
 137
 138             RPM(2) = RMASS(198)
 139
 140           ELSEIF(HQ.EQ.50) THEN
 141
 142             IF(IPRO.EQ.40) THEN
 143
 144               RPM(1) = RMASS(425)
 145
 146               DO I=1,5
 147
 148                 RPM(1) = MIN(RPM(1),RMASS(425+I))
 149
 150               ENDDO
 151
 152               DO I=1,3
 153
 154                 RPM(2) = MIN(RPM(1),RMASS(433+2*I))
 155
 156               ENDDO
 157
 158               RPM(1) = MIN(RPM(1),RPM(2))
 159
 160               RPM(2) = RMASS(203)
 161
 162               DO I=1,2
 163
 164                 RPM(2) = MIN(RPM(2),RMASS(204+I))
 165
 166               ENDDO
 167
 168             ELSE
 169
 170               RPM(1) = RMASS(401)
 171
 172               RPM(2) = RMASS(413)
 173
 174               DO I=1,5
 175
 176                 RPM(1) = MIN(RPM(1),RMASS(401+I))
 177
 178                 RPM(2) = MIN(RPM(2),RMASS(413+I))
 179
 180               ENDDO
 181
 182               RPM(1) = MIN(RPM(1),RPM(2))
 183
 184               RPM(2) = RMASS(203)
 185
 186               DO I=1,2
 187
 188                 RPM(2) = MIN(RPM(2),RMASS(204+I))
 189
 190               ENDDO
 191
 192             ENDIF
 193
 194             RPM(2) = RMASS(203)
 195
 196             DO I=1,2
 197
 198               RPM(2) = MIN(RPM(2),RMASS(204+I))
 199
 200             ENDDO
 201
 202           ELSEIF(HQ.GE.60) THEN
 203
 204             RPM(1) = ZERO
 205
 206           ENDIF
 207
 208           RPM(1) = RPM(1)**2
 209
 210           RPM(2) = RPM(2)**2
 211
 212           S0 = RPM(1)+RPM(2)+TWO*(PTMIN**2+
 213
 214      &         SQRT(RPM(1)*RPM(2)+PTMIN**2*(RPM(1)+RPM(2)+PTMIN**2)))
 215
 216 C--end of mod
 217
 218         ELSEIF (IPRO.EQ.52) THEN
 219
 220           HQ = MOD(IPROC,100)
 221
 222           S0 = (PTMIN+SQRT(PTMIN**2+RMASS(HQ)**2))**2
 223
 224         ELSEIF (IPRO.EQ.60) THEN
 225
 226           HQ = MOD(IPROC,100)
 227
 228           IF (HQ.EQ.0) THEN
 229
 230             S0 = 4.D0*PTMIN**2
 231
 232           ELSE
 233
 234             IF (HQ.GT.6) HQ=2*HQ+107
 235
 236             IF (HQ.EQ.127) HQ=198
 237
 238             S0 = 4.D0*(PTMIN**2+RMASS(HQ)**2)
 239
 240           ENDIF
 241
 242         ELSEIF (IPRO.EQ.80) THEN
 243
 244           S0 = WHMIN**2
 245
 246         ELSEIF (IPRO.EQ.90) THEN
 247
 248           S0 = Q2MIN
 249
 250         ELSEIF (IPRO.EQ.91.OR.IPRO.EQ.92) THEN
 251
 252           S0 = Q2MIN+4.D0*PTMIN**2
 253
 254           HQ = MOD(IPROC,100)
 255
 256           IF (HQ.GT.0) S0 = S0+4.D0*RMASS(HQ)**2
 257
 258           IF (IPRO.EQ.91) S0 = MAX(S0,EMMIN**2)
 259
 260         ELSE
 261
 262           S0 = 0
 263
 264         ENDIF
 265
 266         IF (S0.GT.ZERO) THEN
 267
 268           S0 = (SQRT(S0)+ABS(PHEP(5,IHADIS)))**2-PHEP(5,IHADIS)**2
 269
 270           S0 = MAX(S0,WHMIN**2)
 271
 272           ZMIN = S0 / (SS**2 - PHEP(5,IHEP)**2 - PHEP(5,IHADIS)**2)
 273
 274           ZMAX = ONE
 275
 276         ELSE
 277
 278 C---UNKNOWN PROCESS: USE ENERGY CUTOFF, AND WARN USER
 279
 280           IF (FSTWGT) CALL HWWARN('HWEGAM',1,*999)
 281
 282           ZMIN = EGMIN / PHEP(4,IHEP)
 283
 284           ZMAX = ONE
 285
 286         ENDIF
 287
 288       ELSE
 289
 290         ZMIN=ZMI
 291
 292         ZMAX=ZMA
 293
 294       ENDIF
 295
 296 C---APPLY USER DEFINED CUTS YWWMIN,YWWMAX AND INDIRECT LIMITS ON Z
 297
 298       IF (.NOT.WWA) THEN
 299
 300         ZMIN=MAX(ZMIN,YWWMIN,SQRT(Q2WWMN)/ABS(PHEP(3,IHEP)))
 301
 302         ZMAX=MIN(ZMAX,YWWMAX)
 303
 304         IF (ZMIN.GT.ZMAX) THEN
 305
 306           GAMWT=ZERO
 307
 308           RETURN
 309
 310         ENDIF
 311
 312       ENDIF
 313
 314 C---GENERATE GAMMA MOMENTUM FRACTION
 315
 316       A=HALF
 317
 318  10   IF (HWR().LT.A) THEN
 319
 320         ZGAM=(ZMIN/ZMAX)**HWR()*ZMAX
 321
 322       ELSE
 323
 324         ZGAM=(ZMAX-ZMIN)*HWR()+ZMIN
 325
 326       ENDIF
 327
 328       GAMWT = GAMWT * .5*ALPHEM/PIFAC *
 329
 330      +     (1+(1-ZGAM)**2)/(A/LOG(ZMAX/ZMIN)+(1-A)/(ZMAX-ZMIN)*ZGAM)
 331
 332       IF (WWA) THEN
 333
 334         GAMWT = GAMWT * LOG((ONE-ZGAM)/ZGAM*(SS/PHEP(5,IHEP))**2)
 335
 336       ELSE
 337
 338 C---Q2WWMN AND Q2WWMX ARE USER-DEFINED LIMITS IN THE Q**2 INTEGRATION
 339
 340         QQMAX=MIN(Q2WWMX,(ZGAM*PHEP(3,IHEP))**2)
 341
 342         QQMIN=MAX(Q2WWMN,(PHEP(5,IHEP)*ZGAM)**2/(1.-ZGAM))
 343
 344         IF (QQMIN.GT.QQMAX) CALL HWWARN('HWEGAM',50,*10)
 345
 346         Q2=EXP(HWRUNI(0,LOG(QQMIN),LOG(QQMAX)))
 347
 348         GAMWT = GAMWT * LOG(QQMAX/QQMIN)
 349
 350       ENDIF
 351
 352       IF (GAMWT.LT.ZERO) GAMWT=ZERO
 353
 354 C---FILL PHOTON
 355
 356       NHEP=NHEP+1
 357
 358       IDHW(NHEP)=59
 359
 360       ISTHEP(NHEP)=3
 361
 362       IDHEP(NHEP)=22
 363
 364       JMOHEP(1,NHEP)=IHEP
 365
 366       JMOHEP(2,NHEP)=0
 367
 368       JDAHEP(1,NHEP)=0
 369
 370       JDAHEP(2,NHEP)=0
 371
 372       JDAHEP(1,IHEP)=NHEP
 373
 374       IF (WWA) THEN
 375
 376 C---FOR COLLINEAR KINEMATICS, ZGAM IS THE ENERGY FRACTION
 377
 378         PHEP(4,NHEP)=PHEP(4,IHEP)*ZGAM
 379
 380         PHEP(3,NHEP)=PHEP(3,IHEP)-SIGN(SQRT(
 381
 382      &     (PHEP(4,IHEP)-PHEP(4,NHEP))**2-PHEP(5,IHEP)**2),PHEP(3,IHEP))
 383
 384         PHEP(2,NHEP)=0
 385
 386         PHEP(1,NHEP)=0
 387
 388         CALL HWUMAS(PHEP(1,NHEP))
 389
 390       ELSE
 391
 392 C---FOR EXACT KINEMATICS, ZGAM IS TAKEN TO BE FRACTION OF (E+PZ)
 393
 394         PPL=ZGAM*(ABS(PHEP(3,IHEP))+PHEP(4,IHEP))
 395
 396         QT2=(ONE-ZGAM)*Q2-(ZGAM*PHEP(5,IHEP))**2
 397
 398         PMI=(QT2-Q2)/PPL
 399
 400         PHEP(5,NHEP)=-SQRT(Q2)
 401
 402         PHEP(4,NHEP)=(PPL+PMI)/TWO
 403
 404         PHEP(3,NHEP)=SIGN((PPL-PMI)/TWO,PHEP(3,IHEP))
 405
 406         CALL HWRAZM(SQRT(QT2),PHEP(1,NHEP),PHEP(2,NHEP))
 407
 408       ENDIF
 409
 410 C---UPDATE OVERALL CM FRAME
 411
 412       JMOHEP(IHEP,3)=NHEP
 413
 414       CALL HWVDIF(4,PHEP(1,3),PHEP(1,IHEP),PHEP(1,3))
 415
 416       CALL HWVSUM(4,PHEP(1,NHEP),PHEP(1,3),PHEP(1,3))
 417
 418       CALL HWUMAS(PHEP(1,3))
 419
 420 C---FILL OUTGOING LEPTON
 421
 422       NHEP=NHEP+1
 423
 424       IDHW(NHEP)=IDHW(IHEP)
 425
 426       ISTHEP(NHEP)=1
 427
 428       IDHEP(NHEP)=IDHEP(IHEP)
 429
 430       JMOHEP(1,NHEP)=IHEP
 431
 432       JMOHEP(2,NHEP)=0
 433
 434       JDAHEP(1,NHEP)=0
 435
 436       JDAHEP(2,NHEP)=0
 437
 438       JDAHEP(2,IHEP)=NHEP
 439
 440       CALL HWVDIF(4,PHEP(1,IHEP),PHEP(1,NHEP-1),PHEP(1,NHEP))
 441
 442       PHEP(5,NHEP)=PHEP(5,IHEP)
 443
 444  999  END
 445
 446 CDECK  ID>, HWEINI.
 447
 448 *CMZ :-        -26/04/91  12.42.30  by  Federico Carminati
 449
 450 *-- Author :    Bryan Webber
 451
 452 C-----------------------------------------------------------------------
 453
 454       SUBROUTINE HWEINI
 455
 456 C-----------------------------------------------------------------------
 457
 458 C     INITIALISES ELEMENTARY PROCESS
 459
 460 C-----------------------------------------------------------------------
 461
 462       INCLUDE 'HERWIG61.INC'
 463
 464       DOUBLE PRECISION HWRSET,DUMMY,SAFETY
 465
 466       EXTERNAL HWRSET
 467
 468       PARAMETER (SAFETY=1.001)
 469
 470       INTEGER NBSH,I
 471
 472 C---NO OF WEIGHT GENERATED
 473
 474       NWGTS=0
 475
 476 C---ACCUMULATED WEIGHTS
 477
 478       WGTSUM=0.
 479
 480 C---ACCUMULATED WEIGHT-SQUARED
 481
 482       WSQSUM=0.
 483
 484 C---CURRENT MAX WEIGHT
 485
 486       WBIGST=0.
 487
 488 C---LAST VALUE OF SCALE
 489
 490       EMLST=0.
 491
 492 C---NUMBER OF ERRORS REPORTED
 493
 494       NUMER=0
 495
 496 C---NUMBER OF ERRORS UNREPORTED
 497
 498       NUMERU=0
 499
 500 C---FIND MAXIMUM EVENT WEIGHT IN CASES WHERE THIS IS REQUIRED
 501
 502       IF (NOWGT) THEN
 503
 504         IF (WGTMAX.EQ.ZERO) THEN
 505
 506           NBSH=IBSH
 507
 508           DUMMY = HWRSET(IBRN)
 509
 510           WRITE(6,10) IPROC,IBRN,NBSH
 511
 512    10     FORMAT(/10X,'INITIAL SEARCH FOR MAX WEIGHT'//
 513
 514      &            10X,'PROCESS CODE IPROC = ',I11/
 515
 516      &            10X,'RANDOM NO. SEED 1  = ',I11/
 517
 518      &            10X,'           SEED 2  = ',I11/
 519
 520      &            10X,'NUMBER OF SHOTS    = ',I11)
 521
 522           NEVHEP=0
 523
 524           DO 11 I=1,NBSH
 525
 526           CALL HWEPRO
 527
 528    11     CONTINUE
 529
 530           WRITE(6,20)
 531
 532    20     FORMAT(/10X,'INITIAL SEARCH FINISHED')
 533
 534           IF (WBIGST*NWGTS.LT.SAFETY*WGTSUM)
 535
 536      &                 WGTMAX=SAFETY*WBIGST
 537
 538           CALL HWEFIN
 539
 540           NWGTS=0
 541
 542           WGTSUM=0.
 543
 544           WSQSUM=0.
 545
 546           WBIGST=0.
 547
 548         ELSE
 549
 550           WRITE(6,21) AVWGT,WGTMAX
 551
 552    21     FORMAT(/1P,10X,'INPUT EVT WEIGHT   =',E12.4/
 553
 554      &               10X,'INPUT MAX WEIGHT   =',E12.4)
 555
 556         ENDIF
 557
 558       ENDIF
 559
 560 C---RESET RANDOM NUMBER
 561
 562       DUMMY = HWRSET(NRN)
 563
 564       ISTAT=5
 565
 566   999 END
 567
 568 CDECK  ID>, HWEISR.
 569
 570 *CMZ :-        -01/04/99  19.55.17  by  Mike Seymour
 571
 572 *-- Author :    Mike Seymour
 573
 574 C-----------------------------------------------------------------------
 575
 576       SUBROUTINE HWEISR(IHEP)
 577
 578 C-----------------------------------------------------------------------
 579
 580 C     GENERATES AN ISR PHOTON FROM INCOMING E+, E-, MU+ OR MU-
 581
 582 C-----------------------------------------------------------------------
 583
 584       INCLUDE 'HERWIG61.INC'
 585
 586       DOUBLE PRECISION HWR,QSQMAX,QSQMIN,A,B,B1,B2,B3,B4,B5,B6,B7,B8,
 587
 588      $ R,AA,T0,T1,C1,C2,T,Z(2),QSQ(2),PHI(2),C
 589
 590       INTEGER IHEP,I,J
 591
 592       EXTERNAL HWR
 593
 594       SAVE Z,QSQ,PHI
 595
 596 C---IF ZMXISR IS ZERO, THERE CAN BE NO ISR
 597
 598       IF (ZMXISR.EQ.ZERO .OR. (IPRO.GT.3.AND.IPRO.NE.6)) RETURN
 599
 600 C---CHECK CONSISTENCY OF TMNISR AND ZMXISR
 601
 602       IF (ZMXISR**2.LT.TMNISR) CALL HWWARN('HWEISR',200,*999)
 603
 604 C---CALCULATE VIRTUALITY LIMITS
 605
 606       QSQMAX=4*PHEP(4,IHEP)**2
 607
 608       QSQMIN=PHEP(5,IHEP)**2
 609
 610 C---AND THEREFORE THE Z DEPENDENCE
 611
 612       A=ALPHEM/PIFAC
 613
 614       B=A*(LOG(QSQMAX/QSQMIN)-1)
 615
 616 C---DECIDE HOW MUCH WEIGHT TO GIVE THE Z RESONANCE
 617
 618       IF (IHEP.EQ.1) THEN
 619
 620         IF (IPRO.EQ.1.OR.IPRO.EQ.6) THEN
 621
 622           AA=10
 623
 624         ELSEIF (IPRO.EQ.2) THEN
 625
 626           AA=0
 627
 628         ELSEIF (IPRO.EQ.3) THEN
 629
 630           AA=1
 631
 632         ELSE
 633
 634           RETURN
 635
 636         ENDIF
 637
 638         T0=RMASS(200)**2/QSQMAX
 639
 640         T1=GAMZ*RMASS(200)/QSQMAX
 641
 642         IF (T0.GT.ONE) THEN
 643
 644           T0=0
 645
 646           AA=0
 647
 648         ENDIF
 649
 650         AA=AA*(1-T0)
 651
 652 C---GENERATE A T VALUE BETWEEN TMNISR AND 1 ACCORDING TO:
 653
 654 C   ( b**2*log(zmxisr**2/t)/t + 2*b*(1-(1-zmxisr)**b)*((1-t)**(2*b-1)+1/t
 655
 656 C     +(1-t0)**(2b-1)*aa*t1/((t-t0)**2+t1**2)) ) *theta(zmxisr**2-t)
 657
 658 C  +( 2*b*(1-zmxisr)**b*((1-t)**(b-1)+1/t
 659
 660 C     +(1-t0)**(b-1)*aa*t1/((t-t0)**2+t1**2))  ) *theta(zmxisr-t)
 661
 662 C  +( (1-zmxisr)**(2*b)                        ) *delta(1-t)
 663
 664         B1=(1-ZMXISR)**(2*B)
 665
 666         B2=B1+2*(1-ZMXISR)**B*((1-TMNISR)**B-(1-ZMXISR)**B)
 667
 668         B3=B2+2*B*(1-ZMXISR)**B*LOG(ZMXISR/TMNISR)
 669
 670         B4=B3+2*B*(1-ZMXISR)**B*AA*(1-T0)**(B-1)
 671
 672      $       *(ATAN((ZMXISR-T0)/T1)-ATAN((TMNISR-T0)/T1))
 673
 674         B5=B4+(1-(1-ZMXISR)**B)*((1-TMNISR)**(2*B)-(1-ZMXISR**2)**(2*B))
 675
 676         B6=B5+2*B*(1-(1-ZMXISR)**B)*LOG(ZMXISR**2/TMNISR)
 677
 678         B7=B6+B**2*LOG(ZMXISR**2/TMNISR)**2/2
 679
 680         B8=B7+2*B*(1-(1-ZMXISR)**B)*AA*(1-T0)**(2*B-1)
 681
 682      $       *(ATAN((ZMXISR**2-T0)/T1)-ATAN((TMNISR-T0)/T1))
 683
 684         R=B8*HWR()
 685
 686         IF (R.LE.B1) THEN
 687
 688 C---NEITHER EMITS
 689
 690           T=1
 691
 692           GAMWT=GAMWT*B8/B1
 693
 694           Z(1)=1
 695
 696         ELSEIF (R.LE.B4) THEN
 697
 698 C---ONE EMITS
 699
 700           IF (R.LE.B2) THEN
 701
 702             R=(R-B1)/(B2-B1)
 703
 704             T=1-(1-TMNISR)*(1-R*(1-((1-ZMXISR)/(1-TMNISR))**B))**(1/B)
 705
 706           ELSEIF (R.LE.B3) THEN
 707
 708             R=(R-B2)/(B3-B2)
 709
 710             T=(TMNISR/ZMXISR)**R*ZMXISR
 711
 712           ELSE
 713
 714             R=(R-B3)/(B4-B3)
 715
 716             T=T0+T1*TAN(
 717
 718      $           ATAN((ZMXISR-T0)/T1)*R+ATAN((TMNISR-T0)/T1)*(1-R))
 719
 720           ENDIF
 721
 722           GAMWT=GAMWT*B8/(2*B*(1-ZMXISR)**B*((1-T)**(B-1)+1/T+
 723
 724      $         (1-T0)**(B-1)*AA*T1/((T-T0)**2+T1**2)))
 725
 726           Z(1)=1
 727
 728           IF (HWR().GT.HALF) Z(1)=T
 729
 730           GAMWT=GAMWT*2
 731
 732         ELSE
 733
 734 C---BOTH EMIT
 735
 736           IF (R.LE.B5) THEN
 737
 738             R=(R-B4)/(B5-B4)
 739
 740             T=1-(1-TMNISR)*
 741
 742      $           (1-R*(1-((1-ZMXISR**2)/(1-TMNISR))**(2*B)))**(.5/B)
 743
 744           ELSEIF (R.LE.B6) THEN
 745
 746             R=(R-B5)/(B6-B5)
 747
 748             T=(TMNISR/ZMXISR**2)**R*ZMXISR**2
 749
 750           ELSEIF (R.LE.B7) THEN
 751
 752             R=(R-B6)/(B7-B6)
 753
 754             T=(TMNISR/ZMXISR**2)**SQRT(R)*ZMXISR**2
 755
 756           ELSE
 757
 758             R=(R-B7)/(B8-B7)
 759
 760             T=T0+T1*TAN(
 761
 762      $           ATAN((ZMXISR**2-T0)/T1)*R+ATAN((TMNISR-T0)/T1)*(1-R))
 763
 764           ENDIF
 765
 766           GAMWT=GAMWT*B8/(B**2*LOG(ZMXISR**2/T)/T
 767
 768      $         + 2*B*(1-(1-ZMXISR)**B)*((1-T)**(2*B-1)+1/T+
 769
 770      $         (1-T0)**(B-1)*AA*T1/((T-T0)**2+T1**2)))
 771
 772 C---GENERATE A Z VALUE BETWEEN T/ZMXISR AND ZMXISR ACCORDING TO:
 773
 774 C   1/z+(1-z)**(b-1)+t/z**2*(1-t/z)**(b-1)
 775
 776           C1=LOG(ZMXISR**2/T)
 777
 778           C2=C1+2/B*((1-T/ZMXISR)**B-(1-ZMXISR)**B)
 779
 780           IF (C2.GT.ZERO) THEN
 781
 782             R=C2*HWR()
 783
 784             IF (R.LE.C1) THEN
 785
 786               Z(1)=(T/ZMXISR**2)**HWR()*ZMXISR
 787
 788             ELSE
 789
 790               Z(1)=1-(1-T/ZMXISR)*
 791
 792      $             (1-HWR()*(1-((1-ZMXISR)/(1-T/ZMXISR))**B))**(1/B)
 793
 794               IF (2*R.LE.C2+C1) Z(1)=T/Z(1)
 795
 796             ENDIF
 797
 798           ELSE
 799
 800             Z(1)=SQRT(T)
 801
 802           ENDIF
 803
 804           GAMWT=GAMWT*C2/Z(1)
 805
 806      $         /(1/Z(1)+(1-Z(1))**(B-1)+T/Z(1)**2*(1-T/Z(1))**(B-1))
 807
 808         ENDIF
 809
 810 C---INCLUDE DISTRIBUTION FUNCTIONS
 811
 812         Z(2)=T/Z(1)
 813
 814         DO 10 I=1,2
 815
 816           IF (Z(I).GT.ZMXISR) THEN
 817
 818             Z(I)=1
 819
 820             GAMWT=GAMWT*(1-ZMXISR)**B*EXP(3*B/4)*(1-B**2*PIFAC**2/12)
 821
 822           ELSE
 823
 824             GAMWT=GAMWT*(B*(1-Z(I))**(B-1)*(1+Z(I)**2)/2
 825
 826      $           *EXP(B*Z(I)/2*(1+Z(I)/2))*(1-B**2*PIFAC**2/12)
 827
 828      $           +B**2/8*((1+Z(I))*((1+Z(I))**2+3*LOG(Z(I)))
 829
 830      $           -4*LOG(Z(I))/(1-Z(I))))
 831
 832           ENDIF
 833
 834   10    CONTINUE
 835
 836 C---CHOOSE BOTH QSQ VALUES
 837
 838         DO 30 I=1,2
 839
 840           IF (Z(I).GT.ZMXISR .OR. COLISR) THEN
 841
 842             QSQ(I)=0
 843
 844           ELSE
 845
 846             J=3-I
 847
 848 C---ACCORDING TO 1/(QSQ+QSQMIN) FROM 0 TO (1-Z)*(T/(Z+T))*QSQMAX
 849
 850  20         QSQ(I)=(((1-Z(I))*(T/(Z(I)+T))
 851
 852      $           *QSQMAX/QSQMIN+1)**HWR()-1)*QSQMIN
 853
 854 C---AND REJECT TO QSQ/(QSQ+QSQMIN)**2
 855
 856             IF (HWR()*(QSQ(I)+QSQMIN).GT.QSQ(I)) GOTO 20
 857
 858           ENDIF
 859
 860  30     CONTINUE
 861
 862 C---CHOOSE BOTH AZIMUTHS
 863
 864         PHI(1)=HWR()*2*PIFAC
 865
 866         PHI(2)=HWR()*2*PIFAC
 867
 868 C---USE S-HAT PRESCRIPTION TO MODIFY Z VALUES
 869
 870         I=0
 871
 872         IF ((1-Z(1))*QSQ(1).GT.(1-Z(2))*QSQ(2)) I=1
 873
 874         IF ((1-Z(2))*QSQ(2).GT.(1-Z(1))*QSQ(1)) I=2
 875
 876         IF (I.GT.0) THEN
 877
 878           J=3-I
 879
 880           Z(I)=Z(I)+QSQ(I)/QSQMAX
 881
 882           IF (QSQ(J).GT.ZERO) THEN
 883
 884             Z(J)=((QSQ(I)*QSQMAX+QSQ(J)*QSQMAX
 885
 886      $           -QSQ(I)*QSQ(J))/QSQMAX**2+T)/Z(I)
 887
 888             C=COS(PHI(1)-PHI(2))*SQRT(QSQ(1)*QSQ(2))/QSQMAX
 889
 890             Z(J)=Z(J)+(-2*C**2*(1-Z(I))+2*C*SQRT((1-Z(I))
 891
 892      $           *(C**2*(1-Z(I))+Z(I)**2*(1-Z(J)))))/Z(I)**2
 893
 894           ENDIF
 895
 896         ENDIF
 897
 898       ELSEIF (IHEP.EQ.2) THEN
 899
 900 C---EVERYTHING WAS GENERATED LAST TIME
 901
 902       ELSE
 903
 904 C---ROUTINE CALLED UNEXPECTEDLY
 905
 906         CALL HWWARN('HWEISR',201,*999)
 907
 908       ENDIF
 909
 910 C---IF Z IS TOO LARGE THERE IS NO EMISSION
 911
 912       IF (Z(IHEP).GT.ZMXISR) RETURN
 913
 914 C---PUT NEW LEPTON IN EVENT RECORD
 915
 916       NHEP=NHEP+1
 917
 918       IDHW(NHEP)=IDHW(IHEP)
 919
 920       IDHEP(NHEP)=IDHEP(IHEP)
 921
 922       ISTHEP(NHEP)=3
 923
 924       JMOHEP(1,NHEP)=IHEP
 925
 926       JMOHEP(2,NHEP)=0
 927
 928       JDAHEP(1,NHEP)=0
 929
 930       JDAHEP(2,NHEP)=0
 931
 932       JDAHEP(1,IHEP)=NHEP
 933
 934 C---AND OUTGOING PHOTON
 935
 936       NHEP=NHEP+1
 937
 938       IDHW(NHEP)=59
 939
 940       IDHEP(NHEP)=22
 941
 942       ISTHEP(NHEP)=1
 943
 944       JMOHEP(1,NHEP)=IHEP
 945
 946       JMOHEP(2,NHEP)=0
 947
 948       JDAHEP(1,NHEP)=0
 949
 950       JDAHEP(2,NHEP)=0
 951
 952       JDAHEP(2,IHEP)=NHEP
 953
 954 C---RECONSTRUCT PHOTON KINEMATICS (Z IS LIGHT-CONE MOMENTUM FRACTION)
 955
 956       PHEP(1,NHEP)=SQRT(QSQ(IHEP)*(1-Z(IHEP)))*COS(PHI(IHEP))
 957
 958       PHEP(2,NHEP)=SQRT(QSQ(IHEP)*(1-Z(IHEP)))*SIN(PHI(IHEP))
 959
 960       PHEP(3,NHEP)=(1-Z(IHEP))*PHEP(4,IHEP)-QSQ(IHEP)/(4*PHEP(4,IHEP))
 961
 962       IF (IHEP.EQ.2) PHEP(3,NHEP)=-PHEP(3,NHEP)
 963
 964       PHEP(4,NHEP)=(1-Z(IHEP))*PHEP(4,IHEP)+QSQ(IHEP)/(4*PHEP(4,IHEP))
 965
 966       PHEP(5,NHEP)=0
 967
 968 C---AND LEPTON
 969
 970       CALL HWVDIF(4,PHEP(1,IHEP),PHEP(1,NHEP),PHEP(1,NHEP-1))
 971
 972       CALL HWUMAS(PHEP(1,NHEP-1))
 973
 974 C---UPDATE OVERALL CM FRAME
 975
 976       JMOHEP(IHEP,3)=NHEP-1
 977
 978       CALL HWVDIF(4,PHEP(1,3),PHEP(1,IHEP),PHEP(1,3))
 979
 980       CALL HWVSUM(4,PHEP(1,NHEP-1),PHEP(1,3),PHEP(1,3))
 981
 982       CALL HWUMAS(PHEP(1,3))
 983
 984  999  END
 985
 986 CDECK  ID>, HWEONE.
 987
 988 *CMZ :-        -26/04/91  11.11.55  by  Bryan Webber
 989
 990 *-- Author :    Bryan Webber
 991
 992 C-----------------------------------------------------------------------
 993
 994       SUBROUTINE HWEONE
 995
 996 C-----------------------------------------------------------------------
 997
 998 C     SETS UP 2->1 (COLOUR SINGLET) HARD SUBPROCESS
 999
1000 C-----------------------------------------------------------------------
1001
1002       INCLUDE 'HERWIG61.INC'
1003
1004       DOUBLE PRECISION PA
1005
1006       INTEGER ICMF,I,IBM,IHEP
1007
1008 C---INCOMING LINES
1009
1010       ICMF=NHEP+3
1011
1012       DO 15 I=1,2
1013
1014       IBM=I
1015
1016 C---FIND BEAM AND TARGET
1017
1018       IF (JDAHEP(1,I).NE.0) IBM=JDAHEP(1,I)
1019
1020       IHEP=NHEP+I
1021
1022       IDHW(IHEP)=IDN(I)
1023
1024       IDHEP(IHEP)=IDPDG(IDN(I))
1025
1026       ISTHEP(IHEP)=110+I
1027
1028       JMOHEP(1,IHEP)=ICMF
1029
1030       JMOHEP(I,ICMF)=IHEP
1031
1032       JDAHEP(1,IHEP)=ICMF
1033
1034 C---SPECIAL - IF INCOMING PARTON IS INCOMING BEAM THEN COPY IT
1035
1036       IF (XX(I).EQ.ONE.AND.IDHW(IBM).EQ.IDN(I)) THEN
1037
1038         CALL HWVEQU(5,PHEP(1,IBM),PHEP(1,IHEP))
1039
1040         IF (I.EQ.2) PHEP(3,IHEP)=-PHEP(3,IHEP)
1041
1042       ELSE
1043
1044         PHEP(1,IHEP)=0.
1045
1046         PHEP(2,IHEP)=0.
1047
1048         PHEP(5,IHEP)=RMASS(IDN(I))
1049
1050         PA=XX(I)*(PHEP(4,IBM)+ABS(PHEP(3,IBM)))
1051
1052         PHEP(4,IHEP)=0.5*(PA+PHEP(5,IHEP)**2/PA)
1053
1054         PHEP(3,IHEP)=PA-PHEP(4,IHEP)
1055
1056       ENDIF
1057
1058  15   CONTINUE
1059
1060       PHEP(3,NHEP+2)=-PHEP(3,NHEP+2)
1061
1062 C---HARD CENTRE OF MASS
1063
1064       IDHW(ICMF)=IDCMF
1065
1066       IDHEP(ICMF)=IDPDG(IDCMF)
1067
1068       ISTHEP(ICMF)=110
1069
1070       CALL HWVSUM(4,PHEP(1,NHEP+1),PHEP(1,NHEP+2),PHEP(1,ICMF))
1071
1072       CALL HWUMAS(PHEP(1,ICMF))
1073
1074 C---SET UP COLOUR STRUCTURE LABELS
1075
1076       JMOHEP(2,NHEP+1)=NHEP+2
1077
1078       JDAHEP(2,NHEP+1)=NHEP+2
1079
1080       JMOHEP(2,NHEP+2)=NHEP+1
1081
1082       JDAHEP(2,NHEP+2)=NHEP+1
1083
1084       JDAHEP(1,NHEP+3)=NHEP+3
1085
1086       JDAHEP(2,NHEP+3)=NHEP+3
1087
1088       NHEP=NHEP+3
1089
1090   999 END
1091
1092 CDECK  ID>, HWEPRO.
1093
1094 *CMZ :-        -01/04/99  19.41.18  by  Mike Seymour
1095
1096 *-- Author :    Bryan Webber
1097
1098 C-----------------------------------------------------------------------
1099
1100       SUBROUTINE HWEPRO
1101
1102 C-----------------------------------------------------------------------
1103
1104 C     WHEN NEVHEP=0, CHOOSES X VALUES AND FINDS WEIGHT FOR PROCESS IPROC
1105
1106 C     OTHERWISE, CHOOSES AND LOADS ALL VARIABLES FOR HARD PROCESS
1107
1108 C-----------------------------------------------------------------------
1109
1110       INCLUDE 'HERWIG61.INC'
1111
1112       DOUBLE PRECISION HWR
1113
1114       EXTERNAL HWR
1115
1116       IF (IERROR.NE.0)  RETURN
1117
1118 C---ROUTINE LOOPS BACK TO HERE IF GENERATED WEIGHT WAS NOT ACCEPTED
1119
1120    10 GENEV=.FALSE.
1121
1122 C---FSTWGT IS .TRUE. DURING FIRST CALL TO HARD PROCESS ROUTINE
1123
1124       FSTWGT=NWGTS.EQ.0
1125
1126 C---FSTEVT IS .TRUE. THROUGHOUT THE FIRST EVENT
1127
1128       FSTEVT=NEVHEP.EQ.1
1129
1130 C---SET COLOUR CORRECTION TO FALSE
1131
1132       COLUPD = .FALSE.
1133
1134       HRDCOL(1,1)=0
1135
1136       HRDCOL(1,3)=0
1137
1138 C---SET UP INITIAL STATE
1139
1140       NHEP=1
1141
1142       ISTHEP(NHEP)=101
1143
1144       PHEP(1,NHEP)=0.
1145
1146       PHEP(2,NHEP)=0.
1147
1148       PHEP(3,NHEP)=PBEAM1
1149
1150       PHEP(4,NHEP)=EBEAM1
1151
1152       PHEP(5,NHEP)=RMASS(IPART1)
1153
1154       JMOHEP(1,NHEP)=0
1155
1156       JMOHEP(2,NHEP)=0
1157
1158       JDAHEP(1,NHEP)=0
1159
1160       JDAHEP(2,NHEP)=0
1161
1162       IDHW(NHEP)=IPART1
1163
1164       IDHEP(NHEP)=IDPDG(IPART1)
1165
1166       NHEP=NHEP+1
1167
1168       ISTHEP(NHEP)=102
1169
1170       PHEP(1,NHEP)=0.
1171
1172       PHEP(2,NHEP)=0.
1173
1174       PHEP(3,NHEP)=-PBEAM2
1175
1176       PHEP(4,NHEP)=EBEAM2
1177
1178       PHEP(5,NHEP)=RMASS(IPART2)
1179
1180       JMOHEP(1,NHEP)=0
1181
1182       JMOHEP(2,NHEP)=0
1183
1184       JDAHEP(1,NHEP)=0
1185
1186       JDAHEP(2,NHEP)=0
1187
1188       IDHW(NHEP)=IPART2
1189
1190       IDHEP(NHEP)=IDPDG(IPART2)
1191
1192 C---NEXT ENTRY IS OVERALL CM FRAME
1193
1194       NHEP=NHEP+1
1195
1196       IDHW(NHEP)=14
1197
1198       IDHEP(NHEP)=0
1199
1200       ISTHEP(NHEP)=103
1201
1202       JMOHEP(1,NHEP)=NHEP-2
1203
1204       JMOHEP(2,NHEP)=NHEP-1
1205
1206       JDAHEP(1,NHEP)=0
1207
1208       JDAHEP(2,NHEP)=0
1209
1210       CALL HWVSUM(4,PHEP(1,NHEP-1),PHEP(1,NHEP-2),PHEP(1,NHEP))
1211
1212       CALL HWUMAS(PHEP(1,NHEP))
1213
1214 C Select a primary interaction point
1215
1216       IF (PIPSMR) THEN
1217
1218         CALL HWRPIP
1219
1220       ELSE
1221
1222         CALL HWVZRO(4,VTXPIP)
1223
1224       ENDIF
1225
1226       CALL HWVEQU(3,VTXPIP,VHEP(1,NHEP))
1227
1228       VHEP(4,NHEP)=0.0
1229
1230 C---GENERATE PHOTONS (WEIZSACKER-WILLIAMS APPROX)
1231
1232 C   FOR HADRONIC PROCESSES WITH LEPTON BEAMS
1233
1234       GAMWT=ONE
1235
1236       IF (IPRO.GT.10.AND.IPRO.LT.90) THEN
1237
1238         IF (ABS(IDHEP(1)).EQ.11.OR.ABS(IDHEP(1)).EQ.13)
1239
1240      &       CALL HWEGAM(1,ZERO, ONE,.FALSE.)
1241
1242         IF (ABS(IDHEP(2)).EQ.11.OR.ABS(IDHEP(2)).EQ.13)
1243
1244      &       CALL HWEGAM(2,ZERO, ONE,.FALSE.)
1245
1246       ELSEIF (IPRO.GE.90) THEN
1247
1248         IF (ABS(IDHEP(2)).EQ.11.OR.ABS(IDHEP(2)).EQ.13)
1249
1250      &       CALL HWEGAM(2,ZERO, ONE,.FALSE.)
1251
1252       ENDIF
1253
1254 C---GENERATE ISR PHOTONS FOR LEPTONIC PROCESSES
1255
1256       IF (IPRO.LT.10) THEN
1257
1258         CALL HWEISR(1)
1259
1260         CALL HWEISR(2)
1261
1262       ENDIF
1263
1264 C---IF USER LIMITS WERE TOO TIGHT, MIGHT NOT BE ANY PHASE-SPACE
1265
1266       IF (GAMWT.LE.ZERO) GOTO 30
1267
1268 C---IF CMF HAS ACQUIRED A TRANSVERSE BOOST, OR USER REQUESTS IT ANYWAY,
1269
1270 C   BOOST EVENT RECORD BACK TO CMF
1271
1272       IF (PHEP(1,3)**2+PHEP(2,3)**2.GT.ZERO .OR. USECMF) CALL HWUBST(1)
1273
1274 C---ROUTINE LOOPS BACK TO HERE IF GENERATED WEIGHT WAS ACCEPTED
1275
1276    20 CONTINUE
1277
1278 C---IPRO=MOD(IPROC/100,100)
1279
1280       IF (IPRO.EQ.1) THEN
1281
1282         IF (IPROC.LT.110.OR.IPROC.GE.120) THEN
1283
1284 C--- E+E- -> Q-QBAR OR L-LBAR
1285
1286           CALL HWHEPA
1287
1288         ELSE
1289
1290 C--- E+E- -> Q-QBAR-GLUON
1291
1292           CALL HWHEPG
1293
1294         ENDIF
1295
1296       ELSEIF (IPRO.EQ.2) THEN
1297
1298 C--- E+E- -> W+ W-
1299
1300         CALL HWHEWW
1301
1302       ELSEIF (IPRO.EQ.3) THEN
1303
1304 C---E+E- -> Z H
1305
1306         CALL HWHIGZ
1307
1308       ELSEIF (IPRO.EQ.4) THEN
1309
1310 C---E+E- -> NUEB NUE H
1311
1312         CALL HWHIGW
1313
1314       ELSEIF (IPRO.EQ.5 .AND. IPROC.LT.550) THEN
1315
1316 C---EE -> EE GAMGAM -> EE FFBAR/WW
1317
1318         CALL HWHEGG
1319
1320       ELSEIF (IPRO.EQ.5) THEN
1321
1322 C---EE -> ENU GAMW -> ENU FF'BAR/WZ
1323
1324         CALL HWHEGW
1325
1326       ELSEIF (IPRO.EQ.6) THEN
1327
1328 C---EE -> FOUR JETS
1329
1330         CALL HWH4JT
1331
1332       ELSEIF (IPRO.EQ.13) THEN
1333
1334 C---GAMMA/Z0/Z' DRELL-YAN PROCESS
1335
1336         CALL HWHDYP
1337
1338       ELSEIF (IPRO.EQ.14) THEN
1339
1340 C---W+/- PRODUCTION VIA DRELL-YAN PROCESS
1341
1342         CALL HWHWPR
1343
1344       ELSEIF (IPRO.EQ.15) THEN
1345
1346 C---QCD HARD 2->2 PROCESSES
1347
1348         CALL HWHQCD
1349
1350       ELSEIF (IPRO.EQ.16) THEN
1351
1352 C---HIGGS PRODUCTION VIA GLUON FUSION
1353
1354         CALL HWHIGS
1355
1356       ELSEIF (IPRO.EQ.17) THEN
1357
1358 C---QCD HEAVY FLAVOUR PRODUCTION
1359
1360         CALL HWHHVY
1361
1362       ELSEIF (IPRO.EQ.18) THEN
1363
1364 C---QCD DIRECT PHOTON + JET PRODUCTION
1365
1366         CALL HWHPHO
1367
1368       ELSEIF (IPRO.EQ.19) THEN
1369
1370 C---HIGGS PRODUCTION VIA W FUSION
1371
1372         CALL HWHIGW
1373
1374       ELSEIF (IPRO.EQ.20) THEN
1375
1376 C---TOP PRODUCTION FROM W EXCHANGE
1377
1378         CALL HWHWEX
1379
1380       ELSEIF (IPRO.EQ.21) THEN
1381
1382 C---VECTOR BOSON + JET PRODUCTION
1383
1384         CALL HWHV1J
1385
1386       ELSEIF (IPRO.EQ.22) THEN
1387
1388 C QCD direct photon pair production
1389
1390         CALL HWHPH2
1391
1392       ELSEIF (IPRO.EQ.23) THEN
1393
1394 C QCD Higgs plus jet production
1395
1396         CALL HWHIGJ
1397
1398       ELSEIF (IPRO.EQ.24) THEN
1399
1400 C---COLOUR-SINGLET EXCHANGE
1401
1402         CALL HWHSNG
1403
1404       ELSEIF (IPRO.EQ.30) THEN
1405
1406 C---HADRON-HADRON SUSY PROCESSES
1407
1408         CALL HWHSSP
1409
1410       ELSEIF(IPRO.EQ.40.OR.IPRO.EQ.41) THEN
1411
1412 C---HADRON-HADRON R-PARITY VIOLATING SUSY PROCESSES
1413
1414         CALL HWHRSP
1415
1416       ELSEIF (IPRO.EQ.50) THEN
1417
1418 C Point-like photon two-jet production
1419
1420         CALL HWHPPT
1421
1422       ELSEIF (IPRO.EQ.51) THEN
1423
1424 C Point-like photon/QCD heavy flavour pair production
1425
1426         CALL HWHPPH
1427
1428       ELSEIF (IPRO.EQ.52) THEN
1429
1430 C Point-like photon/QCD heavy flavour single excitation
1431
1432         CALL HWHPPE
1433
1434       ELSEIF (IPRO.EQ.53) THEN
1435
1436 C Compton scattering of point-like photon and (anti)quark
1437
1438         CALL HWHPQS
1439
1440       ELSEIF (IPRO.EQ.55) THEN
1441
1442 C Point-like photon/higher twist meson production
1443
1444         CALL HWHPPM
1445
1446       ELSEIF (IPRO.EQ.60) THEN
1447
1448 C---QPM GAMMA-GAMMA-->QQBAR
1449
1450         CALL HWHQPM
1451
1452       ELSEIF (IPRO.GE.70.AND.IPRO.LE.79) THEN
1453
1454 C---BARYON-NUMBER VIOLATION, AND OTHER MULTI-W PRODUCTION PROCESSES
1455
1456         CALL HVHBVI
1457
1458       ELSEIF (IPRO.EQ.80) THEN
1459
1460 C---MINIMUM-BIAS: NO HARD SUBPROCESS
1461
1462 C   FIND WEIGHT
1463
1464         CALL HWMWGT
1465
1466       ELSEIF (IPRO.EQ.90) THEN
1467
1468 C---DEEP INELASTIC
1469
1470         CALL HWHDIS
1471
1472       ELSEIF(IPRO.EQ.91) THEN
1473
1474 C---BOSON - GLUON(QUARK) FUSION -->  ANTIQUARK(GLUON) + QUARK
1475
1476         CALL HWHBGF
1477
1478       ELSEIF(IPRO.EQ.92) THEN
1479
1480 C---DEEP INELASTIC WITH EXTRA JET: OBSOLETE PROCESS
1481
1482         WRITE (6,40)
1483
1484  40     FORMAT (1X,' IPROC=92** is no longer supported.'
1485
1486      &         /1X,' Please use IPROC=91** instead.')
1487
1488         CALL HWWARN('HWEPRO',500,*999)
1489
1490       ELSEIF(IPRO.EQ.95) THEN
1491
1492 C---HIGGS PRODUCTION VIA W FUSION IN E P
1493
1494         CALL HWHIGW
1495
1496       ELSE
1497
1498 C---UNKNOWN PROCESS
1499
1500         CALL HWWARN('HWEPRO',102,*999)
1501
1502       ENDIF
1503
1504  30   IF (GENEV) THEN
1505
1506         IF (NOWGT) EVWGT=AVWGT
1507
1508         ISTAT=10
1509
1510         RETURN
1511
1512       ELSE
1513
1514 C---IF AN EVENT IS CANCELLED BEFORE IT IS GENERATED, GIVE IT ZERO WEIGHT
1515
1516         IF (IERROR.NE.0) THEN
1517
1518           EVWGT=ZERO
1519
1520           IERROR=0
1521
1522         ENDIF
1523
1524         EVWGT=EVWGT*GAMWT
1525
1526         NWGTS=NWGTS+1
1527
1528         WGTSUM=WGTSUM+EVWGT
1529
1530         WSQSUM=WSQSUM+EVWGT**2
1531
1532         IF (EVWGT.GT.WBIGST) THEN
1533
1534           WBIGST=EVWGT
1535
1536           IF (NOWGT.AND.WBIGST.GT.WGTMAX) THEN
1537
1538             IF (NEVHEP.NE.0) CALL HWWARN('HWEPRO',1,*999)
1539
1540             WGTMAX=WBIGST*1.1
1541
1542             WRITE (6,99) WGTMAX
1543
1544           ENDIF
1545
1546         ELSEIF (EVWGT.LT.ZERO) THEN
1547
1548           IF (EVWGT.LT.-1.D-9) CALL HWWARN('HWEPRO',3,*999)
1549
1550           EVWGT=0.
1551
1552         ENDIF
1553
1554         IF (NEVHEP.NE.0) THEN
1555
1556 C---LOW EFFICIENCY WARNINGS:
1557
1558 C   WARN AT 10*EFFMIN, STOP AT EFFMIN
1559
1560           IF (10*EFFMIN*NWGTS.GT.NEVHEP) THEN
1561
1562             IF (EFFMIN*NWGTS.GT.NEVHEP) CALL HWWARN('HWEPRO',200,*999)
1563
1564             IF (EFFMIN.GT.ZERO) THEN
1565
1566               IF (MOD(NWGTS,INT(10/EFFMIN)).EQ.0) THEN
1567
1568                 CALL HWWARN('HWEPRO',2,*999)
1569
1570                 WRITE (6,98) WGTMAX
1571
1572               ENDIF
1573
1574             ENDIF
1575
1576           ENDIF
1577
1578           IF (NOWGT) THEN
1579
1580             GENEV=EVWGT.GT.WGTMAX*HWR()
1581
1582           ELSE
1583
1584             GENEV=EVWGT.NE.ZERO
1585
1586           ENDIF
1587
1588           IF (GENEV) GOTO 20
1589
1590           GOTO 10
1591
1592         ENDIF
1593
1594       ENDIF
1595
1596  98   FORMAT(10X,'    MAXIMUM WEIGHT =',1PG24.16)
1597
1598  99   FORMAT(10X,'NEW MAXIMUM WEIGHT =',1PG24.16)
1599
1600   999 END
1601
1602 CDECK  ID>, HWETWO.
1603
1604 *CMZ :-        -26/04/91  11.11.55  by  Bryan Webber
1605
1606 *-- Author :    Bryan Webber
1607
1608 C-----------------------------------------------------------------------
1609
1610       SUBROUTINE HWETWO
1611
1612 C-----------------------------------------------------------------------
1613
1614 C     SETS UP 2->2 HARD SUBPROCESS
1615
1616 C-----------------------------------------------------------------------
1617
1618       INCLUDE 'HERWIG61.INC'
1619
1620       DOUBLE PRECISION HWUMBW,HWUPCM,PA,PCM
1621
1622       INTEGER ICMF,IBM,I,J,K,IHEP,NTRY
1623
1624       EXTERNAL HWUPCM
1625
1626 C---INCOMING LINES
1627
1628       ICMF=NHEP+3
1629
1630       DO 15 I=1,2
1631
1632       IBM=I
1633
1634 C---FIND BEAM AND TARGET
1635
1636       IF (JDAHEP(1,I).NE.0) IBM=JDAHEP(1,I)
1637
1638       IHEP=NHEP+I
1639
1640       IDHW(IHEP)=IDN(I)
1641
1642       IDHEP(IHEP)=IDPDG(IDN(I))
1643
1644       ISTHEP(IHEP)=110+I
1645
1646       JMOHEP(1,IHEP)=ICMF
1647
1648       JMOHEP(I,ICMF)=IHEP
1649
1650       JDAHEP(1,IHEP)=ICMF
1651
1652 C---SPECIAL - IF INCOMING PARTON IS INCOMING BEAM THEN COPY IT
1653
1654       IF (XX(I).EQ.ONE.AND.IDHW(IBM).EQ.IDN(I)) THEN
1655
1656         CALL HWVEQU(5,PHEP(1,IBM),PHEP(1,IHEP))
1657
1658         IF (I.EQ.2) PHEP(3,IHEP)=-PHEP(3,IHEP)
1659
1660       ELSE
1661
1662         PHEP(1,IHEP)=0.
1663
1664         PHEP(2,IHEP)=0.
1665
1666         PHEP(5,IHEP)=RMASS(IDN(I))
1667
1668         PA=XX(I)*(PHEP(4,IBM)+ABS(PHEP(3,IBM)))
1669
1670         PHEP(4,IHEP)=0.5*(PA+PHEP(5,IHEP)**2/PA)
1671
1672         PHEP(3,IHEP)=PA-PHEP(4,IHEP)
1673
1674       ENDIF
1675
1676  15   CONTINUE
1677
1678       PHEP(3,NHEP+2)=-PHEP(3,NHEP+2)
1679
1680 C---HARD CENTRE OF MASS
1681
1682       IDHW(ICMF)=IDCMF
1683
1684       IDHEP(ICMF)=IDPDG(IDCMF)
1685
1686       ISTHEP(ICMF)=110
1687
1688       CALL HWVSUM(4,PHEP(1,NHEP+1),PHEP(1,NHEP+2),PHEP(1,ICMF))
1689
1690       CALL HWUMAS(PHEP(1,ICMF))
1691
1692 C---OUTGOING LINES
1693
1694       NTRY=0
1695
1696  19   CONTINUE
1697
1698       DO 20 I=3,4
1699
1700       IHEP=NHEP+I+1
1701
1702       IDHW(IHEP)=IDN(I)
1703
1704       IDHEP(IHEP)=IDPDG(IDN(I))
1705
1706       ISTHEP(IHEP)=110+I
1707
1708       JMOHEP(1,IHEP)=ICMF
1709
1710       JDAHEP(I-2,ICMF)=IHEP
1711
1712    20 PHEP(5,IHEP)=HWUMBW(IDN(I))
1713
1714       PCM=HWUPCM(PHEP(5,NHEP+3),PHEP(5,NHEP+4),PHEP(5,NHEP+5))
1715
1716       IF (PCM.LT.ZERO) THEN
1717
1718         NTRY=NTRY+1
1719
1720         IF (NTRY.LE.NETRY) GO TO 19
1721
1722         CALL HWWARN('HWETWO',103,*999)
1723
1724       ENDIF
1725
1726       IHEP=NHEP+4
1727
1728       PHEP(4,IHEP)=SQRT(PCM**2+PHEP(5,IHEP)**2)
1729
1730       PHEP(3,IHEP)=PCM*COSTH
1731
1732       PHEP(1,IHEP)=SQRT((PCM+PHEP(3,IHEP))*(PCM-PHEP(3,IHEP)))
1733
1734       CALL HWRAZM(PHEP(1,IHEP),PHEP(1,IHEP),PHEP(2,IHEP))
1735
1736       CALL HWULOB(PHEP(1,NHEP+3),PHEP(1,IHEP),PHEP(1,IHEP))
1737
1738       CALL HWVDIF(4,PHEP(1,NHEP+3),PHEP(1,IHEP),PHEP(1,NHEP+5))
1739
1740 C---SET UP COLOUR STRUCTURE LABELS
1741
1742       DO 30 I=1,4
1743
1744       J=I
1745
1746       IF (J.GT.2) J=J+1
1747
1748       K=ICO(I)
1749
1750       IF (K.GT.2) K=K+1
1751
1752       JMOHEP(2,NHEP+J)=NHEP+K
1753
1754    30 JDAHEP(2,NHEP+K)=NHEP+J
1755
1756       NHEP=NHEP+5
1757
1758   999 END
1759
1760 CDECK  ID>, HWH4JT.
1761
1762 *CMZ :-        -01/04/99  19.47.55  by  Mike Seymour
1763
1764 *-- Author :    Ian Knowles
1765
1766 C-----------------------------------------------------------------------
1767
1768       SUBROUTINE HWH4JT
1769
1770 C-----------------------------------------------------------------------
1771
1772 C     Four jet production in e^+e^- annihilation: qqbar+gg & qqbar+qqbar
1773
1774 C     IOP4JT controls the treatment of the colour flow interference term
1775
1776 C     qqbar-gg case:
1777
1778 C     IOP4JT(1)=0 neglect, =1 extreme 2341; =2 extreme 3421
1779
1780 C     qqbar-qqbar (identical quark flavour) case:
1781
1782 C     IOP4JT(2)=0 neglect, =1 extreme 4123; =2 extreme 2143
1783
1784 C
1785
1786 C     Matrix elements based on Ellis Ross & Terrano and Catani & Seymour
1787
1788 C
1789
1790 C     WARNING:  Phase space factor inaccurate for JADE y_cut > 0.14.
1791
1792 C-----------------------------------------------------------------------
1793
1794       INCLUDE 'HERWIG61.INC'
1795
1796       INTEGER LM,LP,IQK,I,J,IDMN,IDMX,ID1,ID2,IST(4)
1797
1798       DOUBLE PRECISION HWR,HWUALF,HWUAEM,HWULDO,HWH4J1,HWH4J2,
1799
1800      & HWH4J4,HWH4J5,HWH4J6,HWH4J7,QNOW,Q2NOW,QLST,SCUT,PSFAC,FACT,X,
1801
1802      & COLA,COLB,COLC,CLF(7,6),P12,P13,P14,P23,P24,P34,FACTR,EP1,EP2,
1803
1804      & EP3,EP4,GG1,GG2,GG12,GG3,GG13,GG23,GGINT,WTGG,QQ,QP,QQINT,QQ1,
1805
1806      & QQ2,WTQQ,WTQP,HCS,WTAB,WTBA,WTOT,RCS
1807
1808      $     ,EF,QF,E(4)
1809
1810       LOGICAL INCLQG(6),INCLQQ(6,6),ORIENT
1811
1812       EXTERNAL HWR,HWUALF,HWUAEM,HWULDO,HWH4J1,HWH4J2,HWH4J4,
1813
1814      & HWH4J5,HWH4J6,HWH4J7
1815
1816       SAVE HCS,QLST,WTQP,WTQQ,WTGG,FACTR,COLA,COLB,COLC,IDMN,IDMX,
1817
1818      & CLF,GG1,GG2,GGINT,INCLQG,INCLQQ,LM,LP,QQ1,QQ2,QQINT,FACT,ORIENT,
1819
1820      & Q2NOW,SCUT
1821
1822       DATA QLST,IST/-1.,113,114,114,114/
1823
1824 C
1825
1826       IF (GENEV) THEN
1827
1828         RCS=HCS*HWR()
1829
1830       ELSE
1831
1832         IF (NHEP+5.GT.NMXHEP) CALL HWWARN('HWH4JT',100,*999)
1833
1834         QNOW=PHEP(5,3)
1835
1836         IF (QNOW.NE.QLST) THEN
1837
1838           QLST=QNOW
1839
1840           Q2NOW=QNOW**2
1841
1842           SCUT=Y4JT*Q2NOW
1843
1844 C Calculate allowed fraction of Phase Space using parameterization
1845
1846           IF (DURHAM) THEN
1847
1848             PSFAC=(1.-6.*Y4JT)**5.50*(1.-173.3*Y4JT*(1.-247.3*Y4JT
1849
1850      &                              *(1.+148.3*Y4JT*(1.+3.913*Y4JT))))
1851
1852      &                              /(1.-8.352*Y4JT*(1.-1102.*Y4JT
1853
1854      &                              *(1.+1603.*Y4JT*(1.+22.99*Y4JT))))
1855
1856           ELSE
1857
1858             PSFAC=(1.-6.*Y4JT)**4.62*(1.-44.72*Y4JT*(1.-176.0*Y4JT
1859
1860      &                              *(1.+102.9*Y4JT*(1.-6.579*Y4JT))))
1861
1862      &                              /(1.-3.392*Y4JT*(1.-946.5*Y4JT
1863
1864      &                              *(1.+423.4*Y4JT*(1.-3.971*Y4JT))))
1865
1866           ENDIF
1867
1868           FACT=GEV2NB*HWUAEM(Q2NOW)**2*CFFAC*FLOAT(NCOLO)*PSFAC
1869
1870      &        /(THREE*16*PIFAC)
1871
1872           COLA=CFFAC
1873
1874           COLB=CFFAC-HALF*CAFAC
1875
1876           COLC=HALF
1877
1878           LM=1
1879
1880           IF (JDAHEP(1,LM).NE.0) LM=JDAHEP(1,LM)
1881
1882           LP=2
1883
1884           IF (JDAHEP(1,LP).NE.0) LP=JDAHEP(1,LP)
1885
1886           IQK=MOD(IPROC,10)
1887
1888           IF (IQK.NE.0) THEN
1889
1890             IDMN=IQK
1891
1892             IDMX=IQK
1893
1894           ELSE
1895
1896             IDMN=1
1897
1898             IDMX=6
1899
1900           ENDIF
1901
1902           DO 10 I=1,6
1903
1904           CALL HWUCFF(11,I,Q2NOW,CLF(1,I))
1905
1906           IF (QNOW.GT.TWO*(RMASS(I)+RMASS(13))) THEN
1907
1908             INCLQG(I)=.TRUE.
1909
1910           ELSE
1911
1912             INCLQG(I)=.FALSE.
1913
1914           ENDIF
1915
1916           DO 10 J=I,6
1917
1918           IF (QNOW.GT.TWO*(RMASS(I)+RMASS(J ))) THEN
1919
1920             INCLQQ(I,J)=.TRUE.
1921
1922             INCLQQ(J,I)=.TRUE.
1923
1924           ELSE
1925
1926             INCLQQ(I,J)=.FALSE.
1927
1928             INCLQQ(J,I)=.FALSE.
1929
1930           ENDIF
1931
1932   10      CONTINUE
1933
1934           IF (MOD(IPROC/10,10).EQ.5) THEN
1935
1936             ORIENT=.FALSE.
1937
1938           ELSE
1939
1940             ORIENT=.TRUE.
1941
1942           ENDIF
1943
1944         ENDIF
1945
1946 C Generate phase space point and check it passes cuts
1947
1948         CALL HWVEQU(5,PHEP(1,3),PHEP(1,NHEP+1))
1949
1950         DO 20 I=2,5
1951
1952   20    PHEP(5,NHEP+I)=0.
1953
1954   30    CALL HWDFOR(PHEP(1,NHEP+1),PHEP(1,NHEP+2),PHEP(1,NHEP+3),
1955
1956      &              PHEP(1,NHEP+4),PHEP(1,NHEP+5))
1957
1958         IF (DURHAM) THEN
1959
1960           P12=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+3))
1961
1962           X=MIN(PHEP(4,NHEP+2)/PHEP(4,NHEP+3),
1963
1964      &          PHEP(4,NHEP+3)/PHEP(4,NHEP+2))*P12
1965
1966           IF (X.GT.SCUT) THEN
1967
1968             P13=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+4))
1969
1970             X=MIN(PHEP(4,NHEP+2)/PHEP(4,NHEP+4),
1971
1972      &            PHEP(4,NHEP+4)/PHEP(4,NHEP+2))*P13
1973
1974             IF (X.GT.SCUT) THEN
1975
1976               P14=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+5))
1977
1978               X=MIN(PHEP(4,NHEP+2)/PHEP(4,NHEP+5),
1979
1980      &              PHEP(4,NHEP+5)/PHEP(4,NHEP+2))*P14
1981
1982               IF (X.GT.SCUT) THEN
1983
1984                 P23=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+4))
1985
1986                 X=MIN(PHEP(4,NHEP+3)/PHEP(4,NHEP+4),
1987
1988      &                PHEP(4,NHEP+4)/PHEP(4,NHEP+3))*P23
1989
1990                 IF (X.GT.SCUT) THEN
1991
1992                   P24=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+5))
1993
1994                   X=MIN(PHEP(4,NHEP+3)/PHEP(4,NHEP+5),
1995
1996      &                  PHEP(4,NHEP+5)/PHEP(4,NHEP+3))*P24
1997
1998                   IF (X.GT.SCUT) THEN
1999
2000                     P34=2*HWULDO(PHEP(1,NHEP+4),PHEP(1,NHEP+5))
2001
2002                     X=MIN(PHEP(4,NHEP+4)/PHEP(4,NHEP+5),
2003
2004      &                    PHEP(4,NHEP+5)/PHEP(4,NHEP+4))*P34
2005
2006                     IF (X.GT.SCUT) GOTO 40
2007
2008                   ENDIF
2009
2010                 ENDIF
2011
2012               ENDIF
2013
2014             ENDIF
2015
2016           ENDIF
2017
2018         ELSE
2019
2020           P12=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+3))
2021
2022           IF (P12.GT.SCUT) THEN
2023
2024             P13=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+4))
2025
2026             IF (P13.GT.SCUT) THEN
2027
2028               P14=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+5))
2029
2030               IF (P14.GT.SCUT) THEN
2031
2032                 P23=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+4))
2033
2034                 IF (P23.GT.SCUT) THEN
2035
2036                   P24=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+5))
2037
2038                   IF (P24.GT.SCUT) THEN
2039
2040                     P34=2*HWULDO(PHEP(1,NHEP+4),PHEP(1,NHEP+5))
2041
2042                     IF (P34.GT.SCUT) GOTO 40
2043
2044                   ENDIF
2045
2046                 ENDIF
2047
2048               ENDIF
2049
2050             ENDIF
2051
2052           ENDIF
2053
2054         ENDIF
2055
2056 C Failed cuts retry
2057
2058         GOTO 30
2059
2060 C Passed cuts: calculate contributions to Matrix Elements
2061
2062   40    EMSCA=SQRT(MIN(P12,P13,P14,P23,P24,P34))
2063
2064         FACTR=FACT*HWUALF(1,EMSCA)**2
2065
2066         IF (ORIENT) THEN
2067
2068           QF=HWULDO(PHEP(1,LP),PHEP(1,3))
2069
2070           EF=Q2NOW/(2*SQRT(QF**2-HWULDO(PHEP(1,LP),PHEP(1,LP))*Q2NOW))
2071
2072           QF=HALF-EF*QF/Q2NOW
2073
2074           DO I=1,4
2075
2076             E(I)=EF*PHEP(I,LP)+QF*PHEP(I,3)
2077
2078           ENDDO
2079
2080           EP1=HWULDO(E,PHEP(1,NHEP+2))
2081
2082           EP2=HWULDO(E,PHEP(1,NHEP+3))
2083
2084           EP3=HWULDO(E,PHEP(1,NHEP+4))
2085
2086           EP4=HWULDO(E,PHEP(1,NHEP+5))
2087
2088         ENDIF
2089
2090 C q-qbar-g-g
2091
2092         GG1=HWH4J1(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2093
2094      &     +HWH4J1(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2095
2096         GG2=HWH4J1(P12,P23,P24,P13,P14,P34,EP2,EP1,EP3,EP4,ORIENT)
2097
2098      &     +HWH4J1(P12,P14,P13,P24,P23,P34,EP1,EP2,EP4,EP3,ORIENT)
2099
2100         GG12=HWH4J2(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2101
2102      &      +HWH4J2(P12,P14,P13,P24,P23,P34,EP1,EP2,EP4,EP3,ORIENT)
2103
2104      &      +HWH4J2(P12,P23,P24,P13,P14,P34,EP2,EP1,EP3,EP4,ORIENT)
2105
2106      &      +HWH4J2(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2107
2108         GG3=HWH4J4(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2109
2110      &     +HWH4J4(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2111
2112         GG13=GG3+HWH4J5(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2113
2114      &          +HWH4J5(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2115
2116         GG23=GG3+HWH4J5(P12,P14,P13,P24,P23,P34,EP1,EP2,EP4,EP3,ORIENT)
2117
2118      &          +HWH4J5(P12,P23,P24,P13,P14,P34,EP2,EP1,EP3,EP4,ORIENT)
2119
2120 C Add up weights
2121
2122         GG1  =COLA*(GG1 +GG13)
2123
2124         GG2  =COLA*(GG2 +GG23)
2125
2126         GGINT=COLB*(GG12-GG13-GG23)
2127
2128         WTGG=FACTR*(GG1+GG2+GGINT)
2129
2130 C q-qbar-q-qbar
2131
2132         QP=HWH4J6(P13,P12,P14,P23,P34,P24,EP1,EP3,EP2,EP4,ORIENT)
2133
2134      &    +HWH4J6(P24,P12,P23,P14,P34,P13,EP2,EP4,EP1,EP3,ORIENT)
2135
2136      &    +HWH4J6(P13,P34,P23,P14,P12,P24,EP3,EP1,EP4,EP2,ORIENT)
2137
2138      &    +HWH4J6(P24,P34,P14,P23,P12,P13,EP4,EP2,EP3,EP1,ORIENT)
2139
2140         QQ=HWH4J6(P13,P23,P34,P12,P14,P24,EP3,EP1,EP2,EP4,ORIENT)
2141
2142      &    +HWH4J6(P24,P23,P12,P34,P14,P13,EP2,EP4,EP3,EP1,ORIENT)
2143
2144      &    +HWH4J6(P13,P14,P12,P34,P23,P24,EP1,EP3,EP4,EP2,ORIENT)
2145
2146      &    +HWH4J6(P24,P14,P34,P12,P23,P13,EP4,EP2,EP1,EP3,ORIENT)
2147
2148         QQINT=HWH4J7(P13,P12,P14,P23,P34,P24,EP1,EP3,EP2,EP4,ORIENT)
2149
2150      &       +HWH4J7(P24,P12,P23,P14,P34,P13,EP2,EP4,EP1,EP3,ORIENT)
2151
2152      &       +HWH4J7(P13,P23,P34,P12,P14,P24,EP3,EP1,EP2,EP4,ORIENT)
2153
2154      &       +HWH4J7(P24,P23,P12,P34,P14,P13,EP2,EP4,EP3,EP1,ORIENT)
2155
2156      &       +HWH4J7(P13,P14,P12,P34,P23,P24,EP1,EP3,EP4,EP2,ORIENT)
2157
2158      &       +HWH4J7(P24,P14,P34,P12,P23,P13,EP4,EP2,EP1,EP3,ORIENT)
2159
2160      &       +HWH4J7(P13,P34,P23,P14,P12,P24,EP3,EP1,EP4,EP2,ORIENT)
2161
2162      &       +HWH4J7(P24,P34,P14,P23,P12,P13,EP4,EP2,EP3,EP1,ORIENT)
2163
2164 C Add up weights
2165
2166         WTQP=FACTR*COLC*QP/TWO
2167
2168         QQ1  =COLC*QP
2169
2170         QQ2  =COLC*QQ
2171
2172         QQINT=COLB*QQINT
2173
2174         WTQQ=FACTR*(QQ1+QQ2+QQINT)/2
2175
2176       ENDIF
2177
2178 C
2179
2180       HCS=0.
2181
2182       DO 60 ID1=IDMN,IDMX
2183
2184       IF (INCLQG(ID1)) THEN
2185
2186 C Gluon channel
2187
2188         HCS=HCS+CLF(1,ID1)*WTGG
2189
2190         IF (GENEV.AND.HCS.GT.RCS) THEN
2191
2192 C Select colour flow
2193
2194           WTAB=GG1
2195
2196           WTBA=GG2
2197
2198           IF (IOP4JT(1).EQ.1) THEN
2199
2200             IF (GGINT.GE.ZERO) THEN
2201
2202               WTAB=WTAB+GGINT
2203
2204             ELSE
2205
2206               WTBA=MAX(WTBA,WTBA+GGINT)
2207
2208             ENDIF
2209
2210           ELSEIF (IOP4JT(1).EQ.2) THEN
2211
2212             IF (GGINT.GE.ZERO) THEN
2213
2214               WTBA=WTBA+GGINT
2215
2216             ELSE
2217
2218               WTAB=MAX(WTAB,WTAB+GGINT)
2219
2220             ENDIF
2221
2222           ELSEIF (IOP4JT(1).NE.0) THEN
2223
2224             CALL HWWARN('HWH4JT',101,*999)
2225
2226           ENDIF
2227
2228           WTOT=WTAB+WTBA
2229
2230           IF (WTAB.GT.HWR()*WTOT) THEN
2231
2232             CALL HWHQCP( 13, 13,3142,91,*99)
2233
2234           ELSE
2235
2236             CALL HWHQCP( 13, 13,4123,92,*99)
2237
2238           ENDIF
2239
2240         ENDIF
2241
2242       ENDIF
2243
2244 C Quark channels
2245
2246       DO 50 ID2=1,6
2247
2248 C Identical quark pairs
2249
2250       IF (ID1.EQ.ID2.AND.INCLQQ(ID1,ID1)) THEN
2251
2252         HCS=HCS+CLF(1,ID1)*WTQQ
2253
2254         IF (GENEV.AND.HCS.GT.RCS) THEN
2255
2256 C Select colour flow
2257
2258           WTAB=QQ1
2259
2260           WTBA=QQ2
2261
2262           IF (IOP4JT(2).EQ.1) THEN
2263
2264             IF (QQINT.GE.ZERO) THEN
2265
2266               WTAB=WTAB+QQINT
2267
2268             ELSE
2269
2270               WTBA=MAX(WTBA,WTBA+QQINT)
2271
2272             ENDIF
2273
2274           ELSEIF (IOP4JT(2).EQ.2) THEN
2275
2276             IF (QQINT.GE.ZERO) THEN
2277
2278               WTBA=WTBA+QQINT
2279
2280             ELSE
2281
2282               WTAB=MAX(WTAB,WTAB+QQINT)
2283
2284             ENDIF
2285
2286           ELSEIF (IOP4JT(2).NE.0) THEN
2287
2288             CALL HWWARN('HWH4JT',102,*999)
2289
2290           ENDIF
2291
2292           WTOT=WTAB+WTBA
2293
2294           IF (WTAB.GT.HWR()*WTOT) THEN
2295
2296             CALL HWHQCP(ID1,ID1+6,4123,93,*99)
2297
2298           ELSE
2299
2300             CALL HWHQCP(ID1,ID1+6,2143,94,*99)
2301
2302           ENDIF
2303
2304         ENDIF
2305
2306 C Unlike quark pairs
2307
2308       ELSEIF (INCLQQ(ID1,ID2)) THEN
2309
2310         HCS=HCS+(CLF(1,ID1)+CLF(1,ID2))*WTQP
2311
2312         IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID2,ID2+6,4123,95,*99)
2313
2314       ENDIF
2315
2316   50  CONTINUE
2317
2318   60  CONTINUE
2319
2320       EVWGT=HCS
2321
2322       RETURN
2323
2324 C Set up labels for selected final state
2325
2326   99  IDN(1)=ID1
2327
2328       IDN(2)=ID1+6
2329
2330       J=NHEP+1
2331
2332       IDHW(J)=200
2333
2334       IDHEP(J)=23
2335
2336       ISTHEP(J)=110
2337
2338       JMOHEP(1,J)=LM
2339
2340       JMOHEP(2,J)=LP
2341
2342       JDAHEP(1,J)=NHEP+2
2343
2344       JDAHEP(2,J)=NHEP+5
2345
2346       DO 100 I=1,4
2347
2348       J=NHEP+1+I
2349
2350       IDHW(J)=IDN(I)
2351
2352       IDHEP(J)=IDPDG(IDN(I))
2353
2354       ISTHEP(J)=IST(I)
2355
2356       JMOHEP(1,J)=NHEP+1
2357
2358   100 JDAHEP(1,J)=0
2359
2360 C And colour structure pointers
2361
2362       DO 110 I=1,4
2363
2364       J=ICO(I)
2365
2366       JMOHEP(2,NHEP+1+I)=NHEP+1+J
2367
2368   110 JDAHEP(2,NHEP+1+J)=NHEP+1+I
2369
2370       NHEP=NHEP+5
2371
2372   999 END
2373
2374 CDECK  ID>, HWH4J1.
2375
2376 *CMZ :-        -01/04/99  19.47.55  by  Mike Seymour
2377
2378 *-- Author :    Ian Knowles
2379
2380 C-----------------------------------------------------------------------
2381
2382       FUNCTION HWH4J1(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2383
2384 C-----------------------------------------------------------------------
2385
2386 C     Evaluate `ERT' functions A, B, C, D, E; S12=(p1+p2)^2 etc.
2387
2388 C-----------------------------------------------------------------------
2389
2390       IMPLICIT NONE
2391
2392       DOUBLE PRECISION HWH4J1,HWH4J2,HWH4J4,HWH4J5,HWH4J6,HWH4J7,
2393
2394      & S12,S13,S14,S23,S24,S34,S123,S124,S134,S234,S,EP1,EP2,EP3,EP4,
2395
2396      & SUM
2397
2398       LOGICAL ORIENT
2399
2400       S123=S12+S13+S23
2401
2402       S124=S12+S14+S24
2403
2404       S134=S13+S14+S34
2405
2406       S234=S23+S24+S34
2407
2408       S=S12+S13+S14+S23+S24+S34
2409
2410       HWH4J1=(S12*((S12+S14+S23+S34)**2+S13*(S12+S14-S24)+S24*(S12+S23))
2411
2412      &       +(S14*S23-S12*S34-S13*S24)*(S14+S23+S34)/2)
2413
2414      &       /(S13*S24*S134*S234)
2415
2416      &      +((S12+S24)*(S13+S34)-S14*S23)/(S13*S134**2)
2417
2418      &      +2*S23*(S-S13)/(S13*S134*S24) + S34/(2*S13*S24)
2419
2420       IF (ORIENT) THEN
2421
2422         HWH4J1=HWH4J1
2423
2424      &  +4*((EP1*EP1*((S-S13)*(S23+S24)-S24*S34)
2425
2426      &      -EP1*EP2*(S12*(S123+S124)+(S+S12)*(S14+S23)+2*S14*S23
2427
2428      &               +S24*S134+S234*(S13+2*S234))
2429
2430      &      +EP1*EP3*(S*(S24-S12)+S12*S13+(S14+2*S234-S34)*S24)
2431
2432      &      -EP1*EP4*(S12*S124+S23*(S+S12+S14))
2433
2434      &      +EP2*EP2*((S-S24)*(S13+S14)+2*(S13+S34)*S234-S13*S34)
2435
2436      &      -EP2*EP3*((S+S23)*(S12+S14)+(S12+2*(S23+S234))*S234)
2437
2438      &      +EP2*EP4*(S12*(S24-S)+S13*(S+S23-S34)+2*(S13+S34-S234)*S234)
2439
2440      &      +EP3*EP3*(S14+2*S234)*S24
2441
2442      &      +EP3*EP4*(-S234*(2*(S12+S23)+S134)+S12*S34-S13*S24-S14*S23)
2443
2444      &      +EP4*EP4*S13*S23)*S134
2445
2446      &      +EP2*(EP1+EP3+EP4)*2*S14*S24*S234)/(S*S13*S24*S134**2*S234)
2447
2448       ELSE
2449
2450         HWH4J1=2*HWH4J1/3
2451
2452       ENDIF
2453
2454       RETURN
2455
2456 C-----------------------------------------------------------------------
2457
2458       ENTRY HWH4J2(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2459
2460 C-----------------------------------------------------------------------
2461
2462       S123=S12+S13+S23
2463
2464       S124=S12+S14+S24
2465
2466       S134=S13+S14+S34
2467
2468       S234=S23+S24+S34
2469
2470       S=S12+S13+S14+S23+S24+S34
2471
2472       HWH4J2=(S12*S14*(S24+S34)+S24*(S12*(S14+S34)+S13*(S14-S24)))
2473
2474      &       /(S14*S23*S13*S134)
2475
2476      &      +S12*(S+S34)*S124/(S24*S234*S14*S134)
2477
2478      &      -(S13*(2*(S12+S24)+S23)+S14**2)/(S134*S13*S14)
2479
2480      &      +S12*S123*S124/(2*S13*S24*S14*S23)
2481
2482       IF (ORIENT) THEN
2483
2484         HWH4J2=HWH4J2
2485
2486      &  +4*((EP1*EP1*(S12*S134*S234-4*S23*S24*S34)
2487
2488      &      +EP1*EP2*(2*(2*S13*S234+S14*S123)*S24-S12*S134*(S+S12+S34))
2489
2490      &      +EP1*EP3*(S12*(4*S24*S34-S134*(S12+S14-S24))
2491
2492      &               -4*(S13*S24-S14*S23)*S24)
2493
2494      &      +EP1*EP4*(4*(S13+S14)*S23*S24-S12*S134*(S12+S13-S23))
2495
2496      &      +EP2*EP2*(S12*S134-4*S13*S24)*S134
2497
2498      &      +EP2*EP3*(4*S13*(S12+S23+S24)*S24-S12*S134*(S12-S14+S24))
2499
2500      &      -EP2*EP4*(4*(S12*(S14+S134)+S13*(S134-S234))*S24
2501
2502      &               +S12*(S12-S13+S23)*S134)
2503
2504      &      -EP3*EP3*4*S12*S14*S24
2505
2506      &      -EP3*EP4*2*S12*(2*S14*S24+S12*S134))*S234
2507
2508      &      +(EP1*(EP1*(S23+S24)+EP2*(S134-2*S))
2509
2510      &       -(EP1+EP2)*(EP3+EP4)*S12+EP2*EP2*(S13+S14))*2*S14*S24*S123)
2511
2512      &    /(2*S*S13*S14*S234*S23*S24*S134)
2513
2514       ELSE
2515
2516         HWH4J2=2*HWH4J2/3
2517
2518       ENDIF
2519
2520       RETURN
2521
2522 C-----------------------------------------------------------------------
2523
2524       ENTRY HWH4J4(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2525
2526 C-----------------------------------------------------------------------
2527
2528       S134=S13+S14+S34
2529
2530       S234=S23+S24+S34
2531
2532       S=S12+S13+S14+S23+S24+S34
2533
2534       HWH4J4=-(S12*(S34*(3*(S+S34)+S12)-S134*S234-2*(S13*S24+S14*S23))
2535
2536      &        +(S14*S23-S13*S24)*(S13-S14+S24-S23))/(2*S134*S234*S34**2)
2537
2538      &      -(S12*(S134**2/2+2*S13*S14+S34*(S13+S14-S34))
2539
2540      &       +S34*((S13+S14)*(S23+S24)+S14*S24+S13*S23)
2541
2542      &       +(S13*S24-S14*S23)*(S14-S13))/(S34*S134)**2
2543
2544       IF (ORIENT) THEN
2545
2546         HWH4J4=HWH4J4
2547
2548      &  +4*((-EP1*EP1*2*(S23+S24)*S34
2549
2550      &       -EP1*EP2*(S13*(S23+3*S24)+S14*(3*S23+S24)-(4*S12-S34)*S34)
2551
2552      &       +EP1*EP3*((2*S12-S24)*S34-(S13-S14)*S24)
2553
2554      &       +EP1*EP4*((2*S12-S23)*S34+(S13-S14)*S23)
2555
2556      &       -EP2*EP2*2*(S13+S14)*S34
2557
2558      &       +EP2*EP3*(2*S12*S34-S14*(S23-S24+S34))
2559
2560      &       +EP2*EP4*(2*S12*S34+S13*(S23-S24-S34))
2561
2562      &       +EP3*EP3*2*S14*S24
2563
2564      &       +EP3*EP4*2*(S12*S34-S13*S24-S14*S23)
2565
2566      &       +EP4*EP4*2*S13*S23)/(S*S134*S234*S34**2)
2567
2568      &      +(EP1*EP2*(S134*(S134+2*S34)+4*(S13*S14-S34**2))
2569
2570      &       +EP2*EP3*2*(2*S13*S34+S14*(S13-S14+S34))
2571
2572      &       +EP2*EP4*2*(2*S14*S34-S13*(S13-S14-S34)))
2573
2574      &  /(S*(S134*S34)**2))
2575
2576       ELSE
2577
2578         HWH4J4=2*HWH4J4/3
2579
2580       ENDIF
2581
2582       RETURN
2583
2584 C-----------------------------------------------------------------------
2585
2586       ENTRY HWH4J5(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2587
2588 C-----------------------------------------------------------------------
2589
2590       S123=S12+S13+S23
2591
2592       S124=S12+S14+S24
2593
2594       S134=S13+S14+S34
2595
2596       S234=S23+S24+S34
2597
2598       S=S12+S13+S14+S23+S24+S34
2599
2600       HWH4J5=(3*S12*S34**2-3*S13*S24*S34+3*S12*S24*S34+3*S14*S23*S34-
2601
2602      $     S13*S24**2-S12*S23*S34+6*S12*S14*S34+2*S12*S13*S34-
2603
2604      $     2*S12**2*S34+S14*S23*S24-3*S13*S23*S24-2*S13*S14*S24+
2605
2606      $     4*S12*S14*S24+2*S12*S13*S24+3*S14*S23**2+2*S14**2*S23+
2607
2608      $     2*S14**2*S12+2*S12**2*S14+6*S12*S14*S23-2*S12*S13**2-
2609
2610      $     2*S12**2*S13)/(2*S13*S134*S234*S34)+
2611
2612      $     (2*S12*S34**2-2*S13*S24*S34+S12*S24*S34+4*S13*S23*S34+
2613
2614      $     4*S12*S14*S34+2*S12*S13*S34+2*S12**2*S34-S13*S24**2+
2615
2616      $     3*S14*S23*S24+4*S13*S23*S24-2*S13*S14*S24+4*S12*S14*S24+
2617
2618      $     2*S12*S13*S24+2*S14*S23**2+4*S13*S23**2+2*S13*S14*S23+
2619
2620      $     2*S12*S14*S23+4*S12*S13*S23+2*S12*S14**2+4*S12**2*S13+
2621
2622      $     4*S12*S13*S14+2*S12**2*S14)/(2*S13*S134*S24*S34)-
2623
2624      $     (S12*S34**2-2*S14*S24*S34-2*S13*S24*S34-S14*S23*S34+
2625
2626      $     S13*S23*S34+S12*S14*S34+2*S12*S13*S34-2*S14**2*S24-
2627
2628      $     4*S13*S14*S24-4*S13**2*S24-S14**2*S23-S13**2*S23+
2629
2630      $     S12*S13*S14-S12*S13**2)/(S13*S34*S134**2)
2631
2632       IF (ORIENT) THEN
2633
2634         SUM=
2635
2636      &    +EP1*EP1*((S13-S14+S23-3*S24)*S34+(S134+S14+2*S34)*S234)
2637
2638      &            *S24*S134
2639
2640      &    +EP1*EP2*((2*(S12-S24)+S34)*S134-S14*(4*S12+S14+3*S23)
2641
2642      &             +S13*(S13+S23)+S24*S34 )*S24*S134
2643
2644      &    -EP1*EP2*(((2*S12*S134+S13*(2*(S12+S14+S23)-S24+S34)
2645
2646      &              +S14*(S14-S23)+(2*S14-S34)*S234)*S234)*S134
2647
2648      &             + 4*S13**2*S24*S234)
2649
2650      &    +EP1*EP3*(S12*(2*S13-S134)+S13*(S24+2*S234)+S14*(3*S24-S234)
2651
2652      &             +S34*(S234-3*S24))*S24*S134
2653
2654      &    +EP1*EP4*((S12*(S13-S14+3*S34)-S23*(S13+3*S14-S34))*S24
2655
2656      &             -(S12*(S13+S134+2*S34)+2*S13*S24
2657
2658      &              +(S13-2*S14)*S23)*S234)*S134
2659
2660      &    +EP2*EP2*(S13*((2*S13+S34)*S234+S24*(S134-2*S34))
2661
2662      &             +2*S14*S134*(S24+S234))*S134
2663
2664         SUM=SUM
2665
2666      &    -EP2*EP3*(((S12*(S13+2*S14-S34)+S14*(S+2*S23-S34))*S24
2667
2668      &              +(S12*(S13+S134)+(S13+S24+2*S234)*S14
2669
2670      &               +2*S13*(2*S23+S34))*S234)*S134
2671
2672      &             +4*S13**2*S24*S234)
2673
2674      &    +EP2*EP4*(((S12*(S13-2*S134)+S13*(S+2*S23-3*S34))*S24
2675
2676      &              -((S-3*S13+S23+2*S24)*S13+2*S12*S14
2677
2678      &                +2*S14*(S23+2*S24))*S234)*S134-4*S13**2*S24*S234)
2679
2680      &    +EP3*EP3*2*(S13*S234+S14*S24)*S24*S134
2681
2682      &    +EP3*EP4*(2*(S12*S34-S13*S24-S14*S23)*S24
2683
2684      &             -(S12*S134+2*S13*S23)*S234)*S134
2685
2686      &    +EP4*EP4*2*(S12*S234+S23*S24)*S13*S134
2687
2688         HWH4J5=HWH4J5+4*SUM/(S*S234*S134**2*S13*S34*S24)
2689
2690       ELSE
2691
2692         HWH4J5=2*HWH4J5/3
2693
2694       ENDIF
2695
2696       RETURN
2697
2698 C-----------------------------------------------------------------------
2699
2700       ENTRY HWH4J6(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2701
2702 C-----------------------------------------------------------------------
2703
2704       S123=S12+S13+S23
2705
2706       S124=S12+S14+S24
2707
2708       S134=S13+S14+S34
2709
2710       S234=S23+S24+S34
2711
2712       S=S12+S13+S14+S23+S24+S34
2713
2714       HWH4J6=(S23*(S123*S234-S*S23)+S12*(S123*S124-S*S12))/(S13*S123)**2
2715
2716      &     -(S12*S34*(S234-2*S23)+S14*S23*(S234-2*S34)
2717
2718      &      -S13*S24*(S234+S13))/(S13**2*S123*S134)
2719
2720       IF (ORIENT) THEN
2721
2722         HWH4J6=HWH4J6
2723
2724      &  +4*(-EP1*EP1*2*S23*S34
2725
2726      &      +EP1*EP2*((S12-S23)*S34-S13*(S24-S34))
2727
2728      &      +(EP1*EP3+EP2*EP4)*2*(S12*S34-S13*S24+S14*S23)
2729
2730      &      -EP1*EP4*(S13*S24-(3*(S13+S14)+S34)*S23)
2731
2732      &      -(EP1+EP2+EP3)*EP4*2
2733
2734      &       *(S12*(S13+S23)+(S12+S13)*S23)*S134/S123
2735
2736      &      +EP2*EP2*S13*(S14+S34)
2737
2738      &      +EP2*EP3*(S13*(S14-S24)-(S12-S23)*S14)
2739
2740      &      -EP3*EP3*2*S12*S14
2741
2742      &      -EP3*EP4*(S13*S24-(3*(S13+S34)+S14)*S12)
2743
2744      &      +EP4*EP4*(S12+S23)*S13)/(S*S134*S123*S13**2)
2745
2746       ELSE
2747
2748         HWH4J6=2*HWH4J6/3
2749
2750       ENDIF
2751
2752       RETURN
2753
2754 C-----------------------------------------------------------------------
2755
2756       ENTRY HWH4J7(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2757
2758 C-----------------------------------------------------------------------
2759
2760       S123=S12+S13+S23
2761
2762       S124=S12+S14+S24
2763
2764       S134=S13+S14+S34
2765
2766       S234=S23+S24+S34
2767
2768       S=S12+S13+S14+S23+S24+S34
2769
2770       HWH4J7=((S12*S34+S13*S24-S14*S23)*(S13+S14+S23+S24)-2*S12*S24*S34)
2771
2772      &      /(S13*S134*S23*S123)
2773
2774      &      -S12*(S12*S-S123*S124)/(S123**2*S13*S23)
2775
2776      &      -(S13+S14)*(S23+S24)*S34/(S13*S134*S23*S234)
2777
2778       IF (ORIENT) THEN
2779
2780         HWH4J7=HWH4J7
2781
2782      &  +4*(+2*(EP1+EP2)*(S23*EP1-S13*EP2)*S34*S134
2783
2784      &      -EP1*EP2*2*S34**2*S123
2785
2786      &      +EP1*EP3*(S123*(S23+S24)*S34+2*S134*(S13*S24-S14*S23))
2787
2788      &      +EP1*EP4*(S123*(S23+S24)*S34+2*S12**2*S134*S234/S123
2789
2790      &               +2*S134*(S24*(S13-S12)-S23*(S12+S14)))
2791
2792      &      +EP2*EP3*(2*(S12*S34+S13*S24-S14*S23)*S134
2793
2794      &               +S123*(S13+S14)*S34)
2795
2796      &      +EP2*EP4*(S123*(S13+S14)*S34+2*S12**2*S234*S134/S123
2797
2798      &               -2*S134*(S12*S234-S13*S24+S14*S23))
2799
2800      &      -EP3*EP3*S12*(2*S24*S134+S123*S34)
2801
2802      &      +EP3*EP4*2*S12*(S134*(S23-S24)-S34*S123+S12*S134*S234/S123)
2803
2804      &      +EP4*EP4*S12*(2*S23*S134-S123*S34))
2805
2806      &     /(S*S13*S23*S123*S134*S234)
2807
2808       ELSE
2809
2810         HWH4J7=2*HWH4J7/3
2811
2812       ENDIF
2813
2814       RETURN
2815
2816       END