HERWIG/src/hwcbct.f

   1
   2 CDECK  ID>, HWCBCT.
   3
   4 *CMZ :-        -20/07/99  10:56:12  by  Peter Richardson
   5
   6 *-- Author :    Peter Richardson
   7
   8 C-----------------------------------------------------------------------
   9
  10       SUBROUTINE HWCBCT(JHEP,KHEP,THEP,PCL,SPLIT)
  11
  12 C-----------------------------------------------------------------------
  13
  14 C  Subroutine to split a baryonic cluster containing two heavy quarks
  15
  16 C  Based on HWCCUT
  17
  18 C-----------------------------------------------------------------------
  19
  20       INCLUDE 'HERWIG61.INC'
  21
  22       DOUBLE PRECISION HWUPCM,HWR,HWVDOT,EMC,QM1,QM2,QM3,QM4,
  23
  24      &                 PXY,PCX,PCY,RCM,PCL(5),AX(5),PA(5),PB(5),PC(5),
  25
  26      &                 VCLUS(4),DQM,EMX,EMY,SKAPPA,RKAPPA,VTMP(4),
  27
  28      &                 DELTM,PDIQUK(5),AY(5)
  29
  30       INTEGER HWRINT,JHEP,KHEP,LHEP,MHEP,THEP,ID1,ID2,ID3,ID4,NTRY,
  31
  32      &        NTRYMX,J,IB
  33
  34       LOGICAL SPLIT
  35
  36       EXTERNAL HWUPCM,HWR,HWVDOT
  37
  38       PARAMETER(SKAPPA=1.,NTRYMX=100)
  39
  40       IF(IERROR.NE.0) RETURN
  41
  42       EMC=PCL(5)
  43
  44       ID1=IDHW(JHEP)
  45
  46       ID2=IDHW(KHEP)
  47
  48       ID3=IDHW(THEP)
  49
  50       QM1=RMASS(ID1)
  51
  52       QM2=RMASS(ID2)
  53
  54       QM3=RMASS(ID3)
  55
  56       SPLIT = .FALSE.
  57
  58       NTRY = 0
  59
  60 C Decide if cluster contains a b-(anti)quark
  61
  62       IF (ID1.EQ.5.OR.ID1.EQ.11.OR.ID2.EQ.5.OR.ID2.EQ.11.OR.
  63
  64      &    ID3.EQ.5.OR.ID3.EQ.11) THEN
  65
  66         IB=2
  67
  68       ELSE
  69
  70         IB=1
  71
  72       ENDIF
  73
  74 C-- Set the positon of the cluster to be that of the heavy quark
  75
  76       CALL HWVEQU(4,VHEP(1,THEP),VCLUS)
  77
  78 C--SPLIT THE BARYONIC CLUSTER INTO A HEAVY FLAVOUR MESON AND A HEAVY
  79
  80 C--FLAVOUR BARYON
  81
  82       PXY=EMC-QM1-QM2-QM3
  83
  84  20   NTRY=NTRY+1
  85
  86       IF(NTRY.GT.NTRYMX) RETURN
  87
  88  30   EMX=QM1+QM2+PXY*HWR()**PSPLT(IB)
  89
  90       EMY=    QM3+PXY*HWR()**PSPLT(IB)
  91
  92       IF(EMX+EMY.GE.EMC) GOTO 30
  93
  94 C--PULL A LIGHT QUARK PAIR OUT OF THE VACUUM
  95
  96  40   ID4=HWRINT(1,3)
  97
  98       IF(QWT(ID4).LT.HWR()) GOTO 40
  99
 100       QM4=RMASS(ID4)
 101
 102 C--Now combine particles 3 & 4 into a diquark
 103
 104 C--If three also heavy this diquark doesn't exist in HERWIG
 105
 106 C--just assume mass is sum of quark masses,as for other diquarks
 107
 108       DQM=QM3+QM4
 109
 110 C--Now obtain the masses for the cluster splitting
 111
 112       PCX=HWUPCM(EMX,QM1,DQM)
 113
 114       IF(PCX.LT.ZERO) GOTO 20
 115
 116       PCY=HWUPCM(EMY,QM2,QM4)
 117
 118       IF(PCY.LT.ZERO) GOTO 20
 119
 120       SPLIT=.TRUE.
 121
 122 C--Now we've decided which light quark to pull out of the vacuum
 123
 124 C--Find the direction of the second heavy quark
 125
 126       CALL HWULOF(PCL,PHEP(1,THEP),AX)
 127
 128       RCM=1./SQRT(HWVDOT(3,AX,AX))
 129
 130       CALL HWVSCA(3,RCM,AX,AX)
 131
 132 C--Construct the new CoM momenta(collinear)
 133
 134       PXY=HWUPCM(EMC,EMX,EMY)
 135
 136       CALL HWVSCA(3,PXY,AX,PC)
 137
 138 C--pc is momenta of Y cluster along 2nd quark dirn in cluster frame
 139
 140       PC(4)=SQRT(PXY**2+EMY**2)
 141
 142       PC(5)=EMY
 143
 144 C--pa is momenta of 2nd quark in Y frame
 145
 146       CALL HWVSCA(3,PCY,AX,PA)
 147
 148       PA(4)=SQRT(PCY**2+QM3**2)
 149
 150       PA(5)=QM3
 151
 152 C--pb is momenta of 2nd quark in cluster frame,pa now momenta of antiquark
 153
 154       CALL HWULOB(PC,PA,PB)
 155
 156       CALL HWVDIF(4,PC,PB,PA)
 157
 158       PA(5)=QM4
 159
 160       LHEP=NHEP+1
 161
 162       MHEP=NHEP+2
 163
 164 C--boost these momenta back to lab frame
 165
 166       CALL HWULOB(PCL,PB,PHEP(1,THEP))
 167
 168       CALL HWULOB(PCL,PA,PHEP(1,MHEP))
 169
 170 C--pc now becomes momenta of X cluster in cluster frame
 171
 172       CALL HWVSCA(3,-ONE,PC,PC)
 173
 174       PC(4)=EMC-PC(4)
 175
 176       PC(5)=EMX
 177
 178 C--find the dirn of the 1st heavy quark in the X frame
 179
 180 C--transform to cluster frame
 181
 182       CALL HWULOF(PCL,PHEP(1,JHEP),AY)
 183
 184 C--transform to X-frame
 185
 186       CALL HWULOF(PC,AY,AY)
 187
 188       RCM=1./SQRT(HWVDOT(3,AY,AY))
 189
 190       CALL HWVSCA(3,RCM,AY,AY)
 191
 192 C--pa now momenta of 1st havy quark along this dirn
 193
 194       CALL HWVSCA(3,PCX,AY,PA)
 195
 196       PA(4)=SQRT(PCX**2+QM1**2)
 197
 198       PA(5)=QM1
 199
 200 C--pb now momenta of 1st heavy quark in cluster frame then to lab
 201
 202       CALL HWULOB(PC,PA,PB)
 203
 204       CALL HWULOB(PCL,PB,PHEP(1,JHEP))
 205
 206 C--now find the diquark momenta by momentum conservation
 207
 208       DO 50 J=1,4
 209
 210  50   PDIQUK(J)=PCL(J)-PHEP(J,THEP)-PHEP(J,MHEP)-PHEP(J,JHEP)
 211
 212       PDIQUK(5)=DQM
 213
 214 C--Now obtain the quark momenta from the diquark
 215
 216       DO 60 J=1,3
 217
 218  60   PA(J) = 0
 219
 220       PA(4) = QM2
 221
 222       PA(5) = QM2
 223
 224       CALL HWULOB(PDIQUK,PA,PHEP(1,KHEP))
 225
 226       CALL HWVDIF(4,PDIQUK,PHEP(1,KHEP),PHEP(1,LHEP))
 227
 228 C--Construct new vertex positions
 229
 230       RKAPPA=GEV2MM/SKAPPA
 231
 232       CALL HWVSCA(3,RKAPPA,AX,AX)
 233
 234       DELTM=(EMX-EMY)*(EMX+EMY)/(TWO*EMC)
 235
 236       CALL HWVSCA(3,DELTM,AX,VTMP)
 237
 238       VTMP(4)=(HALF*EMC-PXY)*RKAPPA
 239
 240       CALL HWULB4(PCL,VTMP,VTMP)
 241
 242       CALL HWVSUM(4,VTMP,VCLUS,VHEP(1,LHEP))
 243
 244       CALL HWVEQU(4,VHEP(1,LHEP),VHEP(1,MHEP))
 245
 246 C--Relabel the colours of the quarks
 247
 248       IDHEP(LHEP) = IDPDG(ID4)
 249
 250       IDHEP(MHEP) = IDPDG(ID4)
 251
 252       IF(IDHEP(JHEP).GT.0) THEN
 253
 254         IDHW(LHEP)  = ID4+6
 255
 256         IDHEP(LHEP) = -IDHEP(LHEP)
 257
 258         IDHW(MHEP)  = ID4
 259
 260         JDAHEP(2,LHEP) = JHEP
 261
 262         JMOHEP(2,LHEP) = MHEP
 263
 264         JMOHEP(2,MHEP) = JMOHEP(2,JHEP)
 265
 266         JDAHEP(2,MHEP) = LHEP
 267
 268         JMOHEP(2,JHEP) = LHEP
 269
 270       ELSE
 271
 272         IDHW(LHEP)  = ID4
 273
 274         IDHW(MHEP)  = ID4+6
 275
 276         IDHEP(MHEP) = -IDHEP(MHEP)
 277
 278         JMOHEP(2,LHEP) = JHEP
 279
 280         JDAHEP(2,MHEP) = JDAHEP(2,JHEP)
 281
 282         JDAHEP(2,LHEP) = MHEP
 283
 284         JMOHEP(2,MHEP) = LHEP
 285
 286         JDAHEP(2,JHEP) = LHEP
 287
 288       ENDIF
 289
 290       ISTHEP(LHEP) = 151
 291
 292       ISTHEP(MHEP) = 151
 293
 294       JMOHEP(1,LHEP) = JMOHEP(1,KHEP)
 295
 296       JDAHEP(1,LHEP) = 0
 297
 298       JMOHEP(1,MHEP) = JMOHEP(1,JHEP)
 299
 300       JDAHEP(1,MHEP) = 0
 301
 302       NHEP = NHEP+2
 303
 304  999  END
 305
 306 CDECK  ID>, HWCBVI.
 307
 308 *CMZ :-
 309
 310 *-- Author :    Mark Gibbs  modified by Peter Richardson
 311
 312 C-----------------------------------------------------------------------
 313
 314       SUBROUTINE HWCBVI
 315
 316 C-----------------------------------------------------------------------
 317
 318 C FINDS UNPAIRED PARTONS AFTER BARYON-NUMBER VIOLATION
 319
 320 C  MODIFIED FOR RPARITY VIOLATING SUSY
 321
 322 C-----------------------------------------------------------------------
 323
 324       INCLUDE 'HERWIG61.INC'
 325
 326       COMMON/HWBVIC/NBV,IBV(18)
 327
 328       DOUBLE PRECISION HWR,PDQ(5)
 329
 330       INTEGER NBV,IBV,JBV,KBV,LBV,IHEP,IP1,IP2,IP3,JP1,JP2,JP3,
 331
 332      & HWCBVT,NBR,MBV,IQ1,IQ2,IQ3,ID1,ID2,IDQ,IDIQK(3,3)
 333
 334       LOGICAL SPLIT,DUNBV(18)
 335
 336       DATA IDIQK/111,110,113,110,109,112,113,112,114/
 337
 338 C---Check for errors
 339
 340       IF (IERROR.NE.0)  RETURN
 341
 342 C---Correct colour connections are gluon splitting
 343
 344       CALL HWCCCC
 345
 346 C---Reset bvi clustering flag
 347
 348       HVFCEN = .FALSE.
 349
 350 C---LIST PARTONS WITH WRONG COLOUR PARTNERS-QUARKS ONLY
 351
 352     5 NBV=0
 353
 354       DO 10 IHEP=1,NHEP
 355
 356       IF (ISTHEP(IHEP).GT.149.AND.ISTHEP(IHEP).LT.155) THEN
 357
 358         IF (QORQQB(IDHW(IHEP))) THEN
 359
 360           IF (.NOT.QORQQB(IDHW(JMOHEP(2,IHEP))).
 361
 362      &        AND.JMOHEP(2,IHEP).GT.6) GOTO 10
 363
 364         ELSE
 365
 366 C---Extra check for Gamma's
 367
 368           IF (IDHW(IHEP).EQ.59) GO TO 10
 369
 370 C---End of bug fix.
 371
 372           IF (QORQQB(IDHW(JDAHEP(2,IHEP)))) GO TO 10
 373
 374           GO TO 10
 375
 376         ENDIF
 377
 378         IF(JMOHEP(2,IHEP).LT.6.AND.
 379
 380      &     .NOT.QBORQQ(IDHW(JMOHEP(2,IHEP)))) GOTO 10
 381
 382 C--new for hard process
 383
 384         NBV=NBV+1
 385
 386         IF (NBV.GT.18) CALL HWWARN('HWCBVI',100,*999)
 387
 388         IBV(NBV)=IHEP
 389
 390         DUNBV(NBV)=.FALSE.
 391
 392       ENDIF
 393
 394    10 CONTINUE
 395
 396 C--NOW FIND THE ANTIQUARKS WITH WRONG COLOUR CONNECTIONS
 397
 398       DO 11 IHEP=1,NHEP
 399
 400       IF(ISTHEP(IHEP).GT.149.AND.ISTHEP(IHEP).LT.155) THEN
 401
 402         IF(QBORQQ(IDHW(IHEP))) THEN
 403
 404           IF(.NOT.QBORQQ(IDHW(JDAHEP(2,IHEP))).AND.
 405
 406      &        JDAHEP(2,IHEP).GT.6) GO TO 11
 407
 408         ELSE
 409
 410 C--Extra check for gamma's
 411
 412           IF(IDHW(IHEP).EQ.59) GO TO 11
 413
 414           IF(QBORQQ(IDHW(JMOHEP(2,IHEP)))) GO TO 11
 415
 416           GO TO 11
 417
 418         ENDIF
 419
 420         IF(JDAHEP(2,IHEP).LT.6.AND.
 421
 422      &    .NOT.QORQQB(IDHW(JDAHEP(2,IHEP)))) GOTO 11
 423
 424         NBV=NBV+1
 425
 426         IF(NBV.GT.18) CALL HWWARN('HWCBVI',100,*999)
 427
 428         IBV(NBV)=IHEP
 429
 430         DUNBV(NBV)=.FALSE.
 431
 432       ENDIF
 433
 434  11   CONTINUE
 435
 436       IF (NBV.EQ.0) RETURN
 437
 438       IF(MOD(NBV,3).NE.0) CALL HWWARN('HWCBVI',101,*999)
 439
 440 C---PROCESS FOUND PARTONS, STARTING AT RANDOM POINT IN LIST
 441
 442       NBR=NBV*HWR()
 443
 444       DO 100 MBV=1,NBV
 445
 446       JBV=MBV+NBR
 447
 448       IF (JBV.GT.NBV) JBV=JBV-NBV
 449
 450       IF (.NOT.DUNBV(JBV)) THEN
 451
 452         DUNBV(JBV)=.TRUE.
 453
 454         IP1=IBV(JBV)
 455
 456         JP1=HWCBVT(IP1)
 457
 458 C---FIND ASSOCIATED PARTONS
 459
 460         DO 20 KBV=1,NBV
 461
 462         IF (.NOT.DUNBV(KBV)) THEN
 463
 464           IP2=IBV(KBV)
 465
 466           JP2=HWCBVT(IP2)
 467
 468           IF (JP2.EQ.JP1) THEN
 469
 470             DUNBV(KBV)=.TRUE.
 471
 472             DO 15 LBV=1,NBV
 473
 474             IF (.NOT.DUNBV(LBV)) THEN
 475
 476               IP3=IBV(LBV)
 477
 478               JP3=HWCBVT(IP3)
 479
 480               IF (JP3.EQ.JP2) THEN
 481
 482                 DUNBV(LBV)=.TRUE.
 483
 484                 GO TO 25
 485
 486               ENDIF
 487
 488             ENDIF
 489
 490    15       CONTINUE
 491
 492           ENDIF
 493
 494         ENDIF
 495
 496    20   CONTINUE
 497
 498         CALL HWWARN('HWCBVI',102,*999)
 499
 500    25   IQ1=0
 501
 502 C---LOOK FOR DIQUARK
 503
 504         IF (ABS(IDHEP(IP1)).GT.100) THEN
 505
 506           IQ1=IP1
 507
 508           IQ2=IP2
 509
 510           IQ3=IP3
 511
 512         ELSEIF (ABS(IDHEP(IP2)).GT.100) THEN
 513
 514           IQ1=IP2
 515
 516           IQ2=IP3
 517
 518           IQ3=IP1
 519
 520         ELSEIF (ABS(IDHEP(IP3)).GT.100) THEN
 521
 522           IQ1=IP3
 523
 524           IQ2=IP1
 525
 526           IQ3=IP2
 527
 528         ENDIF
 529
 530         IF (IQ1.EQ.0) THEN
 531
 532 C---NO DIQUARKS: COMBINE TWO (ANTI)QUARKS
 533
 534           IF (ABS(IDHEP(IP1)).GT.3) THEN
 535
 536             IQ1=IP2
 537
 538             IQ2=IP3
 539
 540             IQ3=IP1
 541
 542           ELSEIF (ABS(IDHEP(IP2)).GT.3) THEN
 543
 544             IQ1=IP3
 545
 546             IQ2=IP1
 547
 548             IQ3=IP2
 549
 550           ELSE
 551
 552             IQ1=IP1
 553
 554             IQ2=IP2
 555
 556             IQ3=IP3
 557
 558           ENDIF
 559
 560           ID1=IDHEP(IQ1)
 561
 562           ID2=IDHEP(IQ2)
 563
 564 C---CHECK FLAVOURS
 565
 566           IF (ID1.GT.0.AND.ID1.LT.4.AND.
 567
 568      &        ID2.GT.0.AND.ID2.LT.4) THEN
 569
 570             IDQ=IDIQK(ID1,ID2)
 571
 572           ELSEIF (ID1.LT.0.AND.ID1.GT.-4.AND.
 573
 574      &            ID1.LT.0.AND.ID2.GT.-4) THEN
 575
 576             IDQ=IDIQK(-ID1,-ID2)+6
 577
 578           ELSE
 579
 580 C---CANT MAKE DIQUARKS WITH HEAVY QUARKS: TRY CLUSTER SPLITTING
 581
 582             CALL HWVSUM(4,PHEP(1,IQ1),PHEP(1,IQ2),PDQ)
 583
 584             CALL HWUMAS(PDQ)
 585
 586 C--Use the original splitting procedure
 587
 588             CALL HWCCUT(IQ1,IQ2,PDQ,.FALSE.,SPLIT)
 589
 590             IF(SPLIT) GOTO 5
 591
 592 C--If it fails try the new procedure
 593
 594             CALL HWVSUM(4,PDQ,PHEP(1,IQ3),PDQ)
 595
 596             CALL HWUMAS(PDQ)
 597
 598             IF(ABS(ID1).GT.3) THEN
 599
 600               CALL HWCBCT(IQ3,IQ2,IQ1,PDQ,SPLIT)
 601
 602             ELSEIF(ABS(ID2).GT.3) THEN
 603
 604               CALL HWCBCT(IQ3,IQ1,IQ2,PDQ,SPLIT)
 605
 606             ELSE
 607
 608               CALL HWWARN('HWCBVI',100,*999)
 609
 610             ENDIF
 611
 612             IF (SPLIT) GO TO 5
 613
 614 C---Unable to form cluster; dispose of event
 615
 616             CALL HWWARN('HWCBVI',-3,*999)
 617
 618           ENDIF
 619
 620 C---OVERWRITE FIRST AND CANCEL SECOND
 621
 622           IDHW(IQ1)=IDQ
 623
 624           IDHEP(IQ1)=IDPDG(IDQ)
 625
 626           CALL HWVSUM(4,PHEP(1,IQ1),PHEP(1,IQ2),PHEP(1,IQ1))
 627
 628           CALL HWUMAS(PHEP(1,IQ1))
 629
 630           ISTHEP(IQ2)=0
 631
 632 C---REMAKE COLOUR CONNECTIONS
 633
 634           IF (QORQQB(IDQ)) THEN
 635
 636             JMOHEP(2,IQ1)=IQ3
 637
 638             JDAHEP(2,IQ3)=IQ1
 639
 640           ELSE
 641
 642             JDAHEP(2,IQ1)=IQ3
 643
 644             JMOHEP(2,IQ3)=IQ1
 645
 646           ENDIF
 647
 648         ELSE
 649
 650 C---SPLIT A DIQUARK
 651
 652           NHEP=NHEP+1
 653
 654           CALL HWVSCA(5,HALF,PHEP(1,IQ1),PHEP(1,IQ1))
 655
 656           CALL HWVEQU(5,PHEP(1,IQ1),PHEP(1,NHEP))
 657
 658           ISTHEP(NHEP)=150
 659
 660           JMOHEP(1,NHEP)=JMOHEP(1,IQ1)
 661
 662           JDAHEP(1,NHEP)=0
 663
 664 C---FIND FLAVOURS
 665
 666           IDQ=IDHW(IQ1)
 667
 668           DO 30 ID2=1,3
 669
 670           DO 30 ID1=1,3
 671
 672           IF (IDIQK(ID1,ID2).EQ.IDQ) THEN
 673
 674             IDHW(IQ1)=ID1
 675
 676             IDHW(NHEP)=ID2
 677
 678 C---REMAKE COLOUR CONNECTIONS (DIQUARK)
 679
 680             JMOHEP(2,IQ1)=IQ2
 681
 682             JMOHEP(2,IQ2)=NHEP
 683
 684             JMOHEP(2,IQ3)=IQ1
 685
 686             JMOHEP(2,NHEP)=IQ3
 687
 688             JDAHEP(2,IQ1)=IQ3
 689
 690             JDAHEP(2,IQ2)=IQ1
 691
 692             JDAHEP(2,IQ3)=NHEP
 693
 694             JDAHEP(2,NHEP)=IQ2
 695
 696             GO TO 35
 697
 698           ELSEIF (IDIQK(ID1,ID2).EQ.IDQ-6) THEN
 699
 700             IDHW(IQ1)=ID1+6
 701
 702             IDHW(NHEP)=ID2+6
 703
 704 C---REMAKE COLOUR CONNECTIONS (ANTIDIQUARK)
 705
 706             JMOHEP(2,IQ1)=IQ3
 707
 708             JMOHEP(2,IQ2)=IQ1
 709
 710             JMOHEP(2,IQ3)=NHEP
 711
 712             JMOHEP(2,NHEP)=IQ2
 713
 714             JDAHEP(2,IQ1)=IQ2
 715
 716             JDAHEP(2,IQ2)=NHEP
 717
 718             JDAHEP(2,IQ3)=IQ1
 719
 720             JDAHEP(2,NHEP)=IQ3
 721
 722             GO TO 35
 723
 724           ENDIF
 725
 726    30     CONTINUE
 727
 728           CALL HWWARN('HWCBVI',104,*999)
 729
 730    35     IDHEP(IQ1)=IDPDG(IDHW(IQ1))
 731
 732           IDHEP(NHEP)=IDPDG(IDHW(NHEP))
 733
 734         ENDIF
 735
 736       ENDIF
 737
 738   100 CONTINUE
 739
 740       RETURN
 741
 742   999 END
 743
 744 CDECK  ID>, HWCBVT.
 745
 746 *CMZ :-
 747
 748 *-- Author :    Peter Richardson
 749
 750 C-----------------------------------------------------------------------
 751
 752       FUNCTION HWCBVT(IP)
 753
 754 C-----------------------------------------------------------------------
 755
 756 C  Function to find the baryon number violating vertex a parton came from
 757
 758 C-----------------------------------------------------------------------
 759
 760       INCLUDE 'HERWIG61.INC'
 761
 762       INTEGER HWCBVT,IP,JP(2),KP,I,J,ID,TYPE,IDM,IDM2,IDM3,IDM4
 763
 764       JP(1) = IP
 765
 766       ID = IDHW(IP)
 767
 768       IF(ID.LE.6.OR.(ID.GE.115.AND.ID.LE.120)) THEN
 769
 770         JP(2) = JMOHEP(2,IP)
 771
 772       ELSE
 773
 774         JP(2) = JDAHEP(2,IP)
 775
 776       ENDIF
 777
 778       DO I=1,2
 779
 780         IDM = JMOHEP(1,JMOHEP(1,JMOHEP(1,JMOHEP(1,JP(I)))))
 781
 782         IF(IDHW(IDM).EQ.6.OR.IDHW(IDM).EQ.12) THEN
 783
 784           JP(I)=IDM
 785
 786         ENDIF
 787
 788       ENDDO
 789
 790       DO J=1,7
 791
 792         DO I=1,2
 793
 794           KP = JMOHEP(1,JP(I))
 795
 796           IDM = IDHW(KP)
 797
 798           IDM2 = IDHW(JDAHEP(1,KP))
 799
 800           IDM3 = IDHW(JDAHEP(2,KP))
 801
 802           IDM4 = IDHW(JDAHEP(1,KP)+1)
 803
 804           IF((ISTHEP(KP).EQ.155.AND.
 805
 806      &      ((IDM.GE.449.AND.IDM.LE.457.AND.IDM2.LE.12.AND.
 807
 808      &       IDM3.LE.12.AND.IDM4.LE.12).OR.
 809
 810      &      (((IDM.GE.411.AND.IDM.LE.424).OR.IDM.EQ.405.OR.IDM.EQ.406)
 811
 812      &      .AND.IDM2.LE.12.AND.IDM3.LE.12)))
 813
 814      &        .OR.(IDM.EQ.15.AND.IDM2.LE.12.AND.
 815
 816      &       IDHW(JMOHEP(1,KP)).LE.12.AND.
 817
 818      &       IDHW(JMOHEP(2,KP)).LE.12.AND.IDM3.GE.449.AND.
 819
 820      &       IDM3.LE.457).OR.
 821
 822      &         (IDM.EQ.15.AND.IDM2.GE.198.AND.IDM2.LE.200.
 823
 824      &          AND.ABS(IDPDG(IDM3)).GT.1000000)) THEN
 825
 826             IF(IDHW(KP).EQ.449.AND.JDAHEP(1,KP).EQ.JP(I)) THEN
 827
 828               KP = JMOHEP(1,KP)
 829
 830             ELSEIF(IDHW(KP).EQ.15) THEN
 831
 832               TYPE=IDHW(JDAHEP(1,KP))
 833
 834               IF(TYPE.GE.7.AND.TYPE.LE.12.AND.
 835
 836      &           JMOHEP(2,JDAHEP(2,KP)).EQ.JP(I)) THEN
 837
 838                 KP=IP
 839
 840               ELSEIF(TYPE.LE.6.AND.
 841
 842      &           JDAHEP(2,JDAHEP(2,KP)).EQ.JP(I)) THEN
 843
 844                 KP=IP
 845
 846               ELSE
 847
 848                 HWCBVT = KP
 849
 850                 RETURN
 851
 852               ENDIF
 853
 854             ELSE
 855
 856               HWCBVT = KP
 857
 858               RETURN
 859
 860             ENDIF
 861
 862           ENDIF
 863
 864           JP(I) =KP
 865
 866         ENDDO
 867
 868       ENDDO
 869
 870       HWCBVT = 0
 871
 872  999  END
 873
 874 CDECK  ID>, HWCCCC.
 875
 876 *CMZ :-
 877
 878 *-- Author :    Peter Richardson
 879
 880 C-----------------------------------------------------------------------
 881
 882       SUBROUTINE HWCCCC
 883
 884 C-----------------------------------------------------------------------
 885
 886 C  Subroutine to correct colour connections after the gluon splitting
 887
 888 C-----------------------------------------------------------------------
 889
 890       INCLUDE 'HERWIG61.INC'
 891
 892       INTEGER IHEP,STFSPT,LHEP,MHEP,RHEP
 893
 894       IF(IERROR.NE.0) RETURN
 895
 896 C--Find the first particle in the event record with status 150
 897
 898       DO IHEP=1,NHEP
 899
 900         IF(ISTHEP(IHEP).GE.150.AND.ISTHEP(IHEP).LE.154) THEN
 901
 902           STFSPT = IHEP
 903
 904           GOTO 10
 905
 906         ENDIF
 907
 908       ENDDO
 909
 910  10   CONTINUE
 911
 912 C--Now find any that are colour connected to earlier particles
 913
 914 C--in the event record
 915
 916       DO IHEP=STFSPT,NHEP
 917
 918 C--First the quarks and antidiquarks
 919
 920         IF(IDHW(IHEP).LT.6.OR.
 921
 922      &     (IDHW(IHEP).GE.115.AND.IDHW(IHEP).LE.120)) THEN
 923
 924           IF(JMOHEP(2,IHEP).LT.STFSPT) THEN
 925
 926             LHEP = IHEP
 927
 928             MHEP = JMOHEP(2,IHEP)
 929
 930             RHEP = MHEP
 931
 932             IF(MHEP.GT.6) RHEP = JDAHEP(1,MHEP)
 933
 934 C--As from Rparity connect to particle not to antiparticle
 935
 936             IF(IDHW(MHEP).NE.13) THEN
 937
 938               JMOHEP(2,LHEP) = RHEP
 939
 940             ELSE
 941
 942               RHEP = RHEP+1
 943
 944               JMOHEP(2,LHEP) = RHEP
 945
 946             ENDIF
 947
 948           ENDIF
 949
 950         ENDIF
 951
 952 C--Now the antiquarks
 953
 954         IF((IDHW(IHEP).GT.6.AND.IDHW(IHEP).LE.12).OR.
 955
 956      &     (IDHW(IHEP).GE.109.AND.IDHW(IHEP).LE.114)) THEN
 957
 958           IF(JDAHEP(2,IHEP).LT.STFSPT) THEN
 959
 960             LHEP = IHEP
 961
 962             MHEP = JDAHEP(2,IHEP)
 963
 964             RHEP = MHEP
 965
 966             IF(MHEP.GT.6) RHEP = JDAHEP(1,MHEP)
 967
 968 C--As from Rparity connect to antiparticle not particle
 969
 970             IF(IDHW(MHEP).NE.13) THEN
 971
 972               JDAHEP(2,LHEP) = RHEP
 973
 974             ELSE
 975
 976               JDAHEP(2,LHEP) = RHEP
 977
 978             ENDIF
 979
 980           ENDIF
 981
 982         ENDIF
 983
 984       ENDDO
 985
 986       END