HERWIG/src/hwhdis.f

   1
   2 CDECK  ID>, HWHDIS.
   3
   4 *CMZ :-        -26/04/91  14.55.44  by  Federico Carminati
   5
   6 *-- Author :    Giovanni Abbiendi & Luca Stanco
   7
   8 C----------------------------------------------------------------------
   9
  10       SUBROUTINE HWHDIS
  11
  12 C----------------------------------------------------------------------
  13
  14 C     DEEP INELASTIC LEPTON-HADRON SCATTERING: MEAN EVWGT = SIGMA IN NB
  15
  16 C----------------------------------------------------------------------
  17
  18       INCLUDE 'HERWIG61.INC'
  19
  20       DOUBLE PRECISION HWR,HWRUNI,HWUPCM,PRAN,PROB,SAMP,SIG,Q2,
  21
  22      & XBJ,Y,W,S,MLEP,MHAD,MLSCAT,LEP,YMIN,YMAX,XXMAX,Q2JAC,XXJAC,
  23
  24      & JACOBI,A1,A2,A3,B1,B2,PCM,PCMEP,PCMLW,PCMEQ,PCMLQ,COSPHI,PA,
  25
  26      & EQ,PZQ,SHAT,PROP,DLEFT,DRGHT,DUP,DWN,FACT,EFACT,OMY2,YPLUS,
  27
  28      & YMNUS,SIGMA,AF(7,12),SMA,Q2SUP,HWUAEM,DCHRG,DNEUT
  29
  30       INTEGER I,IQK,IQKIN,IQKOUT,IDSCAT,IHAD,ILEPT
  31
  32       LOGICAL CHARGD
  33
  34       EXTERNAL HWR,HWRUNI,HWUPCM
  35
  36       SAVE MLEP,MHAD,S,SMA,PCM,MLSCAT,A1,A2,A3,B1,B2,DLEFT,DRGHT,Q2,
  37
  38      & AF,XBJ,Y,YPLUS,YMNUS,OMY2,FACT,EFACT,SIGMA,IDSCAT,CHARGD,
  39
  40      & ILEPT,DCHRG,DNEUT,LEP
  41
  42       IQK=MOD(IPROC,10)
  43
  44       IHAD=2
  45
  46       IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD)
  47
  48       IF (FSTWGT.OR.IHAD.NE.2) THEN
  49
  50 C---INITIALISE PROCESS (MUST BE DONE EVERY TIME IF S VARIES)
  51
  52 C---LEPTON AND HADRON MASSES, INVARIANT MASS, MOMENTUM IN C.M. FRAME
  53
  54         MLEP=PHEP(5,1)
  55
  56         MHAD=PHEP(5,IHAD)
  57
  58         S=PHEP(5,3)**2
  59
  60         SMA=S-MLEP**2-MHAD**2
  61
  62         PCM=HWUPCM(SQRT(S),MLEP,MHAD)
  63
  64 C---LEP = 1 FOR LEPTONS, -1 FOR ANTILEPTONS
  65
  66         IF (IDHW(1).GE.121.AND.IDHW(1).LE.126) THEN
  67
  68           LEP=+ONE
  69
  70         ELSEIF (IDHW(1).GE.127.AND.IDHW(1).LE.132) THEN
  71
  72           LEP=-ONE
  73
  74         ELSE
  75
  76           CALL HWWARN('HWHDIS',500,*999)
  77
  78         ENDIF
  79
  80         DCHRG=FLOAT(MOD(IDHW(1)  ,2))
  81
  82         DNEUT=FLOAT(MOD(IDHW(1)+1,2))
  83
  84         ILEPT=MOD(IDHW(1)-121,6)+11
  85
  86 C---DLEFT,DRIGHT = 1,0 for leptons; = 0,1 for anti-leptons
  87
  88         DLEFT=MAX(+LEP,ZERO)
  89
  90         DRGHT=MAX(-LEP,ZERO)
  91
  92         CHARGD=MOD(IPROC,100)/10.EQ.1
  93
  94 C---Evaluate constant factor in cross section and
  95
  96 C   find and store scattered lepton identity
  97
  98         IF (CHARGD) THEN
  99
 100           IF ((EPOLN(3)-PPOLN(3)).EQ.ONE) THEN
 101
 102              WRITE(6,5)
 103
 104              CALL HWWARN('HWHDIS',501,*999)
 105
 106   5          FORMAT(1X,'WARNING: Cross-section is zero for the',
 107
 108      &                ' specified lepton helicity')
 109
 110           ENDIF
 111
 112           FACT=GEV2NB*(ONE-(EPOLN(3)-PPOLN(3)))*.25D0*PIFAC
 113
 114      &        /(SWEIN*RMASS(198)**2)**2
 115
 116           IDSCAT=IDHW(1)+NINT(DCHRG-DNEUT)
 117
 118         ELSE
 119
 120           FACT=GEV2NB*TWO*PIFAC
 121
 122           IDSCAT=IDHW(1)
 123
 124         ENDIF
 125
 126         MLSCAT=RMASS(IDSCAT)
 127
 128 C---PARAMETERS USED FOR THE WEIGHT GENERATION IN NEUTRAL CURRENT
 129
 130 C   PROCESSES. ASSUME D(SIGMA)/D(Q**2) GOES LIKE A1+A2/Q**2+A3/Q**4
 131
 132 C   AND D(SIGMA)/D(X) LIKE B1+B2/X
 133
 134         A1=0.5
 135
 136         A2=0.5
 137
 138         A3=1.
 139
 140         B1=0.1
 141
 142         B2=1.
 143
 144       ENDIF
 145
 146       IF (GENEV) THEN
 147
 148 C---GENERATE EVENT (KINEMATICAL VARIABLES AND STRUCTURE FUNCTION
 149
 150 C   ALREADY FOUND)
 151
 152         PRAN=SIGMA*HWR()
 153
 154         IF (CHARGD) THEN
 155
 156 C---CHARGED CURRENT PROCESS
 157
 158           IF (IQK.EQ.0) THEN
 159
 160 C---FIND FLAVOUR OF THE STRUCK QUARK (IF NOT SELECTED BY THE USER)
 161
 162             PROB=ZERO
 163
 164             DO 10 I=1,6
 165
 166             DUP=MOD(I+1,2)
 167
 168             DWN=MOD(I  ,2)
 169
 170             PROB=PROB+EFACT*
 171
 172      &          ((DCHRG*(DLEFT*DUP+DRGHT*DWN*OMY2)
 173
 174      &           +DNEUT*(DLEFT*DWN+DRGHT*DUP*OMY2))*DISF(I  ,1)
 175
 176      &          +(DCHRG*(DLEFT*DWN*OMY2+DRGHT*DUP)
 177
 178      &           +DNEUT*(DLEFT*DUP*OMY2+DRGHT*DWN))*DISF(I+6,1))
 179
 180             IF (PROB.GE.PRAN) GOTO 20
 181
 182    10       CONTINUE
 183
 184             I=6
 185
 186    20       IQK=I
 187
 188           ENDIF
 189
 190           DUP=MOD(IQK+1,2)
 191
 192           DWN=MOD(IQK  ,2)
 193
 194           IQKIN=IQK
 195
 196           IF ((LEP.EQ. 1.AND.MOD(IQK+IDHW(1),2).EQ.0)
 197
 198      &    .OR.(LEP.EQ.-1.AND.MOD(IQK+IDHW(1),2).EQ.1)) IQKIN=IQK+6
 199
 200 C---FIND FLAVOUR OF THE OUTGOING QUARK
 201
 202           PRAN=HWR()
 203
 204           PROB=ZERO
 205
 206           IF (DUP.EQ.ONE) THEN
 207
 208             DO 30 I=1,3
 209
 210             PROB=PROB+VCKM(IQK/2,I)
 211
 212             IF (PROB.GE.PRAN) GOTO 40
 213
 214    30       CONTINUE
 215
 216             I=3
 217
 218    40       IQKOUT=2*I-1
 219
 220             IF (IQKIN.GT.6) IQKOUT=IQKOUT+6
 221
 222           ELSE
 223
 224             DO 50 I=1,3
 225
 226             PROB=PROB+VCKM(I,(IQK+1)/2)
 227
 228             IF (PROB.GE.PRAN) GOTO 60
 229
 230    50       CONTINUE
 231
 232             I=3
 233
 234    60       IQKOUT=2*I
 235
 236             IF (IQKIN.GT.6) IQKOUT=IQKOUT+6
 237
 238           ENDIF
 239
 240         ELSE
 241
 242 C---NEUTRAL CURRENT PROCESS
 243
 244           IF (IQK.NE.0) THEN
 245
 246             IQKIN=IQK
 247
 248             PROB=EFACT*(AF(1,IQK)*YPLUS*DISF(IQK,1)+
 249
 250      &              LEP*AF(3,IQK)*YMNUS*DISF(IQK,1))
 251
 252             IF (PROB.LT.PRAN) IQKIN=IQK+6
 253
 254           ELSE
 255
 256 C---FIND FLAVOUR OF THE STRUCK QUARK (IF NOT SELECTED BY THE USER)
 257
 258             PROB=ZERO
 259
 260             SIG=ONE
 261
 262             DO 70 I=1,12
 263
 264             IF (I.GT.6) SIG=-ONE
 265
 266             PROB=PROB+EFACT*(AF(1,I)*YPLUS*DISF(I,1)+
 267
 268      &               LEP*SIG*AF(3,I)*YMNUS*DISF(I,1))
 269
 270             IF (PROB.GE.PRAN) GOTO 80
 271
 272    70       CONTINUE
 273
 274             I=12
 275
 276    80       IQKIN=I
 277
 278           ENDIF
 279
 280           IQKOUT=IQKIN
 281
 282         ENDIF
 283
 284         IDN(1)=IDHW(1)
 285
 286         IDN(2)=IQKIN
 287
 288         IDN(3)=IDSCAT
 289
 290         IDN(4)=IQKOUT
 291
 292         ICO(1)=1
 293
 294         ICO(2)=4
 295
 296         ICO(3)=3
 297
 298         ICO(4)=2
 299
 300         XX(1)=1.
 301
 302         XX(2)=XBJ
 303
 304 C---CHECK PHASE SPACE WITH THE SELECTED FLAVOUR. IF OUTSIDE THE
 305
 306 C   EVENT IS KILLED.
 307
 308         PA=XBJ*(PHEP(4,IHAD)+ABS(PHEP(3,IHAD)))
 309
 310         EQ=HALF*(PA+RMASS(IDN(2))**2/PA)
 311
 312         PZQ=-(PA-EQ)
 313
 314         SHAT=(PHEP(4,1)+EQ)**2-(PHEP(3,1)+PZQ)**2
 315
 316         PCMEQ=HWUPCM(SQRT(SHAT),MLEP,RMASS(IDN(2)))
 317
 318         PCMLQ=HWUPCM(SQRT(SHAT),MLSCAT,RMASS(IDN(4)))
 319
 320         IF (PCMLQ.LT.ZERO) THEN
 321
 322           CALL HWWARN('HWHDIS',101,*999)
 323
 324         ELSEIF (PCMLQ.EQ.ZERO) THEN
 325
 326           COSTH=ZERO
 327
 328         ELSE
 329
 330           COSTH=(TWO*SQRT(PCMEQ**2+MLEP**2)*SQRT(PCMLQ**2+MLSCAT**2)
 331
 332      &         -(Q2+MLEP**2+MLSCAT**2))/(TWO*PCMEQ*PCMLQ)
 333
 334         ENDIF
 335
 336         IF (ABS(COSTH).GT.ONE) CALL HWWARN('HWHDIS',102,*999)
 337
 338         IDCMF=15
 339
 340         CALL HWETWO
 341
 342       ELSE
 343
 344         EVWGT=ZERO
 345
 346         IF (CHARGD) THEN
 347
 348 C---CHOOSE X,Y (CC PROCESS)
 349
 350           YMIN=MAX(YBMIN,Q2MIN/SMA)
 351
 352           YMAX=MIN(YBMAX,ONE)
 353
 354           IF (YMIN.GT.YMAX) GOTO 999
 355
 356           Y=HWRUNI(0,YMIN,YMAX)
 357
 358           XXMIN=Q2MIN/S/Y
 359
 360           XXMAX=MIN(Q2MAX/SMA/Y,ONE)
 361
 362           IF (XXMIN.GT.XXMAX) GOTO 999
 363
 364           XBJ=HWRUNI(0,XXMIN,XXMAX)
 365
 366           Q2=XBJ*Y*(S-MLEP**2-MHAD**2)
 367
 368           JACOBI=(YMAX-YMIN)*(XXMAX-XXMIN)*(S-MLEP**2-MHAD**2)*XBJ
 369
 370         ELSE
 371
 372 C---CHOOSE X,Q**2 (NC PROCESS)
 373
 374           Q2SUP=MIN(Q2MAX,SMA*YBMAX)
 375
 376           IF (Q2MIN.GT.Q2SUP) GOTO 999
 377
 378           SAMP=(A1+A2+A3)*HWR()
 379
 380           IF (SAMP.LE.A1) THEN
 381
 382             Q2=HWRUNI(0,Q2MIN,Q2SUP)
 383
 384           ELSEIF (SAMP.LE.(A1+A2)) THEN
 385
 386             Q2=EXP(HWRUNI(0,LOG(Q2MIN),LOG(Q2SUP)))
 387
 388           ELSE
 389
 390             Q2=-ONE/HWRUNI(0,-ONE/Q2MIN,-ONE/Q2SUP)
 391
 392           ENDIF
 393
 394           Q2JAC=(A1+A2+A3)/
 395
 396      &      (A1/(Q2SUP-Q2MIN)
 397
 398      &      +A2/LOG(Q2SUP/Q2MIN)/Q2
 399
 400      &      +A3*Q2MIN*Q2SUP/(Q2SUP-Q2MIN)/Q2**2)
 401
 402           XXMIN=Q2/SMA/YBMAX
 403
 404           XXMAX=ONE
 405
 406           IF (YBMIN.GT.ZERO) XXMAX=MIN(Q2/SMA/YBMIN,ONE)
 407
 408           IF (XXMIN.GT.XXMAX) GOTO 999
 409
 410           SAMP=(B1+B2)*HWR()
 411
 412           IF (SAMP.LE.B1) THEN
 413
 414             XBJ=HWRUNI(0,XXMIN,XXMAX)
 415
 416           ELSE
 417
 418             XBJ=EXP(HWRUNI(0,LOG(XXMIN),LOG(XXMAX)))
 419
 420           ENDIF
 421
 422           XXJAC=(B1+B2)/(B1/(XXMAX-XXMIN)+B2/LOG(XXMAX/XXMIN)/XBJ)
 423
 424           Y=Q2/(S-MLEP**2-MHAD**2)/XBJ
 425
 426           JACOBI=Q2JAC*XXJAC
 427
 428         ENDIF
 429
 430 C---CHECK IF THE GENERATED POINT IS INSIDE PHASE SPACE. IF NOT
 431
 432 C   RETURN WITH WEIGHT EQUAL TO ZERO.
 433
 434         W=SQRT(MHAD**2+Q2*(ONE-XBJ)/XBJ)
 435
 436         IF (W.LT.WHMIN) RETURN
 437
 438         PCMEP=PCM
 439
 440         PCMLW=HWUPCM(SQRT(S),MLSCAT,W)
 441
 442         IF (PCMLW.LT.ZERO) THEN
 443
 444           EVWGT=ZERO
 445
 446           RETURN
 447
 448         ELSEIF (PCMLW.EQ.ZERO) THEN
 449
 450           COSPHI=ZERO
 451
 452         ELSE
 453
 454           COSPHI=
 455
 456      &    (TWO*SQRT(PCMEP**2+MLEP**2)*SQRT(PCMLW**2+MLSCAT**2)
 457
 458      &    -(Q2+MLEP**2+MLSCAT**2))/(TWO*PCMEP*PCMLW)
 459
 460         ENDIF
 461
 462         IF (ABS(COSPHI).GT.ONE) THEN
 463
 464           EVWGT=ZERO
 465
 466           RETURN
 467
 468         ENDIF
 469
 470 C---SET SCALE EQUAL Q. EVALUATE STRUCTURE FUNCTIONS.
 471
 472         EMSCA=SQRT(Q2)
 473
 474         CALL HWSFUN(XBJ,EMSCA,IDHW(IHAD),NSTRU,DISF,2)
 475
 476 C---SWITCH OFF ANY FLAVOURS THAT ARE BELOW THRESHOLD
 477
 478         DO 90 I=1,12
 479
 480  90     IF (W.LT.2*RMASS(I)) DISF(I,1)=0
 481
 482 C---EVALUATE DIFFERENTIAL CROSS SECTION
 483
 484         IF (CHARGD) THEN
 485
 486           PROP=RMASS(198)**2/(Q2+RMASS(198)**2)
 487
 488           EFACT=FACT*(HWUAEM(-Q2)*PROP)**2/XBJ
 489
 490           OMY2=(ONE-Y)**2
 491
 492           SIGMA=ZERO
 493
 494           DO 100 I=1,6
 495
 496           DUP=MOD(I+1,2)
 497
 498           DWN=MOD(I  ,2)
 499
 500           IF (IQK.NE.0.AND.IQK.NE.I) GOTO 100
 501
 502           SIGMA=SIGMA+EFACT*
 503
 504      &        ((DCHRG*(DLEFT*DUP+DRGHT*DWN*OMY2)
 505
 506      &         +DNEUT*(DLEFT*DWN+DRGHT*DUP*OMY2))*DISF(I  ,1)
 507
 508      &        +(DCHRG*(DLEFT*DWN*OMY2+DRGHT*DUP)
 509
 510      &         +DNEUT*(DLEFT*DUP*OMY2+DRGHT*DWN))*DISF(I+6,1))
 511
 512   100     CONTINUE
 513
 514         ELSE
 515
 516           EFACT=FACT/XBJ*(HWUAEM(-Q2)/Q2)**2
 517
 518           YPLUS=ONE+(ONE-Y)**2
 519
 520           YMNUS=ONE-(ONE-Y)**2
 521
 522           DO 110 I=1,6
 523
 524           CALL HWUCFF(ILEPT,I,-Q2,AF(1,I))
 525
 526           AF(1,I+6)=AF(1,I)
 527
 528           AF(3,I+6)=AF(3,I)
 529
 530   110     CONTINUE
 531
 532           SIGMA=ZERO
 533
 534           DO 200 I=1,6
 535
 536           IF (IQK.NE.0.AND.IQK.NE.I) GOTO 200
 537
 538           SIGMA=SIGMA+EFACT*(AF(1,I)*YPLUS*(DISF(I,1)+DISF(I+6,1))+
 539
 540      &                   LEP*AF(3,I)*YMNUS*(DISF(I,1)-DISF(I+6,1)))
 541
 542   200     CONTINUE
 543
 544         ENDIF
 545
 546 C---FIND WEIGHT: DIFFERENTIAL CROSS SECTION TIME THE JACOBIAN FACTOR
 547
 548         EVWGT=SIGMA*JACOBI
 549
 550         IF (EVWGT.LT.ZERO) EVWGT=ZERO
 551
 552       ENDIF
 553
 554   999 END
 555
 556 CDECK  ID>, HWHDYP.
 557
 558 *CMZ :-        -18/05/99  12.41.07  by  Mike Seymour
 559
 560 *-- Author :    Bryan Webber, Ian Knowles and Mike Seymour
 561
 562 C-----------------------------------------------------------------------
 563
 564       SUBROUTINE HWHDYP
 565
 566 C-----------------------------------------------------------------------
 567
 568 C     Drell-Yan Production of fermion pairs via photon, Z0 & (if ZPRIME)
 569
 570 C     Z' exchange. Lepton universality is assumed for photon and Z, and
 571
 572 C     for Z' if no lepton flavour is specified.
 573
 574 C     MEAN EVWGT = SIGMA IN NB
 575
 576 C-----------------------------------------------------------------------
 577
 578       INCLUDE 'HERWIG61.INC'
 579
 580       DOUBLE PRECISION HWR,HWRUNI,HWUAEM,EPS,C1,C2,C3,EMSQZ,EMGMZ,
 581
 582      & EMSQZP,EMGMZP,CQF(7,6,16),QPOW,RPOW,A01,A1,A02,A2,A03,A3,CRAN,
 583
 584      & EMJ1,EMJ2,EMJ3,EMJAC,FACT,QSQ,HCS,FACTR,RCS,EXTRA,PMAX,PTHETA
 585
 586       INTEGER IMODE,JQMN,JQMX,JQ,JLMN,JLMX,JL,IQ,I,IADD(2,2),ID1,ID2,
 587
 588      & ID3,ID4,JF
 589
 590       EXTERNAL HWR,HWRUNI,HWUAEM
 591
 592       SAVE HCS,JQMN,JQMX,JLMN,JLMX,C1,C2,C3,QPOW,RPOW,EMSQZ,EMGMZ,
 593
 594      & A1,A01,A2,A02,A3,A03,EMSQZP,EMGMZP,FACT,CQF
 595
 596       PARAMETER (EPS=1.D-9)
 597
 598       DATA IADD/0,6,6,0/
 599
 600       IF (GENEV) THEN
 601
 602         RCS=HCS*HWR()
 603
 604       ELSE
 605
 606         IF (FSTWGT) THEN
 607
 608 C Set limits for which particles to include
 609
 610           JLMN=1
 611
 612           JLMX=0
 613
 614           JQMN=1
 615
 616           JQMX=0
 617
 618           IMODE=MOD(IPROC,100)
 619
 620           IF (IMODE.EQ.0) THEN
 621
 622             JQMN=1
 623
 624             JQMX=6
 625
 626           ELSEIF (IMODE.LE.10) THEN
 627
 628             JQMN=IMODE
 629
 630             JQMX=IMODE
 631
 632           ELSEIF (IMODE.EQ.50) THEN
 633
 634             JLMN=11
 635
 636             JLMX=16
 637
 638           ELSEIF (IMODE.GE.50.AND.IMODE.LE.60) THEN
 639
 640             JLMN=IMODE-40
 641
 642             JLMX=IMODE-40
 643
 644           ELSEIF (IMODE.EQ.99) THEN
 645
 646             JQMN=1
 647
 648             JQMX=6
 649
 650             JLMN=11
 651
 652             JLMX=16
 653
 654           ELSE
 655
 656             CALL HWWARN('HWHDYP',500,*999)
 657
 658           ENDIF
 659
 660 C Set up parameters for importance sampling:
 661
 662 C sum of power law and two Breit-Wigners (relative weights C1,C2,C3)
 663
 664           C1=ONE
 665
 666           C2=ONE
 667
 668           C3=ZERO
 669
 670           IF (ZPRIME) C3=ONE
 671
 672           IF (EMPOW.EQ.ONE) CALL HWWARN('HWHDYP',501,*999)
 673
 674           IF (C2.EQ.ZERO) CALL HWWARN('HWHDYP',502,*999)
 675
 676           IF (C3.EQ.ZERO.AND.ZPRIME) CALL HWWARN('HWHDYP',503,*999)
 677
 678           QPOW=-EMPOW+1
 679
 680           RPOW=1/QPOW
 681
 682           EMSQZ=RMASS(200)**2
 683
 684           EMGMZ=RMASS(200)*GAMZ
 685
 686           A01=EMMIN**QPOW
 687
 688           A1=(EMMAX**QPOW-A01)/C1
 689
 690           A02=ATAN((EMMIN**2-EMSQZ)/EMGMZ)
 691
 692           A2=(ATAN((EMMAX**2-EMSQZ)/EMGMZ)-A02)/C2
 693
 694           IF (C3.GT.ZERO) THEN
 695
 696             EMSQZP=RMASS(202)**2
 697
 698             EMGMZP=RMASS(202)*GAMZP
 699
 700             A03=ATAN((EMMIN**2-EMSQZP)/EMGMZP)
 701
 702             A3=(ATAN((EMMAX**2-EMSQZP)/EMGMZP)-A03)/C3
 703
 704           ENDIF
 705
 706         ENDIF
 707
 708         EVWGT=0.
 709
 710 C Select a mass for the produced pair
 711
 712         CRAN=(C1+C2+C3)*HWR()
 713
 714         IF (CRAN.LT.C1) THEN
 715
 716 C Use power law
 717
 718           EMSCA=(A01+A1*CRAN)**RPOW
 719
 720           QSQ=EMSCA**2
 721
 722         ELSEIF (CRAN.LT.C1+C2) THEN
 723
 724 C Use Z Breit-Wigner
 725
 726           CRAN=CRAN-C1
 727
 728           QSQ=EMSQZ+EMGMZ*TAN(A02+A2*CRAN)
 729
 730           EMSCA=SQRT(QSQ)
 731
 732         ELSE
 733
 734 C Use Z' Breit-Wigner
 735
 736           CRAN=CRAN-C1-C2
 737
 738           QSQ=EMSQZP+EMGMZP*TAN(A03+A3*CRAN)
 739
 740           EMSCA=SQRT(QSQ)
 741
 742         ENDIF
 743
 744         EMJ1=EMSCA**EMPOW/(1-EMPOW)*A1
 745
 746         EMJ2=((QSQ-EMSQZ)**2+EMGMZ**2)/(2*EMSCA*EMGMZ)*A2
 747
 748         IF (C3.GT.ZERO) THEN
 749
 750           EMJ3=((QSQ-EMSQZP)**2+EMGMZP**2)/(2*EMSCA*EMGMZP)*A3
 751
 752           EMJAC=(C1+C2+C3)/(1/EMJ1+1/EMJ2+1/EMJ3)
 753
 754         ELSE
 755
 756           EMJAC=(C1+C2)/(1/EMJ1+1/EMJ2)
 757
 758         ENDIF
 759
 760 C Select initial momentum fractions
 761
 762         XXMIN=QSQ/PHEP(5,3)**2
 763
 764         XLMIN=LOG(XXMIN)
 765
 766         CALL HWSGEN(.TRUE.)
 767
 768         FACT=-GEV2NB*HWUAEM(QSQ)**2*PIFAC*8*EMJAC*XLMIN
 769
 770      $       /(3*NCOLO*EMSCA**3)
 771
 772 C Store cross-section coefficients
 773
 774         DO 50 IQ=1,6
 775
 776         DO 30 JQ=JQMN,JQMX
 777
 778         IF (EMSCA.GT.2.*RMASS(JQ)) THEN
 779
 780           CALL HWUCFF(IQ,JQ,QSQ,CQF(1,IQ,JQ))
 781
 782         ELSE
 783
 784           CALL HWVZRO(7,CQF(1,IQ,JQ))
 785
 786         ENDIF
 787
 788   30    CONTINUE
 789
 790         DO 40 JL=JLMN,JLMX
 791
 792         IF (EMSCA.GT.2.*RMASS(JL)) THEN
 793
 794           CALL HWUCFF(IQ,JL,QSQ,CQF(1,IQ,JL))
 795
 796         ELSE
 797
 798           CALL HWVZRO(7,CQF(1,IQ,JL))
 799
 800         ENDIF
 801
 802   40    CONTINUE
 803
 804   50    CONTINUE
 805
 806       ENDIF
 807
 808 C
 809
 810       HCS=0.
 811
 812       DO 90 I=1,2
 813
 814 C I=1 quark first, I=2 anti-quark first
 815
 816       DO 80 IQ=1,6
 817
 818       ID1=IQ+IADD(1,I)
 819
 820       ID2=IQ+IADD(2,I)
 821
 822       IF (DISF(ID1,1).LT.EPS.OR.DISF(ID2,2).LT.EPS) GOTO 80
 823
 824       FACTR=FACT*DISF(ID1,1)*DISF(ID2,2)
 825
 826 C Quark final states
 827
 828       DO 60 JQ=JQMN,JQMX
 829
 830       ID3=JQ
 831
 832       ID4=JQ+6
 833
 834       IF (IQ.EQ.JQ) THEN
 835
 836         HCS=HCS+FACTR*(CQF(1,IQ,JQ)*FLOAT(NCOLO)+3*HALF*QFCH(IQ)**4)
 837
 838         IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID3,ID4,2143,50,*99)
 839
 840       ELSE
 841
 842         HCS=HCS+FACTR*CQF(1,IQ,JQ)*FLOAT(NCOLO)
 843
 844         IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID3,ID4,2143,50,*99)
 845
 846       ENDIF
 847
 848   60  CONTINUE
 849
 850 C Lepton final states
 851
 852       DO 70 JL=JLMN,JLMX
 853
 854       ID3=110+JL
 855
 856       ID4=ID3+6
 857
 858       HCS=HCS+FACTR*CQF(1,IQ,JL)
 859
 860       IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID3,ID4,2134,50,*99)
 861
 862   70  CONTINUE
 863
 864   80  CONTINUE
 865
 866   90  CONTINUE
 867
 868       EVWGT=HCS
 869
 870       RETURN
 871
 872 C Generate event
 873
 874   99  IDN(1)=ID1
 875
 876       IDN(2)=ID2
 877
 878       IDCMF=200
 879
 880       IF (ID3.LE.6) THEN
 881
 882         JF=JQ
 883
 884       ELSE
 885
 886         JF=JL
 887
 888       ENDIF
 889
 890 C Select polar angle from distribution:
 891
 892 C CQF(1,IQ,JF)*(ONE+COSTH**2)+CQF(3,IQ,JF)*COSTH+EXTRA*(ONE+COSTH)
 893
 894       IF (ID1.EQ.ID3.OR.ID2.EQ.ID3) THEN
 895
 896         EXTRA=TWO*QFCH(ID3)**4/NCOLO
 897
 898       ELSE
 899
 900         EXTRA=0
 901
 902       ENDIF
 903
 904       PMAX=2.*(CQF(1,IQ,JF)+EXTRA)+ABS(CQF(3,IQ,JF))
 905
 906   100 COSTH=HWRUNI(0,-ONE,ONE)
 907
 908       PTHETA=CQF(1,IQ,JF)*(ONE+COSTH**2)+CQF(3,IQ,JF)*COSTH
 909
 910      &      +EXTRA*(ONE+COSTH)
 911
 912       IF (PTHETA.LT.PMAX*HWR()) GOTO 100
 913
 914       IF (ID1.GT.ID2) COSTH=-COSTH
 915
 916       IDCMF=200
 917
 918       CALL HWETWO
 919
 920   999 END
 921
 922 CDECK  ID>, HWHEGG.
 923
 924 *CMZ :-        -19/03/92  10.13.56  by  Mike Seymour
 925
 926 *-- Author :    Mike Seymour
 927
 928 C-----------------------------------------------------------------------
 929
 930       SUBROUTINE HWHEGG
 931
 932 C----------------------------------------------------------------------
 933
 934 C     HARD PROCESS: EE --> EEGAMGAM --> EEFFBAR/WW
 935
 936 C     MEAN EVENT WEIGHT = CROSS-SECTION IN NB
 937
 938 C     AFTER CUTS ON PT AND MASS OF CENTRE-OF-MASS SYSTEM
 939
 940 C     AND COS(THETA) IN CENTRE-OF-MASS SYSTEM
 941
 942 C     AND TIMES BRANCHING FRACTION IF WW
 943
 944 C-----------------------------------------------------------------------
 945
 946       INCLUDE 'HERWIG61.INC'
 947
 948       DOUBLE PRECISION HWR,HWULDO,EMSQ,BETA,S,T,U,TMIN,TMAX,TRAT,
 949
 950      & DSDT,PROB,X,Z(2),ZMIN,ZMAX,PCMIN,PCMAX,PCFAC,PLOGMI,PLOGMA,PTCMF,
 951
 952      & Q,PC,BLOG,EMCMIN,EMCMAX,EMLMIN,EMLMAX,WGT(6),RWGT,CV,CA,BR,QT(2),
 953
 954      & QX(2),QY(2),PX,PY,ROOTS,DOT,A,B,C,SHAT,PCF(2),PCM(2),PCMAC,ZZ(2),
 955
 956      & COLFAC
 957
 958       INTEGER I,IGAM,ID,IDL,ID1,ID2,IHEP,JHEP,NADD,NTRY,NQ,JGAM
 959
 960       LOGICAL HWRLOG
 961
 962       EXTERNAL HWR,HWULDO,HWRLOG
 963
 964       SAVE S,BETA,X,ID,NQ,WGT,EMLMIN,EMLMAX,PCFAC,PLOGMA,PLOGMI,SHAT,
 965
 966      &  PCF,PCM,Z,PCMAC,NADD
 967
 968       IF (IERROR.NE.0) RETURN
 969
 970 C---INITIALIZE LOCAL COPIES OF EMMIN,EMMAX
 971
 972       IF (FSTWGT) THEN
 973
 974         EMLMIN=EMMIN
 975
 976         EMLMAX=EMMAX
 977
 978       ENDIF
 979
 980       IF (.NOT.GENEV) THEN
 981
 982 C---CHOOSE Z1,Z2 AND CALCULATE SUB-PROCESS CROSS-SECTION
 983
 984         EVWGT=0
 985
 986 C-----FIND FINAL STATE PARTICLES
 987
 988         IHPRO=MOD(IPROC,100)
 989
 990         IF (IHPRO.EQ.0) THEN
 991
 992           ID=1
 993
 994           NQ=6
 995
 996           COLFAC=FLOAT(NCOLO)
 997
 998           NADD=6
 999
