HERWIG/src/hwdhig.f

   1
   2 CDECK  ID>, HWDHIG.
   3
   4 *CMZ :-        -24/04/92  14.23.44  by  Mike Seymour
   5
   6 *-- Author :    Mike Seymour
   7
   8 C-----------------------------------------------------------------------
   9
  10       SUBROUTINE HWDHIG(GAMINP)
  11
  12 C-----------------------------------------------------------------------
  13
  14 C     HIGGS DECAY ROUTINE
  15
  16 C     A) FOR GAMinp=0 FIND AND DECAY HIGGS
  17
  18 C     B) FOR GAMinp>0 CALCULATE TOTAL HIGGS WIDTH
  19
  20 C                     FOR EMH=GAMINP. STORE RESULT IN GAMINP.
  21
  22 C-----------------------------------------------------------------------
  23
  24       INCLUDE 'HERWIG61.INC'
  25
  26       DOUBLE PRECISION HWDHGF,HWR,HWRUNI,HWUSQR,HWUPCM,GAMINP,EMH,
  27
  28      & EMF,COLFAC,ENF,K1,K0,BET0,BET1,GAM0,GAM1,SCLOG,CFAC,XF,EM,GAMLIM,
  29
  30      & GAM,XW,EMW,XZ,EMZ,YW,YZ,EMI,TAUT,TAUW,WIDHIG,VECDEC,EMB,GAMB,
  31
  32      & TMIN,TMAX1,EM1,TMAX2,EM2,X1,X2,PROB,PCM,SUMR,SUMI,TAUTR,TAUTI,
  33
  34      & TAUWR,TAUWI,GFACTR
  35
  36       INTEGER HWRINT,IHIG,I,IFERM,NLOOK,I1,I2,IPART,IMODE,IDEC,MMAX
  37
  38       LOGICAL HWRLOG
  39
  40       EXTERNAL HWDHGF,HWR,HWRUNI,HWUSQR,HWUPCM,HWRINT,HWRLOG
  41
  42       SAVE GAM,EM,VECDEC
  43
  44       PARAMETER (NLOOK=100)
  45
  46       DIMENSION VECDEC(2,0:NLOOK)
  47
  48       EQUIVALENCE (EMW,RMASS(198)),(EMZ,RMASS(200))
  49
  50       DATA GAMLIM,GAM,EM/10D0,2*0D0/
  51
  52 C---IF DECAY, FIND HIGGS (HWDHAD WILL HAVE GIVEN IT STATUS=1)
  53
  54       IF (GAMINP.EQ.ZERO) THEN
  55
  56         IHIG=0
  57
  58         DO 10 I=1,NHEP
  59
  60  10       IF (IHIG.EQ.0.AND.IDHW(I).EQ.201.AND.ISTHEP(I).EQ.1) IHIG=I
  61
  62         IF (IHIG.EQ.0) CALL HWWARN('HWDHIG',101,*999)
  63
  64         EMH=PHEP(5,IHIG)
  65
  66         IF (EMH.LE.ZERO) CALL HWWARN('HWDHIG',102,*999)
  67
  68         EMSCA=EMH
  69
  70       ELSE
  71
  72         EMH=GAMINP
  73
  74         IF (EMH.LE.ZERO) THEN
  75
  76           GAMINP=0
  77
  78           RETURN
  79
  80         ENDIF
  81
  82       ENDIF
  83
  84 C---CALCULATE BRANCHING FRACTIONS
  85
  86 C---FERMIONS
  87
  88 C---NLL CORRECTION TO QUARK DECAY RATE (HHG eq 2.6-9)
  89
  90       ENF=0
  91
  92       DO 1 I=1,6
  93
  94  1      IF (2*RMASS(I).LT.EMH) ENF=ENF+1
  95
  96       K1=5/PIFAC**2
  97
  98       K0=3/(4*PIFAC**2)
  99
 100       BET0=(11*CAFAC-2*ENF)/3
 101
 102       BET1=(34*CAFAC**2-(10*CAFAC+6*CFFAC)*ENF)/3
 103
 104       GAM0=-8
 105
 106       GAM1=-404./3+40*ENF/9
 107
 108       SCLOG=LOG(EMH**2/QCDLAM**2)
 109
 110       CFAC=1 + ( K1/K0 - 2*GAM0 + GAM0*BET1/BET0**2*LOG(SCLOG)
 111
 112      &       +   (GAM0*BET1-GAM1*BET0)/BET0**2) / (BET0*SCLOG)
 113
 114       DO 100 IFERM=1,9
 115
 116         IF (IFERM.LE.6) THEN
 117
 118           EMF=RMASS(IFERM)
 119
 120           XF=(EMF/EMH)**2
 121
 122           COLFAC=FLOAT(NCOLO)
 123
 124           IF (EMF.GT.QCDLAM)
 125
 126      &      EMF=EMF*(LOG(EMH/QCDLAM)/LOG(EMF/QCDLAM))**(GAM0/(2*BET0))
 127
 128         ELSE
 129
 130           EMF=RMASS(107+IFERM*2)
 131
 132           XF=(EMF/EMH)**2
 133
 134           COLFAC=1
 135
 136           CFAC=1
 137
 138         ENDIF
 139
 140         IF (FOUR*XF.LT.ONE) THEN
 141
 142         GFACTR=ALPHEM/(8.*SWEIN*EMW**2)
 143
 144           BRHIG(IFERM)=COLFAC*GFACTR*EMH*EMF**2 * (1-4*XF)**1.5 * CFAC
 145
 146         ELSE
 147
 148           BRHIG(IFERM)=0
 149
 150         ENDIF
 151
 152  100  CONTINUE
 153
 154 C---W*W*/Z*Z*
 155
 156       IF (ABS(EM-EMH).GE.GAMLIM*GAM) THEN
 157
 158 C---OFF EDGE OF LOOK-UP TABLE
 159
 160         XW=(EMW/EMH)**2
 161
 162         XZ=(EMZ/EMH)**2
 163
 164         YW=EMW*GAMW/EMH**2
 165
 166         YZ=EMZ*GAMZ/EMH**2
 167
 168         BRHIG(10)=.50*GFACTR * EMH**3 * HWDHGF(XW,YW)
 169
 170         BRHIG(11)=.25*GFACTR * EMH**3 * HWDHGF(XZ,YZ)
 171
 172       ELSE
 173
 174 C---LOOK IT UP
 175
 176         EMI=((EMH-EM)/(GAM*GAMLIM)+1)*NLOOK/2.0
 177
 178         I1=INT(EMI)
 179
 180         I2=INT(EMI+1)
 181
 182         BRHIG(10)=.50*GFACTR * EMH**3 * ( VECDEC(1,I1)*(I2-EMI) +
 183
 184      &                                    VECDEC(1,I2)*(EMI-I1) )
 185
 186         BRHIG(11)=.25*GFACTR * EMH**3 * ( VECDEC(2,I1)*(I2-EMI) +
 187
 188      &                                    VECDEC(2,I2)*(EMI-I1) )
 189
 190       ENDIF
 191
 192 C---GAMMAGAMMA
 193
 194       TAUT=(2*RMASS(6)/EMH)**2
 195
 196       TAUW=(2*EMW/EMH)**2
 197
 198       CALL HWDHGC(TAUT,TAUTR,TAUTI)
 199
 200       CALL HWDHGC(TAUW,TAUWR,TAUWI)
 201
 202       SUMR=4./3*(  - 2*TAUT*( 1 + (1-TAUT)*TAUTR ) ) * ENHANC(6)
 203
 204      &         +(2 + 3*TAUW*( 1 + (2-TAUW)*TAUWR ) ) * ENHANC(10)
 205
 206       SUMI=4./3*(  - 2*TAUT*(     (1-TAUT)*TAUTI ) ) * ENHANC(6)
 207
 208      &         +(    3*TAUW*(     (2-TAUW)*TAUWI ) ) * ENHANC(10)
 209
 210       BRHIG(12)=GFACTR*.03125*(ALPHEM/PIFAC)**2
 211
 212      &         *EMH**3 * (SUMR**2 + SUMI**2)
 213
 214       WIDHIG=0
 215
 216       DO 200 IPART=1, 12
 217
 218         IF (IPART.LT.12) BRHIG(IPART)=BRHIG(IPART)*ENHANC(IPART)**2
 219
 220  200    WIDHIG=WIDHIG+BRHIG(IPART)
 221
 222       IF (WIDHIG.EQ.ZERO) CALL HWWARN('HWDHIG',103,*999)
 223
 224       DO 300 IPART=1, 12
 225
 226  300    BRHIG(IPART)=BRHIG(IPART)/WIDHIG
 227
 228       IF (EM.NE.RMASS(201)) THEN
 229
 230 C---SET UP W*W*/Z*Z* LOOKUP TABLES
 231
 232         EM=EMH
 233
 234         GAM=WIDHIG
 235
 236         GAMLIM=MAX(GAMLIM,GAMMAX)
 237
 238         DO 400 I=0,NLOOK
 239
 240           EMH=(I*2.0/NLOOK-1)*GAM*GAMLIM+EM
 241
 242           XW=(EMW/EMH)**2
 243
 244           XZ=(EMZ/EMH)**2
 245
 246           YW=EMW*GAMW/EMH**2
 247
 248           YZ=EMZ*GAMZ/EMH**2
 249
 250           VECDEC(1,I)=HWDHGF(XW,YW)
 251
 252           VECDEC(2,I)=HWDHGF(XZ,YZ)
 253
 254  400    CONTINUE
 255
 256         EMH=EM
 257
 258       ENDIF
 259
 260       IF (GAMINP.GT.ZERO) THEN
 261
 262         GAMINP=WIDHIG
 263
 264         RETURN
 265
 266       ENDIF
 267
 268 C---SEE IF USER SPECIFIED A DECAY MODE
 269
 270       IMODE=MOD(IPROC,100)
 271
 272 C---IF NOT, CHOOSE ONE
 273
 274       IF (IMODE.LT.1.OR.IMODE.GT.12) THEN
 275
 276         MMAX=12
 277
 278         IF (IMODE.LT.1) MMAX=6
 279
 280  500    IMODE=HWRINT(1,MMAX)
 281
 282         IF (BRHIG(IMODE).LT.HWR()) GOTO 500
 283
 284       ENDIF
 285
 286 C---SEE IF SPECIFIED DECAY IS POSSIBLE
 287
 288       IF (BRHIG(IMODE).EQ.ZERO) CALL HWWARN('HWDHIG',104,*999)
 289
 290       IF (IMODE.LE.6) THEN
 291
 292         IDEC=IMODE
 293
 294       ELSEIF (IMODE.LE.9) THEN
 295
 296         IDEC=107+IMODE*2
 297
 298       ELSEIF (IMODE.EQ.10) THEN
 299
 300         IDEC=198
 301
 302       ELSEIF (IMODE.EQ.11) THEN
 303
 304         IDEC=200
 305
 306       ELSEIF (IMODE.EQ.12) THEN
 307
 308         IDEC=59
 309
 310       ENDIF
 311
 312 C---STATUS, IDs AND POINTERS
 313
 314       ISTHEP(IHIG)=195
 315
 316       DO 600 I=1,2
 317
 318         ISTHEP(NHEP+I)=193
 319
 320         IDHW(NHEP+I)=IDEC
 321
 322         IDHEP(NHEP+I)=IDPDG(IDEC)
 323
 324         JDAHEP(I,IHIG)=NHEP+I
 325
 326         JMOHEP(1,NHEP+I)=IHIG
 327
 328         JMOHEP(2,NHEP+I)=NHEP+(3-I)
 329
 330         JDAHEP(2,NHEP+I)=NHEP+(3-I)
 331
 332         PHEP(5,NHEP+I)=RMASS(IDEC)
 333
 334         IDEC=IDEC+6
 335
 336         IF (IDEC.EQ.204) IDEC=199
 337
 338         IF (IDEC.EQ.206) IDEC=200
 339
 340         IF (IDEC.EQ. 65) IDEC= 59
 341
 342  600  CONTINUE
 343
 344 C---ALLOW W/Z TO BE OFF-SHELL
 345
 346       IF (IMODE.EQ.10.OR.IMODE.EQ.11) THEN
 347
 348         IF (IMODE.EQ.10) THEN
 349
 350           EMB=EMW
 351
 352           GAMB=GAMW
 353
 354         ELSE
 355
 356           EMB=EMZ
 357
 358           GAMB=GAMZ
 359
 360         ENDIF
 361
 362 C---STANDARD MASS DISTRIBUTION
 363
 364  700    TMIN=ATAN(-EMB/GAMB)
 365
 366         TMAX1=ATAN((EMH**2/EMB-EMB)/GAMB)
 367
 368         EM1=HWUSQR(EMB*(GAMB*TAN(HWRUNI(0,TMIN,TMAX1))+EMB))
 369
 370         TMAX2=ATAN(((EMH-EM1)**2/EMB-EMB)/GAMB)
 371
 372         EM2=HWUSQR(EMB*(GAMB*TAN(HWRUNI(0,TMIN,TMAX2))+EMB))
 373
 374         X1=(EM1/EMH)**2
 375
 376         X2=(EM2/EMH)**2
 377
 378 C---CORRECT MASS DISTRIBUTION
 379
 380         PROB=HWUSQR(1+X1**2+X2**2-2*X1-2*X2-2*X1*X2)
 381
 382      &        * ((X1+X2-1)**2 + 8*X1*X2)
 383
 384         IF (.NOT.HWRLOG(PROB)) GOTO 700
 385
 386 C---CALCULATE SPIN DENSITY MATRIX
 387
 388         RHOHEP(1,NHEP+1)=4*X1*X2      / (8*X1*X2 + (X1+X2-1)**2)
 389
 390         RHOHEP(2,NHEP+1)=(X1+X2-1)**2 / (8*X1*X2 + (X1+X2-1)**2)
 391
 392         RHOHEP(3,NHEP+1)=RHOHEP(1,NHEP+1)
 393
 394 C---SYMMETRIZE DISTRIBUTIONS IN PARTICLES 1,2
 395
 396         IF (HWRLOG(HALF)) THEN
 397
 398           PHEP(5,NHEP+1)=EM1
 399
 400           PHEP(5,NHEP+2)=EM2
 401
 402         ELSE
 403
 404           PHEP(5,NHEP+1)=EM2
 405
 406           PHEP(5,NHEP+2)=EM1
 407
 408         ENDIF
 409
 410       ENDIF
 411
 412 C---DO DECAY
 413
 414       PCM=HWUPCM(EMH,PHEP(5,NHEP+1),PHEP(5,NHEP+2))
 415
 416       IF (PCM.LT.ZERO) CALL HWWARN('HWDHIG',105,*999)
 417
 418       CALL HWDTWO(PHEP(1,IHIG),PHEP(1,NHEP+1),PHEP(1,NHEP+2),
 419
 420      &            PCM,TWO,.TRUE.)
 421
 422       NHEP=NHEP+2
 423
 424 C---IF QUARK DECAY, HADRONIZE
 425
 426       IF (IMODE.LE.6) THEN
 427
 428         ISTHEP(NHEP-1)=113
 429
 430         ISTHEP(NHEP)=114
 431
 432         CALL HWBGEN
 433
 434         CALL HWDHOB
 435
 436         CALL HWCFOR
 437
 438         CALL HWCDEC
 439
 440       ENDIF
 441
 442   999 END
 443
 444 CDECK  ID>, HWDHOB.
 445
 446 *CMZ :-        -20/10/99  09:46:43  by  Peter Richardson
 447
 448 *-- Author :    Ian Knowles & Bryan Webber
 449
 450 C-----------------------------------------------------------------------
 451
 452       SUBROUTINE HWDHOB
 453
 454 C-----------------------------------------------------------------------
 455
 456 C   Performs decays of heavy objects (heavy quarks & SUSY particles)
 457
 458 C   MODIFIED TO INCLUDE R-PARITY VIOLATING SUSY PR 9/4/99
 459
 460 C-----------------------------------------------------------------------
 461
 462       INCLUDE 'HERWIG61.INC'
 463
 464       DOUBLE PRECISION HWUMBW,HWUPCM,HWR,SDKM,RN,BF,PCM,
 465
 466      & EMMX,EMWSQ,GMWSQ,EMLIM,PW(5),EMTST,HWDPWT,HWDWWT,HWULDO,PDW(5,3)
 467
 468       INTEGER IST(3),IHEP,IS,ID,IM,I,JHEP,KHEP,LHEP,MHEP,NPR,ISM,JCM,
 469
 470      & MTRY,NTRY,IDM,IDM2,THEP,CLSAVE(2),WHEP,RHEP
 471
 472       LOGICAL FOUND
 473
 474       EXTERNAL HWR,HWDPWT,HWDWWT
 475
 476       DATA IST/113,114,114/
 477
 478       IF (IERROR.NE.0) RETURN
 479
 480   10  FOUND=.FALSE.
 481
 482       CLSAVE(1) = 0
 483
 484       CLSAVE(2) = 0
 485
 486       DO 60 IHEP=1,NMXHEP
 487
 488       IS=ISTHEP(IHEP)
 489
 490       ID=IDHW(IHEP)
 491
 492       IF (.NOT.RSTAB(ID).AND.(ID.EQ.6.OR.ID.EQ.12.OR.
 493
 494      & (ID.GE.203.AND.ID.LE.218).OR.ABS(IDPDG(ID)).GT.1000000).AND.
 495
 496      & (IS.EQ.190.OR.(IS.GE.147.AND.IS.LE.151))) THEN
 497
 498         FOUND=.TRUE.
 499
 500         IF(.NOT.RPARTY) THEN
 501
 502           NHEP = NHEP+1
 503
 504           ISTHEP(NHEP) = 3
 505
 506           IDHW(NHEP) = 20
 507
 508           IDHEP(NHEP) = 0
 509
 510           CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP))
 511
 512           CALL HWVEQU(4,VHEP(1,IHEP),VHEP(1,NHEP))
 513
 514           JMOHEP(1,NHEP)=JMOHEP(1,IHEP)
 515
 516           JMOHEP(2,NHEP)=JMOHEP(2,IHEP)
 517
 518           JDAHEP(1,NHEP)=JDAHEP(1,IHEP)
 519
 520           JDAHEP(2,NHEP)=JDAHEP(2,IHEP)
 521
 522         ENDIF
 523
 524 C Make a copy of decaying object
 525
 526         NHEP=NHEP+1
 527
 528         ISTHEP(NHEP)=155
 529
 530         IDHW(NHEP)=IDHW(IHEP)
 531
 532         IDHEP(NHEP)=IDHEP(IHEP)
 533
 534         CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP))
 535
 536         CALL HWVEQU(4,VHEP(1,IHEP),VHEP(1,NHEP))
 537
 538         JMOHEP(1,NHEP)=JMOHEP(1,IHEP)
 539
 540         JMOHEP(2,NHEP)=JMOHEP(2,IHEP)
 541
 542         MTRY=0
 543
 544  15     MTRY=MTRY+1
 545
 546 C Select decay mode
 547
 548         RN=HWR()
 549
 550         BF=0.
 551
 552         IM=LSTRT(ID)
 553
 554         DO 20 I=1,NMODES(ID)
 555
 556         BF=BF+BRFRAC(IM)
 557
 558         IF (BF.GE.RN) GOTO 30
 559
 560   20    IM=LNEXT(IM)
 561
 562         CALL HWWARN('HWDHOB',50,*30)
 563
 564   30    IF (NHEP+5.GT.NMXHEP) CALL HWWARN('HWDHOB',100,*999)
 565
 566         NPR=NPRODS(IM)
 567
 568         JDAHEP(1,NHEP)=NHEP+1
 569
 570         JDAHEP(2,NHEP)=NHEP+NPR
 571
 572 C Reset colour pointers (if set)
 573
 574         JHEP=JMOHEP(2,IHEP)
 575
 576         IF (JHEP.GT.0) THEN
 577
 578           IF (JDAHEP(2,JHEP).EQ.IHEP) JDAHEP(2,JHEP)=NHEP
 579
 580           IF(.NOT.RPARTY.AND.ISTHEP(JHEP).EQ.155
 581
 582      &    .AND.ABS(IDHEP(JHEP)).GT.1000000
 583
 584      &    .AND.JDAHEP(2,JHEP-1).EQ.IHEP) JDAHEP(2,JHEP-1) = NHEP
 585
 586         ENDIF
 587
 588         JHEP=JDAHEP(2,IHEP)
 589
 590         IF (JHEP.GT.0) THEN
 591
 592           IF (JMOHEP(2,JHEP).EQ.IHEP) JMOHEP(2,JHEP)=NHEP
 593
 594           IF(.NOT.RPARTY.AND.ISTHEP(JHEP).EQ.155
 595
 596      &    .AND.ABS(IDHEP(JHEP)).GT.1000000
 597
 598      &    .AND.JMOHEP(2,JHEP-1).EQ.IHEP) JMOHEP(2,JHEP-1) = NHEP
 599
 600         ENDIF
 601
 602 C--Reset colour pointers if baryon number violated
 603
 604         IF(.NOT.RPARTY) THEN
 605
 606           DO JHEP=1,NHEP
 607
 608             IF(ISTHEP(JHEP).EQ.155
 609
 610      &         .AND.ABS(IDHEP(JHEP)).GT.1000000.AND.
 611
 612      &         JDAHEP(2,JHEP-1).EQ.IHEP) JDAHEP(2,JHEP-1)= NHEP
 613
 614             IF(JDAHEP(2,JHEP).EQ.IHEP) JDAHEP(2,JHEP)=NHEP
 615
 616             IF(JMOHEP(2,JHEP).EQ.IHEP) JMOHEP(2,JHEP)=NHEP
 617
 618           ENDDO
 619
 620           IF(HRDCOL(1,1).EQ.IHEP) HRDCOL(1,1)=NHEP
 621
 622         ENDIF
 623
 624 C Relabel original track
 625
 626         ISTHEP(IHEP)=3
 627
 628         JMOHEP(2,IHEP)=JMOHEP(1,IHEP)
 629
 630         JDAHEP(1,IHEP)=NHEP
 631
 632         JDAHEP(2,IHEP)=NHEP
 633
 634 C Label decay products and choose masses
 635
 636         LHEP=NHEP
 637
 638         MHEP=LHEP+1
 639
 640         NTRY=0
 641
 642  35     NTRY=NTRY+1
 643
 644         SDKM=PHEP(5,NHEP)
 645
 646         DO 40 I=1,NPR
 647
 648         NHEP=NHEP+1
 649
 650         IDHW(NHEP)=IDKPRD(I,IM)
 651
 652         IDHEP(NHEP)=IDPDG(IDKPRD(I,IM))
 653
 654         ISTHEP(NHEP)=IST(I)
 655
 656         JMOHEP(1,NHEP)=LHEP
 657
 658         JDAHEP(1,NHEP)=0
 659
 660         PHEP(5,NHEP)=HWUMBW(IDKPRD(I,IM))
 661
 662  40     SDKM=SDKM-PHEP(5,NHEP)
 663
 664         IF (SDKM.LT.ZERO) THEN
 665
 666           NHEP=NHEP-NPR
 667
 668           IF (NTRY.LE.NETRY) GO TO 35
 669
 670           CALL HWWARN('HWDHOB',1,*45)
 671
 672  45       IF (MTRY.LE.NETRY) GO TO 15
 673
 674           CALL HWWARN('HWDHOB',101,*999)
 675
 676         ENDIF
 677
 678 C Assign production vertices to decay products
 679
 680         CALL HWUDKL(ID,PHEP(1,IHEP),VHEP(1,MHEP))
 681
 682         CALL HWVSUM(4,VHEP(1,IHEP),VHEP(1,MHEP),VHEP(1,MHEP))
 683
 684         CALL HWVEQU(4,VHEP(1,MHEP),VHEP(1,NHEP))
 685
 686         IF (NPR.EQ.2) THEN
 687
 688 C Two body decay: LHEP -> MHEP + NHEP
 689
 690           PCM=HWUPCM(PHEP(5,IHEP),PHEP(5,MHEP),PHEP(5,NHEP))
 691
 692           CALL HWDTWO(PHEP(1,IHEP),PHEP(1,MHEP),
 693
 694      &                PHEP(1,NHEP),PCM,TWO,.FALSE.)
 695
 696         ELSEIF (NPR.EQ.3) THEN
 697
 698 C Three body decay: LHEP -> KHEP + MHEP + NHEP
 699
 700           KHEP=MHEP
 701
 702           MHEP=MHEP+1
 703
 704 C Provisional colour self-connection of KHEP
 705
 706           JMOHEP(2,KHEP)=KHEP
 707
 708           JDAHEP(2,KHEP)=KHEP
 709
 710           IF (NME(IM).EQ.100) THEN
 711
 712 C Generate decay momenta using full (V-A)*(V-A) matrix element
 713
 714             EMMX=PHEP(5,IHEP)-PHEP(5,NHEP)
 715
 716             EMWSQ=RMASS(198)**2
 717
 718             GMWSQ=(RMASS(198)*GAMW)**2
 719
 720             EMLIM=GMWSQ
 721
 722             IF (EMMX.LT.RMASS(198)) EMLIM=EMLIM+(EMWSQ-EMMX**2)**2
 723
 724   50        CALL HWDTHR(PHEP(1,IHEP),PHEP(1,MHEP),
 725
 726      &                  PHEP(1,KHEP),PHEP(1,NHEP),HWDWWT)
 727
 728             CALL HWVSUM(4,PHEP(1,KHEP),PHEP(1,MHEP),PW)
 729
 730             PW(5)=HWULDO(PW,PW)
 731
 732             EMTST=(EMWSQ-PW(5))**2
 733
 734             IF ((EMTST+GMWSQ)*HWR().GT.EMLIM) GOTO 50
 735
 736             PW(5)=SQRT(PW(5))
 737
 738 C Assign production vertices to 1 and 2
 739
 740             CALL HWUDKL(198,PW,VHEP(1,KHEP))
 741
 742             CALL HWVSUM(4,VHEP(1,NHEP),VHEP(1,KHEP),VHEP(1,KHEP))
 743
 744           ELSEIF(NME(IM).EQ.300) THEN
 745
 746 C Generate momenta using 3-body RPV matrix element
 747
 748             CALL HWDRME(LHEP,KHEP)
 749
 750           ELSE
 751
 752 C Three body phase space decay
 753
 754             CALL HWDTHR(PHEP(1,IHEP),PHEP(1,MHEP),
 755
 756      &                  PHEP(1,KHEP),PHEP(1,NHEP),HWDPWT)
 757
 758           ENDIF
 759
 760           CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,MHEP))
 761
 762         ELSEIF(NPR.EQ.4) THEN
 763
 764 C Four body decay: LHEP -> KHEP + RHEP + MHEP + NHEP
 765
 766           KHEP = MHEP
 767
 768           RHEP = MHEP+1
 769
 770           MHEP = MHEP+2
 771
 772 C Provisional colour connections of KHEP and RHEP
 773
 774           JMOHEP(2,KHEP)=RHEP
 775
 776           JDAHEP(2,KHEP)=RHEP
 777
 778           JMOHEP(2,RHEP)=KHEP
 779
 780           JDAHEP(2,RHEP)=KHEP
 781
 782 C Four body phase space decay
 783
 784           CALL HWDFOR(PHEP(1,IHEP),PHEP(1,KHEP),PHEP(1,RHEP),
 785
 786      &                  PHEP(1,MHEP),PHEP(1,NHEP))
 787
 788           CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,RHEP))
 789
 790           CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,MHEP))
 791
 792         ELSE
 793
 794           CALL HWWARN('HWDHOB',102,*999)
 795
 796         ENDIF
 797
 798 C Colour connections
 799
 800         IF (ID.EQ.6.OR.ID.EQ.12.OR.(ID.GE.209.AND.ID.LE.212)
 801
 802      &                         .OR.(ID.GE.215.AND.ID.LE.218)) THEN
 803
 804           IF (NPR.EQ.3.AND.NME(IM).EQ.100) THEN
 805
 806 C usual heavy quark decay
 807
 808             JMOHEP(2,KHEP)=MHEP
 809
 810             JDAHEP(2,KHEP)=MHEP
 811
 812             JMOHEP(2,MHEP)=KHEP
 813
 814             JDAHEP(2,MHEP)=KHEP
 815
 816             JMOHEP(2,NHEP)=LHEP
 817
 818             JDAHEP(2,NHEP)=LHEP
 819
 820           ELSEIF (ABS(IDHEP(MHEP)).EQ.37) THEN
 821
 822 C heavy quark to charged Higgs
 823
 824             JMOHEP(2,MHEP)=MHEP
 825
 826             JDAHEP(2,MHEP)=MHEP
 827
 828             JMOHEP(2,NHEP)=LHEP
 829
 830             JDAHEP(2,NHEP)=LHEP
 831
 832           ELSEIF (ABS(IDHEP(NHEP)).EQ.37) THEN
 833
 834             JMOHEP(2,MHEP)=LHEP
 835
 836             JDAHEP(2,MHEP)=LHEP
 837
 838             JMOHEP(2,NHEP)=NHEP
 839
 840             JDAHEP(2,NHEP)=NHEP
 841
 842           ELSE
 843
 844             CALL HWWARN('HWDHOB',103,*999)
 845
 846           ENDIF
 847
 848         ELSE
 849
 850           IF(.NOT.RPARTY.AND.
 851
 852      &       ((NPR.EQ.2.AND.ID.GE.401.AND.ID.LT.448.AND.
 853
 854      &           IDHW(MHEP).LE.132.AND.IDHW(NHEP).LE.132)
 855
 856      &       .OR.(NPR.EQ.3.AND.ID.GE.449.AND.ID.LE.457.AND.
 857
 858      &           IDHW(MHEP).LE.132.AND.IDHW(NHEP).LE.132.AND.
 859
 860      &           IDHW(MHEP-1).LE.132))) THEN
 861
 862 C R-parity violating SUSY decays
 863
 864             IF(NPR.EQ.2) THEN
 865
 866 C--Rparity slepton colour connections
 867
 868               IF(ID.GE.425.AND.ID.LE.448) THEN
 869
 870                 IF(IDHW(MHEP).GT.12) THEN
 871
 872                   JMOHEP(2,MHEP) = MHEP
 873
 874                   JDAHEP(2,MHEP) = MHEP
 875
 876                   JMOHEP(2,NHEP) = NHEP
 877
 878                   JDAHEP(2,NHEP) = NHEP
 879
 880                 ELSE
 881
 882                   JMOHEP(2,MHEP) = NHEP
 883
 884                   JDAHEP(2,MHEP) = NHEP
 885
 886                   JMOHEP(2,NHEP) = MHEP
 887
 888                   JDAHEP(2,NHEP) = MHEP
 889
 890                 ENDIF
 891
 892 C--Rparity squark colour connections
 893
 894               ELSE
 895
 896                 IF(IDHEP(LHEP).GT.0) THEN
 897
 898 C--LQD decay colour connections
 899
 900                   IF(IDHW(MHEP).GT.12) THEN
 901
 902                     JMOHEP(2,MHEP) = MHEP
 903
 904                     JDAHEP(2,MHEP) = MHEP
 905
 906                     JMOHEP(2,NHEP) = LHEP
 907
 908                     JDAHEP(2,NHEP) = LHEP
 909
 910                   ELSE
 911
 912 C--UDD decay colour connections
 913
 914                     HVFCEN = .TRUE.
 915
 916                     CALL HWDRCL(LHEP,MHEP,CLSAVE)
 917
 918                   ENDIF
 919
 920                 ELSE
 921
 922 C--Antisquark connections
 923
 924                   IF(IDHW(MHEP).GT.12) THEN
 925
 926                     JMOHEP(2,MHEP) = MHEP
 927
 928                     JDAHEP(2,MHEP) = MHEP
 929
 930                     JMOHEP(2,NHEP) = LHEP
 931
 932                    JDAHEP(2,NHEP) = LHEP
 933
 934                   ELSE
 935
 936                     HVFCEN = .TRUE.
 937
 938                    CALL HWDRCL(LHEP,MHEP,CLSAVE)
 939
 940                   ENDIF
 941
 942                 ENDIF
 943
 944               ENDIF
 945
 946             ELSE
 947
 948               IF(ID.GE.450.AND.ID.LE.457) THEN
 949
 950 C--Rparity Neutralino/Chargino colour connection
 951
 952                 IF(IDHW(MHEP-1).LE.12.AND.IDHW(MHEP).LE.12.
 953
 954      &                 AND.IDHW(NHEP).LE.12) THEN
 955
 956                   HVFCEN = .TRUE.
 957
 958                   CALL HWDRCL(LHEP,MHEP,CLSAVE)
 959
 960                 ELSE
 961
 962                   JMOHEP(2,MHEP) = NHEP
 963
 964                   JDAHEP(2,MHEP) = NHEP
 965
 966                   JMOHEP(2,NHEP) = MHEP
 967
 968                   JDAHEP(2,NHEP) = MHEP
 969
 970                 ENDIF
 971
 972 C--Rparity gluino colour connections
 973
 974               ELSEIF(ID.EQ.449) THEN
 975
 976                 IF(IDHW(MHEP-1).LE.12.AND.IDHW(MHEP).LE.12.
 977
 978      &                 AND.IDHW(NHEP).LE.12) THEN
 979
 980                   HVFCEN = .TRUE.
 981
 982                   CALL HWDRCL(LHEP,MHEP,CLSAVE)
 983
 984 C--Now the lepton number violating decay
 985
 986                 ELSE
 987
 988                   IF(IDHW(MHEP).LE.6) THEN
 989
 990                     JMOHEP(2,MHEP) = LHEP
 991
 992                     JDAHEP(2,MHEP) = NHEP
 993
 994                     JMOHEP(2,NHEP) = MHEP
 995
 996                     JDAHEP(2,NHEP) = LHEP
 997
 998                   ELSE
 999
