HERWIG/src/hwhbgf.f

   1
   2 CDECK  ID>, HWHBGF.
   3
   4 *CMZ :-        -26/04/91  11.11.55  by  Bryan Webber
   5
   6 *-- Author :    Giovanni Abbiendi & Luca Stanco
   7
   8 C-----------------------------------------------------------------------
   9
  10       SUBROUTINE HWHBGF
  11
  12 C-----------------------------------------------------------------------
  13
  14 C     Order Alpha_s processes in charged lepton-hadron collisions
  15
  16 C
  17
  18 C       Process code IPROC has to be set in the Main Program
  19
  20 C       the following codes IPROC may be selected
  21
  22 C
  23
  24 C                9100 : NC  BOSON-GLUON FUSION
  25
  26 C                9100+IQK (IQK=1,...,6) :  produced flavour is IQK
  27
  28 C                9107 : produced  J/psi + gluon
  29
  30 C
  31
  32 C                9110 : NC  QCD COMPTON
  33
  34 C                9110+IQK (IQK=1,...,12) : struck parton is IQK
  35
  36 C
  37
  38 C                9130 : NC order alpha_s processes (9100+9110)
  39
  40 C
  41
  42 C       Select maximum and minimum generated flavour when IQK=0
  43
  44 C       setting IFLMIN and IFLMAX in the Main Program
  45
  46 C       (allowed values from 1 to 6), default are 1 and 5
  47
  48 C       allowing d,u,s,c,b,dbar,ubar,sbar,cbar,bbar
  49
  50 C
  51
  52 C           CHARGED CURRENT Boson-Gluon Fusion processes
  53
  54 C                9141 : CC  s cbar  (c sbar)
  55
  56 C                9142 : CC  b cbar  (c bbar)
  57
  58 C                9143 : CC  s tbar  (t cbar)
  59
  60 C                9144 : CC  b tbar  (t bbar)
  61
  62 C
  63
  64 C       other inputs : Q2MIN,Q2MAX,YBMIN,YBMAX,PTMIN,EMMIN,EMMAX
  65
  66 C       when IPROC=(1)9107 : as above but Q2WWMN, Q2WWMX substitute
  67
  68 C                            Q2MIN and Q2MAX (EPA is used); ZJMAX cut
  69
  70 C
  71
  72 C      Add 10000 to suppress soft remnant fragmentation
  73
  74 C
  75
  76 C      Mean EVWGT = cross section in nanoBarn
  77
  78 C
  79
  80 C-----------------------------------------------------------------------
  81
  82       INCLUDE 'HERWIG61.INC'
  83
  84       DOUBLE PRECISION HWR,LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,
  85
  86      & ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT,FSIGMA(18),
  87
  88      & SIGSUM,PROB,PRAN,PVRT(4),X
  89
  90       INTEGER IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,IPROO,LEPFIN,ID1,ID2,I,IDD
  91
  92       LOGICAL CHARGD,INCLUD(18),INSIDE(18)
  93
  94       EXTERNAL HWR
  95
  96       SAVE LEPFIN,ID1,ID2,FSIGMA,SIGSUM
  97
  98       COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF,
  99
 100      & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,
 101
 102      & IPROO,CHARGD,INCLUD,INSIDE
 103
 104 C---Initialization
 105
 106       IF (FSTWGT) THEN
 107
 108 C---LEP = 1 FOR LEPTONS, -1 FOR ANTILEPTONS
 109
 110         LEP=ZERO
 111
 112         IF (IDHW(1).GE.121.AND.IDHW(1).LE.126) THEN
 113
 114           LEP=ONE
 115
 116         ELSEIF (IDHW(1).GE.127.AND.IDHW(1).LE.132) THEN
 117
 118           LEP=-ONE
 119
 120         ENDIF
 121
 122         IF (LEP.EQ.ZERO) CALL HWWARN('HWHBGF',500,*999)
 123
 124         IPROO=MOD(IPROC,100)/10
 125
 126         IF (IPROO.EQ.0.OR.IPROO.EQ.4) THEN
 127
 128           IQK=MOD(IPROC,10)
 129
 130           IFL=IQK
 131
 132           IF (IQK.EQ.7) IFL=164
 133
 134           CHARGD=IPROO.EQ.4
 135
 136         ELSEIF (IPROO.EQ.1.OR.IPROO.EQ.2) THEN
 137
 138           IQK=MOD(IPROC,100)-10
 139
 140           IFL=IQK+6
 141
 142           CHARGD=.FALSE.
 143
 144         ELSEIF (IPROO.EQ.3) THEN
 145
 146           IQK=0
 147
 148           IFL=0
 149
 150           CHARGD=.FALSE.
 151
 152         ELSE
 153
 154           CALL HWWARN('HWHBGF',501,*999)
 155
 156         ENDIF
 157
 158 C
 159
 160         LEPFIN = IDHW(1)
 161
 162         IF(CHARGD) THEN
 163
 164           LEPFIN = IDHW(1)+1
 165
 166           IF (IQK.EQ.1) THEN
 167
 168             IFLAVU=4
 169
 170             IFLAVD=3
 171
 172             ID1  = 3
 173
 174             ID2  = 10
 175
 176           ELSEIF (IQK.EQ.2) THEN
 177
 178             IFLAVU=4
 179
 180             IFLAVD=5
 181
 182             ID1  = 5
 183
 184             ID2  = 10
 185
 186           ELSEIF (IQK.EQ.3) THEN
 187
 188             IFLAVU=6
 189
 190             IFLAVD=3
 191
 192             ID1  = 3
 193
 194             ID2  =12
 195
 196           ELSE
 197
 198             IFLAVU=6
 199
 200             IFLAVD=5
 201
 202             ID1  = 5
 203
 204             ID2  =12
 205
 206           ENDIF
 207
 208           IF (LEP.EQ.-1.) THEN
 209
 210             IDD=ID1
 211
 212             ID1=ID2-6
 213
 214             ID2=IDD+6
 215
 216           ENDIF
 217
 218         ENDIF
 219
 220 C
 221
 222         IF (IQK.EQ.0) THEN
 223
 224           DO I=1,18
 225
 226             INCLUD(I)=.TRUE.
 227
 228           ENDDO
 229
 230           IMIN=1
 231
 232           IMAX=18
 233
 234           DO I=1,6
 235
 236             IF (I.LT.IFLMIN.OR.I.GT.IFLMAX) INCLUD(I)=.FALSE.
 237
 238           ENDDO
 239
 240           DO I=7,18
 241
 242             IF (I.LE.12) THEN
 243
 244               IF (I-6.LT.IFLMIN.OR.I-6.GT.IFLMAX) INCLUD(I)=.FALSE.
 245
 246             ELSE
 247
 248               IF (I-12.LT.IFLMIN.OR.I-12.GT.IFLMAX) INCLUD(I)=.FALSE.
 249
 250             ENDIF
 251
 252           ENDDO
 253
 254           IF (IPROO.EQ.0) THEN
 255
 256             DO I=7,18
 257
 258               INCLUD(I)=.FALSE.
 259
 260             ENDDO
 261
 262             IMIN=IFLMIN
 263
 264             IMAX=IFLMAX
 265
 266           ELSEIF (IPROO.EQ.1.OR.IPROO.EQ.2) THEN
 267
 268             DO I=1,6
 269
 270               INCLUD(I)=.FALSE.
 271
 272             ENDDO
 273
 274             IMIN=IFLMIN+6
 275
 276             IMAX=IFLMAX+12
 277
 278           ELSEIF (IPROO.EQ.3) THEN
 279
 280             IMIN=IFLMIN
 281
 282             IMAX=IFLMAX+12
 283
 284           ENDIF
 285
 286         ELSEIF (IQK.NE.0 .AND. (.NOT.CHARGD)) THEN
 287
 288           DO I=1,18
 289
 290             INCLUD(I)=.FALSE.
 291
 292           ENDDO
 293
 294           IF (IFL.LE.18) THEN
 295
 296             INCLUD(IFL)=.TRUE.
 297
 298             IMIN=IFL
 299
 300             IMAX=IFL
 301
 302           ELSEIF (IFL.EQ.164) THEN
 303
 304             INCLUD(7)=.TRUE.
 305
 306             IMIN=7
 307
 308             IMAX=7
 309
 310           ENDIF
 311
 312         ENDIF
 313
 314       ENDIF
 315
 316 C---End of initialization
 317
 318       IF(GENEV) THEN
 319
 320       IF (.NOT.CHARGD) THEN
 321
 322         IF (IQK.EQ.0) THEN
 323
 324           PRAN= SIGSUM * HWR()
 325
 326           PROB=ZERO
 327
 328           DO 10 IFL=IMIN,IMAX
 329
 330             IF (.NOT.INSIDE(IFL)) GOTO 10
 331
 332             PROB=PROB+FSIGMA(IFL)
 333
 334             IF (PROB.GE.PRAN) GOTO 20
 335
 336   10      CONTINUE
 337
 338         ENDIF
 339
 340 C---at this point the subprocess has been selected (IFL)
 341
 342   20    CONTINUE
 343
 344         IF (IFL.LE.6) THEN
 345
 346 C---Boson-Gluon Fusion event
 347
 348           IDHW(NHEP+1)=IDHW(1)
 349
 350           IDHW(NHEP+2)=13
 351
 352           IDHW(NHEP+3)=15
 353
 354           IDHW(NHEP+4)=LEPFIN
 355
 356           IDHW(NHEP+5)=IFL
 357
 358           IDHW(NHEP+6)=IFL+6
 359
 360         ELSEIF (IFL.GE.7.AND.IFL.LE.18) THEN
 361
 362 C---QCD_Compton event
 363
 364           IDHW(NHEP+1)=IDHW(1)
 365
 366           IDHW(NHEP+2)=IFL-6
 367
 368           IDHW(NHEP+3)=15
 369
 370           IDHW(NHEP+4)=LEPFIN
 371
 372           IDHW(NHEP+5)=IFL-6
 373
 374           IDHW(NHEP+6)=13
 375
 376         ELSEIF (IFL.EQ.164) THEN
 377
 378 C---gamma+gluon-->J/Psi+gluon
 379
 380           IDHW(NHEP+1)=IDHW(1)
 381
 382           IDHW(NHEP+2)=13
 383
 384           IDHW(NHEP+3)=15
 385
 386           IDHW(NHEP+4)=LEPFIN
 387
 388           IDHW(NHEP+5)=164
 389
 390           IDHW(NHEP+6)=13
 391
 392         ELSE
 393
 394           CALL HWWARN('HWHBGF',503,*999)
 395
 396         ENDIF
 397
 398       ELSE
 399
 400 C---Charged current event of specified flavours
 401
 402         IDHW(NHEP+1)=IDHW(1)
 403
 404         IDHW(NHEP+2)=13
 405
 406         IDHW(NHEP+3)=15
 407
 408         IDHW(NHEP+4)=LEPFIN
 409
 410         IDHW(NHEP+5)=ID1
 411
 412         IDHW(NHEP+6)=ID2
 413
 414       ENDIF
 415
 416 C
 417
 418       DO 1 I=NHEP+1,NHEP+6
 419
 420     1 IDHEP(I)=IDPDG(IDHW(I))
 421
 422 C
 423
 424 C---Codes common for all processes
 425
 426       ISTHEP(NHEP+1)=111
 427
 428       ISTHEP(NHEP+2)=112
 429
 430       ISTHEP(NHEP+3)=110
 431
 432       ISTHEP(NHEP+4)=113
 433
 434       ISTHEP(NHEP+5)=114
 435
 436       ISTHEP(NHEP+6)=114
 437
 438 C
 439
 440       DO I=NHEP+1,NHEP+6
 441
 442         JMOHEP(1,I)=NHEP+3
 443
 444         JDAHEP(1,I)=0
 445
 446       ENDDO
 447
 448 C---Incoming lepton
 449
 450       JMOHEP(2,NHEP+1)=NHEP+4
 451
 452       JDAHEP(2,NHEP+1)=NHEP+4
 453
 454 C---Hard Process C.M.
 455
 456       JMOHEP(1,NHEP+3)=NHEP+1
 457
 458       JMOHEP(2,NHEP+3)=NHEP+2
 459
 460       JDAHEP(1,NHEP+3)=NHEP+4
 461
 462       JDAHEP(2,NHEP+3)=NHEP+6
 463
 464 C---Outgoing lepton
 465
 466       JMOHEP(2,NHEP+4)=NHEP+1
 467
 468       JDAHEP(2,NHEP+4)=NHEP+1
 469
 470 C
 471
 472       IF (IFL.LE.6 .OR. CHARGD) THEN
 473
 474 C---Codes for boson-gluon fusion processes
 475
 476 C---  Incoming gluon
 477
 478         JMOHEP(2,NHEP+2)=NHEP+6
 479
 480         JDAHEP(2,NHEP+2)=NHEP+5
 481
 482 C---  Outgoing quark
 483
 484         JMOHEP(2,NHEP+5)=NHEP+2
 485
 486         JDAHEP(2,NHEP+5)=NHEP+6
 487
 488 C---  Outgoing antiquark
 489
 490         JMOHEP(2,NHEP+6)=NHEP+5
 491
 492         JDAHEP(2,NHEP+6)=NHEP+2
 493
 494       ELSEIF (IFL.GE.7 .AND. IFL.LE.12) THEN
 495
 496 C---Codes for V+q --> q+g
 497
 498 C---  Incoming quark
 499
 500         JMOHEP(2,NHEP+2)=NHEP+5
 501
 502         JDAHEP(2,NHEP+2)=NHEP+6
 503
 504 C---  Outgoing quark
 505
 506         JMOHEP(2,NHEP+5)=NHEP+6
 507
 508         JDAHEP(2,NHEP+5)=NHEP+2
 509
 510 C---  Outgoing gluon
 511
 512         JMOHEP(2,NHEP+6)=NHEP+2
 513
 514         JDAHEP(2,NHEP+6)=NHEP+5
 515
 516       ELSEIF (IFL.GE.13 .AND. IFL.LE.18) THEN
 517
 518 C---Codes for V+qbar --> qbar+g
 519
 520 C---  Incoming antiquark
 521
 522         JMOHEP(2,NHEP+2)=NHEP+6
 523
 524         JDAHEP(2,NHEP+2)=NHEP+5
 525
 526 C---  Outgoing antiquark
 527
 528         JMOHEP(2,NHEP+5)=NHEP+2
 529
 530         JDAHEP(2,NHEP+5)=NHEP+6
 531
 532 C---  Outgoing gluon
 533
 534         JMOHEP(2,NHEP+6)=NHEP+5
 535
 536         JDAHEP(2,NHEP+6)=NHEP+2
 537
 538       ELSEIF (IFL.EQ.164) THEN
 539
 540 C---Codes for Gamma+gluon --> J/Psi+gluon
 541
 542 C---  Incoming gluon
 543
 544         JMOHEP(2,NHEP+2)=NHEP+6
 545
 546         JDAHEP(2,NHEP+2)=NHEP+6
 547
 548 C---  Outgoing J/Psi
 549
 550         JMOHEP(2,NHEP+5)=NHEP+1
 551
 552         JDAHEP(2,NHEP+5)=NHEP+1
 553
 554 C---  Outgoing gluon
 555
 556         JMOHEP(2,NHEP+6)=NHEP+2
 557
 558         JDAHEP(2,NHEP+6)=NHEP+2
 559
 560       ENDIF
 561
 562 C---Computation of momenta in Laboratory frame of reference
 563
 564       CALL HWHBKI
 565
 566       NHEP=NHEP+6
 567
 568 C Decide which quark radiated and assign production vertices
 569
 570       IF (IFL.LE.6) THEN
 571
 572 C Boson-Gluon fusion case
 573
 574         IF (1-Z.LT.HWR()) THEN
 575
 576 C Gluon splitting to quark
 577
 578           CALL HWVZRO(4,VHEP(1,NHEP-1))
 579
 580           CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP),PVRT)
 581
 582           CALL HWUDKL(IFL,PVRT,VHEP(1,NHEP))
 583
 584           CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-4))
 585
 586         ELSE
 587
 588 C Gluon splitting to antiquark
 589
 590           CALL HWVZRO(4,VHEP(1,NHEP))
 591
 592           CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP-1),PVRT)
 593
 594           CALL HWUDKL(IFL,PVRT,VHEP(1,NHEP-1))
 595
 596           CALL HWVEQU(4,VHEP(1,NHEP-1),VHEP(1,NHEP-4))
 597
 598         ENDIF
 599
 600       ELSEIF (IFL.GE.7.AND.IFL.LE.18) THEN
 601
 602 C QCD Compton case
 603
 604         X=1/(1+SHAT/Q2)
 605
 606         IF (1.LT.HWR()*(1+(1-X-Z)**2+6*X*(1-X)*Z*(1-Z))) THEN
 607
 608 C Incoming quark radiated the gluon
 609
 610           CALL HWVZRO(4,VHEP(1,NHEP-1))
 611
 612           CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP),PVRT)
 613
 614           CALL HWUDKL(IFL-6,PVRT,VHEP(1,NHEP))
 615
 616           CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-4))
 617
 618         ELSE
 619
 620 C Outgoing quark radiated the gluon
 621
 622           CALL HWVZRO(4,VHEP(1,NHEP-4))
 623
 624           CALL HWVSUM(4,PHEP(1,NHEP-1),PHEP(1,NHEP),PVRT)
 625
 626           CALL HWUDKL(IFL-6,PVRT,VHEP(1,NHEP))
 627
 628           CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-1))
 629
 630         ENDIF
 631
 632       ENDIF
 633
 634 C---HERWIG gets confused if lepton momentum is different from beam
 635
 636 C   momentum, which it can be if incoming hadron has negative virtuality
 637
 638 C   As a temporary fix, simply copy the momentum.
 639
 640 C   Momentum conservation somehow gets taken care of HWBGEN!
 641
 642       call hwvequ(5,phep(1,1),phep(1,nhep-5))
 643
 644       ELSE
 645
 646         EVWGT=ZERO
 647
 648 C---generation of the 5 variables Y,Q2,SHAT,Z,PHI and Jacobian computation
 649
 650 C---in the largest phase space avalaible for selected processes and
 651
 652 C---filling of logical vector INSIDE to tag contributing ones
 653
 654         CALL HWHBRN (*999)
 655
 656 C---calculate differential cross section corresponding to the chosen
 657
 658 C---variables and the weight for MC generation
 659
 660         IF (IQK.EQ.0) THEN
 661
 662 C---many subprocesses included
 663
 664           DO I=1,18
 665
 666             FSIGMA(I)=ZERO
 667
 668           ENDDO
 669
 670           SIGSUM=ZERO
 671
 672           DO I=IMIN,IMAX
 673
 674             IF (INSIDE(I)) THEN
 675
 676               IFL=I
 677
 678               DSIGMA=ZERO
 679
 680               CALL HWHBSG
 681
 682               FSIGMA(I)=DSIGMA
 683
 684               SIGSUM=SIGSUM+DSIGMA
 685
 686             ENDIF
 687
 688           ENDDO
 689
 690           EVWGT=SIGSUM * AJACOB
 691
 692         ELSE
 693
 694 C---only one subprocess included
 695
 696           CALL HWHBSG
 697
 698           EVWGT= DSIGMA * AJACOB
 699
 700         ENDIF
 701
 702         IF (EVWGT.LT.ZERO) EVWGT=ZERO
 703
 704       ENDIF
 705
 706   999 END
 707
 708 CDECK  ID>, HWHBKI.
 709
 710 *CMZ :-        -26/04/91  13.19.32  by  Federico Carminati
 711
 712 *-- Author :    Giovanni Abbiendi & Luca Stanco
 713
 714 C----------------------------------------------------------------------
 715
 716       SUBROUTINE HWHBKI
 717
 718 C----------------------------------------------------------------------
 719
 720 C     gives the fourmomenta in the laboratory system for the particles
 721
 722 C     of the hard 2-->3 subprocess, to match with HERWIG routines of
 723
 724 C     jet evolution.
 725
 726 C----------------------------------------------------------------------
 727
 728       INCLUDE 'HERWIG61.INC'
 729
 730       DOUBLE PRECISION HWUECM,HWUPCM,HWUSQR,LEP,Y,Q2,SHAT,Z,PHI,AJACOB,
 731
 732      & DSIGMA,ME,MP,ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT,
 733
 734      & PGAMMA(5),SG,MF1,MF2,EP,PP,EL,PL,E1,E2,Q1,COSBET,SINBET,COSTHE,
 735
 736      & SINTHE,SINAZI,COSAZI,ROTAZI(3,3),EGAM,A,PPROT,MREMIN,PGAM,PEP(5),
 737
 738      & COSPHI,SINPHI,ROT(3,3),EPROT,PROTON(5),MPART
 739
 740       INTEGER IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,IPROO,I,IHAD,J,IS,ICMF
 741
 742       LOGICAL CHARGD,INCLUD(18),INSIDE(18)
 743
 744       EXTERNAL HWUECM,HWUPCM,HWUSQR
 745
 746       COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF,
 747
 748      & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,
 749
 750      & IPROO,CHARGD,INCLUD,INSIDE
 751
 752 C
 753
 754       IHAD=2
 755
 756       IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD)
 757
 758 C---Set masses
 759
 760       IF (CHARGD) THEN
 761
 762         MPART=ZERO
 763
 764         MF1=RMASS(IDHW(NHEP+5))
 765
 766         MF2=RMASS(IDHW(NHEP+6))
 767
 768         MREMIN=MP
 769
 770       ELSE
 771
 772         IS = IFL
 773
 774         IF (IFL.EQ.164) IS=IQK
 775
 776         MPART=ZERO
 777
 778         IF (IFL.GE.7.AND.IFL.LE.18) MPART=RMASS(IFL-6)
 779
 780         MF1=MFIN1(IS)
 781
 782         MF2=MFIN2(IS)
 783
 784         MREMIN = MREMIF(IS)
 785
 786       ENDIF
 787
 788 C---Calculation of kinematical variables for the generated event
 789
 790 C   in the center of mass frame of the incoming boson and parton
 791
 792 C   with parton along +z
 793
 794       EGAM = HWUECM (SHAT, -Q2, MPART**2)
 795
 796       PGAM = SQRT( EGAM**2 + Q2 )
 797
 798       EP = RSHAT-EGAM
 799
 800       PP = PGAM
 801
 802       A = (W2+Q2-MP**2)/TWO
 803
 804       PPROT = (A*PGAM-EGAM*SQRT(A**2+MP**2*Q2))/Q2
 805
 806       IF (PPROT.LT.ZERO) CALL HWWARN('HWHBKI',101,*999)
 807
 808       EPROT = SQRT(PPROT**2+MP**2)
 809
 810       IF ((EPROT+PPROT).LT.(EP+PP)) CALL HWWARN('HWHBKI',102,*999)
 811
 812       EL = ( PGAM / PPROT * SMA - Q2 ) / TWO
 813
 814      +     / (EGAM + PGAM / PPROT * EPROT)
 815
 816       IF (EL.GT.ME) THEN
 817
 818         PL = SQRT ( EL**2 - ME**2 )
 819
 820       ELSE
 821
 822         CALL HWWARN ('HWHBKI',103,*999)
 823
 824       ENDIF
 825
 826       COSBET = (TWO * EPROT * EL - SMA) / (TWO * PPROT * PL)
 827
 828       IF ( ABS(COSBET) .GE. ONE ) THEN
 829
 830         COSBET = SIGN (ONE,COSBET)
 831
 832         SINBET = ZERO
 833
 834       ELSE
 835
 836         SINBET = SQRT (ONE - COSBET**2)
 837
 838       ENDIF
 839
 840       SG = ME**2 + MPART**2 + Q2 + TWO * RSHAT * EL
 841
 842       IF (SG.LE.(RSHAT+ML)**2 .OR. SG.GE.(RS-MREMIN)**2)
 843
 844      +    CALL HWWARN ('HWHBKI',104,*999)
 845
 846       Q1 = HWUPCM( RSHAT, MF1, MF2)
 847
 848       E1 = SQRT(Q1**2+MF1**2)
 849
 850       E2 = SQRT(Q1**2+MF2**2)
 851
 852       IF (Q1 .GT. ZERO) THEN
 853
 854         COSTHE=(TWO*EP*E1 - Z*(SHAT+Q2))/(TWO*PP*Q1)
 855
 856         IF (ABS(COSTHE) .GT. ONE) THEN
 857
 858           COSTHE=SIGN(ONE,COSTHE)
 859
 860           SINTHE=ZERO
 861
 862         ELSE
 863
 864           SINTHE=SQRT(ONE-COSTHE**2)
 865
 866         ENDIF
 867
 868       ELSE
 869
 870         COSTHE=ZERO
 871
 872         SINTHE=ONE
 873
 874       ENDIF
 875
 876 C---Initial lepton
 877
 878       PHEP(1,NHEP+1)=PL*SINBET
 879
 880       PHEP(2,NHEP+1)=ZERO
 881
 882       PHEP(3,NHEP+1)=PL*COSBET
 883
 884       PHEP(4,NHEP+1)=EL
 885
 886       PHEP(5,NHEP+1)=RMASS(IDHW(1))
 887
 888 C---Initial Hadron
 889
 890       PROTON(1)=ZERO
 891
 892       PROTON(2)=ZERO
 893
 894       PROTON(3)=PPROT
 895
 896       PROTON(4)=EPROT
 897
 898       CALL HWUMAS (PROTON)
 899
 900 C---Initial parton
 901
 902       PHEP(1,NHEP+2)=ZERO
 903
 904       PHEP(2,NHEP+2)=ZERO
 905
 906       PHEP(3,NHEP+2)=PP
 907
 908       PHEP(4,NHEP+2)=EP
 909
 910       PHEP(5,NHEP+2)=MPART
 911
 912 C---HARD SUBPROCESS 2-->3 CENTRE OF MASS
 913
 914       PHEP(1,NHEP+3)=PHEP(1,NHEP+1)+PHEP(1,NHEP+2)
 915
 916       PHEP(2,NHEP+3)=PHEP(2,NHEP+1)+PHEP(2,NHEP+2)
 917
 918       PHEP(3,NHEP+3)=PHEP(3,NHEP+1)+PHEP(3,NHEP+2)
 919
 920       PHEP(4,NHEP+3)=PHEP(4,NHEP+1)+PHEP(4,NHEP+2)
 921
 922       CALL HWUMAS  ( PHEP(1,NHEP+3) )
 923
 924 C---Virtual boson
 925
 926       PGAMMA(1)=ZERO
 927
 928       PGAMMA(2)=ZERO
 929
 930       PGAMMA(3)=-PGAM
 931
 932       PGAMMA(4)=EGAM
 933
 934       PGAMMA(5)=HWUSQR(Q2)
 935
 936 C---Scattered lepton
 937
 938       PHEP(1,NHEP+4)=PHEP(1,NHEP+1)-PGAMMA(1)
 939
 940       PHEP(2,NHEP+4)=PHEP(2,NHEP+1)-PGAMMA(2)
 941
 942       PHEP(3,NHEP+4)=PHEP(3,NHEP+1)-PGAMMA(3)
 943
 944       PHEP(4,NHEP+4)=PHEP(4,NHEP+1)-PGAMMA(4)
 945
 946       PHEP(5,NHEP+4)=RMASS(IDHW(1))
 947
 948       IF (CHARGD) PHEP(5,NHEP+4)=ZERO
 949
 950 C---First Final parton:  quark (or J/psi) in Boson-Gluon Fusion
 951
 952 C---                     quark or antiquark in QCD Compton
 953
 954       PHEP(1,NHEP+5)=Q1*SINTHE*COS(PHI)
 955
 956       PHEP(2,NHEP+5)=Q1*SINTHE*SIN(PHI)
 957
 958       PHEP(3,NHEP+5)=Q1*COSTHE
 959
 960       PHEP(4,NHEP+5)=E1
 961
 962       PHEP(5,NHEP+5)=MF1
 963
 964 C---Second Final parton: antiquark in Boson-Gluon Fusion
 965
 966 C---                     gluon in QCD Compton
 967
 968       PHEP(1,NHEP+6)=-PHEP(1,NHEP+5)
 969
 970       PHEP(2,NHEP+6)=-PHEP(2,NHEP+5)
 971
 972       PHEP(3,NHEP+6)=-PHEP(3,NHEP+5)
 973
 974       PHEP(4,NHEP+6)=E2
 975
 976       PHEP(5,NHEP+6)=MF2
 977
 978 C---Boost to lepton-hadron CM frame
 979
 980       PEP(1) = PHEP(1,NHEP+1)
 981
 982       PEP(2) = PHEP(2,NHEP+1)
 983
 984       PEP(3) = PHEP(3,NHEP+1) + PPROT
 985
 986       PEP(4) = PHEP(4,NHEP+1) + EPROT
 987
 988       CALL HWUMAS (PEP)
 989
 990       DO I=1,6
 991
 992         CALL HWULOF (PEP,PHEP(1,NHEP+I),PHEP(1,NHEP+I))
 993
 994       ENDDO
 995
 996       CALL HWULOF (PEP,PROTON,PROTON)
 997
 998       CALL HWULOF (PEP,PGAMMA,PGAMMA)
 999