1000         ELSEIF (IHPRO.LE.6) THEN
1001
1002           ID=IHPRO
1003
1004           NQ=1
1005
1006           COLFAC=FLOAT(NCOLO)
1007
1008           NADD=6
1009
1010           Q=QFCH(ID)
1011
1012         ELSEIF (IHPRO.LE.9) THEN
1013
1014           ID=119+2*(IHPRO-6)
1015
1016           NQ=1
1017
1018           COLFAC=1.
1019
1020           NADD=6
1021
1022           Q=QFCH(ID-110)
1023
1024         ELSEIF (IHPRO.LE.10) THEN
1025
1026           ID=198
1027
1028           NQ=1
1029
1030           NADD=1
1031
1032         ELSE
1033
1034           CALL HWWARN('HWHEGG',200,*999)
1035
1036         ENDIF
1037
1038 C-----SPLIT ELECTRONS TO PHOTONS
1039
1040         NHEP=3
1041
1042         GAMWT=1
1043
1044         S=2*HWULDO(PHEP(1,1),PHEP(1,2))
1045
1046         ROOTS=SQRT(S)
1047
1048         EMCMIN=MAX(EMLMIN,MAX(2*RMASS(ID),PTMIN))
1049
1050         EMCMAX=MIN(EMLMAX,ROOTS)
1051
1052         IF (EMCMIN.GT.EMCMAX) RETURN
1053
1054         ZMIN=EMCMIN**2/S
1055
1056         ZMAX=1-PHEP(5,1)/PHEP(4,1)
1057
1058         IF (ZMIN.GT.ZMAX) RETURN
1059
1060         CALL HWEGAM(1,ZMIN,ZMAX,.TRUE.)
1061
1062         Z(1)=PHEP(4,NHEP-1)/PHEP(4,1)
1063
1064         ZMIN=EMCMIN**2/(Z(1)*S)
1065
1066         ZMAX=MIN(EMCMAX**2/(Z(1)*S), ONE-PHEP(5,2)/PHEP(4,2))
1067
1068         IF (ZMIN.GT.ZMAX) RETURN
1069
1070         CALL HWEGAM(2,ZMIN,ZMAX,.TRUE.)
1071
1072         Z(2)=PHEP(4,NHEP-1)/PHEP(4,2)
1073
1074         EMSCA=PHEP(5,3)
1075
1076         SHAT=EMSCA**2
1077
1078 C-----REMOVE LOG TERMS FROM WEIGHT, CALCULATE NEW ONES FROM PT LIMITS
1079
1080         GAMWT=GAMWT/(0.5*LOG((1-Z(1))*S/(Z(1)*PHEP(5,1)**2))
1081
1082      &              *0.5*LOG((1-Z(2))*Z(1)*S/(Z(2)*PHEP(5,2)**2)))
1083
1084         PCF(1)=Z(1)*PHEP(5,1)
1085
1086         PCF(2)=Z(2)*PHEP(5,2)
1087
1088         PCFAC=SQRT(PCF(1)*PCF(2))
1089
1090         PCM(1)=(1-Z(1))*PHEP(4,1)
1091
1092         PCM(2)=(1-Z(2))*PHEP(4,2)
1093
1094         PCMAC=SQRT(PCM(1)*PCM(2))
1095
1096         PCMIN=MAX(PTMIN,MAX(PCF(1),PCF(2)))
1097
1098         PCMAX=MIN( MIN(PTMAX,PHEP(5,3)) , MIN(PCM(1),PCM(2)) )
1099
1100         IF (PCMIN.GT.PCMAX) RETURN
1101
1102         PLOGMI=(LOG(PCMIN/PCFAC))**2
1103
1104         PLOGMA=(LOG(PCMAX/PCFAC))**2
1105
1106         GAMWT=GAMWT*(PLOGMA-PLOGMI)
1107
1108 C-----CALCULATE CROSS-SECTION
1109
1110         DO 10 IDL=1,NQ
1111
1112           WGT(IDL)=EVWGT
1113
1114           IF (IHPRO.EQ.0) THEN
1115
1116             ID=IDL
1117
1118             Q=QFCH(ID)
1119
1120           ENDIF
1121
1122           EMSQ=RMASS(ID)**2
1123
1124           X=4*EMSQ/SHAT
1125
1126           IF (X.GT.ONE) GOTO 10
1127
1128           BETA=SQRT(1-X)
1129
1130           BLOG=LOG((1+BETA*CTMAX)/(1-BETA*CTMAX))/BETA
1131
1132           IF (IHPRO.LE.9) THEN
1133
1134             EVWGT=EVWGT+GEV2NB*4*PIFAC*COLFAC*Q**4*ALPHEM**2*BETA
1135
1136      &           /SHAT * GAMWT * ( (1+X-0.5*X**2)*BLOG
1137
1138      &                     - CTMAX*(1+X**2/(CTMAX**2*(X-1)+1)) )
1139
1140             WGT(IDL)=EVWGT
1141
1142           ELSE
1143
1144             CALL HWDBOZ(198,ID1,ID2,CV,CA,BR,1)
1145
1146             CALL HWDBOZ(199,ID1,ID2,CV,CA,BR,1)
1147
1148             EVWGT=EVWGT + GEV2NB*6*PIFAC*ALPHEM**2*BETA/SHAT*BR
1149
1150      &        * GAMWT * (-(  X-0.5*X**2)*BLOG
1151
1152      &                     + CTMAX*(1+(X**2+16/3.)/(CTMAX**2*(X-1)+1)) )
1153
1154           ENDIF
1155
1156  10     CONTINUE
1157
1158 C-----GAMWT MUST BE RESET TO ONE, SINCE IT IS REAPPLIED LATER!
1159
1160         GAMWT=ONE
1161
1162       ELSE
1163
1164 C---GENERATE EVENT
1165
1166 C-----CHOOSE PT OF THE CMF
1167
1168         PTCMF=PCFAC*EXP(SQRT(HWR()*(PLOGMA-PLOGMI)+PLOGMI))
1169
1170 C-----CHOOSE WHICH PHOTON USUALLY HAS SMALLER PT
1171
1172         NTRY=0
1173
1174  20     IGAM=1
1175
1176         IF (LOG(PCM(1)/PCF(1)).LT.HWR()*2*LOG(PCMAC/PCFAC)) IGAM=2
1177
1178         JGAM=3-IGAM
1179
1180 C-----CHOOSE ITS PT
1181
1182  30     NTRY=NTRY+1
1183
1184         IF (NTRY.GT.NBTRY) CALL HWWARN('HWHEGG',100,*999)
1185
1186         QT(IGAM)=(PCM(IGAM)/PCF(IGAM))**HWR()
1187
1188         PROB=(QT(IGAM)**2/(QT(IGAM)**2+1))**2
1189
1190         QT(IGAM)=QT(IGAM)*PCF(IGAM)
1191
1192         IF (HWRLOG(1-PROB)) GOTO 30
1193
1194 C-----CHOOSE ITS DIRECTION
1195
1196         CALL HWRAZM(QT(IGAM),QX(IGAM),QY(IGAM))
1197
1198 C-----CALCULATE THE OTHER PHOTON'S PT
1199
1200         QX(JGAM)=PTCMF-QX(IGAM)
1201
1202         QY(JGAM)=     -QY(IGAM)
1203
1204         QT(JGAM)=SQRT(QX(JGAM)**2+QY(JGAM)**2)
1205
1206         IF (QT(JGAM).LT.PCF(JGAM).OR.QT(JGAM).GT.PCM(JGAM)) GOTO 20
1207
1208 C-----APPLY A RANDOM ROTATION AROUND THE BEAM AXIS
1209
1210         CALL HWRAZM(ONE,PX,PY)
1211
1212         IF (PX.EQ.ZERO) PX=1D-20
1213
1214         QX(1)=(QX(1)*PX   -QY(1)*PY)
1215
1216         QY(1)=(QY(1)      +QX(1)*PY)/PX
1217
1218         QX(2)=(QX(2)*PX   -QY(2)*PY)
1219
1220         QY(2)=(QY(2)      +QX(2)*PY)/PX
1221
1222 C-----RECONSTRUCT MOMENTA
1223
1224         IF (QT(IGAM).GT.QT(JGAM)) THEN
1225
1226           IGAM=3-IGAM
1227
1228           JGAM=3-JGAM
1229
1230         ENDIF
1231
1232         DOT=-Z(JGAM)*S+SHAT+2*(QX(1)*QX(2)+QY(1)*QY(2))
1233
1234 C-------SOLVE QUADRATIC IN Z(IGAM) TO FIND ELECTRON ENERGIES
1235
1236         A=S*(S*Z(JGAM)+QT(JGAM)**2)
1237
1238         B=S*DOT*(1+Z(JGAM))
1239
1240         C=DOT**2+S*QT(IGAM)**2*(1-Z(JGAM))**2-4*QT(IGAM)**2*QT(JGAM)**2
1241
1242         IF (B**2.LT.4*A*C) GOTO 20
1243
1244         ZZ(IGAM)=(-B+SQRT(B**2-4*A*C))/(2*A)
1245
1246         IF (ZZ(IGAM).LT.ZERO .OR. ZZ(IGAM).GT.ONE-Z(IGAM)) GOTO 20
1247
1248         ZZ(JGAM)=1-Z(JGAM)
1249
1250 C-------REJECT AGAINST PHOTON DISTRIBUTION FUNCTION
1251
1252         PROB=((1+ZZ(IGAM)**2)/(1-ZZ(IGAM)))/((1+(1-Z(IGAM))**2)/Z(IGAM))
1253
1254      &      *((1+ZZ(JGAM)**2)/(1-ZZ(JGAM)))/((1+(1-Z(JGAM))**2)/Z(JGAM))
1255
1256         IF (HWRLOG(1-PROB)) GOTO 20
1257
1258 C-------RECONSTRUCT ALL OTHER VARIABLES
1259
1260         DO 40 I=1,2
1261
1262           IGAM=2*I+3
1263
1264           PHEP(1,IGAM)=QX(I)
1265
1266           PHEP(2,IGAM)=QY(I)
1267
1268           PHEP(4,IGAM)=ZZ(I)*PHEP(4,I)
1269
1270           PHEP(5,IGAM)=RMASS(IDHW(IGAM))
1271
1272 C---------IF MOMENTUM CANNOT BE CONSERVED TRY AGAIN
1273
1274           IF (PHEP(4,IGAM)**2-PHEP(5,IGAM)**2-QT(I)**2 .LT. 0) GOTO 20
1275
1276           PHEP(3,IGAM)=SIGN(SQRT(PHEP(4,IGAM)**2-PHEP(5,IGAM)**2-
1277
1278      &      QT(I)**2),PHEP(3,IGAM))
1279
1280           CALL HWVDIF(4,PHEP(1,I),PHEP(1,IGAM),PHEP(1,IGAM-1))
1281
1282           CALL HWUMAS(PHEP(1,IGAM-1))
1283
1284  40     CONTINUE
1285
1286 C-----TIDY UP EVENT RECORD
1287
1288         NHEP=NHEP+1
1289
1290         IDHW(NHEP)=IDHW(3)
1291
1292         IDHEP(NHEP)=IDHEP(3)
1293
1294         ISTHEP(NHEP)=110
1295
1296         CALL HWVSUM(4,PHEP(1,4),PHEP(1,6),PHEP(1,NHEP))
1297
1298         CALL HWVSUM(4,PHEP(1,1),PHEP(1,2),PHEP(1,3))
1299
1300         CALL HWUMAS(PHEP(1,NHEP))
1301
1302         CALL HWUMAS(PHEP(1,3))
1303
1304         JMOHEP(1,NHEP)=4
1305
1306         JMOHEP(2,NHEP)=6
1307
1308         JMOHEP(1,3)=0
1309
1310         JMOHEP(2,3)=0
1311
1312 C-----CHOOSE FINAL STATE QUARK
1313
1314         IF (IHPRO.EQ.0) THEN
1315
1316           RWGT=HWR()*EVWGT
1317
1318           ID=1
1319
1320           DO 50 IDL=1,NQ
1321
1322             IF (RWGT.GT.WGT(IDL)) ID=IDL+1
1323
1324  50       CONTINUE
1325
1326           EMSQ=RMASS(ID)**2
1327
1328           X=4*EMSQ/SHAT
1329
1330           BETA=SQRT(1-X)
1331
1332         ENDIF
1333
1334 C-----CHOOSE T (WHERE T = MANDELSTAM_T - EMSQ)
1335
1336         TMIN=-SHAT/2
1337
1338         TMAX=-SHAT/2*(1-BETA*CTMAX)
1339
1340         TRAT=TMAX/TMIN
1341
1342         NTRY=0
1343
1344         IF (IHPRO.LE.9) THEN
1345
1346 C-------FOR FFBAR, CHOOSE T ACCORDING TO -SHAT/T
1347
1348  60       NTRY=NTRY+1
1349
1350           IF (NTRY.GT.NBTRY) CALL HWWARN('HWHEGG',101,*999)
1351
1352           T=TRAT**HWR()*TMIN
1353
1354           U=-T-SHAT
1355
1356 C-------REWEIGHT TO CORRECT DISTRIBUTION
1357
1358           DSDT=(T*U-2*EMSQ*(T+2*EMSQ))/T**2
1359
1360      &        +( 2*EMSQ*(SHAT-4*EMSQ))/(T*U)
1361
1362      &        +(T*U-2*EMSQ*(U+2*EMSQ))/U**2
1363
1364           PROB=-DSDT*T/SHAT / (1 + 2*X - 2*X**2)
1365
1366           IF (HWRLOG(1-PROB)) GOTO 60
1367
1368         ELSE
1369
1370 C-------FOR WW, CHOOSE T ACCORDING TO (SHAT/T)**2
1371
1372  70       NTRY=NTRY+1
1373
1374           IF (NTRY.GT.NBTRY) CALL HWWARN('HWHEGG',102,*999)
1375
1376           T=TMAX/(1-(1-TRAT)*HWR())
1377
1378           U=-T-SHAT
1379
1380 C-------REWEIGHT TO CORRECT DISTRIBUTION
1381
1382           DSDT=( 3*(T*U)**2 - SHAT*T*U*(4*SHAT+6*EMSQ)
1383
1384      &      + SHAT**2*(2*SHAT**2+6*EMSQ**2) ) / (T*U)**2
1385
1386           PROB=DSDT*(T/SHAT)**2 / (4.75 - 1.5*X + 1.5*X**2)
1387
1388           IF (HWRLOG(1-PROB)) GOTO 70
1389
1390         ENDIF
1391
1392 C-----SYMMETRIZE IN T,U
1393
1394         IF (HWRLOG(HALF)) T=U
1395
1396 C-----FILL EVENT RECORD
1397
1398         COSTH=(1+2*T/SHAT)/BETA
1399
1400         PC=0.5*BETA*PHEP(5,NHEP)
1401
1402         PHEP(5,NHEP+1)=RMASS(ID)
1403
1404         PHEP(5,NHEP+2)=RMASS(ID)
1405
1406         CALL HWDTWO(PHEP(1,NHEP),PHEP(1,NHEP+1),PHEP(1,NHEP+2),
1407
1408      &              PC,COSTH,.TRUE.)
1409
1410         DO 80 I=1,2
1411
1412           IHEP=NHEP+I
1413
1414           JHEP=NHEP+3-I
1415
1416           ISTHEP(IHEP)=190
1417
1418           IF (IHPRO.LE.6) ISTHEP(IHEP)=112+I
1419
1420           IDHW(IHEP)=ID+NADD*(I-1)
1421
1422           IDHEP(IHEP)=IDPDG(IDHW(IHEP))
1423
1424           JDAHEP(I,NHEP)=IHEP
1425
1426           JMOHEP(1,IHEP)=NHEP
1427
1428           JMOHEP(2,IHEP)=JHEP
1429
1430           JDAHEP(2,IHEP)=JHEP
1431
1432           IF (IHPRO.EQ.10) THEN
1433
1434             RHOHEP(1,IHEP)=0.3333
1435
1436             RHOHEP(2,IHEP)=0.3333
1437
1438             RHOHEP(3,IHEP)=0.3333
1439
1440           ENDIF
1441
1442  80     CONTINUE
1443
1444         NHEP=NHEP+2
1445
1446       ENDIF
1447
1448  999  END
1449
1450 CDECK  ID>, HWHEGW.
1451
1452 *CMZ :-        -26/04/91  10.18.56  by  Bryan Webber
1453
1454 *-- Author :    Mike Seymour
1455
1456 C-----------------------------------------------------------------------
1457
1458       SUBROUTINE HWHEGW
1459
1460 C----------------------------------------------------------------------
1461
1462 C     W + GAMMA --> FF'BAR :  MEAN EVWGT = CROSS SECTION IN NANOBARN
1463
1464 C     BASED ON BOSON GLUON FUSION OF ABBIENDI AND STANCO
1465
1466 C-----------------------------------------------------------------------
1467
1468       INCLUDE 'HERWIG61.INC'
1469
1470       DOUBLE PRECISION HWR,GMASS,EV(3),RV,LEP,Y,Q2,SHAT,Z,PHI,AJACOB,
1471
1472      & DSIGMA,ME,MP,ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT
1473
1474       INTEGER LEPFIN,ID1,ID2,I,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,IPROO
1475
1476       LOGICAL CHARGD,INCLUD(18),INSIDE(18)
1477
1478       EXTERNAL HWR
1479
1480       SAVE LEPFIN,ID1,ID2
1481
1482       COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF,
1483
1484      & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,
1485
1486      & IPROO,CHARGD,INCLUD,INSIDE
1487
1488       IQK=MOD(IPROC,10)
1489
1490       CHARGD=.TRUE.
1491
1492       IF(GENEV) THEN
1493
1494 C
1495
1496         IDHW(4)=IDHW(1)
1497
1498         IDHW(5)=59
1499
1500         IDHW(6)=15
1501
1502         IDHW(7)=LEPFIN
1503
1504         IDHW(8)=ID1
1505
1506         IDHW(9)=ID2
1507
1508         DO 1 I=4,9
1509
1510     1   IDHEP(I)=IDPDG(IDHW(I))
1511
1512 C
1513
1514         IFLAVD=ID1
1515
1516         IFLAVU=ID2-6
1517
1518 C
1519
1520         ISTHEP(4)=111
1521
1522         ISTHEP(5)=112
1523
1524         ISTHEP(6)=110
1525
1526         ISTHEP(7)=113
1527
1528         ISTHEP(8)=114
1529
1530         ISTHEP(9)=114
1531
1532 C
1533
1534         JMOHEP(1,4)=6
1535
1536         JMOHEP(2,4)=7
1537
1538         JMOHEP(1,5)=6
1539
1540         JMOHEP(2,5)=5
1541
1542         JMOHEP(1,6)=4
1543
1544         JMOHEP(2,6)=5
1545
1546         JMOHEP(1,7)=6
1547
1548         JMOHEP(2,7)=4
1549
1550         JMOHEP(1,8)=6
1551
1552         JMOHEP(2,8)=9
1553
1554         JMOHEP(1,9)=6
1555
1556         JMOHEP(2,9)=8
1557
1558         JDAHEP(1,4)=0
1559
1560         JDAHEP(2,4)=7
1561
1562         JDAHEP(1,5)=0
1563
1564         JDAHEP(2,5)=5
1565
1566         JDAHEP(1,6)=7
1567
1568         JDAHEP(2,6)=9
1569
1570         JDAHEP(1,7)=0
1571
1572         JDAHEP(2,7)=4
1573
1574         JDAHEP(1,8)=0
1575
1576         JDAHEP(2,8)=9
1577
1578         JDAHEP(1,9)=0
1579
1580         JDAHEP(2,9)=8
1581
1582 C---COMPUTATION OF MOMENTA IN LABORATORY FRAME OF REFERENCE
1583
1584 C---Persuade HWHBKI that the gluon is actually a photon...
1585
1586         GMASS=RMASS(13)
1587
1588         RMASS(13)=0
1589
1590         CALL HWHBKI
1591
1592         RMASS(13)=GMASS
1593
1594 C---put the other outgoing lepton in as well
1595
1596         IDHW(10)=IDHW(2)
1597
1598         IDHEP(10)=IDPDG(IDHW(10))
1599
1600         ISTHEP(10)=1
1601
1602         JMOHEP(1,10)=2
1603
1604         JMOHEP(2,10)=0
1605
1606         JDAHEP(1,10)=0
1607
1608         JDAHEP(2,10)=0
1609
1610         JDAHEP(1,2)=5
1611
1612         JDAHEP(2,2)=10
1613
1614         CALL HWVDIF(4,PHEP(1,2),PHEP(1,5),PHEP(1,10))
1615
1616         CALL HWUMAS(PHEP(1,10))
1617
1618         NHEP=10
1619
1620 C
1621
1622 C---if antilepton was first, do charge conjugation
1623
1624         IF (LEP.EQ.ONE) THEN
1625
1626           DO 27 I=7,9
1627
1628             IF (IDHEP(I).NE.0 .AND. ABS(IDHEP(I)).LT.20) THEN
1629
1630               IDHW(I)=IDHW(I) + 6*SIGN(1,IDHEP(I))
1631
1632               IDHEP(I)=-IDHEP(I)
1633
1634             ENDIF
1635
1636  27       CONTINUE
1637
1638         ENDIF
1639
1640 C
1641
1642 C---half the time, do charge conjugation and parity flip
1643
1644         IF (HWR().GT.HALF) THEN
1645
1646           DO 2 I=4,10
1647
1648             IF (IDHEP(I).NE.0 .AND. ABS(IDHEP(I)).LT.20) THEN
1649
1650               IDHW(I)=IDHW(I) + 6*SIGN(1,IDHEP(I))
1651
1652               IDHEP(I)=-IDHEP(I)
1653
1654             ENDIF
1655
1656             PHEP(1,I)=-PHEP(1,I)
1657
1658             PHEP(2,I)=-PHEP(2,I)
1659
1660             PHEP(3,I)=-PHEP(3,I)
1661
1662  2        CONTINUE
1663
1664           JMOHEP(1,10)=3-JMOHEP(1,10)
1665
1666         ENDIF
1667
1668 C
1669
1670       ELSE
1671
1672 C
1673
1674         EVWGT=ZERO
1675
1676 C---LEP = 1 IF TRACK 1 IS A LEPTON, -1 FOR ANTILEPTON
1677
1678         LEP=0.
1679
1680         IF (IDHW(1).GE.121.AND.IDHW(1).LE.126) THEN
1681
1682           LEP=1.
1683
1684         ELSEIF (IDHW(1).GE.127.AND.IDHW(1).LE.132) THEN
1685
1686           LEP=-1.
1687
1688         ENDIF
1689
1690         IF (LEP.EQ.ZERO) CALL HWWARN('HWHEGW',500,*999)
1691
1692 C---program only works if beam and target are charge conjugates
1693
1694         IF (INT(LEP)*(IDHW(2)-IDHW(1)).NE.6)
1695
1696      &   CALL HWWARN('HWHEGW',501,*999)
1697
1698 C---program only works for equal energy beams colliding
1699
1700         IF (PHEP(3,3).NE.ZERO) CALL HWWARN('HWHEGW',503,*999)
1701
1702 C
1703
1704 C---FINAL STATE IS ALWAYS SET UP AS IF PARTICLE IS BEFORE ANTI-PARTICLE
1705
1706 C   AND THEN INVERTED IF NECESSARY
1707
1708         LEPFIN = MIN(IDHW(1),IDHW(2))+1
1709
1710         IF (IQK.LE.2) THEN
1711
1712           IFLAVU=2
1713
1714           IFLAVD=1
1715
1716           ID1  = 1
1717
1718           ID2  = 8
1719
1720         ELSEIF (IQK.LE.4) THEN
1721
1722           IFLAVU=4
1723
1724           IFLAVD=3
1725
1726           ID1  = 3
1727
1728           ID2  =10
1729
1730         ELSEIF (IQK.LE.6) THEN
1731
1732           IFLAVU=6
1733
1734           IFLAVD=5
1735
1736           ID1  = 5
1737
1738           ID2  =12
1739
1740         ELSEIF (IQK.EQ.7) THEN
1741
1742           IFLAVU=122
1743
1744           IFLAVD=121
1745
1746           ID1  = 121
1747
1748           ID2  = 128
1749
1750 C---INTERFERENCE TERMS IN EE -> EE NUE NUEB  NEGLECTED: SIGMA UNRELIABLE
1751
1752           IF (FSTWGT) CALL HWWARN('HWHEGW',1,*999)
1753
1754         ELSEIF (IQK.EQ.8) THEN
1755
1756           IFLAVU=124
1757
1758           IFLAVD=123
1759
1760           ID1  = 123
1761
1762           ID2  = 130
1763
1764         ELSEIF (IQK.EQ.9) THEN
1765
1766           IFLAVU=126
1767
1768           IFLAVD=125
1769
1770           ID1  = 125
1771
1772           ID2  = 132
1773
1774         ELSE
1775
1776           CALL HWWARN('HWHEGW',504,*999)
1777
1778         ENDIF
1779
1780         IF (IQK.GT.0) THEN
1781
1782           IF (IQK.LE.6) IQK=0
1783
1784           CALL HWHBRN(*999)
1785
1786           CALL HWHEGX
1787
1788           EVWGT = 2 * DSIGMA * AJACOB
1789
1790           IF (EVWGT.LT.ZERO) EVWGT=ZERO
1791
1792         ELSE
1793
1794 C---SUM OVER QUARK FLAVOURS
1795
1796           CALL HWHBRN(*999)
1797
1798           DO 3 I=1,3
1799
1800             IF (SHAT.GT.(RMASS(IFLAVD)+RMASS(IFLAVU))**2) THEN
1801
1802               CALL HWHEGX
1803
1804               EV(I) = 2 * DSIGMA * AJACOB
1805
1806               IF (EV(I).LT.ZERO) EV(I)=ZERO
1807
1808             ELSE
1809
1810               EV(I)=ZERO
1811
1812             ENDIF
1813
1814             EVWGT=EVWGT+EV(I)
1815
1816             EV(I)=EVWGT
1817
1818             IFLAVU=IFLAVU+2
1819
1820             IFLAVD=IFLAVD+2
1821
1822  3        CONTINUE
1823
1824 C---CHOOSE QUARK FLAVOUR
1825
1826           RV=EV(3)*HWR()
1827
1828           IF (RV.LT.EV(1)) THEN
1829
1830             ID1 = 1
1831
1832             ID2 = 8
1833
1834           ELSEIF (RV.LT.EV(2)) THEN
1835
1836             ID1 = 3
1837
1838             ID2 =10
1839
1840           ELSE
1841
1842             ID1 = 5
1843
1844             ID2 =12
1845
1846           ENDIF
1847
1848         ENDIF
1849
1850       ENDIF
1851
1852   999 END
1853
1854 CDECK  ID>, HWHEGX.
1855
1856 *CMZ :-        -17/07/92  16.42.56  by  Mike Seymour
1857
1858 *-- Author :    Mike Seymour
1859
1860 C-----------------------------------------------------------------------
1861
1862       SUBROUTINE HWHEGX
1863
1864 C-----------------------------------------------------------------------
1865
1866 C     COMPUTES DIFFERENTIAL CROSS SECTION DSIGMA IN (Y,Q2,ETA,Z,PHI)
1867
1868 C-----------------------------------------------------------------------
1869
1870       INCLUDE 'HERWIG61.INC'
1871
1872       DOUBLE PRECISION TMAX,TMIN,A1,A2,B1,B2,I0,I1,I2,I3,I4,I5,MUSQ,
1873
1874      & MDSQ,ETA,Q1,COSTHE,S,G,T,U,C1,C2,D1,D2,F1,F2,COSBET,WPROP,D(4,4),
1875
1876      & C(4,4),QU,QD,QE,QW,PHOTON,EMWSQ,EMSSQ,CFAC,LEP,Y,Q2,SHAT,Z,PHI,
1877
1878      & AJACOB,DSIGMA,ME,MP,ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,
1879
1880      & RSHAT
1881
1882       INTEGER IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,IPROO,I,J
1883
1884       LOGICAL CHARGD,INCLUD(18),INSIDE(18)
1885
1886       COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF,
1887
1888      & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,
1889
1890      & IPROO,CHARGD,INCLUD,INSIDE
1891
1892 C---INPUT VARIABLES
1893
1894       IF (IERROR.NE.0) RETURN
1895
1896       DSIGMA=0
1897
1898       IF (IFLAVU.LE.12) THEN
1899
1900         QU=QFCH(MOD(IFLAVU-1,6)+1)
1901
1902         QD=QFCH(MOD(IFLAVD-1,6)+1)
1903
1904         CFAC=CAFAC
1905
1906       ELSE
1907
1908         QU=QFCH(MOD(IFLAVU-1,6)+11)
1909
1910         QD=QFCH(MOD(IFLAVD-1,6)+11)
1911
1912         CFAC=1
1913
1914       ENDIF
1915
1916       QE=QFCH(11)
1917
1918       QW=+1
1919
1920       EMWSQ=RMASS(198)**2
1921
1922       EMSCA=PHEP(5,3)
1923
1924       EMSSQ=EMSCA**2
1925
1926       MUSQ=RMASS(IFLAVU)**2
1927
1928       MDSQ=RMASS(IFLAVD)**2
1929
1930       ETA=(SHAT+Q2)/EMSSQ/Y
1931
1932       IF (ETA.GT.ONE) RETURN
1933
1934 C---CALCULATE KINEMATIC TERMS
1935
1936       G=0.5*(ETA*EMSSQ*Y-Q2) -0.5*(MUSQ+MDSQ)
1937
1938       S=0.5*ETA*EMSSQ
1939
1940       T=0.5*ETA*EMSSQ*(1-Y)
1941
1942       U=0.5*Q2
1943
1944       C1=0.5*ETA*EMSSQ*Y*Z
1945
1946       C2=0.5*ETA*EMSSQ*Y*(1-Z)
1947
1948       COSBET=(-ETA*EMSSQ*Y+Q2*(2-Y))/(Y*(ETA*EMSSQ-Q2))
1949
1950       IF (SHAT.LE.(RMASS(IFLAVU)+RMASS(IFLAVD))**2) RETURN
1951
1952       Q1=SQRT((SHAT**2+MUSQ**2+MDSQ**2
1953
1954      &  -2*SHAT*MUSQ-2*SHAT*MDSQ-2*MUSQ*MDSQ)/SHAT**2)
1955
1956       COSTHE=(1+(MDSQ-MUSQ)/SHAT-2*Z)/Q1
1957
1958       IF (ABS(COSTHE).GE.ONE .OR. ABS(COSBET).GE.ONE) RETURN
1959
1960       D1=0.25*(ETA*EMSSQ-Q2)*(1+(MDSQ-MUSQ)/SHAT-Q1*
1961
1962      &     (COSTHE*COSBET+SQRT((1-COSTHE**2)*(1-COSBET**2))*COS(PHI)))
1963
1964       D2=S-U-D1
1965
1966       F1=D1+C1-G            -MDSQ
1967
1968       F2=U+T-F1
1969
1970 C---CALCULATE TRACE TERMS
1971
1972       CALL HWVZRO(16,D)
1973
1974       CALL HWVZRO(16,C)
1975
1976       D(1,1)=2*F1*C2*S
1977
1978       D(2,2)=2*C1*D2*T
1979
1980       D(3,3)=-D1*(2*F2*G-D2*(F1+2*U))
1981
1982      &       -D2*F1*(F2+U-D2+F1)
1983
1984      &       +2*F1*F2*U
1985
1986      &       -G*(-2*D1*(F1+F2+U)-F1*(D2+2*U)+2*D2*(U-F2)+2*U*(F2-U+G))
1987
1988       D(4,4)=2*F1*C2*S
1989
1990       D(1,2)=(D1+U-F2)*(D1*F2-F1*D2)-G*(D1*(F2+U)+U*(U-F2-G)+F1*D2)
1991
1992       D(1,3)=D1*F2*(-2*F1+U-F2+D1)
1993
1994      &      +F1*(F2*(D2-2*U)+F1*D2)
1995
1996      &      +G*(-D1*(2*F1+F2+U)-F1*(D2+2*U)+U*(F2-U+G))
1997
1998       D(1,4)=-2*F1*(D1+U)*(F2+G)
1999
2000       D(2,3)=D1*(D2*(F1+2*(U-F2))+F2*(F2-U-D1))
2001
2002      &      +F1*D2**2
2003
2004      &      +G*(D1*(F2+U)+D2*(F1-2*(U-F2))+U*(U-F2-G))
2005
2006       D(2,4)=-D1*F2*(U-F2+D1)
2007
2008      &       -F1*D2*(U-D1-G-F2)
2009
2010      &       -G*(U*(F2-U+G)-D1*(F2+U))
2011
2012       D(3,4)=D1*(F1*(D2+2*F2)+F2*(F2-U-D1))
2013
2014      &      +F1*(2*F2*U-D2*(U+F1))
2015
2016      &      +G*(D1*(2*F1+F2+U)+U*(2*F1-F2+U-G))
2017
2018 C---REGULATE PROPAGATORS
2019
2020       TMAX=EMSSQ-2*G
2021
2022       TMIN=PHEP(5,2)**2
2023
2024       A1=2*C1+MDSQ*(G+U)/G
2025
2026       A2=2*C2+MUSQ*(G+U)/G
2027
2028       B1=(2*U+MUSQ)/(2*G+2*U)
2029
2030       B2=(2*U+MDSQ)/(2*G+2*U)
2031
2032       I0=LOG(TMAX/TMIN)
2033
2034       I1=1/A1*(I0-LOG((A1+B1*TMAX)/(A1+B1*TMIN)))
2035
2036       I2=1/A2*(I0-LOG((A2+B2*TMAX)/(A2+B2*TMIN)))
2037
2038       I3=(B1*I1-B2*I2)/(B1*A2-B2*A1)
2039
2040       I4=1/A1*(I1+1/(A1+B1*TMAX)-1/(A1+B1*TMIN))
2041
2042       I5=1/A2*(I2+1/(A2+B2*TMAX)-1/(A2+B2*TMIN))
2043
2044       WPROP=1/((2*G-EMWSQ)**2+GAMW**2*EMWSQ)
2045
2046 C---CALCULATE COEFFICIENTS
2047
2048       C(1,1)=    QU**2/(2*U+EMWSQ)**2                       *I5
2049
2050       C(2,2)=    QD**2/(2*U+EMWSQ)**2                       *I4
2051
2052       C(3,3)=    QW**2/(2*U+EMWSQ)**2    *WPROP             *I0
2053
2054       C(4,4)=    QE**2/(2*S)**2          *WPROP             *I0
2055
2056       C(1,2)=  2*QU*QD/(2*U+EMWSQ)**2                       *I3
2057
2058       C(1,3)=  2*QW*QU/(2*U+EMWSQ)**2    *WPROP*(2*G-EMWSQ) *I2
2059
2060       C(1,4)=  2*QU*QE/(2*S*(2*U+EMWSQ)) *WPROP*(2*G-EMWSQ) *I2
2061
2062       C(2,3)=  2*QW*QD/(2*U+EMWSQ)**2    *WPROP*(2*G-EMWSQ) *I1
2063
2064       C(2,4)=  2*QD*QE/(2*S*(2*U+EMWSQ)) *WPROP*(2*G-EMWSQ) *I1
2065
2066       C(3,4)=  2*QW*QE/(2*S*(2*U+EMWSQ)) *WPROP             *I0
2067
2068 C---CALCULATE PHOTON STRUCTURE FUNCTION
2069
2070       PHOTON=ALPHEM * (1+(1-ETA)**2) / (2*PIFAC*ETA)
2071
2072 C---SUM ALL TENSOR CONTRIBUTIONS
2073
2074       DO 10 I=1,4
2075
2076       DO 10 J=1,4
2077
2078  10     DSIGMA=DSIGMA + C(I,J)*D(I,J)
2079
2080 C---CALCULATE TOTAL SUMMED AND AVERAGED MATRIX ELEMENT SQUARED
2081
2082       DSIGMA = DSIGMA * 2*CFAC*(4*PIFAC*ALPHEM)**3/SWEIN**2
2083
2084 C---CALCULATE DIFFERENTIAL CROSS-SECTION
2085
2086       DSIGMA = DSIGMA * GEV2NB*PHOTON/(512*PIFAC**4*ETA*EMSSQ)
2087
2088  999  END
2089
2090 CDECK  ID>, HWHEPA.
2091
2092 *CMZ :-        -26/04/91  14.55.44  by  Federico Carminati
2093
2094 *-- Author :    Bryan Webber and Ian Knowles
2095
2096 C-----------------------------------------------------------------------
2097
2098       SUBROUTINE HWHEPA
2099
2100 C-----------------------------------------------------------------------
2101
2102 C     (Initially polarised) e+e- --> ffbar (f=quark, mu or tau)
2103
2104 C     If IPROC=107: --> gg, distributed as sum of light quarks.
2105
2106 C     If fermion flavour specified mass effects fully included.
2107
2108 C     EVWGT=sig(e+e- --> ffbar) in nb
2109
2110 C-----------------------------------------------------------------------
2111
2112       INCLUDE 'HERWIG61.INC'
2113
2114       DOUBLE PRECISION HWR,HWRUNI,HWUPCM,HWUAEM,Q2NOW,Q2LST,FACTR,
2115
2116      & VF2,VF,CLF(7),PRAN,PQWT,PMAX,PTHETA,SINTH2,CPHI,SPHI,C2PHI,S2PHI,
2117
2118      & PPHI,SINTH,PCM,PP(5),EWGT
2119
2120       INTEGER ID1,ID2,IDF,IQ,IQ1,I
2121
2122       EXTERNAL HWR,HWRUNI,HWUPCM,HWUAEM
2123
2124       SAVE Q2LST,FACTR,ID1,ID2,VF2,VF,CLF,EWGT
2125
2126       DATA Q2LST/0./
2127
2128       IF (GENEV) THEN
2129
2130         IF (ID2.EQ.0) THEN
2131
2132 C Choose quark flavour
2133
2134           PRAN=TQWT*HWR()
2135
2136           PQWT=0.
2137
2138           DO 10 IQ=1,MAXFL
2139
2140           PQWT=PQWT+CLQ(1,IQ)
2141
2142           IF (PQWT.GT.PRAN) GOTO 11
2143
2144    10     CONTINUE
2145
2146           IQ=MAXFL
2147
2148    11     IQ1=MAPQ(IQ)
2149
2150           DO 20 I=1,7
2151
2152    20     CLF(I)=CLQ(I,IQ)
2153
2154         ELSE
2155
2156           IQ1=ID1
2157
2158         ENDIF
2159
2160 C Label particles, assign outgoing particle masses
2161
2162         IDHW(NHEP+1)=200
2163
2164         IDHEP(NHEP+1)=23
2165
2166         ISTHEP(NHEP+1)=110
2167
2168         IF (ID1.EQ.7) THEN
2169
2170           IDHW(NHEP+2)=13
2171
2172           IDHW(NHEP+3)=13
2173
2174           IDHEP(NHEP+2)=21
2175
2176           IDHEP(NHEP+3)=21
2177
2178           PHEP(5,NHEP+2)=RMASS(13)
2179
2180           PHEP(5,NHEP+3)=RMASS(13)
2181
2182         ELSE
2183
2184           IDHW(NHEP+2)=IQ1
2185
2186           IDHW(NHEP+3)=IQ1+6
2187
2188           IDHEP(NHEP+2)=IDPDG(IQ1)
2189
2190           IDHEP(NHEP+3)=-IDHEP(NHEP+2)
2191
2192           PHEP(5,NHEP+2)=RMASS(IQ1)
2193
2194           PHEP(5,NHEP+3)=RMASS(IQ1)
2195
2196         ENDIF
2197
2198         ISTHEP(NHEP+2)=113
2199
2200         ISTHEP(NHEP+3)=114
2201
2202         JMOHEP(1,NHEP+1)=1
2203
2204         IF (JDAHEP(1,1).NE.0) JMOHEP(1,NHEP+1)=JDAHEP(1,1)
2205
2206         JMOHEP(2,NHEP+1)=2
2207
2208         IF (JDAHEP(1,2).NE.0) JMOHEP(1,NHEP+1)=JDAHEP(1,2)
2209
2210         JMOHEP(1,NHEP+2)=NHEP+1
2211
2212         JMOHEP(2,NHEP+2)=NHEP+3
2213
2214         JMOHEP(1,NHEP+3)=NHEP+1
2215
2216         JMOHEP(2,NHEP+3)=NHEP+2
2217
2218         JDAHEP(1,NHEP+1)=NHEP+2
2219
2220         JDAHEP(2,NHEP+1)=NHEP+3
2221
2222         JDAHEP(1,NHEP+2)=0
2223
2224         JDAHEP(2,NHEP+2)=NHEP+3
2225
2226         JDAHEP(1,NHEP+3)=0
2227
2228         JDAHEP(2,NHEP+3)=NHEP+2
2229
2230 C Generate polar and azimuthal angular distributions:
2231
2232 C  CLF(1)*(1+(VF*COSTH)**2)+CLF(2)*(1-VF**2)+CLF(3)*2.*VF*COSTH
2233
2234 C +(VF*SINTH)**2*(CLF(4)*COS(2*PHI-PHI1-PHI2)
2235
2236 C                +CLF(6)*SIN(2*PHI-PHI1-PHI2))
2237
2238         PMAX=CLF(1)*(1.+VF2)+CLF(2)*(1.-VF2)+ABS(CLF(3))*2.*VF
2239
2240   30    COSTH=HWRUNI(0,-ONE, ONE)
2241
2242         PTHETA=CLF(1)*(1.+VF2*COSTH**2)+CLF(2)*(1.-VF2)
2243
2244      &        +CLF(3)*2.*VF*COSTH
2245
2246         IF (PTHETA.LT.PMAX*HWR()) GOTO 30
2247
2248         IF (IDHW(1).GT.IDHW(2)) COSTH=-COSTH
2249
2250         SINTH2=1.-COSTH**2
2251
2252         IF (TPOL) THEN
2253
2254           PMAX=PTHETA+VF2*SINTH2*SQRT(CLF(4)**2+CLF(6)**2)
2255
2256   40      CALL HWRAZM(ONE,CPHI,SPHI)
2257
2258           C2PHI=2.*CPHI**2-1.
2259
2260           S2PHI=2.*CPHI*SPHI
2261
2262           PPHI=PTHETA+(CLF(4)*(C2PHI*COSS+S2PHI*SINS)
2263
2264      &                +CLF(6)*(S2PHI*COSS-C2PHI*SINS))*VF2*SINTH2
2265
2266           IF (PPHI.LT.PMAX*HWR()) GOTO 40
2267
2268         ELSE
2269
2270           CALL HWRAZM(ONE,CPHI,SPHI)
2271
2272         ENDIF
2273
2274 C Construct final state 4-mommenta
2275
2276         CALL HWVEQU(5,PHEP(1,3),PHEP(1,NHEP+1))
2277
2278         PCM=HWUPCM(PHEP(5,NHEP+1),PHEP(5,NHEP+2),PHEP(5,NHEP+3))
2279
2280 C PP is momentum of track NHEP+2 in CoM (track NHEP+1) frame
2281
2282         SINTH=SQRT(SINTH2)
2283
2284         PP(5)=PHEP(5,NHEP+2)
2285
2286         PP(1)=PCM*SINTH*CPHI
2287
2288         PP(2)=PCM*SINTH*SPHI
2289
2290         PP(3)=PCM*COSTH
2291
2292         PP(4)=SQRT(PCM**2+PP(5)**2)
2293
2294         CALL HWULOB(PHEP(1,NHEP+1),PP(1),PHEP(1,NHEP+2))
2295
2296         CALL HWVDIF(4,PHEP(1,NHEP+1),PHEP(1,NHEP+2),PHEP(1,NHEP+3))
2297
2298 C Set production vertices
2299
2300         CALL HWVZRO(4,VHEP(1,NHEP+2))
2301
2302         CALL HWVEQU(4,VHEP(1,NHEP+2),VHEP(1,NHEP+3))
2303
2304         NHEP=NHEP+3
2305
2306       ELSE
2307
2308         EMSCA=PHEP(5,3)
2309
2310         Q2NOW=EMSCA**2
2311
2312         IF (Q2NOW.NE.Q2LST) THEN
2313
2314 C Calculate coefficients for cross-section
2315
2316           EMSCA=PHEP(5,3)
2317
2318           Q2LST=Q2NOW
2319
2320           FACTR=PIFAC*GEV2NB*HWUAEM(Q2NOW)**2/Q2NOW
2321
2322           ID1=MOD(IPROC,10)
2323
2324           ID2=MOD(ID1,7)
2325
2326           IF (ID2.EQ.0) THEN
2327
2328             CALL HWUEEC(1)
2329
2330             VF2=1.
2331
2332             VF=1.
2333
2334             EWGT=FACTR*FLOAT(NCOLO)*TQWT*4./3.
2335
2336           ELSE
2337
2338             IF (IPROC.LT.150) THEN
2339
2340               IDF=ID1
2341
2342               FACTR=FACTR*FLOAT(NCOLO)
2343
2344             ELSE
2345
2346               ID1=2*ID1+119
2347
2348               IDF=ID1-110
2349
2350             ENDIF
2351
2352             IF (EMSCA.LE.2.*RMASS(ID1)) then
2353
2354               EWGT=0.
2355
2356             ELSE
2357
2358               CALL HWUCFF(11,IDF,Q2NOW,CLF(1))
2359
2360               VF2=1.-4.*RMASS(ID1)**2/Q2NOW
2361
2362               VF=SQRT(VF2)
2363
2364               EWGT=FACTR*VF*(CLF(1)*(1.+VF2/3.)+CLF(2)*(1.-VF2))
2365
2366             ENDIF
2367
2368           ENDIF
2369
2370         ENDIF
2371
2372         EVWGT=EWGT
2373
2374       ENDIF
2375
2376   999 END
2377
2378 CDECK  ID>, HWHEPG.
2379
2380 *CMZ :-        -02/05/91  10.57.27  by  Federico Carminati
2381
2382 *-- Author :    Bryan Webber and Ian Knowles
2383
2384 C-----------------------------------------------------------------------
2385
2386       SUBROUTINE HWHEPG
2387
2388 C-----------------------------------------------------------------------
2389
2390 C     (Initially polarised) e-e+ --> qqbar g with parton thrust < THMAX,
2391
2392 C     equivalent to: maximum parton energy < THMAX*EMSCA/2; or a JADE E0
2393
2394 c     scheme, y_cut=1.-THMAX.
2395
2396 C     If flavour specified mass effects fully included.
2397
2398 C     EVWGT=sig(e^-e^+ --> qqbar g) in nb
2399
2400 C-----------------------------------------------------------------------
2401
2402       INCLUDE 'HERWIG61.INC'
2403
2404       DOUBLE PRECISION HWR,HWUALF,HWUAEM,HWULDO,HWDPWT,Q2NOW,Q2LST,
2405
2406      & PHASP,QGMAX,QGMIN,FACTR,QM2,CLF(7),ORDER,PRAN,PQWT,QQG,QBG,SUM,
2407
2408      & RUT,QQLM,QQLP,QBLM,QBLP,DYN1,DYN2,DYN3,DYN4,DYN5,DYN6,XQ2,X2SUM,
2409
2410      & PVRT(4)
2411
2412       INTEGER ID1,IQ,I,LM,LP,IQ1
2413
2414       LOGICAL MASS
2415
2416       EXTERNAL HWR,HWUALF,HWUAEM,HWULDO,HWDPWT
2417
2418       SAVE Q2NOW,Q2LST,QGMAX,QGMIN,FACTR,ORDER,ID1,MASS,QM2,CLF,LM,LP,
2419
2420      & IQ1,QQG,QBG,SUM
2421
2422       DATA Q2LST/0./
2423
2424       IF (GENEV) THEN
2425
2426 C Label produced partons and calculate gluon spin
2427
2428         IDHW(NHEP+1)=200
2429
2430         IDHW(NHEP+2)=IQ1
2431
2432         IDHW(NHEP+3)=13
2433
2434         IDHW(NHEP+4)=IQ1+6
2435
2436         IDHEP(NHEP+1)=23
2437
2438         IDHEP(NHEP+2)=IQ1
2439
2440         IDHEP(NHEP+3)=21
2441
2442         IDHEP(NHEP+4)=-IQ1
2443
2444         ISTHEP(NHEP+1)=110
2445
2446         ISTHEP(NHEP+2)=113
2447
2448         ISTHEP(NHEP+3)=114
2449
2450         ISTHEP(NHEP+4)=114
2451
2452         JMOHEP(1,NHEP+1)=LM
2453
2454         JMOHEP(2,NHEP+1)=LP
2455
2456         JMOHEP(1,NHEP+2)=NHEP+1
2457
2458         JMOHEP(2,NHEP+2)=NHEP+3
2459
2460         JMOHEP(1,NHEP+3)=NHEP+1
2461
2462         JMOHEP(2,NHEP+3)=NHEP+4
2463
2464         JMOHEP(1,NHEP+4)=NHEP+1
2465
2466         JMOHEP(2,NHEP+4)=NHEP+2
2467
2468         JDAHEP(1,NHEP+1)=NHEP+2
2469
2470         JDAHEP(2,NHEP+1)=NHEP+4
2471
2472         JDAHEP(1,NHEP+2)=0
2473
2474         JDAHEP(2,NHEP+2)=NHEP+4
2475
2476         JDAHEP(1,NHEP+3)=0
2477
2478         JDAHEP(2,NHEP+3)=NHEP+2
2479
2480         JDAHEP(1,NHEP+4)=0
2481
2482         JDAHEP(2,NHEP+4)=NHEP+3
2483
2484 C Decide which quark radiated and assign production vertices
2485
2486         XQ2=(Q2NOW-2.*QBG)**2
2487
2488         X2SUM=XQ2+(Q2NOW-2.*QQG)**2
2489
2490         IF (XQ2.LT.HWR()*X2SUM) THEN
2491
2492 C Quark radiated the gluon
2493
2494           CALL HWVZRO(4,VHEP(1,NHEP+4))
2495
2496           CALL HWVSUM(4,PHEP(1,NHEP+2),PHEP(1,NHEP+3),PVRT)
2497
2498           CALL HWUDKL(IQ1,PVRT,VHEP(1,NHEP+3))
2499
2500           CALL HWVEQU(4,VHEP(1,NHEP+3),VHEP(1,NHEP+2))
2501
2502         ELSE
2503
2504 C Anti-quark radiated the gluon
2505
2506           CALL HWVZRO(4,VHEP(1,NHEP+2))
2507
2508           CALL HWVSUM(4,PHEP(1,NHEP+4),PHEP(1,NHEP+3),PVRT)
2509
2510           CALL HWUDKL(IQ1,PVRT,VHEP(1,NHEP+3))
2511
2512           CALL HWVEQU(4,VHEP(1,NHEP+3),VHEP(1,NHEP+4))
2513
2514         ENDIF
2515
2516         IF (AZSPIN) THEN
2517
2518 C  Calculate the transverse polarisation of the gluon
2519
2520 C  Correlation with leptons presently neglected
2521
2522            GPOLN=(QQG**2+QBG**2)/((Q2NOW-2.*SUM)*Q2NOW)
2523
2524            GPOLN=2./(2.+GPOLN)
2525
2526         ENDIF
2527
2528         NHEP=NHEP+4
2529
2530       ELSE
2531
2532         EMSCA=PHEP(5,3)
2533
2534         Q2NOW=EMSCA**2
2535
2536         IF (Q2NOW.NE.Q2LST) THEN
2537
2538           Q2LST=Q2NOW
2539
2540           PHASP=3.*THMAX-2.
2541
2542           IF (PHASP.LE.ZERO) CALL HWWARN('HWHEPG',400,*999)
2543
2544           QGMAX=.5*Q2NOW*THMAX
2545
2546           QGMIN=.5*Q2NOW*(1.-THMAX)
2547
2548           FACTR=GEV2NB*FLOAT(NCOLO)*CFFAC*HWUALF(1,EMSCA)
2549
2550      &         *.5*(HWUAEM(Q2NOW)*PHASP)**2/Q2NOW
2551
2552           LM=1
2553
2554           IF (JDAHEP(1,LM).NE.0) LM=JDAHEP(1,LM)
2555
2556           LP=2
2557
2558           IF (JDAHEP(1,LP).NE.0) LP=JDAHEP(1,LP)
2559
2560           ORDER=1.
2561
2562           IF (IDHW(1).GT.IDHW(2)) ORDER=-ORDER
2563
2564           ID1=MOD(IPROC,10)
2565
2566           IF (ID1.NE.0) THEN
2567
2568              MASS=.TRUE.
2569
2570              QM2=RMASS(ID1)**2
2571
2572              CALL HWUCFF(11,ID1,Q2NOW,CLF(1))
2573
2574              FACTR=FACTR*CLF(1)
2575
2576           ELSE
2577
2578              MASS=.FALSE.
2579
2580              CALL HWUEEC(1)
2581
2582              FACTR=FACTR*TQWT
2583
2584           ENDIF
2585
2586         ENDIF
2587
2588         IF (ID1.EQ.0) THEN
2589
2590 C Select quark flavour
2591
2592           PRAN=TQWT*HWR()
2593
2594           PQWT=0.
2595
2596           DO 10 IQ=1,MAXFL
2597
2598           PQWT=PQWT+CLQ(1,IQ)
2599
2600           IF (PQWT.GT.PRAN) GOTO 11
2601
2602    10     CONTINUE
2603
2604           IQ=MAXFL
2605
2606    11     IQ1=MAPQ(IQ)
2607
2608           DO 20 I=1,7
2609
2610    20     CLF(I)=CLQ(I,IQ)
2611
2612         ELSEIF (Q2NOW.GT.4*QM2/(2*THMAX-1)) THEN
2613
2614           IQ1=ID1
2615
2616         ELSE
2617
2618           EVWGT=0.
2619
2620           RETURN
2621
2622         ENDIF
2623
2624 C Select final state momentum configuration
2625
2626         CALL HWVEQU(5,PHEP(1,3),PHEP(1,NHEP+1))
2627
2628         PHEP(5,NHEP+2)=RMASS(IQ1)
2629
2630         PHEP(5,NHEP+3)=RMASS(13)
2631
2632         PHEP(5,NHEP+4)=RMASS(IQ1)
2633
2634    30   CALL HWDTHR(PHEP(1,NHEP+1),PHEP(1,NHEP+2),
2635
2636      &              PHEP(1,NHEP+3),PHEP(1,NHEP+4),HWDPWT)
2637
2638         QQG=HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+3))
2639
2640         IF (QQG.LT.QGMIN) GOTO 30
2641
2642         QBG=HWULDO(PHEP(1,NHEP+4),PHEP(1,NHEP+3))
2643
2644         SUM=QQG+QBG
2645
2646         IF (QBG.LT.QGMIN.OR.SUM.GT.QGMAX) GOTO 30
2647
2648         QQLM=HWULDO(PHEP(1,NHEP+2),PHEP(1,LM))
2649
2650         QQLP=HWULDO(PHEP(1,NHEP+2),PHEP(1,LP))
2651
2652         QBLM=HWULDO(PHEP(1,NHEP+4),PHEP(1,LM))
2653
2654         QBLP=HWULDO(PHEP(1,NHEP+4),PHEP(1,LP))
2655
2656         DYN1=QQLM**2+QQLP**2+QBLM**2+QBLP**2
2657
2658         DYN2=0.
2659
2660         DYN3=DYN1-2.*(QQLM**2+QBLP**2)
2661
2662         IF (MASS) THEN
2663
2664            RUT=1./QQG+1./QBG
2665
2666            DYN1=DYN1+8.*QM2*(1.-.25*Q2NOW*RUT
2667
2668      &         +QQLM*QQLP/(Q2NOW*QBG)+QBLM*QBLP/(Q2NOW*QQG))
2669
2670            DYN2=QM2*(Q2NOW-SUM*(2.+QM2*RUT)
2671
2672      &         -4.*HWULDO(PHEP(1,NHEP+3),PHEP(1,LM))
2673
2674      &            *HWULDO(PHEP(1,NHEP+3),PHEP(1,LP))/Q2NOW)
2675
2676            DYN3=DYN3+QM2*2.*RUT*(QBG*(QBLP-QBLM)-QQG*(QQLP-QQLM))
2677
2678         ENDIF
2679
2680         EVWGT=CLF(1)*DYN1+CLF(2)*DYN2+ORDER*CLF(3)*DYN3
2681
2682         IF (TPOL) THEN
2683
2684 C Include event plane azimuthal angle
2685
2686            DYN4=.5*Q2NOW
2687
2688            DYN5=DYN4
2689
2690            DYN6=0.
2691
2692            IF (MASS) THEN
2693
2694               DYN4=DYN4-QM2*SUM/QBG
2695
2696               DYN5=DYN5-QM2*SUM/QQG
2697
2698               DYN6=QM2
2699
2700            ENDIF
2701
2702            EVWGT=EVWGT
2703
2704      &     +(CLF(4)*COSS-CLF(6)*SINS)
2705
2706      &      *(DYN4*(PHEP(1,NHEP+2)**2-PHEP(2,NHEP+2)**2)
2707
2708      &       +DYN5*(PHEP(1,NHEP+4)**2-PHEP(2,NHEP+4)**2))
2709
2710      &     +(CLF(4)*SINS+CLF(6)*COSS)*2.
2711
2712      &      *(DYN4*PHEP(1,NHEP+2)*PHEP(2,NHEP+2)
2713
2714      &       +DYN5*PHEP(1,NHEP+4)*PHEP(2,NHEP+4))
2715
2716      &     +(CLF(5)*COSS-CLF(7)*SINS)*DYN6
2717
2718      &      *(PHEP(1,NHEP+3)**2-PHEP(2,NHEP+3)**2)
2719
2720      &     +(CLF(5)*SINS+CLF(7)*COSS)*DYN6*2.
2721
2722      &      *PHEP(1,NHEP+3)*PHEP(2,NHEP+3)
2723
2724         ENDIF
2725
2726 C Assign event weight
2727
2728         EVWGT=EVWGT*FACTR/(QQG*QBG*CLF(1))
2729
2730       ENDIF
2731
2732   999 END
2733
2734 CDECK  ID>, HWHEW0.
2735
2736 *CMZ :-        -26/04/91  11.11.55  by  Bryan Webber
2737
2738 *-- Author :    Zoltan Kunszt, modified by Bryan Webber & Mike Seymour
2739
2740 C-----------------------------------------------------------------------
2741
2742       SUBROUTINE HWHEW0(IP,ETOT,XM,PR,WEIGHT,CR)
2743
2744 C-----------------------------------------------------------------------
2745
2746       INCLUDE 'HERWIG61.INC'
2747
2748       DOUBLE PRECISION HWR,ETOT,XM(2),PR(5,2),WEIGHT,CR,XM1,XM2,S,
2749
2750      & D1,PABS,D,CX,C,E,F,SC,G
2751
2752       INTEGER IP,I
2753
2754       EXTERNAL HWR
2755
2756       WEIGHT=ZERO
2757
2758       XM1=XM(1)**2
2759
2760       XM2=XM(2)**2
2761
2762       S=ETOT*ETOT
2763
2764       D1=S-XM1-XM2
2765
2766       PABS=D1*D1-4.*XM1*XM2
2767
2768       IF (PABS.LE.ZERO) RETURN
2769
2770       PABS=SQRT(PABS)
2771
2772       D=D1/PABS
2773
2774       IF(IP.EQ.2)GOTO3
2775
2776       CX=CR
2777
2778       C=D-(D+CX)*((D-CR)/(D+CX))**HWR()
2779
2780       GOTO 4
2781
2782 3     E=((D+ONE)/(D-ONE))*(TWO*HWR()-ONE)
2783
2784       C=D*((E-ONE)/(E+ONE))
2785
2786 4     F=2D0*PIFAC*HWR()
2787
2788       SC=SQRT(ONE-C*C)
2789
2790       PR(4,1)=(S+XM1-XM2)/(TWO*ETOT)
2791
2792       PR(5,1)=PR(4,1)*PR(4,1)-XM1
2793
2794       IF (PR(5,1).LE.ZERO) RETURN
2795
2796       PR(5,1)=SQRT(PR(5,1))
2797
2798       PR(4,2)=ETOT-PR(4,1)
2799
2800       PR(3,1)=PR(5,1)*C
2801
2802       PR(5,2)=PR(5,1)
2803
2804       PR(2,1)=PR(5,1)*SC*COS(F)
2805
2806       PR(1,1)=PR(5,1)*SC*SIN(F)
2807
2808       DO 7 I=1,3
2809
2810 7     PR(I,2)=-PR(I,1)
2811
2812       G=0.
2813
2814       IF(IP.EQ.1)G=(D-C)*LOG((D+CX)/(D-CR))
2815
2816       IF(IP.EQ.2)G=(D*D-C*C)/D*LOG((D+ONE)/(D-ONE))
2817
2818       WEIGHT=PIFAC*G*PR(5,1)/ETOT*HALF
2819
2820       RETURN
2821
2822       END