1000                     JMOHEP(2,MHEP) = NHEP
1001
1002                     JDAHEP(2,MHEP) = LHEP
1003
1004                     JMOHEP(2,NHEP) = LHEP
1005
1006                     JDAHEP(2,NHEP) = MHEP
1007
1008                   ENDIF
1009
1010                 ENDIF
1011
1012               ELSE
1013
1014                 CALL HWWARN('HWDHOB',104,*999)
1015
1016               ENDIF
1017
1018             ENDIF
1019
1020           ELSE
1021
1022 C Normal SUSY decays
1023
1024             IF (ID.LE.448.AND.ID.GT.207) THEN
1025
1026 C Squark (or slepton)
1027
1028               IF (IDHW(MHEP).EQ.449) THEN
1029
1030                 IF (IDHEP(LHEP).GT.0) THEN
1031
1032                   JMOHEP(2,MHEP)=LHEP
1033
1034                   JDAHEP(2,MHEP)=NHEP
1035
1036                   JMOHEP(2,NHEP)=MHEP
1037
1038                   JDAHEP(2,NHEP)=LHEP
1039
1040                 ELSE
1041
1042                   JMOHEP(2,MHEP)=NHEP
1043
1044                   JDAHEP(2,MHEP)=LHEP
1045
1046                   JMOHEP(2,NHEP)=LHEP
1047
1048                   JDAHEP(2,NHEP)=MHEP
1049
1050                 ENDIF
1051
1052               ELSE
1053
1054                 IF(NPR.EQ.3.AND.IDHW(MHEP).LE.12) THEN
1055
1056                   JMOHEP(2,MHEP)=NHEP
1057
1058                   JDAHEP(2,MHEP)=NHEP
1059
1060                   JMOHEP(2,NHEP)=MHEP
1061
1062                   JDAHEP(2,NHEP)=MHEP
1063
1064                 ELSE
1065
1066                   JMOHEP(2,MHEP)=MHEP
1067
1068                   JDAHEP(2,MHEP)=MHEP
1069
1070                   JMOHEP(2,NHEP)=LHEP
1071
1072                   JDAHEP(2,NHEP)=LHEP
1073
1074                 ENDIF
1075
1076               ENDIF
1077
1078             ELSEIF (ID.EQ.449) THEN
1079
1080 C Gluino
1081
1082               IF (IDHW(NHEP).EQ.13) THEN
1083
1084                 JMOHEP(2,MHEP)=MHEP
1085
1086                 JDAHEP(2,MHEP)=MHEP
1087
1088                 JMOHEP(2,NHEP)=LHEP
1089
1090                 JDAHEP(2,NHEP)=LHEP
1091
1092               ELSEIF (IDHEP(MHEP).GT.0) THEN
1093
1094                 JMOHEP(2,MHEP)=LHEP
1095
1096                 JDAHEP(2,MHEP)=NHEP
1097
1098                 JMOHEP(2,NHEP)=MHEP
1099
1100                 JDAHEP(2,NHEP)=LHEP
1101
1102               ELSE
1103
1104                 JMOHEP(2,MHEP)=NHEP
1105
1106                 JDAHEP(2,MHEP)=LHEP
1107
1108                 JMOHEP(2,NHEP)=LHEP
1109
1110                 JDAHEP(2,NHEP)=MHEP
1111
1112               ENDIF
1113
1114             ELSE
1115
1116 C Gaugino or Higgs
1117
1118               JMOHEP(2,MHEP)=NHEP
1119
1120               JDAHEP(2,MHEP)=NHEP
1121
1122               JMOHEP(2,NHEP)=MHEP
1123
1124               JDAHEP(2,NHEP)=MHEP
1125
1126             ENDIF
1127
1128           ENDIF
1129
1130         ENDIF
1131
1132 C---SPECIAL CASE FOR THREE-BODY TOP DECAYS:
1133
1134 C   RELABEL THEM AS TWO TWO-BODY DECAYS FOR PARTON SHOWERING
1135
1136         IF ((ID.EQ.6.OR.ID.EQ.12).AND.NPR.EQ.3.AND.NME(IM).EQ.100) THEN
1137
1138 C---STORE W DECAY PRODUCTS
1139
1140           CALL HWVEQU(10,PHEP(1,KHEP),PDW)
1141
1142 C---BOOST THEM INTO W REST FRAME
1143
1144           CALL HWULOF(PW,PDW(1,1),PDW(1,3))
1145
1146 C---REPLACE THEM BY W
1147
1148           CALL HWVEQU(5,PW,PHEP(1,KHEP))
1149
1150           WHEP=KHEP
1151
1152           IDHW(KHEP)=198
1153
1154           IF (ID.EQ.12) IDHW(KHEP)=199
1155
1156           IDHEP(KHEP)=IDPDG(IDHW(KHEP))
1157
1158           JMOHEP(2,KHEP)=KHEP
1159
1160           JDAHEP(2,KHEP)=KHEP
1161
1162           CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,KHEP))
1163
1164 C---AND MOVE B UP
1165
1166           CALL HWVEQU(5,PHEP(1,NHEP),PHEP(1,MHEP))
1167
1168           IDHW(MHEP)=IDHW(NHEP)
1169
1170           IDHEP(MHEP)=IDHEP(NHEP)
1171
1172           JDAHEP(2,LHEP)=MHEP
1173
1174           JMOHEP(2,MHEP)=JMOHEP(2,NHEP)
1175
1176           JDAHEP(2,MHEP)=JDAHEP(2,NHEP)
1177
1178           CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,MHEP))
1179
1180           NHEP=MHEP
1181
1182 C---DO PARTON SHOWER
1183
1184           EMSCA=PHEP(5,IHEP)
1185
1186           CALL HWBGEN
1187
1188           IF (IERROR.NE.0) RETURN
1189
1190 C---FIND BOOSTED W MOMENTUM
1191
1192           NTRY=0
1193
1194  41       NTRY=NTRY+1
1195
1196           IF (NTRY.GT.NHEP.OR.WHEP.LE.0.OR.WHEP.GT.NHEP)
1197
1198      $         CALL HWWARN('HWDHOB',101,*999)
1199
1200           WHEP=JDAHEP(1,WHEP)
1201
1202           IF (ISTHEP(WHEP).NE.190) GOTO 41
1203
1204 C---AND HENCE ITS CHILDRENS MOMENTA
1205
1206           CALL HWULOB(PHEP(1,WHEP),PDW(1,3),PHEP(1,NHEP+1))
1207
1208           CALL HWVDIF(4,PHEP(1,WHEP),PHEP(1,NHEP+1),PHEP(1,NHEP+2))
1209
1210           PHEP(5,NHEP+2)=PDW(5,2)
1211
1212 C---LABEL THEM
1213
1214           ISTHEP(WHEP)=195
1215
1216           DO 51 I=1,2
1217
1218             IDHW(NHEP+I)=IDKPRD(I,IM)
1219
1220             IDHEP(NHEP+I)=IDPDG(IDHW(NHEP+I))
1221
1222             ISTHEP(NHEP+I)=112+I
1223
1224             JDAHEP(I,WHEP)=NHEP+I
1225
1226             JMOHEP(1,NHEP+I)=WHEP
1227
1228             JMOHEP(2,NHEP+I)=NHEP+3-I
1229
1230             JDAHEP(2,NHEP+I)=NHEP+3-I
1231
1232  51       CONTINUE
1233
1234           NHEP=NHEP+2
1235
1236 C---ASSIGN PRODUCTION VERTICES TO 1 AND 2
1237
1238           CALL HWUDKL(198,PW,VHEP(1,NHEP))
1239
1240           CALL HWVSUM(4,VHEP(1,WHEP),VHEP(1,NHEP),VHEP(1,NHEP))
1241
1242           CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-1))
1243
1244 C---DO PARTON SHOWERS
1245
1246           EMSCA=PW(5)
1247
1248           CALL HWBGEN
1249
1250           IF (IERROR.NE.0) RETURN
1251
1252         ELSE
1253
1254 C Do parton showers
1255
1256           EMSCA=PHEP(5,IHEP)
1257
1258           CALL HWBGEN
1259
1260           IF (IERROR.NE.0) RETURN
1261
1262         ENDIF
1263
1264       ENDIF
1265
1266 C--New to correct colour connections in Rslash
1267
1268       IF(CLSAVE(1).NE.0) THEN
1269
1270         THEP = MHEP+1
1271
1272         ID   = IDHW(CLSAVE(1))
1273
1274         IDM  = IDHW(JMOHEP(1,CLSAVE(1)))
1275
1276         IDM2 = IDHW(LHEP)
1277
1278         IF(IDM.EQ.15) ID=IDHW(JMOHEP(1,JMOHEP(1,CLSAVE(1))))
1279
1280         IF((ID.LE.6.AND.((IDM.GE.419.AND.IDM.LE.424).OR.IDM.EQ.411.OR.
1281
1282      &      IDM.EQ.412).
1283
1284      &     AND.((IDM2.GE.413.AND.IDM2.LE.418)
1285
1286      &     .OR.IDM2.EQ.449).OR.IDM2.EQ.405.OR.IDM2.EQ.406)
1287
1288      &     .OR.(ID.LE.6.AND.IDM.EQ.449.AND.
1289
1290      &    (((IDM2.GE.413.AND.IDM2.LE.418).OR.IDM2.EQ.405.OR.IDM2.EQ.406)
1291
1292      &     .OR.IDM2.EQ.449)).OR.
1293
1294      &    (IDM.EQ.15.AND.ID.LE.12.AND.ID.GE.7.AND.((IDM2.GE.413.AND.
1295
1296      &     IDM2.LE.418).OR.IDM2.EQ.449.OR.IDM2.
1297
1298      &     EQ.405.OR.IDM2.EQ.406))) THEN
1299
1300           IF(JMOHEP(2,CLSAVE(1)).EQ.MHEP) THEN
1301
1302             IF(IDHW(CLSAVE(1)).NE.13.AND.IDHW(CLSAVE(1)).NE.449)
1303
1304      &                       JMOHEP(2,CLSAVE(2)) = THEP
1305
1306             JDAHEP(2,MHEP) = CLSAVE(1)
1307
1308             JDAHEP(2,THEP) = CLSAVE(2)
1309
1310           ELSE
1311
1312             IF(IDHW(CLSAVE(1)).NE.13.AND.IDHW(CLSAVE(1)).NE.449)
1313
1314      &                       JMOHEP(2,CLSAVE(2)) = MHEP
1315
1316             JDAHEP(2,MHEP) = CLSAVE(2)
1317
1318             JDAHEP(2,THEP) = CLSAVE(1)
1319
1320           ENDIF
1321
1322         ELSEIF((ID.GT.6.AND.ID.LE.12.
1323
1324      &     AND.((IDM.GE.413.AND.IDM.LE.418).OR.IDM.EQ.405.OR.
1325
1326      &     IDM.EQ.406).AND.
1327
1328      &      ((IDM2.GE.419.AND.IDM2.LE.424).OR.IDM2.EQ.449.OR.
1329
1330      &      IDM2.EQ.411.OR.IDM2.EQ.412)).OR.
1331
1332      &        (ID.GT.6.AND.ID.LE.12.AND.IDM.EQ.449.
1333
1334      &   AND.((IDM2.GE.419.AND.IDM2.LE.424).OR.IDM2.EQ.449.OR.
1335
1336      &       IDM2.EQ.411.OR.IDM2.EQ.412)).OR.
1337
1338      &    (IDM.EQ.15.AND.ID.LE.6.AND.((IDM2.GE.419.AND.
1339
1340      &     IDM2.LE.424).OR.IDM2.EQ.449.OR.IDM2.EQ.411.OR.
1341
1342      &     IDM2.EQ.412))) THEN
1343
1344           IF(JDAHEP(2,CLSAVE(1)).EQ.MHEP) THEN
1345
1346             JDAHEP(2,CLSAVE(2))=THEP
1347
1348             JMOHEP(2,MHEP)=CLSAVE(1)
1349
1350             JMOHEP(2,THEP)=CLSAVE(2)
1351
1352           ELSE
1353
1354             JDAHEP(2,CLSAVE(2))=MHEP
1355
1356             JMOHEP(2,MHEP)=CLSAVE(2)
1357
1358             JMOHEP(2,THEP)=CLSAVE(1)
1359
1360           ENDIF
1361
1362         ENDIF
1363
1364         COLUPD = .FALSE.
1365
1366         CALL HWBCON
1367
1368       ENDIF
1369
1370       IF (IHEP.EQ.NHEP) GOTO 70
1371
1372   60  CONTINUE
1373
1374   70  IF (FOUND) THEN
1375
1376 C Fix any SUSY colour disconnections
1377
1378         DO 80 IHEP=1,NHEP
1379
1380           IF (ISTHEP(IHEP).GE.147.AND.ISTHEP(IHEP).LE.151
1381
1382      &      .AND.JDAHEP(2,IHEP).EQ.0) THEN
1383
1384             IM=JMOHEP(1,IHEP)
1385
1386 C Chase connection back through SUSY decays
1387
1388   75        IM=JMOHEP(1,IM)
1389
1390             ISM=ISTHEP(IM)
1391
1392             IF (ISM.EQ.120) GOTO 80
1393
1394             IF (ISM.NE.123.AND.ISM.NE.124.AND.ISM.NE.155) GOTO 75
1395
1396 C Look for unclustered parton to connect
1397
1398             DO JHEP=1,NHEP
1399
1400               IF (ISTHEP(JHEP).GE.147.AND.ISTHEP(JHEP).LE.151) THEN
1401
1402                 JCM=JMOHEP(2,JHEP)
1403
1404                 IF (JCM.EQ.IM) THEN
1405
1406 C Found it: connect
1407
1408                   JMOHEP(2,JHEP)=IHEP
1409
1410                   JDAHEP(2,IHEP)=JHEP
1411
1412                   GOTO 80
1413
1414                 ENDIF
1415
1416               ENDIF
1417
1418             ENDDO
1419
1420 C Not found: need to go further back
1421
1422             GOTO 75
1423
1424           ENDIF
1425
1426    80   CONTINUE
1427
1428 C Go back to check for further heavy decay products
1429
1430         GOTO 10
1431
1432       ENDIF
1433
1434   999 END
1435
1436 CDECK  ID>, HWDHVY.
1437
1438 *CMZ :-        -26/04/91  12.19.24  by  Federico Carminati
1439
1440 *-- Author :    Ian Knowles & Bryan Webber
1441
1442 C-----------------------------------------------------------------------
1443
1444       SUBROUTINE HWDHVY
1445
1446 C-----------------------------------------------------------------------
1447
1448 C     Performs partonic decays of hadrons containing heavy quark(s):
1449
1450 C     either, meson/baryon spectator model weak decays;
1451
1452 C     or, quarkonia -> 2-gluons, q-qbar, 3-gluons, or 2-gluons + photon.
1453
1454 C-----------------------------------------------------------------------
1455
1456       INCLUDE 'HERWIG61.INC'
1457
1458       DOUBLE PRECISION HWULDO,HWR,XS,XB,EMWSQ,GMWSQ,EMLIM,PW(4),
1459
1460      & EMTST,X1,X2,X3,TEST,HWDWWT,HWDPWT
1461
1462       INTEGER IST(3),I,IHEP,IM,ID,IDQ,IQ,IS,J
1463
1464       EXTERNAL HWR,HWDWWT,HWDPWT,HWULDO
1465
1466       DATA IST/113,114,114/
1467
1468       IF (IERROR.NE.0) RETURN
1469
1470       DO 100 I=1,NMXQDK
1471
1472       IF (I.GT.NQDK) THEN
1473
1474         NQDK=0
1475
1476         RETURN
1477
1478       ENDIF
1479
1480       IHEP=LOCQ(I)
1481
1482       IF (ISTHEP(IHEP).EQ.199) GOTO 100
1483
1484       IM=IMQDK(I)
1485
1486       IF (NHEP+NPRODS(IM).GT.NMXHEP) CALL HWWARN('HWDHVY',100,*999)
1487
1488       IF (IDKPRD(4,IM).NE.0) THEN
1489
1490 C Weak decay of meson or baryon
1491
1492 C Idenitify decaying heavy quark and spectator
1493
1494         ID=IDHW(IHEP)
1495
1496         IF (ID.EQ.136.OR.ID.EQ.140.OR.ID.EQ.144.OR.
1497
1498      &      ID.EQ.150.OR.ID.EQ.155.OR.ID.EQ.158.OR.ID.EQ.161.OR.
1499
1500      &     (ID.EQ.254.AND.IDKPRD(4,IM).EQ.11)) THEN
1501
1502 C c hadron or c decay of B_c+
1503
1504           IDQ=4
1505
1506           IQ=NHEP+1
1507
1508           IS=NHEP+2
1509
1510         ELSEIF (ID.EQ.171.OR.ID.EQ.175.OR.ID.EQ.179.OR.
1511
1512      &          ID.EQ.185.OR.ID.EQ.190.OR.ID.EQ.194.OR.ID.EQ.196.OR.
1513
1514      &         (ID.EQ.230.AND.IDKPRD(4,IM).EQ.5)) THEN
1515
1516 C cbar hadron or cbar decay of B_c-
1517
1518           IDQ=10
1519
1520           IS=NHEP+1
1521
1522           IQ=NHEP+2
1523
1524         ELSEIF ((ID.GE.221.AND.ID.LE.229).OR.
1525
1526      &          (ID.EQ.230.AND.IDKPRD(4,IM).EQ.10)) THEN
1527
1528 C b hadron or b decay of B_c-
1529
1530           IDQ=5
1531
1532           IQ=NHEP+1
1533
1534           IS=NHEP+2
1535
1536         ELSEIF ((ID.GE.245.AND.ID.LE.253).OR.
1537
1538      &          (ID.EQ.254.AND.IDKPRD(4,IM).EQ.4)) THEN
1539
1540 C bbar hadron or bbar decay of B_c+
1541
1542           IDQ=11
1543
1544           IS=NHEP+1
1545
1546           IQ=NHEP+2
1547
1548         ELSE
1549
1550 C Decay not recognized
1551
1552           CALL HWWARN('HWDHVY',101,*999)
1553
1554         ENDIF
1555
1556 C Label constituents
1557
1558         IF (NHEP+5.GT.NMXHEP) CALL HWWARN('HWDHVY',102,*999)
1559
1560         ISTHEP(IHEP)=199
1561
1562         JDAHEP(1,IHEP)=NHEP+1
1563
1564         JDAHEP(2,IHEP)=NHEP+2
1565
1566         IDHW(IQ)=IDQ
1567
1568         IDHW(IS)=IDKPRD(4,IM)
1569
1570         IDHEP(IQ)=IDPDG(IDQ)
1571
1572         IDHEP(IS)=IDPDG(IDKPRD(4,IM))
1573
1574         ISTHEP(IQ)=155
1575
1576         ISTHEP(IS)=115
1577
1578         JMOHEP(1,IQ)=IHEP
1579
1580         JMOHEP(2,IQ)=IS
1581
1582         JDAHEP(1,IQ)=NHEP+3
1583
1584         JDAHEP(2,IQ)=NHEP+5
1585
1586         JMOHEP(1,IS)=IHEP
1587
1588         JMOHEP(2,IS)=NHEP+5
1589
1590         JDAHEP(1,IS)=0
1591
1592         JDAHEP(2,IS)=NHEP+5
1593
1594         NHEP=NHEP+2
1595
1596 C and weak decay product jets
1597
1598         DO 10 J=1,3
1599
1600         NHEP=NHEP+1
1601
1602         IDHW(NHEP)=IDKPRD(J,IM)
1603
1604         IDHEP(NHEP)=IDPDG(IDKPRD(J,IM))
1605
1606         ISTHEP(NHEP)=IST(J)
1607
1608         JMOHEP(1,NHEP)=IQ
1609
1610         JDAHEP(1,NHEP)=0
1611
1612   10    PHEP(5,NHEP)=RMASS(IDKPRD(J,IM))
1613
1614         JMOHEP(2,NHEP-2)=NHEP-1
1615
1616         JDAHEP(2,NHEP-2)=NHEP-1
1617
1618         JMOHEP(2,NHEP-1)=NHEP-2
1619
1620         JDAHEP(2,NHEP-1)=NHEP-2
1621
1622         JMOHEP(2,NHEP  )=IQ
1623
1624         JDAHEP(2,NHEP  )=IQ
1625
1626 C Share momenta in ratio of masses, preserving specator mass
1627
1628         XS=RMASS(IDHW(IS))/PHEP(5,IHEP)
1629
1630         XB=ONE-XS
1631
1632         CALL HWVSCA(5,XB,PHEP(1,IHEP),PHEP(1,IQ))
1633
1634         CALL HWVSCA(5,XS,PHEP(1,IHEP),PHEP(1,IS))
1635
1636         IF (NME(IM).EQ.100) THEN
1637
1638 C Generate decay momenta using full (V-A)*(V-A) matrix element
1639
1640           EMWSQ=RMASS(198)**2
1641
1642           GMWSQ=(RMASS(198)*GAMW)**2
1643
1644           EMLIM=GMWSQ+(EMWSQ-(PHEP(5,IQ)-PHEP(5,NHEP))**2)**2
1645
1646   20      CALL HWDTHR(PHEP(1,IQ  ),PHEP(1,NHEP-1),
1647
1648      &                PHEP(1,NHEP-2),PHEP(1,NHEP),HWDWWT)
1649
1650           CALL HWVSUM(4,PHEP(1,NHEP-2),PHEP(1,NHEP-1),PW)
1651
1652           EMTST=(HWULDO(PW,PW)-EMWSQ)**2
1653
1654           IF ((EMTST+GMWSQ)*HWR().GT.EMLIM) GOTO 20
1655
1656         ELSE
1657
1658 C Use phase space
1659
1660           CALL HWDTHR(PHEP(1,IQ  ),PHEP(1,NHEP-2),
1661
1662      &                PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
1663
1664           CALL HWVSUM(4,PHEP(1,NHEP-2),PHEP(1,NHEP-1),PW)
1665
1666         ENDIF
1667
1668 C Set up production vertices
1669
1670         CALL HWVZRO(4,VHEP(1,IQ))
1671
1672         CALL HWVEQU(4,VHEP(1,IQ),VHEP(1,IS))
1673
1674         CALL HWVEQU(4,VHEP(1,IQ),VHEP(1,NHEP))
1675
1676         CALL HWUDKL(198,PW,VHEP(1,NHEP-2))
1677
1678         CALL HWVSUM(4,VHEP(1,IQ),VHEP(1,NHEP-2),VHEP(1,NHEP-2))
1679
1680         CALL HWVEQU(4,VHEP(1,NHEP-2),VHEP(1,NHEP-1))
1681
1682         EMSCA=PHEP(5,IQ)
1683
1684       ELSE
1685
1686 C Quarkonium decay
1687
1688 C Label products
1689
1690         ISTHEP(IHEP)=199
1691
1692         JDAHEP(1,IHEP)=NHEP+1
1693
1694         DO 30 J=1,NPRODS(IM)
1695
1696         NHEP=NHEP+1
1697
1698         IDHW(NHEP)=IDKPRD(J,IM)
1699
1700         IDHEP(NHEP)=IDPDG(IDKPRD(J,IM))
1701
1702         ISTHEP(NHEP)=IST(J)
1703
1704         JMOHEP(1,NHEP)=IHEP
1705
1706         JDAHEP(1,NHEP)=0
1707
1708         PHEP(5,NHEP)=RMASS(IDKPRD(J,IM))
1709
1710   30    CALL HWVZRO(4,VHEP(1,NHEP))
1711
1712         JDAHEP(2,IHEP)=NHEP
1713
1714 C Establish colour connections and select momentum configuration
1715
1716         IF (NPRODS(IM).EQ.3) THEN
1717
1718           IF (IDKPRD(3,IM).EQ.13) THEN
1719
1720 C 3-gluon decay
1721
1722             JMOHEP(2,NHEP-2)=NHEP
1723
1724             JMOHEP(2,NHEP-1)=NHEP-2
1725
1726             JMOHEP(2,NHEP  )=NHEP-1
1727
1728             JDAHEP(2,NHEP-2)=NHEP-1
1729
1730             JDAHEP(2,NHEP-1)=NHEP
1731
1732             JDAHEP(2,NHEP  )=NHEP-2
1733
1734           ELSE
1735
1736 C or 2-gluon + photon decay
1737
1738             JMOHEP(2,NHEP-2)=NHEP-1
1739
1740             JMOHEP(2,NHEP-1)=NHEP-2
1741
1742             JMOHEP(2,NHEP  )=NHEP
1743
1744             JDAHEP(2,NHEP-2)=NHEP-1
1745
1746             JDAHEP(2,NHEP-1)=NHEP-2
1747
1748             JDAHEP(2,NHEP  )=NHEP
1749
1750           ENDIF
1751
1752           IF (NME(IM).EQ.130) THEN
1753
1754 C Use Ore & Powell orthopositronium matrix element
1755
1756   40        CALL HWDTHR(PHEP(1,IHEP),PHEP(1,NHEP-2),
1757
1758      &                               PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
1759
1760             X1=TWO*HWULDO(PHEP(1,IHEP),PHEP(1,NHEP-2))/PHEP(5,IHEP)**2
1761
1762             X2=TWO*HWULDO(PHEP(1,IHEP),PHEP(1,NHEP-1))/PHEP(5,IHEP)**2
1763
1764             X3=TWO-X1-X2
1765
1766             TEST=((X1*(ONE-X1))**2+(X2*(ONE-X2))**2+(X3*(ONE-X3))**2)
1767
1768      &          /(X1*X2*X3)**2
1769
1770             IF (TEST.LT.TWO*HWR()) GOTO 40
1771
1772           ELSE
1773
1774 C Use phase space
1775
1776             CALL HWDTHR(PHEP(1,IHEP),PHEP(1,NHEP-2),
1777
1778      &                               PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
1779
1780           ENDIF
1781
1782         ELSE
1783
1784 C Parapositronium 2-gluon or q-qbar decay
1785
1786           JMOHEP(2,NHEP-1)=NHEP
1787
1788           JMOHEP(2,NHEP  )=NHEP-1
1789
1790           JDAHEP(2,NHEP-1)=NHEP
1791
1792           JDAHEP(2,NHEP  )=NHEP-1
1793
1794           CALL HWDTWO(PHEP(1,IHEP),PHEP(1,NHEP-1),
1795
1796      &                             PHEP(1,NHEP),CMMOM(IM),TWO,.FALSE.)
1797
1798         ENDIF
1799
1800         EMSCA=PHEP(5,IHEP)
1801
1802       ENDIF
1803
1804 C Process this new hard scatter
1805
1806       CALL HWVEQU(4,VTXQDK(1,I),VTXPIP)
1807
1808       CALL HWBGEN
1809
1810       CALL HWCFOR
1811
1812       CALL HWCDEC
1813
1814       CALL HWDHAD
1815
1816   100 CONTINUE
1817
1818       NQDK=0
1819
1820   999 END
1821
1822 CDECK  ID>, HWDRCL.
1823
1824 *CMZ :-        -20/07/99  10:56:12  by  Peter Richardson
1825
1826 *-- Author :    Peter Richardson
1827
1828 C-----------------------------------------------------------------------
1829
1830       SUBROUTINE HWDRCL(IHEP,MHEP,CLSAVE)
1831
1832 C-----------------------------------------------------------------------
1833
1834 C     Sets the colour connections in Baryon number violating decays
1835
1836 C-----------------------------------------------------------------------
1837
1838       INCLUDE 'HERWIG61.INC'
1839
1840       INTEGER IHEP,MHEP,ID,ID2,IDM2,IDM3,COLCON(2,2,3),FLACON(2,3),JHEP,
1841
1842      &        DECAY,COLANT,KHEP,IDM,IDMB,IDMB2,IDMB3,IDMB4,QHEP,IDM4,
1843
1844      &        CLSAVE(2),XHEP,I,HWRINT,THEP
1845
1846       LOGICAL CONBV
1847
1848 C--Colour connections for the decays
1849
1850       DATA COLCON/-1,1,-1,-2,-2,1,-3,-1,-1,1,-2,-1/
1851
1852       DATA FLACON/1,-1,1,-1,-1,0/
1853
1854 C--identify the decay
1855
1856       IF(IERROR.NE.0) RETURN
1857
1858       ID = IDHW(IHEP)
1859
1860       ID2 = IDHW(MHEP)
1861
1862       IF(ID.GE.450.AND.ID.LE.457) THEN
1863
1864         DECAY = 1
1865
1866       ELSEIF(ID.EQ.449) THEN
1867
1868         DECAY = 2
1869
1870       ELSEIF((ID.GE.411.AND.ID.LE.424).OR.ID.EQ.405.OR.ID.EQ.406) THEN
1871
1872         DECAY = 3
1873
1874       ELSE
1875
1876 C--UNKNOWN DECAY
1877
1878         CALL HWWARN('HWDRCL',100,*999)
1879
1880       ENDIF
1881
1882       COLANT = 1
1883
1884 C--identify the colour partner
1885
1886       IF(DECAY.GT.1.AND.ID2.LE.6) THEN
1887
1888 C--colour partner
1889
1890         COLANT = 2
1891
1892         KHEP = JDAHEP(2,IHEP-1)
1893
1894       ELSEIF(DECAY.GT.1.AND.ID2.GE.7) THEN
1895
1896 C--anticolour partner
1897
1898         COLANT = 3
1899
1900         KHEP = JMOHEP(2,IHEP)
1901
1902       ELSE
1903
1904         KHEP=IHEP
1905
1906       ENDIF
1907
1908       IDM   = IDHW(JMOHEP(1,KHEP))
1909
1910       IF(ABS(IDPDG(IDM)).GT.1000000.OR.IDM.EQ.15) THEN
1911
1912         IDM2  = IDHW(JDAHEP(1,JMOHEP(1,KHEP)))
1913
1914         IDM3  = IDHW(JDAHEP(2,JMOHEP(1,KHEP)))
1915
1916         IDM4  = IDHW(JDAHEP(2,JMOHEP(1,KHEP))-1)
1917
1918         QHEP  = JMOHEP(1,KHEP)
1919
1920         IDMB  = IDHW(JMOHEP(1,QHEP))
1921
1922         IDMB2 = IDHW(JMOHEP(2,QHEP))
1923
1924         IDMB3 = IDHW(JDAHEP(1,QHEP))
1925
1926         IDMB4 = IDHW(JDAHEP(2,QHEP))
1927
1928       ENDIF
1929
1930 C--Now decide if the colour partner decayed via BV
1931
1932       IF(COLANT.EQ.2.AND.((((IDM.GE.413.AND.IDM.LE.418).OR.
1933
1934      &                     IDM.EQ.449.OR.IDM.EQ.405.OR.IDM.EQ.406).AND.
1935
1936      &                       (IDM2.GE.7.AND.IDM2.LE.12.AND.
1937
1938      &                       IDM3.GE.7.AND.IDM3.LE.12.AND.
1939
1940      &                       IDM4.GE.7.AND.IDM4.LE.12)).OR.
1941
1942      &             (IDM.EQ.15.AND.IDMB.LE.6.AND.IDMB2.LE.6.AND.
1943
1944      &              ((IDMB3.GE.7.AND.IDMB4.GE.12.AND.IDMB4.EQ.449).OR.
1945
1946      &               (IDMB3.GE.198.AND.IDMB3.LE.207.AND.
1947
1948      &                ABS(IDPDG(IDMB4)).GT.1000000))))) THEN
1949
1950         CONBV = .TRUE.
1951
1952         COLUPD = .TRUE.
1953
1954         HVFCEN = .FALSE.
1955
1956         XHEP = JMOHEP(2,JDAHEP(2,JMOHEP(1,KHEP)))
1957
1958       ELSEIF(COLANT.EQ.3.AND.((((IDM.GE.419.AND.IDM.LE.424).OR.
1959
1960      &                   IDM.EQ.449.OR.IDM.EQ.411.OR.IDM.EQ.412).AND.
1961
1962      &                    (IDM2.LE.6.AND.IDM3.LE.6.AND.IDM4.LE.6)).OR.
1963
1964      &               (IDM.EQ.15.AND.IDMB.GE.7.AND.IDMB.LE.12.AND.
1965
1966      &                IDMB2.GE.7.AND.IDMB2.LE.12.AND.((IDMB3.LE.6.AND.
1967
1968      &                IDMB4.EQ.449).OR.(ABS(IDPDG(IDMB4)).GT.1000000
1969
1970      &                .AND.IDMB3.GE.198.AND.IDMB3.LE.207))))) THEN
1971
1972         CONBV = .TRUE.
1973
1974         COLUPD = .TRUE.
1975
1976         HVFCEN = .FALSE.
1977
1978         XHEP = JDAHEP(2,JDAHEP(2,JMOHEP(1,KHEP)))
1979
1980       ELSE
1981
1982         CONBV = .FALSE.
1983
1984         COLUPD = .FALSE.
1985
1986         XHEP = 0
1987
1988       ENDIF
1989
1990       IF(CONBV) THEN
1991
1992         IF(IDM.NE.15) THEN
1993
1994           CLSAVE(1) = JDAHEP(2,JMOHEP(1,KHEP))-1
1995
1996           CLSAVE(2) = CLSAVE(1)+1
1997
1998         ELSE
1999
2000           IF(IDMB4.EQ.449) THEN
2001
2002             DO I=1,2
2003
2004               CLSAVE(I) = JMOHEP(I,JMOHEP(1,KHEP))
2005
2006               IF(CLSAVE(I).EQ.XHEP) CLSAVE(I)=JDAHEP(1,JMOHEP(1,KHEP))
2007
2008             ENDDO
2009
2010           ELSE
2011
2012             CLSAVE(1) = JMOHEP(1,JMOHEP(1,KHEP))
2013
2014             CLSAVE(2) = JMOHEP(2,JMOHEP(1,KHEP))
2015
2016           ENDIF
2017
2018         ENDIF
2019
2020       ELSE
2021
2022         CLSAVE(1)=0
2023
2024         CLSAVE(2)=0
2025
2026       ENDIF
2027
2028 C--Now set the colours for angular ordering
2029
2030       THEP = MHEP-1
2031
2032       IF(DECAY.EQ.1) THEN
2033
2034         IF(ID2.LE.6) THEN
2035
2036           JMOHEP(2,THEP) = THEP+HWRINT(1,2)
2037
2038           JDAHEP(2,THEP) = THEP
2039
2040         ELSE
2041
2042           JMOHEP(2,THEP) = THEP
2043
2044           JDAHEP(2,THEP) = THEP+HWRINT(1,2)
2045
2046         ENDIF
2047
2048       ELSEIF(DECAY.EQ.2) THEN
2049
2050         IF(ID2.LE.6) THEN
2051
2052           JMOHEP(2,THEP) = IHEP
2053
2054           JDAHEP(2,THEP) = THEP
2055
2056         ELSE
2057
2058           JMOHEP(2,THEP) = THEP
2059
2060           JDAHEP(2,THEP) = IHEP
2061
2062         ENDIF
2063
2064       ENDIF
2065
2066 C--Colour of the second two
2067
2068       DO JHEP=1,2
2069
2070         IF(ID2.LE.6) THEN
2071
2072           JMOHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+
2073
2074      &                           COLCON(HWRINT(1,2),JHEP,DECAY)
2075
2076           JDAHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+FLACON(JHEP,DECAY)
2077
2078         ELSE
2079
2080           JDAHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+
2081
2082      &                           COLCON(HWRINT(1,2),JHEP,DECAY)
2083
2084           JMOHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+FLACON(JHEP,DECAY)
2085
2086         ENDIF
2087
2088       ENDDO
2089
2090 C--Now set the colours of the colour partner
2091
2092       IF(DECAY.GT.1.AND..NOT.CONBV) THEN
2093
2094         IF(ID2.LE.6) JMOHEP(2,KHEP) = MHEP+HWRINT(0,1)
2095
2096         IF(ID2.GE.7) JDAHEP(2,KHEP) = MHEP+HWRINT(0,1)
2097
2098       ELSEIF(CONBV) THEN
2099
2100         IF(ID2.GT.6) THEN
2101
2102           JMOHEP(2,CLSAVE(1)) = MHEP+HWRINT(0,1)
2103
2104           IF(JMOHEP(2,CLSAVE(1)).EQ.MHEP) THEN
2105
2106             JMOHEP(2,CLSAVE(2)) = MHEP+1
2107
2108           ELSE
2109
2110             JMOHEP(2,CLSAVE(2)) = MHEP
2111
2112           ENDIF
2113
2114         ELSE
2115
2116           JDAHEP(2,CLSAVE(1)) = MHEP+HWRINT(0,1)
2117
2118           IF(JDAHEP(2,CLSAVE(1)).EQ.MHEP) THEN
2119
2120             JDAHEP(2,CLSAVE(2)) = MHEP+1
2121
2122           ELSE
2123
2124             JDAHEP(2,CLSAVE(2)) = MHEP
2125
2126           ENDIF
2127
2128         ENDIF
2129
2130       ENDIF
2131
2132  999  END
2133
2134 CDECK  ID>, HWDRME.
2135
2136 *CMZ :-        -20/07/99  10:56:12  by  Peter Richardson
2137
2138 *-- Author :    Peter Richardson
2139
2140 C-----------------------------------------------------------------------
2141
2142       SUBROUTINE HWDRME(LHEP,MHEP)
2143
2144 C-----------------------------------------------------------------------
2145
2146 C     SUBROUTINE TO IMPLEMENT ALL RPARITY DECAY MATRIX ELEMENTS
2147
2148 C-----------------------------------------------------------------------
2149
2150       INCLUDE 'HERWIG61.INC'
2151
2152       DOUBLE PRECISION SM(6),SW(6),HWULDO,INFCOL,AM, M12SQ,M23SQ,MSGN,
2153
2154      &                 M13SQ,A(6),B(6),SWEAK,MW,DECMOM(5),TEST(4),EPS,
2155
2156      &                 M12SQT(6),M23SQT(6),M13SQT(6),LIMIT,M(4),RAND,
2157
2158      &                 MC(2),MX2(6),MX(6),HWDPWT,HWR,HWDRM1,LAMD(3)
2159
2160       EXTERNAL         HWDRM1,HWULDO,HWDPWT,HWR
2161
2162       INTEGER K,SN(3),LHEP,CSP,I,SB(3),J,ND,LTRY,MHEP,NSP,ID(3),IG,
2163
2164      &        IDHWTP,IDHPTP,MTRY
2165
2166       PARAMETER(EPS=1D-20)
2167
2168       IF(IERROR.NE.0) RETURN
2169
2170 C--Electroweak parameters, etc
2171
2172       SWEAK = SQRT(SWEIN)
2173
2174       MW    = RMASS(198)
2175
2176       M(4)  = PHEP(5,LHEP)
2177
2178       IG    = IDHW(LHEP)
2179
2180 C--Find the masses of the final state and zero parameters
2181
2182       DO K=1,3
2183
2184         ID(K) = IDHW(MHEP+K-1)
2185
2186         IF(ID(K).LE.12) THEN
2187
2188           SN(K)=ID(K)
2189
2190         ELSE
2191
2192           SN(K)=ID(K)-120
2193
2194         ENDIF
2195
2196         IF(SN(K).GT.6) SN(K)=SN(K)-6
2197
2198         M(K) = PHEP(5,LHEP+K)
2199
2200         SB(K)=SN(K)
2201
2202         LAMD(K) = ZERO
2203
2204       ENDDO
2205
2206       DO J=1,6
2207
2208         MX2(J) = ZERO
2209
2210         MX(J)  = ZERO
2211
2212         M13SQT(J) = ZERO
2213
2214         M23SQT(J) = ZERO
2215
2216         M12SQT(J) = ZERO
2217
2218       ENDDO
2219
2220 C--Evaluate the coefficents for the mode we want
2221
2222       IF(IG.GE.450.AND.IG.LE.453) THEN
2223
2224 C--NEUTRALINO
2225
2226         NSP = IG-449
2227
2228         AM = RMASS(IG)
2229
2230         MSGN = ZSGNSS(NSP)
2231
2232         MC(1) =  ZMIXSS(NSP,3)/(2*MW*COSB*SWEAK)
2233
2234         MC(2) =  ZMIXSS(NSP,4)/(2*MW*SINB*SWEAK)
2235
2236 C--Calculate the combinations of couplings needed
2237
2238         IF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
2239
2240 C--first for the UDD modes
2241
2242           DO J=1,2
2243
2244             A(J) = M(1)*MC(2)*QMIXSS(SN(1),2,J)
2245
2246      &             +SLFCH(SN(1),NSP)*QMIXSS(SN(1),1,J)
2247
2248             B(J) = MSGN*(M(1)*MC(2)*QMIXSS(SN(1),1,J)
2249
2250      &             +SRFCH(SN(1),NSP)*QMIXSS(SN(1),2,J))
2251
2252             MX2(J) = QMIXSS(SN(1),2,J)
2253
2254             A(J+2) = M(2)*MC(1)*QMIXSS(SN(2),2,J)
2255
2256      &               +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
2257
2258             B(J+2) = MSGN*(M(2)*MC(1)*QMIXSS(SN(2),1,J)
2259
2260      &               +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
2261
2262             MX2(J+2) = QMIXSS(SN(2),2,J)
2263
2264             A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
2265
2266      &              +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
2267
2268             B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
2269
2270      &              +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
2271
2272             MX2(J+2) = QMIXSS(SN(3),2,J)
2273
2274           ENDDO
2275
2276           DO K=1,3
2277
2278             SN(K) = SN(K)+400
2279
2280             SB(K) = SB(K)+412
2281
2282           ENDDO
2283
2284         ELSEIF(ID(1).GE.121.AND.ID(2).GE.121.AND.ID(3).GE.121) THEN
2285
2286 C--Now for the LLE modes
2287
2288           DO J=1,2
2289
2290             A(J)  = MSGN*(M(1)*MC(1)*LMIXSS(SN(1),1,J)
2291
2292      &              +SRFCH(10+SN(1),NSP)*LMIXSS(SN(1),2,J))
2293
2294             B(J)  = M(1)*MC(1)*LMIXSS(SN(1),2,J)
2295
2296      &              +SLFCH(10+SN(1),NSP)*LMIXSS(SN(2),1,J)
2297
2298             MX2(J)= LMIXSS(SN(1),1,J)
2299
2300             A(J+2) = ZERO
2301
2302             B(J+2) = SLFCH(10+SN(2),NSP)*LMIXSS(SN(2),1,J)
2303
2304             MX2(J+2) =  LMIXSS(SN(2),1,J)
2305
2306             A(J+4) = M(3)*MC(1)*LMIXSS(SN(3),2,J)
2307
2308      &      +SLFCH(10+SN(3),NSP)*LMIXSS(SN(3),1,J)
2309
2310             B(J+4) = MSGN*(M(3)*MC(1)*LMIXSS(SN(3),1,J)
2311
2312      &      +SRFCH(10+SN(3),NSP)*LMIXSS(SN(3),2,J))
2313
2314             MX2(4+J) = LMIXSS(SN(3),2,J)
2315
2316           ENDDO
2317
2318           DO J=1,3
2319
2320             SN(J) = SN(J) + 424
2321
2322             SB(J) = SB(J) + 436
2323
2324           ENDDO
2325
2326         ELSE
2327
2328 C--Now for both types of LQD modes
2329
2330           IF(MOD(SN(1),2).EQ.0) THEN
2331
2332 C--First the neutrino,down,antidown mode
2333
2334             DO J=1,2
2335
2336               A(J) = ZERO
2337
2338               B(J) = SLFCH(10+SN(1),NSP)*
2339
2340      &               LMIXSS(SN(1),1,J)
2341
2342               MX2(J) = LMIXSS(SN(1),1,J)
2343
2344               A(J+2) = MSGN*(M(2)*MC(1)*QMIXSS(SN(2),1,J)
2345
2346      &        +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
2347
2348               B(J+2) = M(2)*MC(1)*QMIXSS(SN(2),2,J)
2349
2350      &        +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
2351
2352               MX2(2+J) = QMIXSS(SN(2),1,J)
2353
2354               A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
2355
2356      &        +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
2357
2358               B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
2359
2360      &        +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
2361
2362               MX2(J+4) = QMIXSS(SN(3),2,J)
2363
2364             ENDDO
2365
2366           ELSE
2367
2368 C--Now the charged lepton, antiup,down modes
2369
2370             DO J=1,2
2371
2372               A(J) = MSGN*(M(1)*MC(1)*LMIXSS(SN(1),1,J)
2373
2374      &        +SRFCH(10+SN(1),NSP)*LMIXSS(SN(1),2,J))
2375
2376               B(J) = M(1)*MC(1)*LMIXSS(SN(1),2,J)
2377
2378      &        +SLFCH(10+SN(1),NSP)*LMIXSS(SN(1),1,J)
2379
2380               MX2(J) = LMIXSS(SN(1),1,J)
2381
2382               A(J+2) =MSGN*(M(2)*MC(2)*QMIXSS(SN(2),1,J)
2383
2384      &        +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
2385
2386               B(J+2) = M(2)*MC(2)*QMIXSS(SN(2),2,J)
2387
2388      &        +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
2389
2390               MX2(2+J) = QMIXSS(SN(2),1,J)
2391
2392               A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
2393
2394      &        +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
2395
2396               B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
2397
2398      &        +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
2399
2400               MX2(J+4) = QMIXSS(SN(3),2,J)
2401
2402             ENDDO
2403
2404           ENDIF
2405
2406           SN(1) = SN(1) + 424
2407
2408           SB(1) = SB(1) + 436
2409
2410           DO J=2,3
2411
2412             SN(J) = SN(J) + 400
2413
2414             SB(J) = SB(J) + 412
2415
2416           ENDDO
2417
2418         ENDIF
2419
2420         DO K=1,3
2421
2422           SM(2*K-1) = RMASS(SN(K))
2423
2424           SM(2*K)   = RMASS(SB(K))
2425
2426           SW(2*K-1) = HBAR/RLTIM(SN(K))
2427
2428           SW(2*K)   = HBAR/RLTIM(SB(K))
2429
2430         ENDDO
2431
2432         ND = 3
2433
2434         DO K=1,3
2435
2436           LAMD(K) = ONE
2437
2438         ENDDO
2439
2440         INFCOL = ONE
2441
2442       ELSEIF(IG.EQ.449) THEN
2443
2444 C--GLUINO
2445
2446 C--First obtian the masses and widths needed
2447
2448         AM  = RMASS(IG)
2449
2450         ND = 3
2451
2452 C--Calculate the combinations of couplings needed
2453
2454         IF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
2455
2456 C--first for the UDD modes
2457
2458           INFCOL = -0.5D0
2459
2460 C--Couplings
2461
2462           DO I=1,3
2463
2464             DO J=1,2
2465
2466               A(2*I-2+J)  = -QMIXSS(SN(I),1,J)
2467
2468               B(2*I-2+J)  =  QMIXSS(SN(I),2,J)
2469
2470               MX2(2*I-2+J) =  QMIXSS(SN(I),2,J)
2471
2472             ENDDO
2473
2474             SN(I) = SN(I)+400
2475
2476             SB(I) = SB(I)+412
2477
2478           ENDDO
2479
2480         ELSE
2481
2482           INFCOL = ONE
2483
2484 C--Now for both types of LQD modes
2485
2486           IF(MOD(SN(1),2).EQ.0) THEN
2487
2488 C--First the neutrino,down,antidown mode
2489
2490             DO J=1,2
2491
2492               A(J)  = ZERO
2493
2494               B(J)  = ZERO
2495
2496               MX2(J) = ZERO
2497
2498               A(J+2)   =  QMIXSS(SN(2),2,J)
2499
2500               B(J+2)   = -QMIXSS(SN(2),1,J)
2501
2502               MX2(J+2) =  QMIXSS(SN(2),1,J)
2503
2504               A(J+4)   = -QMIXSS(SN(3),1,J)
2505
2506               B(J+4)   =  QMIXSS(SN(3),2,J)
2507
2508               MX2(4+J) =  QMIXSS(SN(3),2,J)
2509
2510             ENDDO
2511
2512           ELSEIF(MOD(SN(1),2).EQ.1) THEN
2513
2514 C--Now the charged lepton, antiup,down modes
2515
2516             DO J=1,2
2517
2518               A(J)  = ZERO
2519
2520               B(J)  = ZERO
2521
2522               MX2(J) = ZERO
2523
2524               A(J+2)   =  QMIXSS(SN(2),2,J)
2525
2526               B(J+2)   = -QMIXSS(SN(2),1,J)
2527
2528               MX2(J+2) =  QMIXSS(SN(2),1,J)
2529
2530               A(J+4)     = -QMIXSS(SN(3),1,J)
2531
2532               B(J+4)   =  QMIXSS(SN(3),2,J)
2533
2534               MX2(J+4) =  QMIXSS(SN(3),2,J)
2535
2536             ENDDO
2537
2538           ENDIF
2539
2540           SN(1) = SN(1) + 424
2541
2542           SB(1) = SB(1) + 436
2543
2544           DO K=2,3
2545
2546             SN(K) = SN(K) + 400
2547
2548             SB(K) = SB(K) + 412
2549
2550           ENDDO
2551
2552         ENDIF
2553
2554         DO K=1,3
2555
2556           SM(2*K-1) = RMASS(SN(K))
2557
2558           SM(2*K)   = RMASS(SB(K))
2559
2560           SW(2*K-1) = HBAR/RLTIM(SN(K))
2561
2562           SW(2*K)   = HBAR/RLTIM(SB(K))
2563
2564         ENDDO
2565
2566         DO K=1,3
2567
2568           LAMD(K) = ONE
2569
2570         ENDDO
2571
2572       ELSEIF(IG.GE.454.AND.IG.LE.457) THEN
2573
2574 C--CHARGINO
2575
2576         CSP = IG-453
2577
2578         IF(CSP.GT.2) CSP = CSP-2
2579
2580         AM  = RMASS(IG)
2581
2582         INFCOL = -ONE
2583
2584         MSGN = WSGNSS(CSP)
2585
2586         MC(1) =  ONE/(SQRT(2.0D0)*MW*COSB)
2587
2588         MC(2) =  ONE/(SQRT(2.0D0)*MW*SINB)
2589
2590 C--Calculate the combinations of the couplings needed
2591
2592         IF(ID(1).GT.120.AND.ID(2).GT.120.AND.ID(3).GT.120) THEN
2593
2594 C--first for the LLE modes, three modes
2595
2596           IF(MOD(SN(1),2).EQ.0.AND.MOD(SN(3),2).EQ.0) THEN
2597
2598 C--the one diagram mode nubar,positron, nu
2599
2600             DO J=1,2
2601
2602               A(J+4) = LMIXSS(SN(3)-1,1,J)*WMXUSS(CSP,1)
2603
2604      & -RMASS(SN(3)+119)*MC(1)*LMIXSS(SN(3)-1,2,J)*WMXUSS(CSP,2)
2605
2606               B(J+4) = ZERO
2607
2608               MX2(J+4) = LMIXSS(SN(3)-1,2,J)
2609
2610             ENDDO
2611
2612             ND = 1
2613
2614             SN(3) = SN(3)+423
2615
2616             SB(3) = SB(3)+435
2617
2618           ELSEIF(MOD(SN(1),2).EQ.0.AND.MOD(SN(2),2).EQ.0) THEN
2619
2620 C--the first two diagram mode nu, nu, positron
2621
2622             DO J=1,2
2623
2624               A(J)   = ZERO
2625
2626               B(J)   = LMIXSS(SN(1)-1,1,J)*WMXUSS(CSP,1)
2627
2628      & -RMASS(SN(1)+119)*MC(1)*LMIXSS(SN(1)-1,2,J)*WMXUSS(CSP,2)
2629
2630               A(J+2) = ZERO
2631
2632               B(J+2) = LMIXSS(SN(2)-1,1,J)*WMXUSS(CSP,1)
2633
2634      & -RMASS(SN(2)+119)*MC(1)*LMIXSS(SN(2)-1,2,J)*WMXUSS(CSP,2)
2635
2636               MX2(J)   = LMIXSS(SN(1)-1,1,J)
2637
2638               MX2(J+2) = LMIXSS(SN(2)-1,1,J)
2639
2640             ENDDO
2641
2642             ND = 2
2643
2644             DO J=1,2
2645
2646               SN(J) = SN(J)+423
2647
2648               SB(J) = SB(J)+435
2649
2650             ENDDO
2651
2652           ELSE
2653
2654 C--the second two diagram mode positron, positron, electron
2655
2656             DO J=1,2
2657
2658               A(J)   = -M(1)*WMXUSS(CSP,2)*MC(1)*LMIXSS(SN(1)+1,1,J)
2659
2660               B(J)   = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(1)+1,1,J)
2661
2662               A(J+2) = -M(2)*WMXUSS(CSP,2)*MC(1)*LMIXSS(SN(2)+1,1,J)
2663
2664               B(J+2) = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(2)+1,1,J)
2665
2666               MX2(J)   = LMIXSS(SN(1)+1,1,J)
2667
2668               MX2(J+2) = LMIXSS(SN(2)+1,1,J)
2669
2670             ENDDO
2671
2672             DO J=1,2
2673
2674               SN(J) = SN(J)+425
2675
2676               SB(J) = SB(J)+437
2677
2678             ENDDO
2679
2680             ND = 2
2681
2682           ENDIF
2683
2684           DO K=1,3
2685
2686             LAMD(K) = ONE
2687
2688           ENDDO
2689
2690         ELSEIF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
2691
2692 C--now for the UDD
2693
2694           IF(MOD(SN(1),2).EQ.0) THEN
2695
2696 C--two diagram mode
2697
2698             LAMD(1) = LAMDA3(SN(2)/2,SN(1)/2,(SN(3)+1)/2)
2699
2700             LAMD(2) = LAMDA3(SN(1)/2,SN(2)/2,(SN(3)+1)/2)
2701
2702             DO J=1,2
2703
2704               A(J)   = WMXUSS(CSP,1)*QMIXSS(SN(1)-1,1,J)
2705
2706      & -RMASS(SN(1)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(1)-1,2,J)
2707
2708               B(J)   = -MSGN*M(2)*WMXVSS(CSP,2)*QMIXSS(SN(1)-1,1,J)
2709
2710               A(J+2) = WMXUSS(CSP,1)*QMIXSS(SN(2)-1,1,J)
2711
2712      & -RMASS(SN(2)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(2)-1,2,J)
2713
2714               B(J+2) = -MSGN*M(2)*WMXVSS(CSP,2)*QMIXSS(SN(2)-1,1,J)
2715
2716               MX2(J)   = QMIXSS(SN(1)-1,2,J)
2717
2718               MX2(J+2) = QMIXSS(SN(2)-1,2,J)
2719
2720             ENDDO
2721
2722             DO J=1,2
2723
2724               SN(J) = SN(J) + 399
2725
2726               SB(J) = SB(J) + 411
2727
2728             ENDDO
2729
2730             ND = 2
2731
2732           ELSE
2733
2734 C--three diagram mode
2735
2736             LAMD(1) = LAMDA3((SN(1)+1)/2,(SN(2)+1)/2,(SN(3)+1)/2)
2737
2738             LAMD(2) = LAMDA3((SN(2)+1)/2,(SN(1)+1)/2,(SN(3)+1)/2)
2739
2740             LAMD(3) = LAMDA3((SN(3)+1)/2,(SN(2)+1)/2,(SN(1)+1)/2)
2741
2742             DO I=1,3
2743
2744               DO J=1,2
2745
2746                 A(J+2*I-2) = MSGN*(WMXVSS(CSP,1)*QMIXSS(SN(I)+1,1,J)
2747
2748      & -RMASS(SN(I)+1)*MC(2)*WMXVSS(CSP,2)*QMIXSS(SN(I)+1,2,J))
2749
2750                 B(J+2*I-2) = -M(I)*MC(1)*WMXUSS(CSP,2)
2751
2752      &                       *QMIXSS(SN(I)+1,1,J)
2753
2754                 MX2(J+2*I-2)   = QMIXSS(SN(I)+1,2,J)
2755
2756               ENDDO
2757
2758               SN(I) = SN(I) + 401
2759
2760               SB(I) = SB(I) + 413
2761
2762             ENDDO
2763
2764             ND = 3
2765
2766           ENDIF
2767
2768         ELSE
2769
2770 C--now for the LQD modes
2771
2772           IF(MOD(SN(2),2).EQ.1.AND.MOD(SN(3),2).EQ.0) THEN
2773
2774 C--first one diagram mode nubar, dbar, up
2775
2776             DO J=1,2
2777
2778               A(J+4) = -MSGN*M(3)*WMXVSS(CSP,2)*MC(2)*
2779
2780      &                  QMIXSS(SN(3)-1,1,J)
2781
2782               B(J+4) = WMXUSS(CSP,1)*QMIXSS(SN(3)-1,1,J)
2783
2784      &        -RMASS(SN(3)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(3)-1,2,1)
2785
2786               MX2(J+4)   = QMIXSS(SN(3)-1,2,J)
2787
2788             ENDDO
2789
2790             SN(3) = SN(3) + 399
2791
2792             SB(3) = SB(3) + 411
2793
2794             ND = 1
2795
2796           ELSEIF(MOD(SN(2),2).EQ.0.AND.MOD(SN(3),2).EQ.0) THEN
2797
2798 C--second one diagram mode positron, ubar, up
2799
2800             DO J=1,2
2801
2802               A(J+4) = -MSGN*M(3)*WMXVSS(CSP,2)*MC(2)*
2803
2804      &                  QMIXSS(SN(3)-1,1,J)
2805
2806               B(J+4) = WMXUSS(CSP,1)*QMIXSS(SN(3)-1,1,J)
2807
2808      &   -RMASS(SN(3)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(3)-1,2,1)
2809
2810               MX2(J+4)   = QMIXSS(SN(3)-1,2,J)
2811
2812             ENDDO
2813
2814             SN(3) = SN(3) + 399
2815
2816             SB(3) = SB(3) + 411
2817
2818             ND = 1
2819
2820           ELSEIF(MOD(SN(2),2).EQ.1.AND.MOD(SN(3),2).EQ.1) THEN
2821
2822 C--first two diagram mode positron, dbar, down
2823
2824             DO J=1,2
2825
2826               A(J)   = -M(1)*MC(1)*WMXUSS(CSP,2)*LMIXSS(SN(1)+1,1,J)
2827
2828               B(J)   = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(2)+1,1,J)
2829
2830               A(J+2) = -M(2)*WMXUSS(CSP,2)*MC(1)*QMIXSS(SN(2)+1,1,J)
2831
2832               B(J+2) = MSGN*(WMXVSS(CSP,1)*QMIXSS(SN(2)+1,1,J)
2833
2834      &   -RMASS(SN(2)+1)*MC(2)*WMXVSS(CSP,2)*QMIXSS(SN(2)+1,2,J))
2835
2836               MX2(J)   = LMIXSS(SN(1)+1,1,J)
2837
2838               MX2(J+2) = QMIXSS(SN(2)+1,1,J)
2839
2840             ENDDO
2841
2842             SN(1) = SN(1) + 425
2843
2844             SB(1) = SB(1) + 437
2845
2846             SN(2) = SN(2) + 401
2847
2848             SB(2) = SB(2) + 413
2849
2850             ND = 2
2851
2852           ELSE
2853
2854 C--second two diagram mode nu, up, dbar
2855
2856             DO J=1,2
2857
2858               A(J)   = ZERO
2859
2860               B(J)   = WMXUSS(CSP,1)*LMIXSS(SN(1)-1,1,J)
2861
2862      &   -RMASS(119+SN(1))*MC(1)*WMXUSS(CSP,2)*LMIXSS(SN(1)-1,2,J)
2863
2864               A(J+2) = -MSGN*M(2)*MC(2)*WMXVSS(CSP,2)*
2865
2866      &                 QMIXSS(SN(2)-1,1,J)
2867
2868               B(J+2) = WMXUSS(CSP,1)*QMIXSS(SN(2)-1,1,J)
2869
2870      &   -RMASS(SN(2)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(2)-1,2,J)
2871
2872               MX2(J)   = LMIXSS(SN(1)-1,1,J)
2873
2874               MX2(J+2) = QMIXSS(SN(2)-1,1,J)
2875
2876             ENDDO
2877
2878             SN(1) = SN(1) + 423
2879
2880             SB(1) = SB(1) + 435
2881
2882             SN(2) = SN(2) + 399
2883
2884             SB(2) = SB(2) + 411
2885
2886             ND = 2
2887
2888           ENDIF
2889
2890           DO K=1,3
2891
2892             LAMD(K) = ONE
2893
2894           ENDDO
2895
2896         ENDIF
2897
2898         IF(ND.EQ.1) THEN
2899
2900           DO K=1,2
2901
2902             SM(2*K-1) = 0.0D0
2903
2904             SM(2*K)   = 0.0D0
2905
2906             SW(2*K-1) = 0.0D0
2907
2908             SW(2*K)   = 0.0D0
2909
2910           ENDDO
2911
2912           SM(5) = RMASS(SN(3))
2913
2914           SM(6)   = RMASS(SB(3))
2915
2916           SW(5) = HBAR/RLTIM(SN(3))
2917
2918           SW(6)   = HBAR/RLTIM(SB(3))
2919
2920         ELSE
2921
2922           DO K=1,2
2923
2924             SM(2*K-1) = RMASS(SN(K))
2925
2926             SM(2*K)   = RMASS(SB(K))
2927
2928             SW(2*K-1) = HBAR/RLTIM(SN(K))
2929
2930             SW(2*K)   = HBAR/RLTIM(SB(K))
2931
2932             SM(4+K)   = ZERO
2933
2934             SW(4+K)   = ZERO
2935
2936           ENDDO
2937
2938         ENDIF
2939
2940       ELSE
2941
2942 C--UNKNOWN
2943
2944         CALL HWWARN('HWDRME',500,*999)
2945
2946       ENDIF
2947
2948 C--Set mixing to zero if diagram not available
2949
2950       IF((AM.LT.(ABS(SM(1))+M(1)).OR.ABS(SM(1)).LT.(M(2)+M(3)))
2951
2952      &   .AND.ABS(MX2(1)).GT.ZERO.AND.ND.NE.1) MX(1) = MX2(1)*LAMD(1)
2953
2954         IF((AM.LT.(ABS(SM(2))+M(1)).OR.ABS(SM(2)).LT.(M(2)+M(3)))
2955
2956      &   .AND.ABS(MX2(2)).GT.ZERO.AND.ND.NE.1) MX(2) = MX2(2)*LAMD(1)
2957
2958         IF((AM.LT.(ABS(SM(3))+M(2)).OR.ABS(SM(3)).LT.(M(1)+M(3)))
2959
2960      &   .AND.ABS(MX2(3)).GT.ZERO.AND.ND.NE.1) MX(3) = MX2(3)*LAMD(2)
2961
2962         IF((AM.LT.(ABS(SM(4))+M(2)).OR.ABS(SM(4)).LT.(M(1)+M(3)))
2963
2964      &   .AND.ABS(MX2(4)).GT.ZERO.AND.ND.NE.1) MX(4) = MX2(4)*LAMD(2)
2965
2966         IF((AM.LT.(ABS(SM(5))+M(3)).OR.ABS(SM(5)).LT.(M(1)+M(2)))
2967
2968      &   .AND.ABS(MX2(5)).GT.ZERO.AND.ND.NE.2) MX(5) = MX2(5)*LAMD(3)
2969
2970         IF((AM.LT.(ABS(SM(6))+M(3)).OR.ABS(SM(6)).LT.(M(1)+M(2)))
2971
2972      &   .AND.ABS(MX2(6)).GT.ZERO.AND.ND.NE.2) MX(6) = MX2(6)*LAMD(3)
2973
2974 C--Calculate the limiting points
2975
2976       DO J=1,2
2977
2978         IF(ND.NE.1) THEN
2979
2980           IF(ABS(MX(J)).GT.EPS) CALL HWDRM5(M23SQT(J),M13SQT(J),
2981
2982      &      M12SQT(J),A(J),B(J),M(2),M(3),M(1),M(4),SM(J),SW(J))
2983
2984           IF(ABS(MX(J+2)).GT.EPS) CALL HWDRM5(M13SQT(2+J),M23SQT(2+J),
2985
2986      &    M12SQT(2+J),A(2+J),B(2+J),M(1),M(3),M(2),M(4),SM(2+J),SW(2+J))
2987
2988         ENDIF
2989
2990         IF(ND.NE.2) THEN
2991
2992           IF(ABS(MX(J+4)).GT.EPS) CALL HWDRM5(M12SQT(4+J),M23SQT(4+J),
2993
2994      &    M13SQT(4+J),A(4+J),B(4+J),M(1),M(2),M(3),M(4),SM(4+J),SW(4+J))
2995
2996         ENDIF
2997
2998       ENDDO
2999
3000 C--Now evaluate the limit using these points
3001
3002       LIMIT = ZERO
3003
3004       DO 100 I=1,6
3005
3006         IF(ABS(MX(I)).LT.EPS) GOTO 100
3007
3008         LIMIT = LIMIT+HWDRM1(TEST,M12SQT(I),M13SQT(I),M23SQT(I),A,B,MX,
3009
3010      &                       M,SM,SW,INFCOL,AM,0,ND)
3011
3012  100  CONTINUE
3013
3014       LIMIT = TWO*LIMIT
3015
3016 C--Now evaluate at a random point
3017
3018       MTRY = 0
3019
3020  25   MTRY = MTRY+1
3021
3022       LTRY = 0
3023
3024  35   LTRY = LTRY+1
3025
3026       CALL HWDTHR(PHEP(1,LHEP),PHEP(1,MHEP),
3027
3028      &                  PHEP(1,MHEP+1),PHEP(1,MHEP+2),HWDPWT)
3029
3030 C--Now calculate the m12sq etc for the actual point
3031
3032       M12SQ = M(1)**2+M(2)**2+2*HWULDO(PHEP(1,MHEP),PHEP(1,MHEP+1))
3033
3034       M13SQ = M(1)**2+M(3)**2+2*HWULDO(PHEP(1,MHEP),PHEP(1,MHEP+2))
3035
3036       M23SQ = M(2)**2+M(3)**2+2*HWULDO(PHEP(1,MHEP+1),PHEP(1,MHEP+2))
3037
3038 C--Now calulate the matrix element
3039
3040       TEST(4) = HWDRM1(TEST,M12SQ,M13SQ,M23SQ,A,B,MX,
3041
3042      &                       M,SM,SW,INFCOL,AM,1,ND)
3043
3044 C--Now test the value againest the limit
3045
3046       RAND = HWR()*LIMIT
3047
3048       IF(TEST(4).GT.LIMIT) THEN
3049
3050         LIMIT = 1.1D0*TEST(4)
3051
3052         CALL HWWARN('HWDRME',51,*150)
3053
3054       ENDIF
3055
3056  150  IF(TEST(4).LT.RAND.AND.LTRY.LT.NETRY) THEN
3057
3058         GOTO 35
3059
3060       ELSEIF(LTRY.GE.NETRY) THEN
3061
3062         IF(MTRY.LE.NETRY) THEN
3063
3064           LIMIT = LIMIT*0.9D0
3065
3066           CALL HWWARN('HWDRME',52,*25)
3067
3068         ELSE
3069
3070           CALL HWWARN('HWDRME',100,*999)
3071
3072         ENDIF
3073
3074       ENDIF
3075
3076 C--Reorder the particles in gluino decay to get angular ordering right
3077
3078       IF(IG.EQ.449.AND.ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
3079
3080         DO LTRY=1,3
3081
3082           IF(TEST(LTRY).GT.RAND) THEN
3083
3084             IF(LTRY.EQ.2) THEN
3085
3086               IDHWTP        = IDHW(MHEP)
3087
3088               IDHW(MHEP)    = IDHW(MHEP+1)
3089
3090               IDHW(MHEP+1)  = IDHWTP
3091
3092               IDHPTP        = IDHEP(MHEP)
3093
3094               IDHEP(MHEP)   = IDHEP(MHEP+1)
3095
3096               IDHEP(MHEP+1) = IDHPTP
3097
3098               CALL HWVEQU(5,PHEP(1,MHEP),DECMOM)
3099
3100               CALL HWVEQU(5,PHEP(1,MHEP+1),PHEP(1,MHEP))
3101
3102               CALL HWVEQU(5,DECMOM,PHEP(1,MHEP+1))
3103
3104             ELSEIF(LTRY.EQ.3) THEN
3105
3106               IDHWTP        = IDHW(MHEP)
3107
3108               IDHW(MHEP)    = IDHW(MHEP+2)
3109
3110               IDHW(MHEP+2)    = IDHWTP
3111
3112               IDHPTP        = IDHEP(MHEP)
3113
3114               IDHEP(MHEP)   = IDHEP(MHEP+2)
3115
3116               IDHEP(MHEP+2)   = IDHPTP
3117
3118               DO I=1,5
3119
3120               CALL HWVEQU(5,PHEP(1,MHEP),DECMOM)
3121
3122               CALL HWVEQU(5,PHEP(1,MHEP+2),PHEP(1,MHEP))
3123
3124               CALL HWVEQU(5,DECMOM,PHEP(1,MHEP+2))
3125
3126               ENDDO
3127
3128             ENDIF
3129
3130             GOTO 52
3131
3132           ENDIF
3133
3134           RAND=RAND-TEST(LTRY)
3135
3136         ENDDO
3137
3138       ENDIF
3139
3140  52   CONTINUE
3141
3142  999  END
3143
3144 CDECK  ID>, HWDRM1.
3145
3146 *CMZ :-        -20/07/99  10:56:12  by  Peter Richardson
3147
3148 *-- Author :    Peter Richardson
3149
3150 C-----------------------------------------------------------------------
3151
3152       FUNCTION HWDRM1(TEST,M12SQ,M13SQ,M23SQ,A,B,MX,M,SM,SW
3153
3154      &                ,INFCOL,AM,LM,ND)
3155
3156 C-----------------------------------------------------------------------
3157
3158 C     FUNCTION TO GIVE THE R-PARITY VIOLATING MATRIX ELEMENT AT A GIVEN
3159
3160 C     PHASE-SPACE POINT
3161
3162 C-----------------------------------------------------------------------
3163
3164       IMPLICIT NONE
3165
3166       DOUBLE PRECISION M12SQ,M13SQ,M23SQ,MX(6),A(6),B(6),SM(6),SW(6),
3167
3168      &                 INFCOL,AM,TERM(21),TEST(4),PLN,NPLN,ZERO,
3169
3170      &                 M(4),HWDRM1,HWDRM2,HWDRM3,HWDRM4
3171
3172       PARAMETER (ZERO=0)
3173
3174       EXTERNAL HWDRM2,HWDRM3,HWDRM4
3175
3176       INTEGER LM,K,ND
3177
3178 C--Zero the array
3179
3180         DO K=1,21
3181
3182           TERM(K) = 0.0D0
3183
3184         ENDDO
3185
3186         HWDRM1 = 0.0D0
3187
3188 C--The amplitude
3189
3190       IF(ABS(MX(1)).GT.ZERO.AND.ND.NE.1) THEN
3191
3192         TERM(1) = MX(1)**2*HWDRM2(M23SQ,M(2),M(3),M(1),M(4),SM(1),
3193
3194      &            SW(1),A(1),B(1))
3195
3196         IF(ABS(MX(2)).GT.ZERO) TERM(7)= MX(1)*MX(2)*HWDRM3(M23SQ,M(2),
3197
3198      &   M(3),M(1),M(4),SM(1),SM(2),SW(1),SW(2),A(1),A(2),B(1),B(2))
3199
3200         IF(ABS(MX(3)).GT.ZERO) TERM(10)=-MX(1)*MX(3)*HWDRM4(M13SQ,M23SQ,
3201
3202      &  M(1),M(3),M(2),M(4),SM(3),SM(1),SW(3),SW(1),A(1),A(3),B(1),B(3))
3203
3204         IF(ABS(MX(4)).GT.ZERO) TERM(11)=-MX(1)*MX(4)*HWDRM4(M13SQ,M23SQ,
3205
3206      &  M(1),M(3),M(2),M(4),SM(4),SM(1),SW(4),SW(1),A(1),A(4),B(1),B(4))
3207
3208         IF(ABS(MX(5)).GT.ZERO) TERM(12)=-MX(1)*MX(5)*HWDRM4(M23SQ,M12SQ,
3209
3210      &  M(3),M(2),M(1),M(4),SM(1),SM(5),SW(1),SW(5),A(5),A(1),B(5),B(1))
3211
3212         IF(ABS(MX(6)).GT.ZERO) TERM(13)=-MX(1)*MX(6)*HWDRM4(M23SQ,M12SQ,
3213
3214      &  M(3),M(2),M(1),M(4),SM(1),SM(6),SW(1),SW(6),A(6),A(1),B(6),B(1))
3215
3216       ENDIF
3217
3218       IF(ABS(MX(2)).GT.ZERO.AND.ND.NE.1) THEN
3219
3220         TERM(2) = MX(2)**2*HWDRM2(M23SQ,M(2),M(3),M(1),M(4),SM(2),
3221
3222      &            SW(2),A(2),B(2))
3223
3224         IF(ABS(MX(3)).GT.ZERO) TERM(14)=-MX(2)*MX(3)*HWDRM4(M13SQ,M23SQ,
3225
3226      &  M(1),M(3),M(2),M(4),SM(3),SM(2),SW(3),SW(2),A(2),A(3),B(2),B(3))
3227
3228         IF(ABS(MX(4)).GT.ZERO) TERM(15)=-MX(2)*MX(4)*HWDRM4(M13SQ,M23SQ,
3229
3230      &  M(1),M(3),M(2),M(4),SM(4),SM(2),SW(4),SW(2),A(2),A(4),B(2),B(4))
3231
3232         IF(ABS(MX(5)).GT.ZERO) TERM(16)=-MX(2)*MX(5)*HWDRM4(M23SQ,M12SQ,
3233
3234      &  M(3),M(2),M(1),M(4),SM(2),SM(5),SW(2),SW(5),A(5),A(2),B(5),B(2))
3235
3236         IF(ABS(MX(6)).GT.ZERO) TERM(17)=-MX(2)*MX(6)*HWDRM4(M23SQ,M12SQ,
3237
3238      &  M(3),M(2),M(1),M(4),SM(2),SM(6),SW(2),SW(6),A(6),A(2),B(6),B(2))
3239
3240       ENDIF
3241
3242       IF(ABS(MX(3)).GT.ZERO.AND.ND.NE.1) THEN
3243
3244         TERM(3) = MX(3)**2*HWDRM2(M13SQ,M(1),M(3),M(2),M(4),SM(3),
3245
3246      &            SW(3),A(3),B(3))
3247
3248         IF(ABS(MX(4)).GT.ZERO) TERM(8)= MX(3)*MX(4)*HWDRM3(M13SQ,M(1),
3249
3250      &   M(3),M(2),M(4),SM(3),SM(4),SW(3),SW(4),A(3),A(4),B(3),B(4))
3251
3252         IF(ABS(MX(5)).GT.ZERO) TERM(18)=-MX(3)*MX(5)*HWDRM4(M12SQ,M13SQ,
3253
3254      &  M(2),M(1),M(3),M(4),SM(5),SM(3),SW(5),SW(3),A(3),A(5),B(3),B(5))
3255
3256         IF(ABS(MX(6)).GT.ZERO) TERM(19)=-MX(3)*MX(6)*HWDRM4(M12SQ,M13SQ,
3257
3258      &  M(2),M(1),M(3),M(4),SM(6),SM(3),SW(6),SW(3),A(3),A(6),B(3),B(6))
3259
3260       ENDIF
3261
3262       IF(ABS(MX(4)).GT.ZERO.AND.ND.NE.1) THEN
3263
3264         TERM(4) = MX(4)**2*HWDRM2(M13SQ,M(1),M(3),M(2),M(4),SM(4),
3265
3266      &            SW(4),A(4),B(4))
3267
3268         IF(ABS(MX(5)).GT.ZERO) TERM(20)=-MX(4)*MX(5)*HWDRM4(M12SQ,M13SQ,
3269
3270      &  M(2),M(1),M(3),M(4),SM(5),SM(4),SW(5),SW(4),A(4),A(5),B(4),B(5))
3271
3272         IF(ABS(MX(6)).GT.ZERO) TERM(21)=-MX(4)*MX(6)*HWDRM4(M12SQ,M13SQ,
3273
3274      &  M(2),M(1),M(3),M(4),SM(6),SM(4),SW(6),SW(4),A(4),A(6),B(4),B(6))
3275
3276       ENDIF
3277
3278       IF(ABS(MX(5)).GT.ZERO.AND.ND.NE.2) THEN
3279
3280         TERM(5) = MX(5)**2*HWDRM2(M12SQ,M(1),M(2),M(3),M(4),SM(5),
3281
3282      &            SW(5),A(5),B(5))
3283
3284         IF(ABS(MX(6)).GT.ZERO) TERM(9)= MX(5)*MX(6)*HWDRM3(M12SQ,M(1),
3285
3286      &     M(2),M(3),M(4),SM(5),SM(6),SW(5),SW(6),A(5),A(6),B(5),B(6))
3287
3288       ENDIF
3289
3290       IF(ABS(MX(6)).GT.ZERO.AND.ND.NE.2) TERM(6) = MX(6)**2*
3291
3292      &    HWDRM2(M12SQ,M(1),M(2),M(3),M(4),SM(6),SW(6),A(6),B(6))
3293
3294       DO K=10,21
3295
3296         TERM(K)=TERM(K)*INFCOL
3297
3298       ENDDO
3299
3300 C--Add them up
3301
3302       DO K=1,21
3303
3304         HWDRM1 = HWDRM1+TERM(K)
3305
3306       ENDDO
3307
3308 C--Different colour flows for the gluino
3309
3310       IF(LM.NE.0) THEN
3311
3312         NPLN = 0.0D0
3313
3314         PLN = 0.0D0
3315
3316         DO K=1,9
3317
3318           PLN = PLN+TERM(K)
3319
3320         ENDDO
3321
3322         DO K=10,21
3323
3324           NPLN= NPLN+TERM(K)
3325
3326         ENDDO
3327
3328         DO K=1,3
3329
3330           TEST(K) = (TERM(2*K-1)+TERM(2*K)+TERM(6+K))*(1+NPLN/PLN)
3331
3332         ENDDO
3333
3334       ELSE
3335
3336         DO K=1,3
3337
3338           TEST(K) = 0.0D0
3339
3340         ENDDO
3341
3342       ENDIF
3343
3344       IF(TEST(4).LT.ZERO) CALL HWWARN('HWDRM1',50,*999)
3345
3346  999  END
3347
3348 CDECK  ID>, HWDRM2.
3349
3350 *CMZ :-        -20/07/99  10:56:12  by  Peter Richardson
3351
3352 *-- Author :    Peter Richardson
3353
3354 C-----------------------------------------------------------------------
3355
3356       FUNCTION HWDRM2(X,MA,MB,MC,MD,MR1,GAM1,A,B)
3357
3358 C-----------------------------------------------------------------------
3359
3360 C     Function to compute the matrix element squared part of a 3-body
3361
3362 C     R-parity decay
3363
3364 C-----------------------------------------------------------------------
3365
3366       IMPLICIT NONE
3367
3368       DOUBLE PRECISION X,MA,MB,MC,MD,A,B,HWDRM2,MR1,GAM1
3369
3370       HWDRM2  = (X - MA**2 - MB**2)*(4*A*B*MC*MD +
3371
3372      &    (A**2 + B**2)*(-X + MC**2 + MD**2))/
3373
3374      &     ((X-MR1**2)**2+GAM1**2*MR1**2)
3375
3376       END