1000 C---Rotation around y-axis to align lepton beam with z-axis
1001
1002       COSPHI = PHEP(3,NHEP+1) /
1003
1004      &           SQRT( PHEP(1,NHEP+1)**2 + PHEP(3,NHEP+1)**2 )
1005
1006       SINPHI = PHEP(1,NHEP+1) /
1007
1008      &           SQRT( PHEP(1,NHEP+1)**2 + PHEP(3,NHEP+1)**2 )
1009
1010       DO I=1,3
1011
1012       DO J=1,3
1013
1014         ROT(I,J)=ZERO
1015
1016       ENDDO
1017
1018       ENDDO
1019
1020         ROT(1,1) = COSPHI
1021
1022         ROT(1,3) = -SINPHI
1023
1024         ROT(2,2) = ONE
1025
1026         ROT(3,1) = SINPHI
1027
1028         ROT(3,3) = COSPHI
1029
1030       DO I=1,6
1031
1032         CALL HWUROF (ROT,PHEP(1,NHEP+I),PHEP(1,NHEP+I))
1033
1034       ENDDO
1035
1036       CALL HWUROF (ROT,PROTON,PROTON)
1037
1038       CALL HWUROF (ROT,PGAMMA,PGAMMA)
1039
1040 C---Boost to the LAB frame
1041
1042       ICMF=3
1043
1044       DO I=1,6
1045
1046         CALL HWULOB (PHEP(1,ICMF),PHEP(1,NHEP+I),PHEP(1,NHEP+I))
1047
1048       ENDDO
1049
1050       CALL HWULOB (PHEP(1,ICMF),PROTON,PROTON)
1051
1052       CALL HWULOB (PHEP(1,ICMF),PGAMMA,PGAMMA)
1053
1054 C---Random azimuthal rotation
1055
1056       CALL HWRAZM (ONE,COSAZI,SINAZI)
1057
1058       DO I=1,3
1059
1060       DO J=1,3
1061
1062         ROTAZI(I,J)=ZERO
1063
1064       ENDDO
1065
1066       ENDDO
1067
1068         ROTAZI(1,1) = COSAZI
1069
1070         ROTAZI(1,2) = SINAZI
1071
1072         ROTAZI(2,1) = -SINAZI
1073
1074         ROTAZI(2,2) = COSAZI
1075
1076         ROTAZI(3,3) = ONE
1077
1078       DO I=1,6
1079
1080         CALL HWUROF (ROTAZI,PHEP(1,NHEP+I),PHEP(1,NHEP+I))
1081
1082       ENDDO
1083
1084       CALL HWUROF (ROTAZI,PROTON,PROTON)
1085
1086       CALL HWUROF (ROTAZI,PGAMMA,PGAMMA)
1087
1088   999 END
1089
1090 CDECK  ID>, HWHBRN.
1091
1092 *CMZ :-        -03/07/95  19.02.12  by  Giovanni Abbiendi
1093
1094 *-- Author :    Giovanni Abbiendi & Luca Stanco
1095
1096 C-----------------------------------------------------------------------
1097
1098       SUBROUTINE HWHBRN (*)
1099
1100 C----------------------------------------------------------------------
1101
1102 C     Returns a point in the phase space (Y,Q2,SHAT,Z,PHI) and the
1103
1104 C     corresponding Jacobian factor AJACOB
1105
1106 C     Fill the logical vector INSIDE to tag contributing subprocesses
1107
1108 C     to the cross-section
1109
1110 C-----------------------------------------------------------------------
1111
1112       INCLUDE 'HERWIG61.INC'
1113
1114       DOUBLE PRECISION HWRUNI,HWR,HWUPCM,LEP,Y,Q2,SHAT,Z,PHI,AJACOB,
1115
1116      & DSIGMA,ME,MP,ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT,
1117
1118      & MF1,MF2,YMIN,YMAX,YJAC,Q2INF,Q2SUP,Q2JAC,EMW2,ZMIN,ZMAX,ZJAC,
1119
1120      & GAMMA2,LAMBDA,PHIJAC,ZINT,ZLMIN,ZL,EMW,TMIN,TMAX,EMLMIN,EMLMAX,
1121
1122      & SHMIN,EMMIF(18),EMMAF(18),WMIF(18),WMIN,MREMIN,YMIF(18),Q1CM(18),
1123
1124      & Q2MAF(18),EMMAWF(18),ZMIF(18),ZMAF(18),PLMAX,PINC,SHINF,SHSUP,
1125
1126      & SHJAC,CTHLIM,Q1,DETDSH,SRY,SRY0,SRY1
1127
1128       INTEGER IQK,IFLAVU,IFLAVD,I,IMIN,IMAX,IFL,IPROO,IHAD,NTRY,DEBUG
1129
1130       LOGICAL CHARGD,INCLUD(18),INSIDE(18)
1131
1132       EXTERNAL HWRUNI,HWR,HWUPCM
1133
1134       SAVE EMLMIN,EMLMAX,EMMIF,EMMAF,MREMIN,MF1,MF2,YMIF,
1135
1136      &     YMIN,YMAX,WMIN,WMIF
1137
1138       COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF,
1139
1140      & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,
1141
1142      & IPROO,CHARGD,INCLUD,INSIDE
1143
1144       EQUIVALENCE (EMW,RMASS(198))
1145
1146 C
1147
1148       IHAD=2
1149
1150       IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD)
1151
1152 C---Initialization
1153
1154       IF (FSTWGT.OR.IHAD.NE.2) THEN
1155
1156         ME = RMASS(IDHW(1))
1157
1158         MP = RMASS(IDHW(IHAD))
1159
1160         RS = PHEP(5,3)
1161
1162         SMA = RS**2-ME**2-MP**2
1163
1164         PINC = HWUPCM(RS,ME,MP)
1165
1166 C---Charged current
1167
1168         IF (CHARGD) THEN
1169
1170           ML=RMASS(IDHW(1)+1)
1171
1172           YMAX = ONE - TWO*ML*MP / SMA
1173
1174           YMAX = MIN(YMAX,YBMAX)
1175
1176           MREMIN=MP
1177
1178           IF (LEP.EQ.ONE) THEN
1179
1180             MF1=RMASS(IFLAVD)
1181
1182             MF2=RMASS(IFLAVU)
1183
1184           ELSE
1185
1186             MF1=RMASS(IFLAVU)
1187
1188             MF2=RMASS(IFLAVD)
1189
1190           ENDIF
1191
1192           SHMIN = MF1**2+MF2**2 + TWO * PTMIN**2 +
1193
1194      +            TWO * SQRT(PTMIN**2+MF1**2) * SQRT(PTMIN**2+MF2**2)
1195
1196           EMLMIN=MAX(EMMIN,SQRT(SHMIN))
1197
1198           EMLMAX=MIN(EMMAX,RS-ML-MREMIN)
1199
1200           DEBUG=1
1201
1202           IF (EMLMIN.GT.EMLMAX) GOTO 888
1203
1204           WMIN=EMLMIN+MREMIN
1205
1206           PLMAX=HWUPCM(RS,ML,WMIN)
1207
1208           YMIN = ONE-TWO*(SQRT(PINC**2+MP**2)*SQRT(PLMAX**2+ML**2)+
1209
1210      +                    PINC*PLMAX)/SMA
1211
1212           YMIN = MAX(YMIN,YBMIN)
1213
1214           DEBUG=2
1215
1216           IF (YMIN.GT.YMAX) GOTO 888
1217
1218         ELSE
1219
1220 C---Neutral current
1221
1222           ML = ME
1223
1224           YMAX = ONE - TWO*ML*MP / SMA
1225
1226           YMAX = MIN(YMAX,YBMAX)
1227
1228           DO I=1,18
1229
1230             YMIF(I)=ZERO
1231
1232             EMMIF(I)=ZERO
1233
1234             EMMAF(I)=ZERO
1235
1236             WMIF(I)=ZERO
1237
1238             IF (I.LE.8) THEN
1239
1240 C---Boson-Gluon Fusion (also J/Psi) and QCD Compton with struck u or d
1241
1242               MREMIF(I)=MP
1243
1244               IF (I.LE.6) THEN
1245
1246                 MFIN1(I)=RMASS(I)
1247
1248                 MFIN2(I)=RMASS(I+6)
1249
1250               ELSE
1251
1252                 MFIN1(I)=RMASS(I-6)
1253
1254                 MFIN2(I)=ZERO
1255
1256               ENDIF
1257
1258             ELSE
1259
1260 C---QCD Compton with struck non-valence parton
1261
1262               MREMIF(I)=MP+RMASS(I-6)
1263
1264               MFIN1(I)=RMASS(I-6)
1265
1266               MFIN2(I)=ZERO
1267
1268             ENDIF
1269
1270           ENDDO
1271
1272           IF (IFL.EQ.164) THEN
1273
1274 C---J/Psi
1275
1276             MFIN1(7)=RMASS(164)
1277
1278             MFIN2(7)=ZERO
1279
1280           ENDIF
1281
1282 C---y boundaries for different flavours and processes
1283
1284           DO 100 I=IMIN,IMAX
1285
1286             IF (INCLUD(I)) THEN
1287
1288               MF1=MFIN1(I)
1289
1290               MF2=MFIN2(I)
1291
1292               MREMIN=MREMIF(I)
1293
1294               SHMIN = MF1**2+MF2**2 + TWO * PTMIN**2 +
1295
1296      +              TWO * SQRT(PTMIN**2+MF1**2) * SQRT(PTMIN**2+MF2**2)
1297
1298               EMMIF(I) = MAX(EMMIN,SQRT(SHMIN))
1299
1300               EMMAF(I) = MIN(EMMAX,RS-ML-MREMIN)
1301
1302               IF (EMMIF(I).GT.EMMAF(I)) THEN
1303
1304                 INCLUD(I)=.FALSE.
1305
1306                 CALL HWWARN('HWHBRN',3,*999)
1307
1308                 GOTO 100
1309
1310               ENDIF
1311
1312               WMIF(I) = EMMIF(I)+MREMIF(I)
1313
1314               WMIN = WMIF(I)
1315
1316               PLMAX = HWUPCM(RS,ML,WMIN)
1317
1318               YMIF(I)=ONE-TWO*(SQRT(PINC**2+MP**2)*SQRT(PLMAX**2+ML**2)+
1319
1320      +                         PINC*PLMAX)/SMA
1321
1322               IF (YMIF(I).GT.YMAX) THEN
1323
1324                 INCLUD(I)=.FALSE.
1325
1326                 CALL HWWARN('HWHBRN',4,*999)
1327
1328                 GOTO 100
1329
1330               ENDIF
1331
1332             ENDIF
1333
1334  100      CONTINUE
1335
1336 C---considering the largest boundaries
1337
1338           EMLMIN=EMMIF(IMIN)
1339
1340           EMLMAX=EMMAF(IMIN)
1341
1342           IF (IPROO.EQ.3) THEN
1343
1344             EMLMIN=MIN(EMMIF(IMIN),EMMIF(IMIN+6))
1345
1346             EMLMAX=MAX(EMMAF(IMIN),EMMAF(IMIN+6))
1347
1348           ENDIF
1349
1350           DEBUG=3
1351
1352           IF (EMLMIN.GT.EMLMAX) GOTO 888
1353
1354           YMIN=YMIF(IMIN)
1355
1356           IF (IPROO.EQ.3) YMIN=MIN(YMIF(IMIN),YMIF(IMIN+6))
1357
1358           YMIN = MAX(YMIN,YBMIN)
1359
1360           DEBUG=4
1361
1362           IF (YMIN.GT.YMAX) GOTO 888
1363
1364           WMIN = WMIF(IMIN)
1365
1366           MREMIN = MREMIF(IMIN)
1367
1368           MF1=MFIN1(IMIN)
1369
1370           MF2=MFIN2(IMIN)
1371
1372           IF (IPROO.EQ.3) THEN
1373
1374             WMIN = MIN(WMIF(IMIN),WMIF(IMIN+6))
1375
1376             MREMIN = MIN(MREMIF(IMIN),MREMIF(IMIN+6))
1377
1378           ENDIF
1379
1380         ENDIF
1381
1382       ENDIF
1383
1384 C---Random generation in largest phase space
1385
1386       Y=ZERO
1387
1388       Q2=ZERO
1389
1390       SHAT=ZERO
1391
1392       Z=ZERO
1393
1394       PHI=ZERO
1395
1396       AJACOB=ZERO
1397
1398 C---y generation
1399
1400       IF (.NOT.CHARGD) THEN
1401
1402         IF (IFL.LE.5.OR.(IFL.GE.7.AND.IFL.LE.18)) THEN
1403
1404           SRY0 = SQRT(YMIN)
1405
1406           SRY1 = SQRT(YMAX)
1407
1408           SRY = HWRUNI(0,SRY0,SRY1)
1409
1410           Y = SRY**2
1411
1412           YJAC = TWO*SRY*(SRY1-SRY0)
1413
1414         ELSEIF (IFL.EQ.6) THEN
1415
1416           Y = SQRT(HWRUNI(0,YMIN**2,YMAX**2))
1417
1418           YJAC = HALF * (YMAX**2-YMIN**2) / Y
1419
1420         ELSEIF (IFL.EQ.164) THEN
1421
1422 C---in J/psi photoproduction Y and Q2 are given by the Equivalent Photon
1423
1424 C   Approximation
1425
1426    10     NTRY=0
1427
1428    20     NTRY=NTRY+1
1429
1430           IF (NTRY.GT.NETRY) CALL HWWARN('HWHBRN',50,*10)
1431
1432           Y = (YMIN/YMAX)**HWR()*YMAX
1433
1434           IF (ONE+(ONE-Y)**2.LT.TWO*HWR()) GOTO 20
1435
1436           YJAC=(TWO*LOG(YMAX/YMIN)-TWO*(YMAX-YMIN)
1437
1438      &                            +HALF*(YMAX**2-YMIN**2))
1439
1440         ENDIF
1441
1442       ELSE
1443
1444         IF (IPRO.EQ.5) THEN
1445
1446           Y = EXP(HWRUNI(0,LOG(YMIN),LOG(YMAX)))
1447
1448           YJAC = Y * LOG(YMAX/YMIN)
1449
1450         ELSE
1451
1452           Y = HWRUNI(0,YMIN,YMAX)
1453
1454           YJAC = YMAX - YMIN
1455
1456         ENDIF
1457
1458       ENDIF
1459
1460 C---Q**2 generation
1461
1462       Q2INF = ME**2*Y**2 / (ONE-Y)
1463
1464       Q2SUP = MP**2 + SMA*Y - WMIN**2
1465
1466       IF (IFL.EQ.164) THEN
1467
1468         Q2INF = MAX(Q2INF,Q2WWMN)
1469
1470         Q2SUP = MIN(Q2SUP,Q2WWMX)
1471
1472       ELSE
1473
1474         Q2INF = MAX(Q2INF,Q2MIN)
1475
1476         Q2SUP = MIN(Q2SUP,Q2MAX)
1477
1478       ENDIF
1479
1480       DEBUG=5
1481
1482       IF (Q2INF .GT. Q2SUP) GOTO 888
1483
1484 C
1485
1486       IF (.NOT.CHARGD) THEN
1487
1488         IF (IFL.EQ.164) THEN
1489
1490           Q2 = EXP(HWRUNI(0,LOG(Q2INF),LOG(Q2SUP)))
1491
1492           Q2JAC = LOG(Q2SUP/Q2INF)
1493
1494         ELSEIF (Q2INF.LT.RMASS(4)**2) THEN
1495
1496           Q2 = EXP(HWRUNI(0,LOG(Q2INF),LOG(Q2SUP)))
1497
1498           Q2JAC = Q2 * LOG(Q2SUP/Q2INF)
1499
1500         ELSE
1501
1502           Q2 = Q2INF*Q2SUP/HWRUNI(0,Q2INF,Q2SUP)
1503
1504           Q2JAC = Q2**2 * (Q2SUP-Q2INF)/(Q2SUP*Q2INF)
1505
1506         ENDIF
1507
1508       ELSE
1509
1510         EMW2=EMW**2
1511
1512         Q2=(Q2INF+EMW2)*(Q2SUP+EMW2)/(HWRUNI(0,Q2INF,Q2SUP)+EMW2)-EMW2
1513
1514         Q2JAC=(Q2+EMW2)**2*(Q2SUP-Q2INF)/((Q2SUP+EMW2)*(Q2INF+EMW2))
1515
1516       ENDIF
1517
1518       W2 = MP**2 + SMA*Y - Q2
1519
1520 C---s_hat generation
1521
1522       SHINF = EMLMIN **2
1523
1524       SHSUP = (MIN(SQRT(W2)-MREMIN,EMLMAX))**2
1525
1526       DEBUG=6
1527
1528       IF (SHINF .GT. SHSUP) GOTO 888
1529
1530 C
1531
1532       IF (IPRO.EQ.91) THEN
1533
1534         IF (.NOT.CHARGD) THEN
1535
1536           SHAT = SHINF*SHSUP/HWRUNI(0,SHINF,SHSUP)
1537
1538           SHJAC = SHAT**2 * (SHSUP-SHINF)/(SHSUP*SHINF)
1539
1540         ELSE
1541
1542           SHAT = EXP(HWRUNI(0,LOG(SHINF),LOG(SHSUP)))
1543
1544           SHJAC = SHAT*(LOG(SHSUP/SHINF))
1545
1546         ENDIF
1547
1548       ELSE
1549
1550         EMW2=EMW**2
1551
1552         IF (SHINF.GT.EMW2+10*GAMW*EMW) THEN
1553
1554           SHAT = SHINF*SHSUP/HWRUNI(0,SHINF,SHSUP)
1555
1556           SHJAC = SHAT**2 * (SHSUP-SHINF)/(SHSUP*SHINF)
1557
1558         ELSEIF (SHSUP.LT.EMW2-10*EMW*GAMW) THEN
1559
1560           SHAT = HWRUNI(0,SHINF,SHSUP)
1561
1562           SHJAC = SHSUP-SHINF
1563
1564         ELSE
1565
1566           TMIN=ATAN((SHINF-EMW2)/(GAMW*EMW))
1567
1568           TMAX=ATAN((SHSUP-EMW2)/(GAMW*EMW))
1569
1570           SHAT = GAMW*EMW*TAN(HWRUNI(0,TMIN,TMAX))+EMW2
1571
1572           SHJAC=((SHAT-EMW2)**2+(GAMW*EMW)**2)/(GAMW*EMW)*(TMAX-TMIN)
1573
1574         ENDIF
1575
1576       ENDIF
1577
1578       DETDSH = ONE/SMA/Y
1579
1580       SHJAC=SHJAC*DETDSH
1581
1582       RSHAT = SQRT (SHAT)
1583
1584 C--- z generation
1585
1586       ZMIN = 10E10
1587
1588       ZMAX = -ONE
1589
1590       IF (.NOT.CHARGD) THEN
1591
1592         DO I=1,18
1593
1594           Q1CM(I) = ZERO
1595
1596           ZMIF(I) = ZERO
1597
1598           ZMAF(I) = ZERO
1599
1600         ENDDO
1601
1602         DO 150 I=IMIN,IMAX
1603
1604           IF (INCLUD(I)) THEN
1605
1606             Q1CM(I) = HWUPCM( RSHAT, MFIN1(I), MFIN2(I) )
1607
1608             IF (Q1CM(I) .LT. PTMIN) THEN
1609
1610               ZMAF(I)=-ONE
1611
1612               GOTO 150
1613
1614             ENDIF
1615
1616             CTHLIM = SQRT(ONE - (PTMIN / Q1CM(I))**2)
1617
1618             GAMMA2 = SHAT + MFIN1(I)**2 - MFIN2(I)**2
1619
1620             LAMBDA = (SHAT-MFIN1(I)**2-MFIN2(I)**2)**2 -
1621
1622      +                4.D0*MFIN1(I)**2*MFIN2(I)**2
1623
1624             ZMIF(I) = (GAMMA2 - SQRT(LAMBDA)*CTHLIM)/TWO/SHAT
1625
1626             ZMIF(I) = MAX(ZMIF(I),ZERO)
1627
1628             ZMAF(I) = (GAMMA2 + SQRT(LAMBDA)*CTHLIM)/TWO/SHAT
1629
1630             ZMAF(I) = MIN(ZMAF(I),ONE)
1631
1632             ZMIN = MIN( ZMIN, ZMIF(I) )
1633
1634             ZMAX = MAX( ZMAX, ZMAF(I) )
1635
1636           ENDIF
1637
1638  150    CONTINUE
1639
1640         IF (IFL.EQ.164) ZMAX=MIN(ZMAX,ZJMAX)
1641
1642       ELSE
1643
1644         Q1 = HWUPCM(RSHAT,MF1,MF2)
1645
1646         DEBUG=7
1647
1648         IF (Q1.LT.PTMIN) GOTO 888
1649
1650         CTHLIM = SQRT(ONE-(PTMIN/Q1)**2)
1651
1652         GAMMA2 = SHAT+MF1**2-MF2**2
1653
1654         LAMBDA = (SHAT-MF1**2-MF2**2)**2-4.D0*MF1**2*MF2**2
1655
1656         ZMIN = (GAMMA2-SQRT(LAMBDA)*CTHLIM)/TWO/SHAT
1657
1658         ZMIN = MAX(ZMIN,1D-6)
1659
1660         ZMAX = (GAMMA2+SQRT(LAMBDA)*CTHLIM)/TWO/SHAT
1661
1662         ZMAX = MIN(ZMAX,ONE-1D-6)
1663
1664       ENDIF
1665
1666       DEBUG=8
1667
1668       IF (ZMIN .GT. ZMAX) GOTO 888
1669
1670       ZLMIN = LOG(ZMIN/(ONE-ZMIN))
1671
1672       ZINT = LOG(ZMAX/(ONE-ZMAX)) - LOG(ZMIN/(ONE-ZMIN))
1673
1674       ZL = ZLMIN+HWR()*ZINT
1675
1676       Z = EXP(ZL)/(ONE+EXP(ZL))
1677
1678       ZJAC = Z*(ONE-Z)*ZINT
1679
1680 C
1681
1682       DEBUG=9
1683
1684       IF ((Y.LT.YMIN.OR.Y.GT.YMAX).OR.(Q2.LT.Q2INF.OR.Q2.GT.Q2SUP).OR.
1685
1686      +   (SHAT.LT.SHINF.OR.SHAT.GT.SHSUP).OR.(Z.LT.ZMIN.OR.Z.GT.ZMAX))
1687
1688      +     GOTO 888
1689
1690 C---Phi generation
1691
1692       PHI = HWRUNI(0,ZERO,2*PIFAC)
1693
1694       PHIJAC = 2 * PIFAC
1695
1696       IF (IFL.EQ.164) PHIJAC=ONE
1697
1698 C
1699
1700       AJACOB = YJAC * Q2JAC * SHJAC * ZJAC * PHIJAC
1701
1702 C
1703
1704       IF (IQK.NE.0.OR.IPRO.EQ.5) GOTO 999
1705
1706 C---contributing subprocesses: filling of logical vector INSIDE
1707
1708       DO I=1,18
1709
1710         INSIDE(I)=.FALSE.
1711
1712         Q2MAF(I)=ZERO
1713
1714         EMMAWF(I)=ZERO
1715
1716       ENDDO
1717
1718       DO 200 I=IMIN,IMAX
1719
1720       IF (INCLUD(I)) THEN
1721
1722       IF ( Y.LT.YMIF(I) ) GOTO 200
1723
1724 C
1725
1726       Q2MAF(I) = MP**2 + SMA*Y - WMIF(I)**2
1727
1728       Q2MAF(I) = MIN( Q2MAF(I), Q2MAX)
1729
1730       IF (Q2INF .GT. Q2MAF(I)) GOTO 200
1731
1732       IF (Q2.LT.Q2INF .OR. Q2.GT.Q2MAF(I)) GOTO 200
1733
1734 C
1735
1736       EMMAWF(I) = SQRT(W2) - MREMIF(I)
1737
1738       EMMAWF(I) = MIN( EMMAWF(I), EMLMAX )
1739
1740 C
1741
1742       IF (EMMIF(I) .GT. EMMAWF(I)) GOTO 200
1743
1744       IF (SHAT.LT.EMMIF(I)**2.OR.SHAT.GT.EMMAWF(I)**2) GOTO 200
1745
1746 C
1747
1748       IF (ZMIF(I) .GT. ZMAF(I)) GOTO 200
1749
1750       IF (Z.LT.ZMIF(I) .OR. Z.GT.ZMAF(I)) GOTO 200
1751
1752       INSIDE(I)=.TRUE.
1753
1754       ENDIF
1755
1756  200  CONTINUE
1757
1758  999  RETURN
1759
1760  888  EVWGT=ZERO
1761
1762 C---UNCOMMENT THIS LINE TO GET A DEBUGGING WARNING FOR NO PHASE-SPACE
1763
1764 C      CALL HWWARN('HWHBRN',DEBUG,*777)
1765
1766  777  RETURN 1
1767
1768       END