4 *CMZ :- -24/04/92 14.23.44 by Mike Seymour
6 *-- Author : Mike Seymour
8 C-----------------------------------------------------------------------
10 SUBROUTINE HWDHIG(GAMINP)
12 C-----------------------------------------------------------------------
16 C A) FOR GAMinp=0 FIND AND DECAY HIGGS
18 C B) FOR GAMinp>0 CALCULATE TOTAL HIGGS WIDTH
20 C FOR EMH=GAMINP. STORE RESULT IN GAMINP.
22 C-----------------------------------------------------------------------
24 INCLUDE 'HERWIG61.INC'
26 DOUBLE PRECISION HWDHGF,HWR,HWRUNI,HWUSQR,HWUPCM,GAMINP,EMH,
28 & EMF,COLFAC,ENF,K1,K0,BET0,BET1,GAM0,GAM1,SCLOG,CFAC,XF,EM,GAMLIM,
30 & GAM,XW,EMW,XZ,EMZ,YW,YZ,EMI,TAUT,TAUW,WIDHIG,VECDEC,EMB,GAMB,
32 & TMIN,TMAX1,EM1,TMAX2,EM2,X1,X2,PROB,PCM,SUMR,SUMI,TAUTR,TAUTI,
36 INTEGER HWRINT,IHIG,I,IFERM,NLOOK,I1,I2,IPART,IMODE,IDEC,MMAX
40 EXTERNAL HWDHGF,HWR,HWRUNI,HWUSQR,HWUPCM,HWRINT,HWRLOG
46 DIMENSION VECDEC(2,0:NLOOK)
48 EQUIVALENCE (EMW,RMASS(198)),(EMZ,RMASS(200))
50 DATA GAMLIM,GAM,EM/10D0,2*0D0/
52 C---IF DECAY, FIND HIGGS (HWDHAD WILL HAVE GIVEN IT STATUS=1)
54 IF (GAMINP.EQ.ZERO) THEN
60 10 IF (IHIG.EQ.0.AND.IDHW(I).EQ.201.AND.ISTHEP(I).EQ.1) IHIG=I
62 IF (IHIG.EQ.0) CALL HWWARN('HWDHIG',101,*999)
66 IF (EMH.LE.ZERO) CALL HWWARN('HWDHIG',102,*999)
84 C---CALCULATE BRANCHING FRACTIONS
88 C---NLL CORRECTION TO QUARK DECAY RATE (HHG eq 2.6-9)
94 1 IF (2*RMASS(I).LT.EMH) ENF=ENF+1
100 BET0=(11*CAFAC-2*ENF)/3
102 BET1=(34*CAFAC**2-(10*CAFAC+6*CFFAC)*ENF)/3
106 GAM1=-404./3+40*ENF/9
108 SCLOG=LOG(EMH**2/QCDLAM**2)
110 CFAC=1 + ( K1/K0 - 2*GAM0 + GAM0*BET1/BET0**2*LOG(SCLOG)
112 & + (GAM0*BET1-GAM1*BET0)/BET0**2) / (BET0*SCLOG)
126 & EMF=EMF*(LOG(EMH/QCDLAM)/LOG(EMF/QCDLAM))**(GAM0/(2*BET0))
130 EMF=RMASS(107+IFERM*2)
140 IF (FOUR*XF.LT.ONE) THEN
142 GFACTR=ALPHEM/(8.*SWEIN*EMW**2)
144 BRHIG(IFERM)=COLFAC*GFACTR*EMH*EMF**2 * (1-4*XF)**1.5 * CFAC
156 IF (ABS(EM-EMH).GE.GAMLIM*GAM) THEN
158 C---OFF EDGE OF LOOK-UP TABLE
168 BRHIG(10)=.50*GFACTR * EMH**3 * HWDHGF(XW,YW)
170 BRHIG(11)=.25*GFACTR * EMH**3 * HWDHGF(XZ,YZ)
176 EMI=((EMH-EM)/(GAM*GAMLIM)+1)*NLOOK/2.0
182 BRHIG(10)=.50*GFACTR * EMH**3 * ( VECDEC(1,I1)*(I2-EMI) +
184 & VECDEC(1,I2)*(EMI-I1) )
186 BRHIG(11)=.25*GFACTR * EMH**3 * ( VECDEC(2,I1)*(I2-EMI) +
188 & VECDEC(2,I2)*(EMI-I1) )
194 TAUT=(2*RMASS(6)/EMH)**2
198 CALL HWDHGC(TAUT,TAUTR,TAUTI)
200 CALL HWDHGC(TAUW,TAUWR,TAUWI)
202 SUMR=4./3*( - 2*TAUT*( 1 + (1-TAUT)*TAUTR ) ) * ENHANC(6)
204 & +(2 + 3*TAUW*( 1 + (2-TAUW)*TAUWR ) ) * ENHANC(10)
206 SUMI=4./3*( - 2*TAUT*( (1-TAUT)*TAUTI ) ) * ENHANC(6)
208 & +( 3*TAUW*( (2-TAUW)*TAUWI ) ) * ENHANC(10)
210 BRHIG(12)=GFACTR*.03125*(ALPHEM/PIFAC)**2
212 & *EMH**3 * (SUMR**2 + SUMI**2)
218 IF (IPART.LT.12) BRHIG(IPART)=BRHIG(IPART)*ENHANC(IPART)**2
220 200 WIDHIG=WIDHIG+BRHIG(IPART)
222 IF (WIDHIG.EQ.ZERO) CALL HWWARN('HWDHIG',103,*999)
226 300 BRHIG(IPART)=BRHIG(IPART)/WIDHIG
228 IF (EM.NE.RMASS(201)) THEN
230 C---SET UP W*W*/Z*Z* LOOKUP TABLES
236 GAMLIM=MAX(GAMLIM,GAMMAX)
240 EMH=(I*2.0/NLOOK-1)*GAM*GAMLIM+EM
250 VECDEC(1,I)=HWDHGF(XW,YW)
252 VECDEC(2,I)=HWDHGF(XZ,YZ)
260 IF (GAMINP.GT.ZERO) THEN
268 C---SEE IF USER SPECIFIED A DECAY MODE
272 C---IF NOT, CHOOSE ONE
274 IF (IMODE.LT.1.OR.IMODE.GT.12) THEN
278 IF (IMODE.LT.1) MMAX=6
280 500 IMODE=HWRINT(1,MMAX)
282 IF (BRHIG(IMODE).LT.HWR()) GOTO 500
286 C---SEE IF SPECIFIED DECAY IS POSSIBLE
288 IF (BRHIG(IMODE).EQ.ZERO) CALL HWWARN('HWDHIG',104,*999)
294 ELSEIF (IMODE.LE.9) THEN
298 ELSEIF (IMODE.EQ.10) THEN
302 ELSEIF (IMODE.EQ.11) THEN
306 ELSEIF (IMODE.EQ.12) THEN
312 C---STATUS, IDs AND POINTERS
322 IDHEP(NHEP+I)=IDPDG(IDEC)
324 JDAHEP(I,IHIG)=NHEP+I
326 JMOHEP(1,NHEP+I)=IHIG
328 JMOHEP(2,NHEP+I)=NHEP+(3-I)
330 JDAHEP(2,NHEP+I)=NHEP+(3-I)
332 PHEP(5,NHEP+I)=RMASS(IDEC)
336 IF (IDEC.EQ.204) IDEC=199
338 IF (IDEC.EQ.206) IDEC=200
340 IF (IDEC.EQ. 65) IDEC= 59
344 C---ALLOW W/Z TO BE OFF-SHELL
346 IF (IMODE.EQ.10.OR.IMODE.EQ.11) THEN
348 IF (IMODE.EQ.10) THEN
362 C---STANDARD MASS DISTRIBUTION
364 700 TMIN=ATAN(-EMB/GAMB)
366 TMAX1=ATAN((EMH**2/EMB-EMB)/GAMB)
368 EM1=HWUSQR(EMB*(GAMB*TAN(HWRUNI(0,TMIN,TMAX1))+EMB))
370 TMAX2=ATAN(((EMH-EM1)**2/EMB-EMB)/GAMB)
372 EM2=HWUSQR(EMB*(GAMB*TAN(HWRUNI(0,TMIN,TMAX2))+EMB))
378 C---CORRECT MASS DISTRIBUTION
380 PROB=HWUSQR(1+X1**2+X2**2-2*X1-2*X2-2*X1*X2)
382 & * ((X1+X2-1)**2 + 8*X1*X2)
384 IF (.NOT.HWRLOG(PROB)) GOTO 700
386 C---CALCULATE SPIN DENSITY MATRIX
388 RHOHEP(1,NHEP+1)=4*X1*X2 / (8*X1*X2 + (X1+X2-1)**2)
390 RHOHEP(2,NHEP+1)=(X1+X2-1)**2 / (8*X1*X2 + (X1+X2-1)**2)
392 RHOHEP(3,NHEP+1)=RHOHEP(1,NHEP+1)
394 C---SYMMETRIZE DISTRIBUTIONS IN PARTICLES 1,2
396 IF (HWRLOG(HALF)) THEN
414 PCM=HWUPCM(EMH,PHEP(5,NHEP+1),PHEP(5,NHEP+2))
416 IF (PCM.LT.ZERO) CALL HWWARN('HWDHIG',105,*999)
418 CALL HWDTWO(PHEP(1,IHIG),PHEP(1,NHEP+1),PHEP(1,NHEP+2),
424 C---IF QUARK DECAY, HADRONIZE
446 *CMZ :- -20/10/99 09:46:43 by Peter Richardson
448 *-- Author : Ian Knowles & Bryan Webber
450 C-----------------------------------------------------------------------
454 C-----------------------------------------------------------------------
456 C Performs decays of heavy objects (heavy quarks & SUSY particles)
458 C MODIFIED TO INCLUDE R-PARITY VIOLATING SUSY PR 9/4/99
460 C-----------------------------------------------------------------------
462 INCLUDE 'HERWIG61.INC'
464 DOUBLE PRECISION HWUMBW,HWUPCM,HWR,SDKM,RN,BF,PCM,
466 & EMMX,EMWSQ,GMWSQ,EMLIM,PW(5),EMTST,HWDPWT,HWDWWT,HWULDO,PDW(5,3)
468 INTEGER IST(3),IHEP,IS,ID,IM,I,JHEP,KHEP,LHEP,MHEP,NPR,ISM,JCM,
470 & MTRY,NTRY,IDM,IDM2,THEP,CLSAVE(2),WHEP,RHEP
474 EXTERNAL HWR,HWDPWT,HWDWWT
476 DATA IST/113,114,114/
478 IF (IERROR.NE.0) RETURN
492 IF (.NOT.RSTAB(ID).AND.(ID.EQ.6.OR.ID.EQ.12.OR.
494 & (ID.GE.203.AND.ID.LE.218).OR.ABS(IDPDG(ID)).GT.1000000).AND.
496 & (IS.EQ.190.OR.(IS.GE.147.AND.IS.LE.151))) THEN
510 CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP))
512 CALL HWVEQU(4,VHEP(1,IHEP),VHEP(1,NHEP))
514 JMOHEP(1,NHEP)=JMOHEP(1,IHEP)
516 JMOHEP(2,NHEP)=JMOHEP(2,IHEP)
518 JDAHEP(1,NHEP)=JDAHEP(1,IHEP)
520 JDAHEP(2,NHEP)=JDAHEP(2,IHEP)
524 C Make a copy of decaying object
530 IDHW(NHEP)=IDHW(IHEP)
532 IDHEP(NHEP)=IDHEP(IHEP)
534 CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP))
536 CALL HWVEQU(4,VHEP(1,IHEP),VHEP(1,NHEP))
538 JMOHEP(1,NHEP)=JMOHEP(1,IHEP)
540 JMOHEP(2,NHEP)=JMOHEP(2,IHEP)
558 IF (BF.GE.RN) GOTO 30
562 CALL HWWARN('HWDHOB',50,*30)
564 30 IF (NHEP+5.GT.NMXHEP) CALL HWWARN('HWDHOB',100,*999)
568 JDAHEP(1,NHEP)=NHEP+1
570 JDAHEP(2,NHEP)=NHEP+NPR
572 C Reset colour pointers (if set)
578 IF (JDAHEP(2,JHEP).EQ.IHEP) JDAHEP(2,JHEP)=NHEP
580 IF(.NOT.RPARTY.AND.ISTHEP(JHEP).EQ.155
582 & .AND.ABS(IDHEP(JHEP)).GT.1000000
584 & .AND.JDAHEP(2,JHEP-1).EQ.IHEP) JDAHEP(2,JHEP-1) = NHEP
592 IF (JMOHEP(2,JHEP).EQ.IHEP) JMOHEP(2,JHEP)=NHEP
594 IF(.NOT.RPARTY.AND.ISTHEP(JHEP).EQ.155
596 & .AND.ABS(IDHEP(JHEP)).GT.1000000
598 & .AND.JMOHEP(2,JHEP-1).EQ.IHEP) JMOHEP(2,JHEP-1) = NHEP
602 C--Reset colour pointers if baryon number violated
608 IF(ISTHEP(JHEP).EQ.155
610 & .AND.ABS(IDHEP(JHEP)).GT.1000000.AND.
612 & JDAHEP(2,JHEP-1).EQ.IHEP) JDAHEP(2,JHEP-1)= NHEP
614 IF(JDAHEP(2,JHEP).EQ.IHEP) JDAHEP(2,JHEP)=NHEP
616 IF(JMOHEP(2,JHEP).EQ.IHEP) JMOHEP(2,JHEP)=NHEP
620 IF(HRDCOL(1,1).EQ.IHEP) HRDCOL(1,1)=NHEP
624 C Relabel original track
628 JMOHEP(2,IHEP)=JMOHEP(1,IHEP)
634 C Label decay products and choose masses
650 IDHW(NHEP)=IDKPRD(I,IM)
652 IDHEP(NHEP)=IDPDG(IDKPRD(I,IM))
660 PHEP(5,NHEP)=HWUMBW(IDKPRD(I,IM))
662 40 SDKM=SDKM-PHEP(5,NHEP)
664 IF (SDKM.LT.ZERO) THEN
668 IF (NTRY.LE.NETRY) GO TO 35
670 CALL HWWARN('HWDHOB',1,*45)
672 45 IF (MTRY.LE.NETRY) GO TO 15
674 CALL HWWARN('HWDHOB',101,*999)
678 C Assign production vertices to decay products
680 CALL HWUDKL(ID,PHEP(1,IHEP),VHEP(1,MHEP))
682 CALL HWVSUM(4,VHEP(1,IHEP),VHEP(1,MHEP),VHEP(1,MHEP))
684 CALL HWVEQU(4,VHEP(1,MHEP),VHEP(1,NHEP))
688 C Two body decay: LHEP -> MHEP + NHEP
690 PCM=HWUPCM(PHEP(5,IHEP),PHEP(5,MHEP),PHEP(5,NHEP))
692 CALL HWDTWO(PHEP(1,IHEP),PHEP(1,MHEP),
694 & PHEP(1,NHEP),PCM,TWO,.FALSE.)
696 ELSEIF (NPR.EQ.3) THEN
698 C Three body decay: LHEP -> KHEP + MHEP + NHEP
704 C Provisional colour self-connection of KHEP
710 IF (NME(IM).EQ.100) THEN
712 C Generate decay momenta using full (V-A)*(V-A) matrix element
714 EMMX=PHEP(5,IHEP)-PHEP(5,NHEP)
718 GMWSQ=(RMASS(198)*GAMW)**2
722 IF (EMMX.LT.RMASS(198)) EMLIM=EMLIM+(EMWSQ-EMMX**2)**2
724 50 CALL HWDTHR(PHEP(1,IHEP),PHEP(1,MHEP),
726 & PHEP(1,KHEP),PHEP(1,NHEP),HWDWWT)
728 CALL HWVSUM(4,PHEP(1,KHEP),PHEP(1,MHEP),PW)
732 EMTST=(EMWSQ-PW(5))**2
734 IF ((EMTST+GMWSQ)*HWR().GT.EMLIM) GOTO 50
738 C Assign production vertices to 1 and 2
740 CALL HWUDKL(198,PW,VHEP(1,KHEP))
742 CALL HWVSUM(4,VHEP(1,NHEP),VHEP(1,KHEP),VHEP(1,KHEP))
744 ELSEIF(NME(IM).EQ.300) THEN
746 C Generate momenta using 3-body RPV matrix element
748 CALL HWDRME(LHEP,KHEP)
752 C Three body phase space decay
754 CALL HWDTHR(PHEP(1,IHEP),PHEP(1,MHEP),
756 & PHEP(1,KHEP),PHEP(1,NHEP),HWDPWT)
760 CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,MHEP))
762 ELSEIF(NPR.EQ.4) THEN
764 C Four body decay: LHEP -> KHEP + RHEP + MHEP + NHEP
772 C Provisional colour connections of KHEP and RHEP
782 C Four body phase space decay
784 CALL HWDFOR(PHEP(1,IHEP),PHEP(1,KHEP),PHEP(1,RHEP),
786 & PHEP(1,MHEP),PHEP(1,NHEP))
788 CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,RHEP))
790 CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,MHEP))
794 CALL HWWARN('HWDHOB',102,*999)
800 IF (ID.EQ.6.OR.ID.EQ.12.OR.(ID.GE.209.AND.ID.LE.212)
802 & .OR.(ID.GE.215.AND.ID.LE.218)) THEN
804 IF (NPR.EQ.3.AND.NME(IM).EQ.100) THEN
806 C usual heavy quark decay
820 ELSEIF (ABS(IDHEP(MHEP)).EQ.37) THEN
822 C heavy quark to charged Higgs
832 ELSEIF (ABS(IDHEP(NHEP)).EQ.37) THEN
844 CALL HWWARN('HWDHOB',103,*999)
852 & ((NPR.EQ.2.AND.ID.GE.401.AND.ID.LT.448.AND.
854 & IDHW(MHEP).LE.132.AND.IDHW(NHEP).LE.132)
856 & .OR.(NPR.EQ.3.AND.ID.GE.449.AND.ID.LE.457.AND.
858 & IDHW(MHEP).LE.132.AND.IDHW(NHEP).LE.132.AND.
860 & IDHW(MHEP-1).LE.132))) THEN
862 C R-parity violating SUSY decays
866 C--Rparity slepton colour connections
868 IF(ID.GE.425.AND.ID.LE.448) THEN
870 IF(IDHW(MHEP).GT.12) THEN
872 JMOHEP(2,MHEP) = MHEP
874 JDAHEP(2,MHEP) = MHEP
876 JMOHEP(2,NHEP) = NHEP
878 JDAHEP(2,NHEP) = NHEP
882 JMOHEP(2,MHEP) = NHEP
884 JDAHEP(2,MHEP) = NHEP
886 JMOHEP(2,NHEP) = MHEP
888 JDAHEP(2,NHEP) = MHEP
892 C--Rparity squark colour connections
896 IF(IDHEP(LHEP).GT.0) THEN
898 C--LQD decay colour connections
900 IF(IDHW(MHEP).GT.12) THEN
902 JMOHEP(2,MHEP) = MHEP
904 JDAHEP(2,MHEP) = MHEP
906 JMOHEP(2,NHEP) = LHEP
908 JDAHEP(2,NHEP) = LHEP
912 C--UDD decay colour connections
916 CALL HWDRCL(LHEP,MHEP,CLSAVE)
922 C--Antisquark connections
924 IF(IDHW(MHEP).GT.12) THEN
926 JMOHEP(2,MHEP) = MHEP
928 JDAHEP(2,MHEP) = MHEP
930 JMOHEP(2,NHEP) = LHEP
932 JDAHEP(2,NHEP) = LHEP
938 CALL HWDRCL(LHEP,MHEP,CLSAVE)
948 IF(ID.GE.450.AND.ID.LE.457) THEN
950 C--Rparity Neutralino/Chargino colour connection
952 IF(IDHW(MHEP-1).LE.12.AND.IDHW(MHEP).LE.12.
954 & AND.IDHW(NHEP).LE.12) THEN
958 CALL HWDRCL(LHEP,MHEP,CLSAVE)
962 JMOHEP(2,MHEP) = NHEP
964 JDAHEP(2,MHEP) = NHEP
966 JMOHEP(2,NHEP) = MHEP
968 JDAHEP(2,NHEP) = MHEP
972 C--Rparity gluino colour connections
974 ELSEIF(ID.EQ.449) THEN
976 IF(IDHW(MHEP-1).LE.12.AND.IDHW(MHEP).LE.12.
978 & AND.IDHW(NHEP).LE.12) THEN
982 CALL HWDRCL(LHEP,MHEP,CLSAVE)
984 C--Now the lepton number violating decay
988 IF(IDHW(MHEP).LE.6) THEN
990 JMOHEP(2,MHEP) = LHEP
992 JDAHEP(2,MHEP) = NHEP
994 JMOHEP(2,NHEP) = MHEP
996 JDAHEP(2,NHEP) = LHEP
1000 JMOHEP(2,MHEP) = NHEP
1002 JDAHEP(2,MHEP) = LHEP
1004 JMOHEP(2,NHEP) = LHEP
1006 JDAHEP(2,NHEP) = MHEP
1014 CALL HWWARN('HWDHOB',104,*999)
1022 C Normal SUSY decays
1024 IF (ID.LE.448.AND.ID.GT.207) THEN
1026 C Squark (or slepton)
1028 IF (IDHW(MHEP).EQ.449) THEN
1030 IF (IDHEP(LHEP).GT.0) THEN
1054 IF(NPR.EQ.3.AND.IDHW(MHEP).LE.12) THEN
1078 ELSEIF (ID.EQ.449) THEN
1082 IF (IDHW(NHEP).EQ.13) THEN
1092 ELSEIF (IDHEP(MHEP).GT.0) THEN
1132 C---SPECIAL CASE FOR THREE-BODY TOP DECAYS:
1134 C RELABEL THEM AS TWO TWO-BODY DECAYS FOR PARTON SHOWERING
1136 IF ((ID.EQ.6.OR.ID.EQ.12).AND.NPR.EQ.3.AND.NME(IM).EQ.100) THEN
1138 C---STORE W DECAY PRODUCTS
1140 CALL HWVEQU(10,PHEP(1,KHEP),PDW)
1142 C---BOOST THEM INTO W REST FRAME
1144 CALL HWULOF(PW,PDW(1,1),PDW(1,3))
1146 C---REPLACE THEM BY W
1148 CALL HWVEQU(5,PW,PHEP(1,KHEP))
1154 IF (ID.EQ.12) IDHW(KHEP)=199
1156 IDHEP(KHEP)=IDPDG(IDHW(KHEP))
1162 CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,KHEP))
1166 CALL HWVEQU(5,PHEP(1,NHEP),PHEP(1,MHEP))
1168 IDHW(MHEP)=IDHW(NHEP)
1170 IDHEP(MHEP)=IDHEP(NHEP)
1174 JMOHEP(2,MHEP)=JMOHEP(2,NHEP)
1176 JDAHEP(2,MHEP)=JDAHEP(2,NHEP)
1178 CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,MHEP))
1182 C---DO PARTON SHOWER
1188 IF (IERROR.NE.0) RETURN
1190 C---FIND BOOSTED W MOMENTUM
1196 IF (NTRY.GT.NHEP.OR.WHEP.LE.0.OR.WHEP.GT.NHEP)
1198 $ CALL HWWARN('HWDHOB',101,*999)
1202 IF (ISTHEP(WHEP).NE.190) GOTO 41
1204 C---AND HENCE ITS CHILDRENS MOMENTA
1206 CALL HWULOB(PHEP(1,WHEP),PDW(1,3),PHEP(1,NHEP+1))
1208 CALL HWVDIF(4,PHEP(1,WHEP),PHEP(1,NHEP+1),PHEP(1,NHEP+2))
1210 PHEP(5,NHEP+2)=PDW(5,2)
1218 IDHW(NHEP+I)=IDKPRD(I,IM)
1220 IDHEP(NHEP+I)=IDPDG(IDHW(NHEP+I))
1222 ISTHEP(NHEP+I)=112+I
1224 JDAHEP(I,WHEP)=NHEP+I
1226 JMOHEP(1,NHEP+I)=WHEP
1228 JMOHEP(2,NHEP+I)=NHEP+3-I
1230 JDAHEP(2,NHEP+I)=NHEP+3-I
1236 C---ASSIGN PRODUCTION VERTICES TO 1 AND 2
1238 CALL HWUDKL(198,PW,VHEP(1,NHEP))
1240 CALL HWVSUM(4,VHEP(1,WHEP),VHEP(1,NHEP),VHEP(1,NHEP))
1242 CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-1))
1244 C---DO PARTON SHOWERS
1250 IF (IERROR.NE.0) RETURN
1260 IF (IERROR.NE.0) RETURN
1266 C--New to correct colour connections in Rslash
1268 IF(CLSAVE(1).NE.0) THEN
1272 ID = IDHW(CLSAVE(1))
1274 IDM = IDHW(JMOHEP(1,CLSAVE(1)))
1278 IF(IDM.EQ.15) ID=IDHW(JMOHEP(1,JMOHEP(1,CLSAVE(1))))
1280 IF((ID.LE.6.AND.((IDM.GE.419.AND.IDM.LE.424).OR.IDM.EQ.411.OR.
1284 & AND.((IDM2.GE.413.AND.IDM2.LE.418)
1286 & .OR.IDM2.EQ.449).OR.IDM2.EQ.405.OR.IDM2.EQ.406)
1288 & .OR.(ID.LE.6.AND.IDM.EQ.449.AND.
1290 & (((IDM2.GE.413.AND.IDM2.LE.418).OR.IDM2.EQ.405.OR.IDM2.EQ.406)
1292 & .OR.IDM2.EQ.449)).OR.
1294 & (IDM.EQ.15.AND.ID.LE.12.AND.ID.GE.7.AND.((IDM2.GE.413.AND.
1296 & IDM2.LE.418).OR.IDM2.EQ.449.OR.IDM2.
1298 & EQ.405.OR.IDM2.EQ.406))) THEN
1300 IF(JMOHEP(2,CLSAVE(1)).EQ.MHEP) THEN
1302 IF(IDHW(CLSAVE(1)).NE.13.AND.IDHW(CLSAVE(1)).NE.449)
1304 & JMOHEP(2,CLSAVE(2)) = THEP
1306 JDAHEP(2,MHEP) = CLSAVE(1)
1308 JDAHEP(2,THEP) = CLSAVE(2)
1312 IF(IDHW(CLSAVE(1)).NE.13.AND.IDHW(CLSAVE(1)).NE.449)
1314 & JMOHEP(2,CLSAVE(2)) = MHEP
1316 JDAHEP(2,MHEP) = CLSAVE(2)
1318 JDAHEP(2,THEP) = CLSAVE(1)
1322 ELSEIF((ID.GT.6.AND.ID.LE.12.
1324 & AND.((IDM.GE.413.AND.IDM.LE.418).OR.IDM.EQ.405.OR.
1328 & ((IDM2.GE.419.AND.IDM2.LE.424).OR.IDM2.EQ.449.OR.
1330 & IDM2.EQ.411.OR.IDM2.EQ.412)).OR.
1332 & (ID.GT.6.AND.ID.LE.12.AND.IDM.EQ.449.
1334 & AND.((IDM2.GE.419.AND.IDM2.LE.424).OR.IDM2.EQ.449.OR.
1336 & IDM2.EQ.411.OR.IDM2.EQ.412)).OR.
1338 & (IDM.EQ.15.AND.ID.LE.6.AND.((IDM2.GE.419.AND.
1340 & IDM2.LE.424).OR.IDM2.EQ.449.OR.IDM2.EQ.411.OR.
1342 & IDM2.EQ.412))) THEN
1344 IF(JDAHEP(2,CLSAVE(1)).EQ.MHEP) THEN
1346 JDAHEP(2,CLSAVE(2))=THEP
1348 JMOHEP(2,MHEP)=CLSAVE(1)
1350 JMOHEP(2,THEP)=CLSAVE(2)
1354 JDAHEP(2,CLSAVE(2))=MHEP
1356 JMOHEP(2,MHEP)=CLSAVE(2)
1358 JMOHEP(2,THEP)=CLSAVE(1)
1370 IF (IHEP.EQ.NHEP) GOTO 70
1376 C Fix any SUSY colour disconnections
1380 IF (ISTHEP(IHEP).GE.147.AND.ISTHEP(IHEP).LE.151
1382 & .AND.JDAHEP(2,IHEP).EQ.0) THEN
1386 C Chase connection back through SUSY decays
1392 IF (ISM.EQ.120) GOTO 80
1394 IF (ISM.NE.123.AND.ISM.NE.124.AND.ISM.NE.155) GOTO 75
1396 C Look for unclustered parton to connect
1400 IF (ISTHEP(JHEP).GE.147.AND.ISTHEP(JHEP).LE.151) THEN
1420 C Not found: need to go further back
1428 C Go back to check for further heavy decay products
1438 *CMZ :- -26/04/91 12.19.24 by Federico Carminati
1440 *-- Author : Ian Knowles & Bryan Webber
1442 C-----------------------------------------------------------------------
1446 C-----------------------------------------------------------------------
1448 C Performs partonic decays of hadrons containing heavy quark(s):
1450 C either, meson/baryon spectator model weak decays;
1452 C or, quarkonia -> 2-gluons, q-qbar, 3-gluons, or 2-gluons + photon.
1454 C-----------------------------------------------------------------------
1456 INCLUDE 'HERWIG61.INC'
1458 DOUBLE PRECISION HWULDO,HWR,XS,XB,EMWSQ,GMWSQ,EMLIM,PW(4),
1460 & EMTST,X1,X2,X3,TEST,HWDWWT,HWDPWT
1462 INTEGER IST(3),I,IHEP,IM,ID,IDQ,IQ,IS,J
1464 EXTERNAL HWR,HWDWWT,HWDPWT,HWULDO
1466 DATA IST/113,114,114/
1468 IF (IERROR.NE.0) RETURN
1482 IF (ISTHEP(IHEP).EQ.199) GOTO 100
1486 IF (NHEP+NPRODS(IM).GT.NMXHEP) CALL HWWARN('HWDHVY',100,*999)
1488 IF (IDKPRD(4,IM).NE.0) THEN
1490 C Weak decay of meson or baryon
1492 C Idenitify decaying heavy quark and spectator
1496 IF (ID.EQ.136.OR.ID.EQ.140.OR.ID.EQ.144.OR.
1498 & ID.EQ.150.OR.ID.EQ.155.OR.ID.EQ.158.OR.ID.EQ.161.OR.
1500 & (ID.EQ.254.AND.IDKPRD(4,IM).EQ.11)) THEN
1502 C c hadron or c decay of B_c+
1510 ELSEIF (ID.EQ.171.OR.ID.EQ.175.OR.ID.EQ.179.OR.
1512 & ID.EQ.185.OR.ID.EQ.190.OR.ID.EQ.194.OR.ID.EQ.196.OR.
1514 & (ID.EQ.230.AND.IDKPRD(4,IM).EQ.5)) THEN
1516 C cbar hadron or cbar decay of B_c-
1524 ELSEIF ((ID.GE.221.AND.ID.LE.229).OR.
1526 & (ID.EQ.230.AND.IDKPRD(4,IM).EQ.10)) THEN
1528 C b hadron or b decay of B_c-
1536 ELSEIF ((ID.GE.245.AND.ID.LE.253).OR.
1538 & (ID.EQ.254.AND.IDKPRD(4,IM).EQ.4)) THEN
1540 C bbar hadron or bbar decay of B_c+
1550 C Decay not recognized
1552 CALL HWWARN('HWDHVY',101,*999)
1556 C Label constituents
1558 IF (NHEP+5.GT.NMXHEP) CALL HWWARN('HWDHVY',102,*999)
1562 JDAHEP(1,IHEP)=NHEP+1
1564 JDAHEP(2,IHEP)=NHEP+2
1568 IDHW(IS)=IDKPRD(4,IM)
1570 IDHEP(IQ)=IDPDG(IDQ)
1572 IDHEP(IS)=IDPDG(IDKPRD(4,IM))
1596 C and weak decay product jets
1602 IDHW(NHEP)=IDKPRD(J,IM)
1604 IDHEP(NHEP)=IDPDG(IDKPRD(J,IM))
1612 10 PHEP(5,NHEP)=RMASS(IDKPRD(J,IM))
1614 JMOHEP(2,NHEP-2)=NHEP-1
1616 JDAHEP(2,NHEP-2)=NHEP-1
1618 JMOHEP(2,NHEP-1)=NHEP-2
1620 JDAHEP(2,NHEP-1)=NHEP-2
1626 C Share momenta in ratio of masses, preserving specator mass
1628 XS=RMASS(IDHW(IS))/PHEP(5,IHEP)
1632 CALL HWVSCA(5,XB,PHEP(1,IHEP),PHEP(1,IQ))
1634 CALL HWVSCA(5,XS,PHEP(1,IHEP),PHEP(1,IS))
1636 IF (NME(IM).EQ.100) THEN
1638 C Generate decay momenta using full (V-A)*(V-A) matrix element
1642 GMWSQ=(RMASS(198)*GAMW)**2
1644 EMLIM=GMWSQ+(EMWSQ-(PHEP(5,IQ)-PHEP(5,NHEP))**2)**2
1646 20 CALL HWDTHR(PHEP(1,IQ ),PHEP(1,NHEP-1),
1648 & PHEP(1,NHEP-2),PHEP(1,NHEP),HWDWWT)
1650 CALL HWVSUM(4,PHEP(1,NHEP-2),PHEP(1,NHEP-1),PW)
1652 EMTST=(HWULDO(PW,PW)-EMWSQ)**2
1654 IF ((EMTST+GMWSQ)*HWR().GT.EMLIM) GOTO 20
1660 CALL HWDTHR(PHEP(1,IQ ),PHEP(1,NHEP-2),
1662 & PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
1664 CALL HWVSUM(4,PHEP(1,NHEP-2),PHEP(1,NHEP-1),PW)
1668 C Set up production vertices
1670 CALL HWVZRO(4,VHEP(1,IQ))
1672 CALL HWVEQU(4,VHEP(1,IQ),VHEP(1,IS))
1674 CALL HWVEQU(4,VHEP(1,IQ),VHEP(1,NHEP))
1676 CALL HWUDKL(198,PW,VHEP(1,NHEP-2))
1678 CALL HWVSUM(4,VHEP(1,IQ),VHEP(1,NHEP-2),VHEP(1,NHEP-2))
1680 CALL HWVEQU(4,VHEP(1,NHEP-2),VHEP(1,NHEP-1))
1692 JDAHEP(1,IHEP)=NHEP+1
1694 DO 30 J=1,NPRODS(IM)
1698 IDHW(NHEP)=IDKPRD(J,IM)
1700 IDHEP(NHEP)=IDPDG(IDKPRD(J,IM))
1708 PHEP(5,NHEP)=RMASS(IDKPRD(J,IM))
1710 30 CALL HWVZRO(4,VHEP(1,NHEP))
1714 C Establish colour connections and select momentum configuration
1716 IF (NPRODS(IM).EQ.3) THEN
1718 IF (IDKPRD(3,IM).EQ.13) THEN
1722 JMOHEP(2,NHEP-2)=NHEP
1724 JMOHEP(2,NHEP-1)=NHEP-2
1726 JMOHEP(2,NHEP )=NHEP-1
1728 JDAHEP(2,NHEP-2)=NHEP-1
1730 JDAHEP(2,NHEP-1)=NHEP
1732 JDAHEP(2,NHEP )=NHEP-2
1736 C or 2-gluon + photon decay
1738 JMOHEP(2,NHEP-2)=NHEP-1
1740 JMOHEP(2,NHEP-1)=NHEP-2
1742 JMOHEP(2,NHEP )=NHEP
1744 JDAHEP(2,NHEP-2)=NHEP-1
1746 JDAHEP(2,NHEP-1)=NHEP-2
1748 JDAHEP(2,NHEP )=NHEP
1752 IF (NME(IM).EQ.130) THEN
1754 C Use Ore & Powell orthopositronium matrix element
1756 40 CALL HWDTHR(PHEP(1,IHEP),PHEP(1,NHEP-2),
1758 & PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
1760 X1=TWO*HWULDO(PHEP(1,IHEP),PHEP(1,NHEP-2))/PHEP(5,IHEP)**2
1762 X2=TWO*HWULDO(PHEP(1,IHEP),PHEP(1,NHEP-1))/PHEP(5,IHEP)**2
1766 TEST=((X1*(ONE-X1))**2+(X2*(ONE-X2))**2+(X3*(ONE-X3))**2)
1770 IF (TEST.LT.TWO*HWR()) GOTO 40
1776 CALL HWDTHR(PHEP(1,IHEP),PHEP(1,NHEP-2),
1778 & PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
1784 C Parapositronium 2-gluon or q-qbar decay
1786 JMOHEP(2,NHEP-1)=NHEP
1788 JMOHEP(2,NHEP )=NHEP-1
1790 JDAHEP(2,NHEP-1)=NHEP
1792 JDAHEP(2,NHEP )=NHEP-1
1794 CALL HWDTWO(PHEP(1,IHEP),PHEP(1,NHEP-1),
1796 & PHEP(1,NHEP),CMMOM(IM),TWO,.FALSE.)
1804 C Process this new hard scatter
1806 CALL HWVEQU(4,VTXQDK(1,I),VTXPIP)
1824 *CMZ :- -20/07/99 10:56:12 by Peter Richardson
1826 *-- Author : Peter Richardson
1828 C-----------------------------------------------------------------------
1830 SUBROUTINE HWDRCL(IHEP,MHEP,CLSAVE)
1832 C-----------------------------------------------------------------------
1834 C Sets the colour connections in Baryon number violating decays
1836 C-----------------------------------------------------------------------
1838 INCLUDE 'HERWIG61.INC'
1840 INTEGER IHEP,MHEP,ID,ID2,IDM2,IDM3,COLCON(2,2,3),FLACON(2,3),JHEP,
1842 & DECAY,COLANT,KHEP,IDM,IDMB,IDMB2,IDMB3,IDMB4,QHEP,IDM4,
1844 & CLSAVE(2),XHEP,I,HWRINT,THEP
1848 C--Colour connections for the decays
1850 DATA COLCON/-1,1,-1,-2,-2,1,-3,-1,-1,1,-2,-1/
1852 DATA FLACON/1,-1,1,-1,-1,0/
1854 C--identify the decay
1856 IF(IERROR.NE.0) RETURN
1862 IF(ID.GE.450.AND.ID.LE.457) THEN
1866 ELSEIF(ID.EQ.449) THEN
1870 ELSEIF((ID.GE.411.AND.ID.LE.424).OR.ID.EQ.405.OR.ID.EQ.406) THEN
1878 CALL HWWARN('HWDRCL',100,*999)
1884 C--identify the colour partner
1886 IF(DECAY.GT.1.AND.ID2.LE.6) THEN
1892 KHEP = JDAHEP(2,IHEP-1)
1894 ELSEIF(DECAY.GT.1.AND.ID2.GE.7) THEN
1896 C--anticolour partner
1900 KHEP = JMOHEP(2,IHEP)
1908 IDM = IDHW(JMOHEP(1,KHEP))
1910 IF(ABS(IDPDG(IDM)).GT.1000000.OR.IDM.EQ.15) THEN
1912 IDM2 = IDHW(JDAHEP(1,JMOHEP(1,KHEP)))
1914 IDM3 = IDHW(JDAHEP(2,JMOHEP(1,KHEP)))
1916 IDM4 = IDHW(JDAHEP(2,JMOHEP(1,KHEP))-1)
1918 QHEP = JMOHEP(1,KHEP)
1920 IDMB = IDHW(JMOHEP(1,QHEP))
1922 IDMB2 = IDHW(JMOHEP(2,QHEP))
1924 IDMB3 = IDHW(JDAHEP(1,QHEP))
1926 IDMB4 = IDHW(JDAHEP(2,QHEP))
1930 C--Now decide if the colour partner decayed via BV
1932 IF(COLANT.EQ.2.AND.((((IDM.GE.413.AND.IDM.LE.418).OR.
1934 & IDM.EQ.449.OR.IDM.EQ.405.OR.IDM.EQ.406).AND.
1936 & (IDM2.GE.7.AND.IDM2.LE.12.AND.
1938 & IDM3.GE.7.AND.IDM3.LE.12.AND.
1940 & IDM4.GE.7.AND.IDM4.LE.12)).OR.
1942 & (IDM.EQ.15.AND.IDMB.LE.6.AND.IDMB2.LE.6.AND.
1944 & ((IDMB3.GE.7.AND.IDMB4.GE.12.AND.IDMB4.EQ.449).OR.
1946 & (IDMB3.GE.198.AND.IDMB3.LE.207.AND.
1948 & ABS(IDPDG(IDMB4)).GT.1000000))))) THEN
1956 XHEP = JMOHEP(2,JDAHEP(2,JMOHEP(1,KHEP)))
1958 ELSEIF(COLANT.EQ.3.AND.((((IDM.GE.419.AND.IDM.LE.424).OR.
1960 & IDM.EQ.449.OR.IDM.EQ.411.OR.IDM.EQ.412).AND.
1962 & (IDM2.LE.6.AND.IDM3.LE.6.AND.IDM4.LE.6)).OR.
1964 & (IDM.EQ.15.AND.IDMB.GE.7.AND.IDMB.LE.12.AND.
1966 & IDMB2.GE.7.AND.IDMB2.LE.12.AND.((IDMB3.LE.6.AND.
1968 & IDMB4.EQ.449).OR.(ABS(IDPDG(IDMB4)).GT.1000000
1970 & .AND.IDMB3.GE.198.AND.IDMB3.LE.207))))) THEN
1978 XHEP = JDAHEP(2,JDAHEP(2,JMOHEP(1,KHEP)))
1994 CLSAVE(1) = JDAHEP(2,JMOHEP(1,KHEP))-1
1996 CLSAVE(2) = CLSAVE(1)+1
2000 IF(IDMB4.EQ.449) THEN
2004 CLSAVE(I) = JMOHEP(I,JMOHEP(1,KHEP))
2006 IF(CLSAVE(I).EQ.XHEP) CLSAVE(I)=JDAHEP(1,JMOHEP(1,KHEP))
2012 CLSAVE(1) = JMOHEP(1,JMOHEP(1,KHEP))
2014 CLSAVE(2) = JMOHEP(2,JMOHEP(1,KHEP))
2028 C--Now set the colours for angular ordering
2036 JMOHEP(2,THEP) = THEP+HWRINT(1,2)
2038 JDAHEP(2,THEP) = THEP
2042 JMOHEP(2,THEP) = THEP
2044 JDAHEP(2,THEP) = THEP+HWRINT(1,2)
2048 ELSEIF(DECAY.EQ.2) THEN
2052 JMOHEP(2,THEP) = IHEP
2054 JDAHEP(2,THEP) = THEP
2058 JMOHEP(2,THEP) = THEP
2060 JDAHEP(2,THEP) = IHEP
2066 C--Colour of the second two
2072 JMOHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+
2074 & COLCON(HWRINT(1,2),JHEP,DECAY)
2076 JDAHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+FLACON(JHEP,DECAY)
2080 JDAHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+
2082 & COLCON(HWRINT(1,2),JHEP,DECAY)
2084 JMOHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+FLACON(JHEP,DECAY)
2090 C--Now set the colours of the colour partner
2092 IF(DECAY.GT.1.AND..NOT.CONBV) THEN
2094 IF(ID2.LE.6) JMOHEP(2,KHEP) = MHEP+HWRINT(0,1)
2096 IF(ID2.GE.7) JDAHEP(2,KHEP) = MHEP+HWRINT(0,1)
2102 JMOHEP(2,CLSAVE(1)) = MHEP+HWRINT(0,1)
2104 IF(JMOHEP(2,CLSAVE(1)).EQ.MHEP) THEN
2106 JMOHEP(2,CLSAVE(2)) = MHEP+1
2110 JMOHEP(2,CLSAVE(2)) = MHEP
2116 JDAHEP(2,CLSAVE(1)) = MHEP+HWRINT(0,1)
2118 IF(JDAHEP(2,CLSAVE(1)).EQ.MHEP) THEN
2120 JDAHEP(2,CLSAVE(2)) = MHEP+1
2124 JDAHEP(2,CLSAVE(2)) = MHEP
2136 *CMZ :- -20/07/99 10:56:12 by Peter Richardson
2138 *-- Author : Peter Richardson
2140 C-----------------------------------------------------------------------
2142 SUBROUTINE HWDRME(LHEP,MHEP)
2144 C-----------------------------------------------------------------------
2146 C SUBROUTINE TO IMPLEMENT ALL RPARITY DECAY MATRIX ELEMENTS
2148 C-----------------------------------------------------------------------
2150 INCLUDE 'HERWIG61.INC'
2152 DOUBLE PRECISION SM(6),SW(6),HWULDO,INFCOL,AM, M12SQ,M23SQ,MSGN,
2154 & M13SQ,A(6),B(6),SWEAK,MW,DECMOM(5),TEST(4),EPS,
2156 & M12SQT(6),M23SQT(6),M13SQT(6),LIMIT,M(4),RAND,
2158 & MC(2),MX2(6),MX(6),HWDPWT,HWR,HWDRM1,LAMD(3)
2160 EXTERNAL HWDRM1,HWULDO,HWDPWT,HWR
2162 INTEGER K,SN(3),LHEP,CSP,I,SB(3),J,ND,LTRY,MHEP,NSP,ID(3),IG,
2164 & IDHWTP,IDHPTP,MTRY
2166 PARAMETER(EPS=1D-20)
2168 IF(IERROR.NE.0) RETURN
2170 C--Electroweak parameters, etc
2180 C--Find the masses of the final state and zero parameters
2184 ID(K) = IDHW(MHEP+K-1)
2186 IF(ID(K).LE.12) THEN
2196 IF(SN(K).GT.6) SN(K)=SN(K)-6
2198 M(K) = PHEP(5,LHEP+K)
2220 C--Evaluate the coefficents for the mode we want
2222 IF(IG.GE.450.AND.IG.LE.453) THEN
2232 MC(1) = ZMIXSS(NSP,3)/(2*MW*COSB*SWEAK)
2234 MC(2) = ZMIXSS(NSP,4)/(2*MW*SINB*SWEAK)
2236 C--Calculate the combinations of couplings needed
2238 IF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
2240 C--first for the UDD modes
2244 A(J) = M(1)*MC(2)*QMIXSS(SN(1),2,J)
2246 & +SLFCH(SN(1),NSP)*QMIXSS(SN(1),1,J)
2248 B(J) = MSGN*(M(1)*MC(2)*QMIXSS(SN(1),1,J)
2250 & +SRFCH(SN(1),NSP)*QMIXSS(SN(1),2,J))
2252 MX2(J) = QMIXSS(SN(1),2,J)
2254 A(J+2) = M(2)*MC(1)*QMIXSS(SN(2),2,J)
2256 & +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
2258 B(J+2) = MSGN*(M(2)*MC(1)*QMIXSS(SN(2),1,J)
2260 & +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
2262 MX2(J+2) = QMIXSS(SN(2),2,J)
2264 A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
2266 & +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
2268 B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
2270 & +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
2272 MX2(J+2) = QMIXSS(SN(3),2,J)
2284 ELSEIF(ID(1).GE.121.AND.ID(2).GE.121.AND.ID(3).GE.121) THEN
2286 C--Now for the LLE modes
2290 A(J) = MSGN*(M(1)*MC(1)*LMIXSS(SN(1),1,J)
2292 & +SRFCH(10+SN(1),NSP)*LMIXSS(SN(1),2,J))
2294 B(J) = M(1)*MC(1)*LMIXSS(SN(1),2,J)
2296 & +SLFCH(10+SN(1),NSP)*LMIXSS(SN(2),1,J)
2298 MX2(J)= LMIXSS(SN(1),1,J)
2302 B(J+2) = SLFCH(10+SN(2),NSP)*LMIXSS(SN(2),1,J)
2304 MX2(J+2) = LMIXSS(SN(2),1,J)
2306 A(J+4) = M(3)*MC(1)*LMIXSS(SN(3),2,J)
2308 & +SLFCH(10+SN(3),NSP)*LMIXSS(SN(3),1,J)
2310 B(J+4) = MSGN*(M(3)*MC(1)*LMIXSS(SN(3),1,J)
2312 & +SRFCH(10+SN(3),NSP)*LMIXSS(SN(3),2,J))
2314 MX2(4+J) = LMIXSS(SN(3),2,J)
2328 C--Now for both types of LQD modes
2330 IF(MOD(SN(1),2).EQ.0) THEN
2332 C--First the neutrino,down,antidown mode
2338 B(J) = SLFCH(10+SN(1),NSP)*
2342 MX2(J) = LMIXSS(SN(1),1,J)
2344 A(J+2) = MSGN*(M(2)*MC(1)*QMIXSS(SN(2),1,J)
2346 & +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
2348 B(J+2) = M(2)*MC(1)*QMIXSS(SN(2),2,J)
2350 & +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
2352 MX2(2+J) = QMIXSS(SN(2),1,J)
2354 A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
2356 & +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
2358 B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
2360 & +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
2362 MX2(J+4) = QMIXSS(SN(3),2,J)
2368 C--Now the charged lepton, antiup,down modes
2372 A(J) = MSGN*(M(1)*MC(1)*LMIXSS(SN(1),1,J)
2374 & +SRFCH(10+SN(1),NSP)*LMIXSS(SN(1),2,J))
2376 B(J) = M(1)*MC(1)*LMIXSS(SN(1),2,J)
2378 & +SLFCH(10+SN(1),NSP)*LMIXSS(SN(1),1,J)
2380 MX2(J) = LMIXSS(SN(1),1,J)
2382 A(J+2) =MSGN*(M(2)*MC(2)*QMIXSS(SN(2),1,J)
2384 & +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
2386 B(J+2) = M(2)*MC(2)*QMIXSS(SN(2),2,J)
2388 & +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
2390 MX2(2+J) = QMIXSS(SN(2),1,J)
2392 A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
2394 & +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
2396 B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
2398 & +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
2400 MX2(J+4) = QMIXSS(SN(3),2,J)
2422 SM(2*K-1) = RMASS(SN(K))
2424 SM(2*K) = RMASS(SB(K))
2426 SW(2*K-1) = HBAR/RLTIM(SN(K))
2428 SW(2*K) = HBAR/RLTIM(SB(K))
2442 ELSEIF(IG.EQ.449) THEN
2446 C--First obtian the masses and widths needed
2452 C--Calculate the combinations of couplings needed
2454 IF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
2456 C--first for the UDD modes
2466 A(2*I-2+J) = -QMIXSS(SN(I),1,J)
2468 B(2*I-2+J) = QMIXSS(SN(I),2,J)
2470 MX2(2*I-2+J) = QMIXSS(SN(I),2,J)
2484 C--Now for both types of LQD modes
2486 IF(MOD(SN(1),2).EQ.0) THEN
2488 C--First the neutrino,down,antidown mode
2498 A(J+2) = QMIXSS(SN(2),2,J)
2500 B(J+2) = -QMIXSS(SN(2),1,J)
2502 MX2(J+2) = QMIXSS(SN(2),1,J)
2504 A(J+4) = -QMIXSS(SN(3),1,J)
2506 B(J+4) = QMIXSS(SN(3),2,J)
2508 MX2(4+J) = QMIXSS(SN(3),2,J)
2512 ELSEIF(MOD(SN(1),2).EQ.1) THEN
2514 C--Now the charged lepton, antiup,down modes
2524 A(J+2) = QMIXSS(SN(2),2,J)
2526 B(J+2) = -QMIXSS(SN(2),1,J)
2528 MX2(J+2) = QMIXSS(SN(2),1,J)
2530 A(J+4) = -QMIXSS(SN(3),1,J)
2532 B(J+4) = QMIXSS(SN(3),2,J)
2534 MX2(J+4) = QMIXSS(SN(3),2,J)
2556 SM(2*K-1) = RMASS(SN(K))
2558 SM(2*K) = RMASS(SB(K))
2560 SW(2*K-1) = HBAR/RLTIM(SN(K))
2562 SW(2*K) = HBAR/RLTIM(SB(K))
2572 ELSEIF(IG.GE.454.AND.IG.LE.457) THEN
2578 IF(CSP.GT.2) CSP = CSP-2
2586 MC(1) = ONE/(SQRT(2.0D0)*MW*COSB)
2588 MC(2) = ONE/(SQRT(2.0D0)*MW*SINB)
2590 C--Calculate the combinations of the couplings needed
2592 IF(ID(1).GT.120.AND.ID(2).GT.120.AND.ID(3).GT.120) THEN
2594 C--first for the LLE modes, three modes
2596 IF(MOD(SN(1),2).EQ.0.AND.MOD(SN(3),2).EQ.0) THEN
2598 C--the one diagram mode nubar,positron, nu
2602 A(J+4) = LMIXSS(SN(3)-1,1,J)*WMXUSS(CSP,1)
2604 & -RMASS(SN(3)+119)*MC(1)*LMIXSS(SN(3)-1,2,J)*WMXUSS(CSP,2)
2608 MX2(J+4) = LMIXSS(SN(3)-1,2,J)
2618 ELSEIF(MOD(SN(1),2).EQ.0.AND.MOD(SN(2),2).EQ.0) THEN
2620 C--the first two diagram mode nu, nu, positron
2626 B(J) = LMIXSS(SN(1)-1,1,J)*WMXUSS(CSP,1)
2628 & -RMASS(SN(1)+119)*MC(1)*LMIXSS(SN(1)-1,2,J)*WMXUSS(CSP,2)
2632 B(J+2) = LMIXSS(SN(2)-1,1,J)*WMXUSS(CSP,1)
2634 & -RMASS(SN(2)+119)*MC(1)*LMIXSS(SN(2)-1,2,J)*WMXUSS(CSP,2)
2636 MX2(J) = LMIXSS(SN(1)-1,1,J)
2638 MX2(J+2) = LMIXSS(SN(2)-1,1,J)
2654 C--the second two diagram mode positron, positron, electron
2658 A(J) = -M(1)*WMXUSS(CSP,2)*MC(1)*LMIXSS(SN(1)+1,1,J)
2660 B(J) = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(1)+1,1,J)
2662 A(J+2) = -M(2)*WMXUSS(CSP,2)*MC(1)*LMIXSS(SN(2)+1,1,J)
2664 B(J+2) = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(2)+1,1,J)
2666 MX2(J) = LMIXSS(SN(1)+1,1,J)
2668 MX2(J+2) = LMIXSS(SN(2)+1,1,J)
2690 ELSEIF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
2694 IF(MOD(SN(1),2).EQ.0) THEN
2698 LAMD(1) = LAMDA3(SN(2)/2,SN(1)/2,(SN(3)+1)/2)
2700 LAMD(2) = LAMDA3(SN(1)/2,SN(2)/2,(SN(3)+1)/2)
2704 A(J) = WMXUSS(CSP,1)*QMIXSS(SN(1)-1,1,J)
2706 & -RMASS(SN(1)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(1)-1,2,J)
2708 B(J) = -MSGN*M(2)*WMXVSS(CSP,2)*QMIXSS(SN(1)-1,1,J)
2710 A(J+2) = WMXUSS(CSP,1)*QMIXSS(SN(2)-1,1,J)
2712 & -RMASS(SN(2)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(2)-1,2,J)
2714 B(J+2) = -MSGN*M(2)*WMXVSS(CSP,2)*QMIXSS(SN(2)-1,1,J)
2716 MX2(J) = QMIXSS(SN(1)-1,2,J)
2718 MX2(J+2) = QMIXSS(SN(2)-1,2,J)
2734 C--three diagram mode
2736 LAMD(1) = LAMDA3((SN(1)+1)/2,(SN(2)+1)/2,(SN(3)+1)/2)
2738 LAMD(2) = LAMDA3((SN(2)+1)/2,(SN(1)+1)/2,(SN(3)+1)/2)
2740 LAMD(3) = LAMDA3((SN(3)+1)/2,(SN(2)+1)/2,(SN(1)+1)/2)
2746 A(J+2*I-2) = MSGN*(WMXVSS(CSP,1)*QMIXSS(SN(I)+1,1,J)
2748 & -RMASS(SN(I)+1)*MC(2)*WMXVSS(CSP,2)*QMIXSS(SN(I)+1,2,J))
2750 B(J+2*I-2) = -M(I)*MC(1)*WMXUSS(CSP,2)
2752 & *QMIXSS(SN(I)+1,1,J)
2754 MX2(J+2*I-2) = QMIXSS(SN(I)+1,2,J)
2770 C--now for the LQD modes
2772 IF(MOD(SN(2),2).EQ.1.AND.MOD(SN(3),2).EQ.0) THEN
2774 C--first one diagram mode nubar, dbar, up
2778 A(J+4) = -MSGN*M(3)*WMXVSS(CSP,2)*MC(2)*
2780 & QMIXSS(SN(3)-1,1,J)
2782 B(J+4) = WMXUSS(CSP,1)*QMIXSS(SN(3)-1,1,J)
2784 & -RMASS(SN(3)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(3)-1,2,1)
2786 MX2(J+4) = QMIXSS(SN(3)-1,2,J)
2796 ELSEIF(MOD(SN(2),2).EQ.0.AND.MOD(SN(3),2).EQ.0) THEN
2798 C--second one diagram mode positron, ubar, up
2802 A(J+4) = -MSGN*M(3)*WMXVSS(CSP,2)*MC(2)*
2804 & QMIXSS(SN(3)-1,1,J)
2806 B(J+4) = WMXUSS(CSP,1)*QMIXSS(SN(3)-1,1,J)
2808 & -RMASS(SN(3)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(3)-1,2,1)
2810 MX2(J+4) = QMIXSS(SN(3)-1,2,J)
2820 ELSEIF(MOD(SN(2),2).EQ.1.AND.MOD(SN(3),2).EQ.1) THEN
2822 C--first two diagram mode positron, dbar, down
2826 A(J) = -M(1)*MC(1)*WMXUSS(CSP,2)*LMIXSS(SN(1)+1,1,J)
2828 B(J) = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(2)+1,1,J)
2830 A(J+2) = -M(2)*WMXUSS(CSP,2)*MC(1)*QMIXSS(SN(2)+1,1,J)
2832 B(J+2) = MSGN*(WMXVSS(CSP,1)*QMIXSS(SN(2)+1,1,J)
2834 & -RMASS(SN(2)+1)*MC(2)*WMXVSS(CSP,2)*QMIXSS(SN(2)+1,2,J))
2836 MX2(J) = LMIXSS(SN(1)+1,1,J)
2838 MX2(J+2) = QMIXSS(SN(2)+1,1,J)
2854 C--second two diagram mode nu, up, dbar
2860 B(J) = WMXUSS(CSP,1)*LMIXSS(SN(1)-1,1,J)
2862 & -RMASS(119+SN(1))*MC(1)*WMXUSS(CSP,2)*LMIXSS(SN(1)-1,2,J)
2864 A(J+2) = -MSGN*M(2)*MC(2)*WMXVSS(CSP,2)*
2866 & QMIXSS(SN(2)-1,1,J)
2868 B(J+2) = WMXUSS(CSP,1)*QMIXSS(SN(2)-1,1,J)
2870 & -RMASS(SN(2)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(2)-1,2,J)
2872 MX2(J) = LMIXSS(SN(1)-1,1,J)
2874 MX2(J+2) = QMIXSS(SN(2)-1,1,J)
2912 SM(5) = RMASS(SN(3))
2914 SM(6) = RMASS(SB(3))
2916 SW(5) = HBAR/RLTIM(SN(3))
2918 SW(6) = HBAR/RLTIM(SB(3))
2924 SM(2*K-1) = RMASS(SN(K))
2926 SM(2*K) = RMASS(SB(K))
2928 SW(2*K-1) = HBAR/RLTIM(SN(K))
2930 SW(2*K) = HBAR/RLTIM(SB(K))
2944 CALL HWWARN('HWDRME',500,*999)
2948 C--Set mixing to zero if diagram not available
2950 IF((AM.LT.(ABS(SM(1))+M(1)).OR.ABS(SM(1)).LT.(M(2)+M(3)))
2952 & .AND.ABS(MX2(1)).GT.ZERO.AND.ND.NE.1) MX(1) = MX2(1)*LAMD(1)
2954 IF((AM.LT.(ABS(SM(2))+M(1)).OR.ABS(SM(2)).LT.(M(2)+M(3)))
2956 & .AND.ABS(MX2(2)).GT.ZERO.AND.ND.NE.1) MX(2) = MX2(2)*LAMD(1)
2958 IF((AM.LT.(ABS(SM(3))+M(2)).OR.ABS(SM(3)).LT.(M(1)+M(3)))
2960 & .AND.ABS(MX2(3)).GT.ZERO.AND.ND.NE.1) MX(3) = MX2(3)*LAMD(2)
2962 IF((AM.LT.(ABS(SM(4))+M(2)).OR.ABS(SM(4)).LT.(M(1)+M(3)))
2964 & .AND.ABS(MX2(4)).GT.ZERO.AND.ND.NE.1) MX(4) = MX2(4)*LAMD(2)
2966 IF((AM.LT.(ABS(SM(5))+M(3)).OR.ABS(SM(5)).LT.(M(1)+M(2)))
2968 & .AND.ABS(MX2(5)).GT.ZERO.AND.ND.NE.2) MX(5) = MX2(5)*LAMD(3)
2970 IF((AM.LT.(ABS(SM(6))+M(3)).OR.ABS(SM(6)).LT.(M(1)+M(2)))
2972 & .AND.ABS(MX2(6)).GT.ZERO.AND.ND.NE.2) MX(6) = MX2(6)*LAMD(3)
2974 C--Calculate the limiting points
2980 IF(ABS(MX(J)).GT.EPS) CALL HWDRM5(M23SQT(J),M13SQT(J),
2982 & M12SQT(J),A(J),B(J),M(2),M(3),M(1),M(4),SM(J),SW(J))
2984 IF(ABS(MX(J+2)).GT.EPS) CALL HWDRM5(M13SQT(2+J),M23SQT(2+J),
2986 & M12SQT(2+J),A(2+J),B(2+J),M(1),M(3),M(2),M(4),SM(2+J),SW(2+J))
2992 IF(ABS(MX(J+4)).GT.EPS) CALL HWDRM5(M12SQT(4+J),M23SQT(4+J),
2994 & M13SQT(4+J),A(4+J),B(4+J),M(1),M(2),M(3),M(4),SM(4+J),SW(4+J))
3000 C--Now evaluate the limit using these points
3006 IF(ABS(MX(I)).LT.EPS) GOTO 100
3008 LIMIT = LIMIT+HWDRM1(TEST,M12SQT(I),M13SQT(I),M23SQT(I),A,B,MX,
3010 & M,SM,SW,INFCOL,AM,0,ND)
3016 C--Now evaluate at a random point
3026 CALL HWDTHR(PHEP(1,LHEP),PHEP(1,MHEP),
3028 & PHEP(1,MHEP+1),PHEP(1,MHEP+2),HWDPWT)
3030 C--Now calculate the m12sq etc for the actual point
3032 M12SQ = M(1)**2+M(2)**2+2*HWULDO(PHEP(1,MHEP),PHEP(1,MHEP+1))
3034 M13SQ = M(1)**2+M(3)**2+2*HWULDO(PHEP(1,MHEP),PHEP(1,MHEP+2))
3036 M23SQ = M(2)**2+M(3)**2+2*HWULDO(PHEP(1,MHEP+1),PHEP(1,MHEP+2))
3038 C--Now calulate the matrix element
3040 TEST(4) = HWDRM1(TEST,M12SQ,M13SQ,M23SQ,A,B,MX,
3042 & M,SM,SW,INFCOL,AM,1,ND)
3044 C--Now test the value againest the limit
3048 IF(TEST(4).GT.LIMIT) THEN
3050 LIMIT = 1.1D0*TEST(4)
3052 CALL HWWARN('HWDRME',51,*150)
3056 150 IF(TEST(4).LT.RAND.AND.LTRY.LT.NETRY) THEN
3060 ELSEIF(LTRY.GE.NETRY) THEN
3062 IF(MTRY.LE.NETRY) THEN
3066 CALL HWWARN('HWDRME',52,*25)
3070 CALL HWWARN('HWDRME',100,*999)
3076 C--Reorder the particles in gluino decay to get angular ordering right
3078 IF(IG.EQ.449.AND.ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
3082 IF(TEST(LTRY).GT.RAND) THEN
3088 IDHW(MHEP) = IDHW(MHEP+1)
3090 IDHW(MHEP+1) = IDHWTP
3092 IDHPTP = IDHEP(MHEP)
3094 IDHEP(MHEP) = IDHEP(MHEP+1)
3096 IDHEP(MHEP+1) = IDHPTP
3098 CALL HWVEQU(5,PHEP(1,MHEP),DECMOM)
3100 CALL HWVEQU(5,PHEP(1,MHEP+1),PHEP(1,MHEP))
3102 CALL HWVEQU(5,DECMOM,PHEP(1,MHEP+1))
3104 ELSEIF(LTRY.EQ.3) THEN
3108 IDHW(MHEP) = IDHW(MHEP+2)
3110 IDHW(MHEP+2) = IDHWTP
3112 IDHPTP = IDHEP(MHEP)
3114 IDHEP(MHEP) = IDHEP(MHEP+2)
3116 IDHEP(MHEP+2) = IDHPTP
3120 CALL HWVEQU(5,PHEP(1,MHEP),DECMOM)
3122 CALL HWVEQU(5,PHEP(1,MHEP+2),PHEP(1,MHEP))
3124 CALL HWVEQU(5,DECMOM,PHEP(1,MHEP+2))
3134 RAND=RAND-TEST(LTRY)
3146 *CMZ :- -20/07/99 10:56:12 by Peter Richardson
3148 *-- Author : Peter Richardson
3150 C-----------------------------------------------------------------------
3152 FUNCTION HWDRM1(TEST,M12SQ,M13SQ,M23SQ,A,B,MX,M,SM,SW
3156 C-----------------------------------------------------------------------
3158 C FUNCTION TO GIVE THE R-PARITY VIOLATING MATRIX ELEMENT AT A GIVEN
3162 C-----------------------------------------------------------------------
3166 DOUBLE PRECISION M12SQ,M13SQ,M23SQ,MX(6),A(6),B(6),SM(6),SW(6),
3168 & INFCOL,AM,TERM(21),TEST(4),PLN,NPLN,ZERO,
3170 & M(4),HWDRM1,HWDRM2,HWDRM3,HWDRM4
3174 EXTERNAL HWDRM2,HWDRM3,HWDRM4
3190 IF(ABS(MX(1)).GT.ZERO.AND.ND.NE.1) THEN
3192 TERM(1) = MX(1)**2*HWDRM2(M23SQ,M(2),M(3),M(1),M(4),SM(1),
3196 IF(ABS(MX(2)).GT.ZERO) TERM(7)= MX(1)*MX(2)*HWDRM3(M23SQ,M(2),
3198 & M(3),M(1),M(4),SM(1),SM(2),SW(1),SW(2),A(1),A(2),B(1),B(2))
3200 IF(ABS(MX(3)).GT.ZERO) TERM(10)=-MX(1)*MX(3)*HWDRM4(M13SQ,M23SQ,
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))
3204 IF(ABS(MX(4)).GT.ZERO) TERM(11)=-MX(1)*MX(4)*HWDRM4(M13SQ,M23SQ,
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))
3208 IF(ABS(MX(5)).GT.ZERO) TERM(12)=-MX(1)*MX(5)*HWDRM4(M23SQ,M12SQ,
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))
3212 IF(ABS(MX(6)).GT.ZERO) TERM(13)=-MX(1)*MX(6)*HWDRM4(M23SQ,M12SQ,
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))
3218 IF(ABS(MX(2)).GT.ZERO.AND.ND.NE.1) THEN
3220 TERM(2) = MX(2)**2*HWDRM2(M23SQ,M(2),M(3),M(1),M(4),SM(2),
3224 IF(ABS(MX(3)).GT.ZERO) TERM(14)=-MX(2)*MX(3)*HWDRM4(M13SQ,M23SQ,
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))
3228 IF(ABS(MX(4)).GT.ZERO) TERM(15)=-MX(2)*MX(4)*HWDRM4(M13SQ,M23SQ,
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))
3232 IF(ABS(MX(5)).GT.ZERO) TERM(16)=-MX(2)*MX(5)*HWDRM4(M23SQ,M12SQ,
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))
3236 IF(ABS(MX(6)).GT.ZERO) TERM(17)=-MX(2)*MX(6)*HWDRM4(M23SQ,M12SQ,
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))
3242 IF(ABS(MX(3)).GT.ZERO.AND.ND.NE.1) THEN
3244 TERM(3) = MX(3)**2*HWDRM2(M13SQ,M(1),M(3),M(2),M(4),SM(3),
3248 IF(ABS(MX(4)).GT.ZERO) TERM(8)= MX(3)*MX(4)*HWDRM3(M13SQ,M(1),
3250 & M(3),M(2),M(4),SM(3),SM(4),SW(3),SW(4),A(3),A(4),B(3),B(4))
3252 IF(ABS(MX(5)).GT.ZERO) TERM(18)=-MX(3)*MX(5)*HWDRM4(M12SQ,M13SQ,
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))
3256 IF(ABS(MX(6)).GT.ZERO) TERM(19)=-MX(3)*MX(6)*HWDRM4(M12SQ,M13SQ,
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))
3262 IF(ABS(MX(4)).GT.ZERO.AND.ND.NE.1) THEN
3264 TERM(4) = MX(4)**2*HWDRM2(M13SQ,M(1),M(3),M(2),M(4),SM(4),
3268 IF(ABS(MX(5)).GT.ZERO) TERM(20)=-MX(4)*MX(5)*HWDRM4(M12SQ,M13SQ,
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))
3272 IF(ABS(MX(6)).GT.ZERO) TERM(21)=-MX(4)*MX(6)*HWDRM4(M12SQ,M13SQ,
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))
3278 IF(ABS(MX(5)).GT.ZERO.AND.ND.NE.2) THEN
3280 TERM(5) = MX(5)**2*HWDRM2(M12SQ,M(1),M(2),M(3),M(4),SM(5),
3284 IF(ABS(MX(6)).GT.ZERO) TERM(9)= MX(5)*MX(6)*HWDRM3(M12SQ,M(1),
3286 & M(2),M(3),M(4),SM(5),SM(6),SW(5),SW(6),A(5),A(6),B(5),B(6))
3290 IF(ABS(MX(6)).GT.ZERO.AND.ND.NE.2) TERM(6) = MX(6)**2*
3292 & HWDRM2(M12SQ,M(1),M(2),M(3),M(4),SM(6),SW(6),A(6),B(6))
3296 TERM(K)=TERM(K)*INFCOL
3304 HWDRM1 = HWDRM1+TERM(K)
3308 C--Different colour flows for the gluino
3330 TEST(K) = (TERM(2*K-1)+TERM(2*K)+TERM(6+K))*(1+NPLN/PLN)
3344 IF(TEST(4).LT.ZERO) CALL HWWARN('HWDRM1',50,*999)
3350 *CMZ :- -20/07/99 10:56:12 by Peter Richardson
3352 *-- Author : Peter Richardson
3354 C-----------------------------------------------------------------------
3356 FUNCTION HWDRM2(X,MA,MB,MC,MD,MR1,GAM1,A,B)
3358 C-----------------------------------------------------------------------
3360 C Function to compute the matrix element squared part of a 3-body
3364 C-----------------------------------------------------------------------
3368 DOUBLE PRECISION X,MA,MB,MC,MD,A,B,HWDRM2,MR1,GAM1
3370 HWDRM2 = (X - MA**2 - MB**2)*(4*A*B*MC*MD +
3372 & (A**2 + B**2)*(-X + MC**2 + MD**2))/
3374 & ((X-MR1**2)**2+GAM1**2*MR1**2)