HERWIG/src/hwbfin.f

   1
   2 CDECK  ID>, HWBFIN.
   3
   4 *CMZ :-        -26/04/91  10.18.56  by  Bryan Webber
   5
   6 *-- Author :    Bryan Webber
   7
   8 C-----------------------------------------------------------------------
   9
  10       SUBROUTINE HWBFIN(IHEP)
  11
  12 C-----------------------------------------------------------------------
  13
  14 C     DELETES INTERNAL LINES FROM SHOWER, MAKES COLOUR CONNECTION INDEX
  15
  16 C     AND COPIES INTO /HEPEVT/ IN COLOUR ORDER.
  17
  18 C-----------------------------------------------------------------------
  19
  20       INCLUDE 'HERWIG61.INC'
  21
  22       INTEGER IHEP,ID,IJET,KHEP,IPAR,JPAR,NXPAR,IP,JP
  23
  24       IF (IERROR.NE.0) RETURN
  25
  26 C---SAVE VIRTUAL PARTON DATA
  27
  28       NHEP=NHEP+1
  29
  30       IF(NHEP.GT.NMXHEP) CALL HWWARN('HWBFIN',100,*999)
  31
  32       ID=IDPAR(2)
  33
  34       IDHW(NHEP)=ID
  35
  36       IDHEP(NHEP)=IDPDG(ID)
  37
  38       ISTHEP(NHEP)=ISTHEP(IHEP)+20
  39
  40       JMOHEP(1,NHEP)=IHEP
  41
  42       JMOHEP(2,NHEP)=JMOHEP(1,IHEP)
  43
  44       JDAHEP(1,IHEP)=NHEP
  45
  46       JDAHEP(1,NHEP)=0
  47
  48       JDAHEP(2,NHEP)=0
  49
  50       CALL HWVEQU(5,PPAR(1,2),PHEP(1,NHEP))
  51
  52       CALL HWVEQU(4,VPAR(1,2),VHEP(1,NHEP))
  53
  54 C---FINISHED FOR SPECTATOR OR NON-PARTON JETS
  55
  56       IF (ISTHEP(NHEP).GT.136) RETURN
  57
  58       IF (ID.GT.13.AND.ID.LT.209 .AND. ID.NE.59) RETURN
  59
  60       IF (ID.GT.220.AND.ABS(IDPDG(ID)).LT.1000000) RETURN
  61
  62       IF (ID.GT.424.AND.ID.NE.449) RETURN
  63
  64       IF (.NOT.TMPAR(2).AND.ID.EQ.59) RETURN
  65
  66       IDHEP(NHEP)=94
  67
  68       IJET=NHEP
  69
  70       IF (NPAR.GT.2) THEN
  71
  72 C---SAVE CONE DATA
  73
  74         NHEP=NHEP+1
  75
  76         IF(NHEP.GT.NMXHEP) CALL HWWARN('HWBFIN',101,*999)
  77
  78         IDHW(NHEP)=IDPAR(1)
  79
  80         IDHEP(NHEP)=0
  81
  82         ISTHEP(NHEP)=100
  83
  84         JMOHEP(1,NHEP)=IHEP
  85
  86         JMOHEP(2,NHEP)=JCOPAR(1,1)
  87
  88         JDAHEP(1,NHEP)=0
  89
  90         JDAHEP(2,NHEP)=0
  91
  92         CALL HWVEQU(5,PPAR,PHEP(1,NHEP))
  93
  94         CALL HWVEQU(4,VPAR(1,2),VHEP(1,NHEP))
  95
  96       ENDIF
  97
  98       KHEP=NHEP
  99
 100 C---START WITH ANTICOLOUR DAUGHTER OF HARDEST PARTON
 101
 102       IPAR=2
 103
 104       JPAR=JCOPAR(4,IPAR)
 105
 106       NXPAR=NPAR/2
 107
 108       DO 20 IP=1,NXPAR
 109
 110       DO 10 JP=1,NXPAR
 111
 112       IF (JPAR.EQ.0) GOTO 15
 113
 114       IF (JCOPAR(2,JPAR).EQ.IPAR) THEN
 115
 116         IPAR=JPAR
 117
 118         JPAR=JCOPAR(4,IPAR)
 119
 120       ELSE
 121
 122         IPAR=JPAR
 123
 124         JPAR=JCOPAR(1,IPAR)
 125
 126       ENDIF
 127
 128    10 CONTINUE
 129
 130 C---COULDN'T FIND COLOUR PARTNER
 131
 132       CALL HWWARN('HWBFIN',1,*999)
 133
 134    15 JPAR=JCOPAR(1,IPAR)
 135
 136       KHEP=KHEP+1
 137
 138       IF(KHEP.GT.NMXHEP) CALL HWWARN('HWBFIN',102,*999)
 139
 140       ID=IDPAR(IPAR)
 141
 142       IF (TMPAR(IPAR)) THEN
 143
 144         IF (ID.LT.14) THEN
 145
 146           ISTHEP(KHEP)=139
 147
 148         ELSEIF (ID.EQ.59) THEN
 149
 150           ISTHEP(KHEP)=139
 151
 152         ELSEIF (ID.LT.109) THEN
 153
 154           ISTHEP(KHEP)=130
 155
 156         ELSEIF (ID.LT.120) THEN
 157
 158           ISTHEP(KHEP)=139
 159
 160         ELSEIF (ABS(IDPDG(ID)).LT.1000000) THEN
 161
 162           ISTHEP(KHEP)=130
 163
 164         ELSEIF (ID.LT.425) THEN
 165
 166           ISTHEP(KHEP)=139
 167
 168         ELSEIF (ID.EQ.449) THEN
 169
 170           ISTHEP(KHEP)=139
 171
 172         ELSE
 173
 174           ISTHEP(KHEP)=130
 175
 176         ENDIF
 177
 178       ELSE
 179
 180         ISTHEP(KHEP)=ISTHEP(IHEP)+24
 181
 182       ENDIF
 183
 184       IDHW(KHEP)=ID
 185
 186       IDHEP(KHEP)=IDPDG(ID)
 187
 188       CALL HWVEQU(5,PPAR(1,IPAR),PHEP(1,KHEP))
 189
 190       CALL HWVEQU(4,VPAR(1,IPAR),VHEP(1,KHEP))
 191
 192       JMOHEP(1,KHEP)=IJET
 193
 194       JMOHEP(2,KHEP)=KHEP+1
 195
 196       JDAHEP(1,KHEP)=0
 197
 198       JDAHEP(2,KHEP)=KHEP-1
 199
 200    20 CONTINUE
 201
 202       JMOHEP(2,KHEP)=0
 203
 204       JDAHEP(2,NHEP+1)=0
 205
 206       JDAHEP(1,IJET)=NHEP+1
 207
 208       JDAHEP(2,IJET)=KHEP
 209
 210       NHEP=KHEP
 211
 212   999 END
 213
 214 CDECK  ID>, HWBGEN.
 215
 216 *CMZ :-        -14/10/99  18.04.56  by  Mike Seymour
 217
 218 *-- Author :    Bryan Webber
 219
 220 C-----------------------------------------------------------------------
 221
 222       SUBROUTINE HWBGEN
 223
 224 C-----------------------------------------------------------------------
 225
 226 C     BRANCHING GENERATOR WITH INTERFERING GLUONS
 227
 228 C     HWBGEN EVOLVES QCD JETS ACCORDING TO THE METHOD OF
 229
 230 C     G.MARCHESINI & B.R.WEBBER, NUCL. PHYS. B238(1984)1
 231
 232 C-----------------------------------------------------------------------
 233
 234       INCLUDE 'HERWIG61.INC'
 235
 236       DOUBLE PRECISION HWULDO,HWRGAU,EINHEP,ERTXI,RTXI,XF
 237
 238       INTEGER NTRY,LASHEP,IHEP,NRHEP,ID,IST,JHEP,KPAR,I,J,IRHEP(NMXJET),
 239
 240      & IRST(NMXJET)
 241
 242       LOGICAL HWRLOG
 243
 244       EXTERNAL HWULDO,HWRGAU
 245
 246       IF (IERROR.NE.0) RETURN
 247
 248       IF (IPRO.EQ.80) RETURN
 249
 250 C---CHECK THAT EMSCA IS SET
 251
 252       IF (EMSCA.LE.ZERO) CALL HWWARN('HWBGEN',200,*999)
 253
 254       IF (HARDME) THEN
 255
 256 C---FORCE A BRANCH INTO THE `DEAD ZONE' IN E+E-
 257
 258         IF (IPROC/10.EQ.10) CALL HWBDED(1)
 259
 260 C---FORCE A BRANCH INTO THE `DEAD ZONE' IN DIS
 261
 262         IF (IPRO.EQ.90) CALL HWBDIS(1)
 263
 264 C---FORCE A BRANCH INTO THE `DEAD ZONE' IN DRELL-YAN PROCESSES
 265
 266         IF (IPRO.EQ.13.OR.IPRO.EQ.14) CALL HWBDYP(1)
 267
 268 C---FORCE A BRANCH INTO THE `DEAD ZONE' IN TOP DECAYS
 269
 270         CALL HWBTOP
 271
 272       ENDIF
 273
 274 C---GENERATE INTRINSIC PT ONCE AND FOR ALL
 275
 276       DO 5 JNHAD=1,2
 277
 278         IF (PTRMS.NE.0.) THEN
 279
 280           PTINT(1,JNHAD)=HWRGAU(1,ZERO,PXRMS)
 281
 282           PTINT(2,JNHAD)=HWRGAU(2,ZERO,PXRMS)
 283
 284           PTINT(3,JNHAD)=PTINT(1,JNHAD)**2+PTINT(2,JNHAD)**2
 285
 286         ELSE
 287
 288           CALL HWVZRO(3,PTINT(1,JNHAD))
 289
 290         ENDIF
 291
 292  5    CONTINUE
 293
 294       NTRY=0
 295
 296       LASHEP=NHEP
 297
 298  10   NTRY=NTRY+1
 299
 300       IF (NTRY.GT.NETRY) CALL HWWARN('HWBGEN',ISLENT*100,*999)
 301
 302       NRHEP=0
 303
 304       NHEP=LASHEP
 305
 306       FROST=.FALSE.
 307
 308       DO 100 IHEP=1,LASHEP
 309
 310       IST=ISTHEP(IHEP)
 311
 312       IF (IST.GE.111.AND.IST.LE.115) THEN
 313
 314        NRHEP=NRHEP+1
 315
 316        IRHEP(NRHEP)=IHEP
 317
 318        IRST(NRHEP)=IST
 319
 320        ID=IDHW(IHEP)
 321
 322        IF (IST.NE.115) THEN
 323
 324 C---FOUND A PARTON TO EVOLVE
 325
 326         NEVPAR=IHEP
 327
 328         NPAR=2
 329
 330         IDPAR(1)=17
 331
 332         IDPAR(2)=ID
 333
 334         TMPAR(1)=.TRUE.
 335
 336         PPAR(2,1)=0.
 337
 338         PPAR(4,1)=1.
 339
 340         DO 15 J=1,2
 341
 342         DO 15 I=1,2
 343
 344         JMOPAR(I,J)=0
 345
 346  15     JCOPAR(I,J)=0
 347
 348 C---SET UP EVOLUTION SCALE AND FRAME
 349
 350         JHEP=JMOHEP(2,IHEP)
 351
 352         IF (ID.EQ.13) THEN
 353
 354           IF (HWRLOG(HALF)) JHEP=JDAHEP(2,IHEP)
 355
 356         ELSEIF (IST.GT.112) THEN
 357
 358           IF ((ID.GT.6.AND.ID.LT.13).OR.
 359
 360      &        (ID.GT.214.AND.ID.LT.221)) JHEP=JDAHEP(2,IHEP)
 361
 362         ELSE
 363
 364           IF (ID.LT.7.OR.(ID.GT.208.AND.ID.LT.215)) JHEP=JDAHEP(2,IHEP)
 365
 366         ENDIF
 367
 368         IF (JHEP.LE.0.OR.JHEP.GT.NHEP) THEN
 369
 370           CALL HWWARN('HWBGEN',1,*999)
 371
 372           JHEP=IHEP
 373
 374         ENDIF
 375
 376         JCOPAR(1,1)=JHEP
 377
 378         EINHEP=PHEP(4,IHEP)
 379
 380         ERTXI=HWULDO(PHEP(1,IHEP),PHEP(1,JHEP))
 381
 382         IF (ERTXI.LT.ZERO) ERTXI=0.
 383
 384         IF (IST.LE.112.AND.IHEP.EQ.JHEP) ERTXI=0.
 385
 386         IF (ISTHEP(JHEP).EQ.155) THEN
 387
 388           ERTXI=ERTXI/PHEP(5,JHEP)
 389
 390           RTXI=1.
 391
 392         ELSE
 393
 394           ERTXI=SQRT(ERTXI)
 395
 396           RTXI=ERTXI/EINHEP
 397
 398         ENDIF
 399
 400         IF (RTXI.EQ.ZERO) THEN
 401
 402           XF=1.
 403
 404           PPAR(1,1)=0.
 405
 406           PPAR(3,1)=1.
 407
 408           PPAR(1,2)=EINHEP
 409
 410           PPAR(2,2)=0.
 411
 412           PPAR(4,2)=EINHEP
 413
 414         ELSE
 415
 416           XF=1./RTXI
 417
 418           PPAR(1,1)=1.
 419
 420           PPAR(3,1)=0.
 421
 422           PPAR(1,2)=ERTXI
 423
 424           PPAR(2,2)=1.
 425
 426           PPAR(4,2)=ERTXI
 427
 428         ENDIF
 429
 430         IF (PPAR(4,2).LT.PHEP(5,IHEP)) PPAR(4,2)=PHEP(5,IHEP)
 431
 432 C---STORE MASS
 433
 434         PPAR(5,2)=PHEP(5,IHEP)
 435
 436         CALL HWVZRO(4,VPAR(1,1))
 437
 438         CALL HWVZRO(4,VPAR(1,2))
 439
 440         IF (IST.GT.112) THEN
 441
 442           TMPAR(2)=.TRUE.
 443
 444           INHAD=0
 445
 446           JNHAD=0
 447
 448           XFACT=0.
 449
 450         ELSE
 451
 452           TMPAR(2)=.FALSE.
 453
 454           JNHAD=IST-110
 455
 456           INHAD=JNHAD
 457
 458           IF (JDAHEP(1,JNHAD).NE.0) INHAD=JDAHEP(1,JNHAD)
 459
 460           XFACT=XF/PHEP(4,INHAD)
 461
 462           ANOMSC(1,JNHAD)=ZERO
 463
 464           ANOMSC(2,JNHAD)=ZERO
 465
 466         ENDIF
 467
 468 C---FOR QUARKS IN A COLOUR SINGLET, ALLOW SOFT MATRIX-ELEMENT CORRECTION
 469
 470         HARDST=PPAR(4,2)
 471
 472         IF (SOFTME.AND.IDHW(IHEP).LT.13.AND.
 473
 474      $       ((JMOHEP(2,JHEP).EQ.IHEP.AND.JDAHEP(2,JHEP).EQ.IHEP).OR.
 475
 476      $       ISTHEP(JHEP).EQ.155)) HARDST=0
 477
 478 C---CREATE BRANCHES AND COMPUTE ENERGIES
 479
 480         DO 20 KPAR=2,NMXPAR
 481
 482         IF (TMPAR(KPAR)) THEN
 483
 484           CALL HWBRAN(KPAR)
 485
 486         ELSE
 487
 488           CALL HWSBRN(KPAR)
 489
 490         ENDIF
 491
 492         IF (IERROR.NE.0) RETURN
 493
 494         IF (FROST) GOTO 100
 495
 496         IF (KPAR.EQ.NPAR) GOTO 30
 497
 498  20     CONTINUE
 499
 500 C---COMPUTE MASSES AND 3-MOMENTA
 501
 502  30     CONTINUE
 503
 504         CALL HWBMAS
 505
 506         IF (AZSPIN) CALL HWBSPN
 507
 508         IF (TMPAR(2)) THEN
 509
 510            CALL HWBTIM(2,1)
 511
 512         ELSE
 513
 514            CALL HWBSPA
 515
 516         ENDIF
 517
 518 C---ENTER PARTON JET IN /HEPEVT/
 519
 520         CALL HWBFIN(IHEP)
 521
 522        ELSE
 523
 524 C---COPY SPECTATOR
 525
 526         NHEP=NHEP+1
 527
 528         IF (ID.GT.120.AND.ID.LT.133 .OR. ID.GE.198.AND.ID.LE.201) THEN
 529
 530           ISTHEP(NHEP)=190
 531
 532         ELSE
 533
 534           ISTHEP(NHEP)=152
 535
 536         ENDIF
 537
 538         IDHW(NHEP)=ID
 539
 540         IDHEP(NHEP)=IDPDG(ID)
 541
 542         JMOHEP(1,NHEP)=IHEP
 543
 544         JMOHEP(2,NHEP)=0
 545
 546         JDAHEP(2,NHEP)=0
 547
 548         JDAHEP(1,IHEP)=NHEP
 549
 550         CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP))
 551
 552        ENDIF
 553
 554        ISTHEP(IHEP)=ISTHEP(IHEP)+10
 555
 556       ENDIF
 557
 558  100  CONTINUE
 559
 560       IF (.NOT.FROST) THEN
 561
 562 C---COMBINE JETS
 563
 564         ISTAT=20
 565
 566         CALL HWBJCO
 567
 568       ENDIF
 569
 570       IF (.NOT.FROST) THEN
 571
 572 C---ATTACH SPECTATORS
 573
 574         ISTAT=30
 575
 576         CALL HWSSPC
 577
 578       ENDIF
 579
 580       IF (FROST) THEN
 581
 582 C---BAD JET: RESTORE PARTONS AND RE-EVOLVE
 583
 584          DO 120 I=1,NRHEP
 585
 586  120     ISTHEP(IRHEP(I))=IRST(I)
 587
 588          GOTO 10
 589
 590       ENDIF
 591
 592 C---CONNECT COLOURS
 593
 594       CALL HWBCON
 595
 596       ISTAT=40
 597
 598       LASHEP=NHEP
 599
 600       IF (HARDME) THEN
 601
 602 C---CLEAN UP IF THERE WAS A BRANCH IN THE `DEAD ZONE' IN E+E-
 603
 604         IF (IPROC/10.EQ.10) CALL HWBDED(2)
 605
 606 C---CLEAN UP IF THERE WAS A BRANCH IN THE `DEAD ZONE' IN DIS
 607
 608         IF (IPRO.EQ.90) CALL HWBDIS(2)
 609
 610 C---CLEAN UP IF THERE WAS A BRANCH IN THE `DEAD ZONE' IN DRELL-YAN PROC
 611
 612         IF (IPRO.EQ.13.OR.IPRO.EQ.14) CALL HWBDYP(2)
 613
 614       ENDIF
 615
 616 C---IF THE CLEAN-UP OPERATION ADDED ANY PARTONS TO THE EVENT RECORD
 617
 618 C   IT MIGHT NEED RESHOWERING
 619
 620       IF (NHEP.GT.LASHEP) THEN
 621
 622         LASHEP=NHEP
 623
 624         GOTO 10
 625
 626       ENDIF
 627
 628  999  END
 629
 630 CDECK  ID>, HWBJCO.
 631
 632 *CMZ :-        -26/04/91  14.25.31  by  Federico Carminati
 633
 634 *-- Author :    Bryan Webber
 635
 636 C-----------------------------------------------------------------------
 637
 638       SUBROUTINE HWBJCO
 639
 640 C-----------------------------------------------------------------------
 641
 642 C     COMBINES JETS WITH REQUIRED KINEMATICS
 643
 644 C-----------------------------------------------------------------------
 645
 646       INCLUDE 'HERWIG61.INC'
 647
 648       DOUBLE PRECISION HWULDO,EPS,PTX,PTY,PF,PTINF,PTCON,CN,CP,SP,PP0,
 649
 650      & PM0,ET0,DET,ECM,EMJ,EMP,EMS,DMS,ES,DPF,ALF,AL(2),ET(2),PP(2),
 651
 652      & PT(3),PB(5),PC(5),PQ(5),PR(5),PS(5),RR(3,3),RS(3,3),ETC,
 653
 654      & PJ(NMXJET),PM(NMXJET),PBR(5),RBR(3,3),DISP(4)
 655
 656       INTEGER LJET,IJ1,IST,IP,ICM,IP1,IP2,NP,IHEP,MHEP,JP,KP,LP,KHEP,
 657
 658      & JHEP,NE,IJT,IEND(2),IJET(NMXJET),IPAR(NMXJET)
 659
 660       LOGICAL AZCOR,JETRAD,DISPRO,DISLOW
 661
 662       EXTERNAL HWULDO
 663
 664       PARAMETER (EPS=1.D-4)
 665
 666       IF (IERROR.NE.0) RETURN
 667
 668       AZCOR=AZSOFT.OR.AZSPIN
 669
 670 C---FIRST LOOK FOR SPACELIKE JETS
 671
 672       LJET=131
 673
 674   10  IJET(1)=1
 675
 676   20  IJ1=IJET(1)
 677
 678       DO 40 IHEP=IJ1,NHEP
 679
 680       IST=ISTHEP(IHEP)
 681
 682       IF (IST.EQ.137.OR.IST.EQ.138) IST=133
 683
 684       IF (IST.EQ.LJET) THEN
 685
 686 C---FOUND AN UNBOOSTED JET - FIND PARTNERS
 687
 688         IP=JMOHEP(1,IHEP)
 689
 690         ICM=JMOHEP(1,IP)
 691
 692         DISPRO=IPRO/10.EQ.9.AND.IDHW(ICM).EQ.15
 693
 694         DISLOW=DISPRO.AND.JDAHEP(1,ICM).EQ.JDAHEP(2,ICM)-1
 695
 696         IF (IST.EQ.131) THEN
 697
 698           IP1=JMOHEP(1,ICM)
 699
 700           IP2=JMOHEP(2,ICM)
 701
 702         ELSE
 703
 704           IP1=JDAHEP(1,ICM)
 705
 706           IP2=JDAHEP(2,ICM)
 707
 708         ENDIF
 709
 710         IF (IP1.NE.IP) CALL HWWARN('HWBJCO',100,*999)
 711
 712         NP=0
 713
 714         DO 30 JHEP=IP1,IP2
 715
 716         NP=NP+1
 717
 718         IPAR(NP)=JHEP
 719
 720   30    IJET(NP)=JDAHEP(1,JHEP)
 721
 722         GOTO 50
 723
 724       ENDIF
 725
 726   40  CONTINUE
 727
 728 C---NO MORE JETS?
 729
 730       IF (LJET.EQ.131) THEN
 731
 732         LJET=133
 733
 734         GOTO 10
 735
 736       ENDIF
 737
 738       RETURN
 739
 740   50  IF (LJET.EQ.131) THEN
 741
 742 C---SPACELIKE JETS: FIND SPACELIKE PARTONS
 743
 744         IF (NP.NE.2) CALL HWWARN('HWBJCO',103,*999)
 745
 746 C---special for DIS: FIND BOOST AND ROTATION FROM LAB TO BREIT FRAME
 747
 748         IF (DISPRO.AND.BREIT) THEN
 749
 750           IP=2
 751
 752           IF (JDAHEP(1,IP).NE.0) IP=JDAHEP(1,IP)
 753
 754           CALL HWVDIF(4,PHEP(1,JMOHEP(1,ICM)),PHEP(1,JDAHEP(1,ICM)),PB)
 755
 756           CALL HWUMAS(PB)
 757
 758 C---IF Q**2<10**-2, SOMETHING MUST HAVE ALREADY GONE WRONG
 759
 760           IF (PB(5)**2.LT.1.D-2) CALL HWWARN('HWBJCO',102,*999)
 761
 762           CALL HWVSCA(4,PB(5)**2/HWULDO(PHEP(1,IP),PB),PHEP(1,IP),PBR)
 763
 764           CALL HWVSUM(4,PB,PBR,PBR)
 765
 766           CALL HWUMAS(PBR)
 767
 768           CALL HWULOF(PBR,PB,PB)
 769
 770           CALL HWUROT(PB,ONE,ZERO,RBR)
 771
 772         ENDIF
 773
 774         PTX=0.
 775
 776         PTY=0.
 777
 778         PF=1.D0
 779
 780         DO 90 IP=1,2
 781
 782         MHEP=IJET(IP)
 783
 784         IF (JDAHEP(1,MHEP).EQ.0) THEN
 785
 786 C---SPECIAL FOR NON-PARTON JETS
 787
 788           IHEP=MHEP
 789
 790           GOTO 70
 791
 792         ELSE
 793
 794           IST=134+IP
 795
 796           DO 60 IHEP=MHEP,NHEP
 797
 798   60      IF (ISTHEP(IHEP).EQ.IST) GOTO 70
 799
 800 C---COULDN'T FIND SPACELIKE PARTON
 801
 802           CALL HWWARN('HWBJCO',101,*999)
 803
 804         ENDIF
 805
 806   70    CALL HWVSCA(3,PF,PHEP(1,IHEP),PS)
 807
 808         IF (PTINT(3,IP).GT.ZERO) THEN
 809
 810 C---ADD INTRINSIC PT
 811
 812           PT(1)=PTINT(1,IP)
 813
 814           PT(2)=PTINT(2,IP)
 815
 816           PT(3)=0.
 817
 818           CALL HWUROT(PS, ONE,ZERO,RS)
 819
 820           CALL HWUROB(RS,PT,PT)
 821
 822           CALL HWVSUM(3,PS,PT,PS)
 823
 824         ENDIF
 825
 826         JP=IJET(IP)+1
 827
 828         IF (AZCOR.AND.JP.LE.NHEP.AND.IDHW(JP).EQ.17) THEN
 829
 830 C---ALIGN CONE WITH INTERFERING PARTON
 831
 832           CALL HWUROT(PS, ONE,ZERO,RS)
 833
 834           CALL HWUROF(RS,PHEP(1,JP),PR)
 835
 836           PTCON=PR(1)**2+PR(2)**2
 837
 838           KP=JMOHEP(2,JP)
 839
 840           IF (KP.EQ.0) THEN
 841
 842             CALL HWWARN('HWBJCO',1,*999)
 843
 844             PTINF=0.
 845
 846           ELSE
 847
 848             CALL HWVEQU(4,PHEP(1,KP),PB)
 849
 850             IF (DISPRO.AND.BREIT) THEN
 851
 852               CALL HWULOF(PBR,PB,PB)
 853
 854               CALL HWUROF(RBR,PB,PB)
 855
 856             ENDIF
 857
 858             PTINF=PB(1)**2+PB(2)**2
 859
 860             IF (PTINF.LT.EPS) THEN
 861
 862 C---COLLINEAR JETS: ALIGN CONES
 863
 864               KP=JDAHEP(1,KP)+1
 865
 866               IF (ISTHEP(KP).EQ.100.AND.ISTHEP(KP-1)/10.EQ.14) THEN
 867
 868                 CALL HWVEQU(4,PHEP(1,KP),PB)
 869
 870                 IF (DISPRO.AND.BREIT) THEN
 871
 872                   CALL HWULOF(PBR,PB,PB)
 873
 874                   CALL HWUROF(RBR,PB,PB)
 875
 876                 ENDIF
 877
 878                 PTINF=PB(1)**2+PB(2)**2
 879
 880               ELSE
 881
 882                 PTINF=0.
 883
 884               ENDIF
 885
 886             ENDIF
 887
 888           ENDIF
 889
 890           IF (PTCON.NE.ZERO.AND.PTINF.NE.ZERO) THEN
 891
 892             CN=1./SQRT(PTINF*PTCON)
 893
 894             CP=CN*(PR(1)*PB(1)+PR(2)*PB(2))
 895
 896             SP=CN*(PR(1)*PB(2)-PR(2)*PB(1))
 897
 898           ELSE
 899
 900             CALL HWRAZM( ONE,CP,SP)
 901
 902           ENDIF
 903
 904         ELSE
 905
 906           CALL HWRAZM( ONE,CP,SP)
 907
 908         ENDIF
 909
 910 C---ROTATE SO SPACELIKE IS ALONG AXIS (APART FROM INTRINSIC PT)
 911
 912         CALL HWUROT(PS,CP,SP,RS)
 913
 914         IHEP=IJET(IP)
 915
 916         KHEP=JDAHEP(2,IHEP)
 917
 918         IF (KHEP.LT.IHEP) KHEP=IHEP
 919
 920         IEND(IP)=KHEP
 921
 922         DO 80 JHEP=IHEP,KHEP
 923
 924         CALL HWUROF(RS,PHEP(1,JHEP),PHEP(1,JHEP))
 925
 926   80    CALL HWUROF(RS,VHEP(1,JHEP),VHEP(1,JHEP))
 927
 928         PP(IP)=PHEP(4,IHEP)+PF*PHEP(3,IHEP)
 929
 930         ET(IP)=PHEP(1,IHEP)**2+PHEP(2,IHEP)**2-PHEP(5,IHEP)**2
 931
 932 C---REDEFINE HARD CM
 933
 934         PTX=PTX+PHEP(1,IHEP)
 935
 936         PTY=PTY+PHEP(2,IHEP)
 937
 938   90    PF=-PF
 939
 940         PHEP(1,ICM)=PTX
 941
 942         PHEP(2,ICM)=PTY
 943
 944 C---special for DIS: keep lepton momenta fixed
 945
 946         IF (DISPRO) THEN
 947
 948           IP1=JMOHEP(1,ICM)
 949
 950           IP2=JDAHEP(1,ICM)
 951
 952           IJT=IJET(1)
 953
 954 C---IJT will be used to store lepton momentum transfer
 955
 956           CALL HWVDIF(4,PHEP(1,IP1),PHEP(1,IP2),PHEP(1,IJT))
 957
 958           CALL HWUMAS(PHEP(1,IJT))
 959
 960           IF (IDHEP(IP1).EQ.IDHEP(IP2)) THEN
 961
 962             IDHW(IJT)=200
 963
 964           ELSEIF (IDHEP(IP1).LT.IDHEP(IP2)) THEN
 965
 966             IDHW(IJT)=199
 967
 968           ELSE
 969
 970             IDHW(IJT)=198
 971
 972           ENDIF
 973
 974           IDHEP(IJT)=IDPDG(IDHW(IJT))
 975
 976           ISTHEP(IJT)=3
 977
 978 C---calculate boost for struck parton
 979
 980 C   PC is momentum of outgoing parton(s)
 981
 982           IP2=JDAHEP(2,ICM)
 983
 984           IF (.NOT.DISLOW) THEN
 985
 986 C---FOR heavy QQbar PQ and PC are old and new QQbar momenta
 987
 988             CALL HWVSUM(4,PHEP(1,IP2-1),PHEP(1,IP2),PQ)
 989
 990             CALL HWUMAS(PQ)
 991
 992             PC(5)=PQ(5)
 993
 994           ELSE
 995
 996             PC(5)=PHEP(5,JDAHEP(1,IP2))
 997
 998           ENDIF
 999
1000           CALL HWVSUM(2,PHEP(1,IJT),PHEP(1,IJET(2)),PC)
1001
1002           ET(1)=ET(2)
1003
1004 C---USE BREIT FRAME BOSON MOMENTUM IF NECESSARY
1005
1006           IF (BREIT) THEN
1007
1008             ET(2)=ET(1)+PC(5)**2+PHEP(5,IJET(2))**2
1009
1010             PM0=PHEP(5,IJT)
1011
1012             PP0=-PM0
1013
1014           ELSE
1015
1016             ET(2)=PC(1)**2+PC(2)**2+PC(5)**2
1017
1018             PP0=PHEP(4,IJT)+PHEP(3,IJT)
1019
1020             PM0=PHEP(4,IJT)-PHEP(3,IJT)
1021
1022           ENDIF
1023
1024           ET0=(PP0*PM0)+ET(1)-ET(2)
1025
1026           DET=ET0**2-4.*(PP0*PM0)*ET(1)
1027
1028           IF (DET.LT.ZERO) THEN
1029
1030             FROST=.TRUE.
1031
1032             RETURN
1033
1034           ENDIF
1035
1036           ALF=(SQRT(DET)-ET0)/(2.*PP0*PP(2))
1037
1038           PB(1)=0.
1039
1040           PB(2)=0.
1041
1042           PB(5)=2.D0
1043
1044           PB(3)=ALF-(1./ALF)
1045
1046           PB(4)=ALF+(1./ALF)
1047
1048           DO 100 IHEP=IJET(2),IEND(2)
1049
1050           CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP))
1051
1052           CALL HWULF4(PB,VHEP(1,IHEP),VHEP(1,IHEP))
1053
1054 C---BOOST FROM BREIT FRAME IF NECESSARY
1055
1056           IF (BREIT) THEN
1057
1058             CALL HWUROB(RBR,PHEP(1,IHEP),PHEP(1,IHEP))
1059
1060             CALL HWULOB(PBR,PHEP(1,IHEP),PHEP(1,IHEP))
1061
1062             CALL HWUROB(RBR,VHEP(1,IHEP),VHEP(1,IHEP))
1063
1064             CALL HWULB4(PBR,VHEP(1,IHEP),VHEP(1,IHEP))
1065
1066           ENDIF
1067
1068   100     ISTHEP(IHEP)=ISTHEP(IHEP)+10
1069
1070           CALL HWVDIF(4,VHEP(1,IPAR(2)),VHEP(1,IJET(2)),DISP)
1071
1072           DO 110 IHEP=IJET(2),IEND(2)
1073
1074   110     CALL HWVSUM(4,DISP,VHEP(1,IHEP),VHEP(1,IHEP))
1075
1076           IF (IEND(2).GT.IJET(2)+1) ISTHEP(IJET(2)+1)=100
1077
1078           CALL HWVSUM(4,PHEP(1,IJT),PHEP(1,IJET(2)),PC)
1079
1080           CALL HWVSUM(4,PHEP(1,IP1),PHEP(1,IJET(2)),PHEP(1,ICM))
1081
1082           CALL HWUMAS(PHEP(1,ICM))
1083
1084         ELSEIF (IPRO/10.EQ.5) THEN
1085
1086 C Special to preserve photon momentum
1087
1088            ETC=PTX**2+PTY**2+PHEP(5,ICM)**2
1089
1090            ET0=ETC+ET(1)-ET(2)
1091
1092            DET=ET0**2-4.*ETC*ET(1)
1093
1094            IF (DET.LT.ZERO) THEN
1095
1096               FROST=.TRUE.
1097
1098               RETURN
1099
1100            ENDIF
1101
1102            ALF=(SQRT(DET)+ET0-2.*ET(1))/(2.*PP(1)*PP(2))
1103
1104            PB(1)=0.
1105
1106            PB(2)=0.
1107
1108            PB(3)=ALF-1./ALF
1109
1110            PB(4)=ALF+1./ALF
1111
1112            PB(5)=2.
1113
1114            IJT=IJET(2)
1115
1116            DO 120 IHEP=IJT,IEND(2)
1117
1118            CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP))
1119
1120            CALL HWULF4(PB,VHEP(1,IHEP),VHEP(1,IHEP))
1121
1122   120      ISTHEP(IHEP)=ISTHEP(IHEP)+10
1123
1124            CALL HWVDIF(4,VHEP(1,IPAR(2)),VHEP(1,IJT),DISP)
1125
1126            DO 130 IHEP=IJT,IEND(2)
1127
1128   130      CALL HWVSUM(4,DISP,VHEP(1,IHEP),VHEP(1,IHEP))
1129
1130            IF (IEND(2).GT.IJT+1) ISTHEP(IJT+1)=100
1131
1132            ISTHEP(IJET(1))=ISTHEP(IJET(1))+10
1133
1134            CALL HWVSUM(2,PHEP(3,IPAR(1)),PHEP(3,IJT),PHEP(3,ICM))
1135
1136         ELSE
1137
1138           PHEP(4,ICM)=SQRT(PTX**2+PTY**2+PHEP(3,ICM)**2+PHEP(5,ICM)**2)
1139
1140 C---NOW BOOST TO REQUIRED Q**2 AND X-F
1141
1142           PP0=PHEP(4,ICM)+PHEP(3,ICM)
1143
1144           PM0=PHEP(4,ICM)-PHEP(3,ICM)
1145
1146           ET0=(PP0*PM0)+ET(1)-ET(2)
1147
1148           DET=ET0**2-4.*(PP0*PM0)*ET(1)
1149
1150           IF (DET.LT.ZERO) THEN
1151
1152             FROST=.TRUE.
1153
1154             RETURN
1155
1156           ENDIF
1157
1158           DET=SQRT(DET)+ET0
1159
1160           AL(1)= 2.*PM0*PP(1)/DET
1161
1162           AL(2)=(PM0/PP(2))*(1.-2.*ET(1)/DET)
1163
1164           PB(1)=0.
1165
1166           PB(2)=0.
1167
1168           PB(5)=2.
1169
1170           DO 160 IP=1,2
1171
1172           PB(3)=AL(IP)-(1./AL(IP))
1173
1174           PB(4)=AL(IP)+(1./AL(IP))
1175
1176           IJT=IJET(IP)
1177
1178           DO 140 IHEP=IJT,IEND(IP)
1179
1180           CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP))
1181
1182           CALL HWULF4(PB,VHEP(1,IHEP),VHEP(1,IHEP))
1183
1184   140     ISTHEP(IHEP)=ISTHEP(IHEP)+10
1185
1186           CALL HWVDIF(4,VHEP(1,IPAR(IP)),VHEP(1,IJT),DISP)
1187
1188           DO 150 IHEP=IJT,IEND(IP)
1189
1190   150     CALL HWVSUM(4,DISP,VHEP(1,IHEP),VHEP(1,IHEP))
1191
1192           IF (IEND(IP).GT.IJT+1) THEN
1193
1194             ISTHEP(IJT+1)=100
1195
1196           ELSEIF (IEND(IP).EQ.IJT) THEN
1197
1198 C---NON-PARTON JET
1199
1200             ISTHEP(IJT)=3
1201
1202           ENDIF
1203
1204   160     CONTINUE
1205
1206         ENDIF
1207
1208         ISTHEP(ICM)=120
1209
1210       ELSE
1211
1212 C---TIMELIKE JETS
1213
1214 C   special for DIS: preserve outgoing lepton momentum
1215
1216         IF (DISPRO) THEN
1217
1218           CALL HWVEQU(5,PHEP(1,IPAR(1)),PHEP(1,IJET(1)))
1219
1220           ISTHEP(IJET(1))=1
1221
1222           LP=2
1223
1224         ELSE
1225
1226           CALL HWVEQU(5,PHEP(1,ICM),PC)
1227
1228 C--- PQ AND PC ARE OLD AND NEW PARTON CM
1229
1230           CALL HWVSUM(4,PHEP(1,IPAR(1)),PHEP(1,IPAR(2)),PQ)
1231
1232           PQ(5)=PHEP(5,ICM)
1233
1234           IF (NP.GT.2) THEN
1235
1236             DO 170 KP=3,NP
1237
1238   170       CALL HWVSUM(4,PHEP(1,IPAR(KP)),PQ,PQ)
1239
1240           ENDIF
1241
1242           LP=1
1243
1244         ENDIF
1245
1246         IF (.NOT.DISLOW) THEN
1247
1248 C---FIND JET CM MOMENTA
1249
1250           ECM=PQ(5)
1251
1252           EMS=0.
1253
1254           JETRAD=.FALSE.
1255
1256           DO 180 KP=LP,NP
1257
1258           EMJ=PHEP(5,IJET(KP))
1259
1260           EMP=PHEP(5,IPAR(KP))
1261
1262           JETRAD=JETRAD.OR.EMJ.NE.EMP
1263
1264           EMS=EMS+EMJ
1265
1266           PM(KP)= EMJ**2
1267
1268 C---N.B. ROUNDING ERRORS HERE AT HIGH ENERGIES
1269
1270           PJ(KP)=(HWULDO(PHEP(1,IPAR(KP)),PQ)/ECM)**2-EMP**2
1271
1272           IF (PJ(KP).LE.ZERO) CALL HWWARN('HWBJCO',104,*999)
1273
1274   180     CONTINUE
1275
1276           PF=1.
1277
1278           IF (JETRAD) THEN
1279
1280 C---JETS DID RADIATE
1281
1282             IF (EMS.GE.ECM) THEN
1283
1284               FROST=.TRUE.
1285
1286               RETURN
1287
1288             ENDIF
1289
1290             DO 200 NE=1,NETRY
1291
1292             EMS=-ECM
1293
1294             DMS=0.
1295
1296             DO 190 KP=LP,NP
1297
1298             ES=SQRT(PF*PJ(KP)+PM(KP))
1299
1300             EMS=EMS+ES
1301
1302   190       DMS=DMS+PJ(KP)/ES
1303
1304             DPF=2.*EMS/DMS
1305
1306             IF (DPF.GT.PF) DPF=0.9*PF
1307
1308             PF=PF-DPF
1309
1310   200       IF (ABS(DPF).LT.EPS) GOTO 210
1311
1312             CALL HWWARN('HWBJCO',105,*999)
1313
1314           ENDIF
1315
1316   210     CONTINUE
1317
1318         ENDIF
1319
1320 C---BOOST PC AND PQ TO BREIT FRAME IF NECESSARY
1321
1322         IF (DISPRO.AND.BREIT) THEN
1323
1324           CALL HWULOF(PBR,PC,PC)
1325
1326           CALL HWUROF(RBR,PC,PC)
1327
1328           IF (.NOT.DISLOW) THEN
1329
1330             CALL HWULOF(PBR,PQ,PQ)
1331
1332             CALL HWUROF(RBR,PQ,PQ)
1333
1334           ENDIF
1335
1336         ENDIF
1337
1338         DO 230 IP=LP,NP
1339
1340 C---FIND CM ROTATION FOR JET IP
1341
1342         IF (.NOT.DISLOW) THEN
1343
1344           CALL HWVEQU(4,PHEP(1,IPAR(IP)),PR)
1345
1346           IF (DISPRO.AND.BREIT) THEN
1347
1348             CALL HWULOF(PBR,PR,PR)
1349
1350             CALL HWUROF(RBR,PR,PR)
1351
1352           ENDIF
1353
1354           CALL HWULOF(PQ,PR,PR)
1355
1356           CALL HWUROT(PR, ONE,ZERO,RR)
1357
1358           PR(1)=0.
1359
1360           PR(2)=0.
1361
1362           PR(3)=SQRT(PF*PJ(IP))
1363
1364           PR(4)=SQRT(PF*PJ(IP)+PM(IP))
1365
1366           PR(5)=PHEP(5,IJET(IP))
1367
1368           CALL HWUROB(RR,PR,PR)
1369
1370           CALL HWULOB(PC,PR,PR)
1371
1372         ELSE
1373
1374           CALL HWVEQU(5,PC,PR)
1375
1376         ENDIF
1377
1378 C---NOW PR IS LAB/BREIT MOMENTUM OF JET IP
1379
1380         KP=IJET(IP)+1
1381
1382         IF (AZCOR.AND.KP.LE.NHEP.AND.IDHW(KP).EQ.17) THEN
1383
1384 C---ALIGN CONE WITH INTERFERING PARTON
1385
1386           CALL HWUROT(PR, ONE,ZERO,RS)
1387
1388           JP=JMOHEP(2,KP)
1389
1390           IF (JP.EQ.0) THEN
1391
1392             CALL HWWARN('HWBJCO',2,*999)
1393
1394             PTINF=0.
1395
1396           ELSE
1397
1398             CALL HWVEQU(4,PHEP(1,JP),PS)
1399
1400             IF (DISPRO.AND.BREIT) THEN
1401
1402               CALL HWULOF(PBR,PS,PS)
1403
1404               CALL HWUROF(RBR,PS,PS)
1405
1406             ENDIF
1407
1408             CALL HWUROF(RS,PS,PS)
1409
1410             PTINF=PS(1)**2+PS(2)**2
1411
1412             IF (PTINF.LT.EPS) THEN
1413
1414 C---COLLINEAR JETS: ALIGN CONES
1415
1416               JP=JDAHEP(1,JP)+1
1417
1418               IF (ISTHEP(JP).EQ.100.AND.ISTHEP(JP-1)/10.EQ.14) THEN
1419
1420                 CALL HWVEQU(4,PHEP(1,JP),PS)
1421
1422                 IF (DISPRO.AND.BREIT) THEN
1423
1424                   CALL HWULOF(PBR,PS,PS)
1425
1426                   CALL HWUROF(RBR,PS,PS)
1427
1428                 ENDIF
1429
1430                 CALL HWUROF(RS,PS,PS)
1431
1432                 PTINF=PS(1)**2+PS(2)**2
1433
1434               ELSE
1435
1436                 PTINF=0.
1437
1438               ENDIF
1439
1440             ENDIF
1441
1442           ENDIF
1443
1444           CALL HWVEQU(4,PHEP(1,KP),PB)
1445
1446           IF (DISPRO.AND.BREIT) THEN
1447
1448             CALL HWULOF(PBR,PB,PB)
1449
1450             CALL HWUROF(RBR,PB,PB)
1451
1452           ENDIF
1453
1454           PTCON=PB(1)**2+PB(2)**2
1455
1456           IF (PTCON.NE.ZERO.AND.PTINF.NE.ZERO) THEN
1457
1458             CN=1./SQRT(PTINF*PTCON)
1459
1460             CP=CN*(PS(1)*PB(1)+PS(2)*PB(2))
1461
1462             SP=CN*(PS(1)*PB(2)-PS(2)*PB(1))
1463
1464           ELSE
1465
1466             CALL HWRAZM( ONE,CP,SP)
1467
1468           ENDIF
1469
1470         ELSE
1471
1472           CALL HWRAZM( ONE,CP,SP)
1473
1474         ENDIF
1475
1476         CALL HWUROT(PR,CP,SP,RS)
1477
1478 C---FIND BOOST FOR JET IP
1479
1480         ALF=(PHEP(3,IJET(IP))+PHEP(4,IJET(IP)))/
1481
1482      &      (PR(4)+SQRT((PR(4)+PR(5))*(PR(4)-PR(5))))
1483
1484         PB(1)=0.
1485
1486         PB(2)=0.
1487
1488         PB(3)=ALF-(1./ALF)
1489
1490         PB(4)=ALF+(1./ALF)
1491
1492         PB(5)=2.
1493
1494         IHEP=IJET(IP)
1495
1496         KHEP=JDAHEP(2,IHEP)
1497
1498         IF (KHEP.LT.IHEP) KHEP=IHEP
1499
1500         DO 220 JHEP=IHEP,KHEP
1501
1502         CALL HWULOF(PB,PHEP(1,JHEP),PHEP(1,JHEP))
1503
1504         CALL HWUROB(RS,PHEP(1,JHEP),PHEP(1,JHEP))
1505
1506         CALL HWULF4(PB,VHEP(1,JHEP),VHEP(1,JHEP))
1507
1508         CALL HWUROB(RS,VHEP(1,JHEP),VHEP(1,JHEP))
1509
1510 C---BOOST FROM BREIT FRAME IF NECESSARY
1511
1512         IF (DISPRO.AND.BREIT) THEN
1513
1514           CALL HWUROB(RBR,PHEP(1,JHEP),PHEP(1,JHEP))
1515
1516           CALL HWULOB(PBR,PHEP(1,JHEP),PHEP(1,JHEP))
1517
1518           CALL HWUROB(RBR,VHEP(1,JHEP),VHEP(1,JHEP))
1519
1520           CALL HWULB4(PBR,VHEP(1,JHEP),VHEP(1,JHEP))
1521
1522         ENDIF
1523
1524         CALL HWVSUM(4,VHEP(1,JHEP),VHEP(1,IPAR(IP)),VHEP(1,JHEP))
1525
1526   220   ISTHEP(JHEP)=ISTHEP(JHEP)+10
1527
1528         IF (KHEP.GT.IHEP+1) THEN
1529
1530           ISTHEP(IHEP+1)=100
1531
1532         ELSEIF (KHEP.EQ.IHEP) THEN
1533
1534 C---NON-PARTON JET
1535
1536           ISTHEP(IHEP)=190
1537
1538         ENDIF
1539
1540   230   CONTINUE
1541
1542         IF (ISTHEP(ICM).EQ.110) ISTHEP(ICM)=120
1543
1544       ENDIF
1545
1546       GOTO 20
1547
1548   999 END
1549
1550 CDECK  ID>, HWBMAS.
1551
1552 *CMZ :-        -26/04/91  11.11.54  by  Bryan Webber
1553
1554 *-- Author :    Bryan Webber
1555
1556 C-----------------------------------------------------------------------
1557
1558       SUBROUTINE HWBMAS
1559
1560 C-----------------------------------------------------------------------
1561
1562 C     Passes  backwards through a  jet cascade  calculating the masses
1563
1564 C     and magnitudes of the longitudinal and transverse three momenta.
1565
1566 C     Components given relative to direction of parent for a time-like
1567
1568 C     vertex and with respect to z-axis for space-like vertices.
1569
1570 C
1571
1572 C     On input PPAR(1-5,*) contains:
1573
1574 C     (E*sqrt(Xi),Xi,3-mom (if external),E,M-sq (if external))
1575
1576 C
1577
1578 C     On output PPAR(1-5,*) (if TMPAR(*)), containts:
1579
1580 C     (P-trans,Xi or Xilast,P-long,E,M)
1581
1582 C-----------------------------------------------------------------------
1583
1584       INCLUDE 'HERWIG61.INC'
1585
1586       DOUBLE PRECISION HWUSQR,EXI,PISQ,PJPK,EJEK,PTSQ,Z,ZMIN,ZMAX,
1587
1588      $     EMI,EMJ,EMK,C,NQ,HWBVMC,RHO,POLD,PNEW,EOLD,ENEW,A,B
1589
1590       INTEGER IPAR,JPAR,KPAR,MPAR,I,J,K
1591
1592       EXTERNAL HWUSQR
1593
1594       IF (IERROR.NE.0) RETURN
1595
1596       IF (NPAR.GT.2) THEN
1597
1598         DO 30 MPAR=NPAR-1,3,-2
1599
1600          JPAR=MPAR
1601
1602 C Find parent and partner of this branch
1603
1604          IPAR=JMOPAR(1,JPAR)
1605
1606          KPAR=JPAR+1
1607
1608 C Determine type of branching
1609
1610          IF (TMPAR(IPAR)) THEN
1611
1612 C Time-like branching
1613
1614 C           Compute mass of parent
1615
1616             EXI=PPAR(1,JPAR)*PPAR(1,KPAR)
1617
1618             PPAR(5,IPAR)=PPAR(5,JPAR)+PPAR(5,KPAR)+2.*EXI
1619
1620 C           Compute three momentum of parent
1621
1622             PISQ=PPAR(4,IPAR)*PPAR(4,IPAR)-PPAR(5,IPAR)
1623
1624             PPAR(3,IPAR)=HWUSQR(PISQ)
1625
1626 C---SPECIAL FOR G-->QQBAR: READJUST ANGULAR DISTRIBUTION
1627
1628             IF (IDPAR(IPAR).EQ.13 .AND. IDPAR(JPAR).LT.13) THEN
1629
1630               Z=PPAR(4,JPAR)/PPAR(4,IPAR)
1631
1632               ZMIN=HWBVMC(IDPAR(JPAR))/PPAR(1,JPAR)*Z
1633
1634               RHO=(Z*(3-Z*(3-2*Z))-ZMIN*(3-ZMIN*(3-2*ZMIN)))
1635
1636      $             /(2*(1-2*ZMIN)*(1-ZMIN*(1-ZMIN)))
1637
1638               NQ=PPAR(3,IPAR)*(PPAR(3,IPAR)+PPAR(4,IPAR))
1639
1640               EMI=PPAR(5,IPAR)
1641
1642               EMJ=PPAR(5,JPAR)
1643
1644               EMK=PPAR(5,KPAR)
1645
1646               ZMIN=MAX((EMI+EMJ-EMK)/(2*(EMI+NQ)),
1647
1648      $           (EMI+EMJ-EMK-SQRT((EMI-EMJ-EMK)**2-4*EMJ*EMK))/(2*EMI))
1649
1650               ZMAX=1-MAX((EMI-EMJ+EMK)/(2*(EMI+NQ)),
1651
1652      $           (EMI-EMJ+EMK-SQRT((EMI-EMJ-EMK)**2-4*EMJ*EMK))/(2*EMI))
1653
1654               C=2*RMASS(IDPAR(JPAR))**2/EMI
1655
1656               Z=(4*ZMIN*(1.5*(1+C-ZMIN)+ZMIN**2)*(1-RHO)
1657
1658      $          +4*ZMAX*(1.5*(1+C-ZMAX)+ZMAX**2)*RHO-2-3*C)/(1+2*C)**1.5
1659
1660               Z=SQRT(1+2*C)*SINH(LOG(Z+SQRT(Z**2+1))/3)+0.5
1661
1662               Z=(Z*NQ+(EMI+EMJ-EMK)/2)/(NQ+EMI)
1663
1664               PPAR(4,JPAR)=Z*PPAR(4,IPAR)
1665
1666               PPAR(4,KPAR)=PPAR(4,IPAR)-PPAR(4,JPAR)
1667
1668               PPAR(3,JPAR)=HWUSQR(PPAR(4,JPAR)**2-EMJ)
1669
1670               PPAR(3,KPAR)=HWUSQR(PPAR(4,KPAR)**2-EMK)
1671
1672               PPAR(2,JPAR)=EXI/(PPAR(4,JPAR)*PPAR(4,KPAR))
1673
1674               IF(JDAPAR(2,JPAR).NE.0)PPAR(2,JDAPAR(2,JPAR))=PPAR(2,JPAR)
1675
1676               IF(JDAPAR(2,KPAR).NE.0)PPAR(2,JDAPAR(2,KPAR))=PPAR(2,JPAR)
1677
1678 C---FIND DESCENDENTS OF THIS SPLITTING AND READJUST THEIR MOMENTA TOO
1679
1680               DO 20 J=JPAR+2,NPAR-1,2
1681
1682                 I=J
1683
1684  10             I=JMOPAR(1,I)
1685
1686                 IF (I.GT.IPAR) GOTO 10
1687
1688                 IF (I.EQ.IPAR) THEN
1689
1690                   I=JMOPAR(1,J)
1691
1692                   K=J+1
1693
1694                   POLD=PPAR(3,J)+PPAR(3,K)
1695
1696                   EOLD=PPAR(4,J)+PPAR(4,K)
1697
1698                   PNEW=HWUSQR(PPAR(4,I)**2-PPAR(5,I))
1699
1700                   ENEW=PPAR(4,I)
1701
1702                   A=(ENEW*EOLD-PNEW*POLD)/PPAR(5,I)
1703
1704                   B=(PNEW*EOLD-ENEW*POLD)/PPAR(5,I)
1705
1706                   PPAR(3,J)=A*PPAR(3,J)+B*PPAR(4,J)
1707
1708                   PPAR(4,J)=(PPAR(4,J)+B*PPAR(3,J))/A
1709
1710                   PPAR(3,K)=PNEW-PPAR(3,J)
1711
1712                   PPAR(4,K)=ENEW-PPAR(4,J)
1713
1714                   PPAR(2,J)=1-(PPAR(3,J)*PPAR(3,K)+PPAR(1,J)*PPAR(1,K))
1715
1716      $                 /(PPAR(4,J)*PPAR(4,K))
1717
1718                   IF (JDAPAR(2,J).NE.0) PPAR(2,JDAPAR(2,J))=PPAR(2,J)
1719
1720                   IF (JDAPAR(2,K).NE.0) PPAR(2,JDAPAR(2,K))=PPAR(2,J)
1721
1722                 ENDIF
1723
1724  20           CONTINUE
1725
1726             ENDIF
1727
1728 C           Compute daughter' transverse and longitudinal momenta
1729
1730             PJPK=PPAR(3,JPAR)*PPAR(3,KPAR)
1731
1732             EJEK=PPAR(4,JPAR)*PPAR(4,KPAR)-EXI
1733
1734             PTSQ=(PJPK+EJEK)*(PJPK-EJEK)/PISQ
1735
1736             PPAR(1,JPAR)=HWUSQR(PTSQ)
1737
1738             PPAR(3,JPAR)=HWUSQR(PPAR(3,JPAR)*PPAR(3,JPAR)-PTSQ)
1739
1740             PPAR(1,KPAR)=-PPAR(1,JPAR)
1741
1742             PPAR(3,KPAR)= PPAR(3,IPAR)-PPAR(3,JPAR)
1743
1744          ELSE
1745
1746 C Space-like branching
1747
1748 C           Re-arrange such that JPAR is time-like
1749
1750             IF (TMPAR(KPAR)) THEN
1751
1752                KPAR=JPAR
1753
1754                JPAR=JPAR+1
1755
1756             ENDIF
1757
1758 C           Compute time-like branch
1759
1760             PTSQ=(2.-PPAR(2,JPAR))*PPAR(1,JPAR)*PPAR(1,JPAR)
1761
1762      &          -PPAR(5,JPAR)
1763
1764             PPAR(1,JPAR)=HWUSQR(PTSQ)
1765
1766             PPAR(3,JPAR)=(1.-PPAR(2,JPAR))*PPAR(4,JPAR)
1767
1768             PPAR(3,IPAR)=PPAR(3,KPAR)-PPAR(3,JPAR)
1769
1770             PPAR(5,IPAR)=0.
1771
1772             PPAR(1,KPAR)=0.
1773
1774          ENDIF
1775
1776 C Reset Xi to Xilast
1777
1778          PPAR(2,KPAR)=PPAR(2,IPAR)
1779
1780  30    CONTINUE
1781
1782       ENDIF
1783
1784       DO 40 IPAR=2,NPAR
1785
1786  40   PPAR(5,IPAR)=HWUSQR(PPAR(5,IPAR))
1787
1788       PPAR(1,2)=0.
1789
1790       PPAR(2,2)=0.
1791
1792       END