+++ /dev/null
-
-CDECK ID>, HWDHIG.
-
-*CMZ :- -24/04/92 14.23.44 by Mike Seymour
-
-*-- Author : Mike Seymour
-
-C-----------------------------------------------------------------------
-
- SUBROUTINE HWDHIG(GAMINP)
-
-C-----------------------------------------------------------------------
-
-C HIGGS DECAY ROUTINE
-
-C A) FOR GAMinp=0 FIND AND DECAY HIGGS
-
-C B) FOR GAMinp>0 CALCULATE TOTAL HIGGS WIDTH
-
-C FOR EMH=GAMINP. STORE RESULT IN GAMINP.
-
-C-----------------------------------------------------------------------
-
- INCLUDE 'HERWIG61.INC'
-
- DOUBLE PRECISION HWDHGF,HWR,HWRUNI,HWUSQR,HWUPCM,GAMINP,EMH,
-
- & EMF,COLFAC,ENF,K1,K0,BET0,BET1,GAM0,GAM1,SCLOG,CFAC,XF,EM,GAMLIM,
-
- & GAM,XW,EMW,XZ,EMZ,YW,YZ,EMI,TAUT,TAUW,WIDHIG,VECDEC,EMB,GAMB,
-
- & TMIN,TMAX1,EM1,TMAX2,EM2,X1,X2,PROB,PCM,SUMR,SUMI,TAUTR,TAUTI,
-
- & TAUWR,TAUWI,GFACTR
-
- INTEGER HWRINT,IHIG,I,IFERM,NLOOK,I1,I2,IPART,IMODE,IDEC,MMAX
-
- LOGICAL HWRLOG
-
- EXTERNAL HWDHGF,HWR,HWRUNI,HWUSQR,HWUPCM,HWRINT,HWRLOG
-
- SAVE GAM,EM,VECDEC
-
- PARAMETER (NLOOK=100)
-
- DIMENSION VECDEC(2,0:NLOOK)
-
- EQUIVALENCE (EMW,RMASS(198)),(EMZ,RMASS(200))
-
- DATA GAMLIM,GAM,EM/10D0,2*0D0/
-
-C---IF DECAY, FIND HIGGS (HWDHAD WILL HAVE GIVEN IT STATUS=1)
-
- IF (GAMINP.EQ.ZERO) THEN
-
- IHIG=0
-
- DO 10 I=1,NHEP
-
- 10 IF (IHIG.EQ.0.AND.IDHW(I).EQ.201.AND.ISTHEP(I).EQ.1) IHIG=I
-
- IF (IHIG.EQ.0) CALL HWWARN('HWDHIG',101,*999)
-
- EMH=PHEP(5,IHIG)
-
- IF (EMH.LE.ZERO) CALL HWWARN('HWDHIG',102,*999)
-
- EMSCA=EMH
-
- ELSE
-
- EMH=GAMINP
-
- IF (EMH.LE.ZERO) THEN
-
- GAMINP=0
-
- RETURN
-
- ENDIF
-
- ENDIF
-
-C---CALCULATE BRANCHING FRACTIONS
-
-C---FERMIONS
-
-C---NLL CORRECTION TO QUARK DECAY RATE (HHG eq 2.6-9)
-
- ENF=0
-
- DO 1 I=1,6
-
- 1 IF (2*RMASS(I).LT.EMH) ENF=ENF+1
-
- K1=5/PIFAC**2
-
- K0=3/(4*PIFAC**2)
-
- BET0=(11*CAFAC-2*ENF)/3
-
- BET1=(34*CAFAC**2-(10*CAFAC+6*CFFAC)*ENF)/3
-
- GAM0=-8
-
- GAM1=-404./3+40*ENF/9
-
- SCLOG=LOG(EMH**2/QCDLAM**2)
-
- CFAC=1 + ( K1/K0 - 2*GAM0 + GAM0*BET1/BET0**2*LOG(SCLOG)
-
- & + (GAM0*BET1-GAM1*BET0)/BET0**2) / (BET0*SCLOG)
-
- DO 100 IFERM=1,9
-
- IF (IFERM.LE.6) THEN
-
- EMF=RMASS(IFERM)
-
- XF=(EMF/EMH)**2
-
- COLFAC=FLOAT(NCOLO)
-
- IF (EMF.GT.QCDLAM)
-
- & EMF=EMF*(LOG(EMH/QCDLAM)/LOG(EMF/QCDLAM))**(GAM0/(2*BET0))
-
- ELSE
-
- EMF=RMASS(107+IFERM*2)
-
- XF=(EMF/EMH)**2
-
- COLFAC=1
-
- CFAC=1
-
- ENDIF
-
- IF (FOUR*XF.LT.ONE) THEN
-
- GFACTR=ALPHEM/(8.*SWEIN*EMW**2)
-
- BRHIG(IFERM)=COLFAC*GFACTR*EMH*EMF**2 * (1-4*XF)**1.5 * CFAC
-
- ELSE
-
- BRHIG(IFERM)=0
-
- ENDIF
-
- 100 CONTINUE
-
-C---W*W*/Z*Z*
-
- IF (ABS(EM-EMH).GE.GAMLIM*GAM) THEN
-
-C---OFF EDGE OF LOOK-UP TABLE
-
- XW=(EMW/EMH)**2
-
- XZ=(EMZ/EMH)**2
-
- YW=EMW*GAMW/EMH**2
-
- YZ=EMZ*GAMZ/EMH**2
-
- BRHIG(10)=.50*GFACTR * EMH**3 * HWDHGF(XW,YW)
-
- BRHIG(11)=.25*GFACTR * EMH**3 * HWDHGF(XZ,YZ)
-
- ELSE
-
-C---LOOK IT UP
-
- EMI=((EMH-EM)/(GAM*GAMLIM)+1)*NLOOK/2.0
-
- I1=INT(EMI)
-
- I2=INT(EMI+1)
-
- BRHIG(10)=.50*GFACTR * EMH**3 * ( VECDEC(1,I1)*(I2-EMI) +
-
- & VECDEC(1,I2)*(EMI-I1) )
-
- BRHIG(11)=.25*GFACTR * EMH**3 * ( VECDEC(2,I1)*(I2-EMI) +
-
- & VECDEC(2,I2)*(EMI-I1) )
-
- ENDIF
-
-C---GAMMAGAMMA
-
- TAUT=(2*RMASS(6)/EMH)**2
-
- TAUW=(2*EMW/EMH)**2
-
- CALL HWDHGC(TAUT,TAUTR,TAUTI)
-
- CALL HWDHGC(TAUW,TAUWR,TAUWI)
-
- SUMR=4./3*( - 2*TAUT*( 1 + (1-TAUT)*TAUTR ) ) * ENHANC(6)
-
- & +(2 + 3*TAUW*( 1 + (2-TAUW)*TAUWR ) ) * ENHANC(10)
-
- SUMI=4./3*( - 2*TAUT*( (1-TAUT)*TAUTI ) ) * ENHANC(6)
-
- & +( 3*TAUW*( (2-TAUW)*TAUWI ) ) * ENHANC(10)
-
- BRHIG(12)=GFACTR*.03125*(ALPHEM/PIFAC)**2
-
- & *EMH**3 * (SUMR**2 + SUMI**2)
-
- WIDHIG=0
-
- DO 200 IPART=1, 12
-
- IF (IPART.LT.12) BRHIG(IPART)=BRHIG(IPART)*ENHANC(IPART)**2
-
- 200 WIDHIG=WIDHIG+BRHIG(IPART)
-
- IF (WIDHIG.EQ.ZERO) CALL HWWARN('HWDHIG',103,*999)
-
- DO 300 IPART=1, 12
-
- 300 BRHIG(IPART)=BRHIG(IPART)/WIDHIG
-
- IF (EM.NE.RMASS(201)) THEN
-
-C---SET UP W*W*/Z*Z* LOOKUP TABLES
-
- EM=EMH
-
- GAM=WIDHIG
-
- GAMLIM=MAX(GAMLIM,GAMMAX)
-
- DO 400 I=0,NLOOK
-
- EMH=(I*2.0/NLOOK-1)*GAM*GAMLIM+EM
-
- XW=(EMW/EMH)**2
-
- XZ=(EMZ/EMH)**2
-
- YW=EMW*GAMW/EMH**2
-
- YZ=EMZ*GAMZ/EMH**2
-
- VECDEC(1,I)=HWDHGF(XW,YW)
-
- VECDEC(2,I)=HWDHGF(XZ,YZ)
-
- 400 CONTINUE
-
- EMH=EM
-
- ENDIF
-
- IF (GAMINP.GT.ZERO) THEN
-
- GAMINP=WIDHIG
-
- RETURN
-
- ENDIF
-
-C---SEE IF USER SPECIFIED A DECAY MODE
-
- IMODE=MOD(IPROC,100)
-
-C---IF NOT, CHOOSE ONE
-
- IF (IMODE.LT.1.OR.IMODE.GT.12) THEN
-
- MMAX=12
-
- IF (IMODE.LT.1) MMAX=6
-
- 500 IMODE=HWRINT(1,MMAX)
-
- IF (BRHIG(IMODE).LT.HWR()) GOTO 500
-
- ENDIF
-
-C---SEE IF SPECIFIED DECAY IS POSSIBLE
-
- IF (BRHIG(IMODE).EQ.ZERO) CALL HWWARN('HWDHIG',104,*999)
-
- IF (IMODE.LE.6) THEN
-
- IDEC=IMODE
-
- ELSEIF (IMODE.LE.9) THEN
-
- IDEC=107+IMODE*2
-
- ELSEIF (IMODE.EQ.10) THEN
-
- IDEC=198
-
- ELSEIF (IMODE.EQ.11) THEN
-
- IDEC=200
-
- ELSEIF (IMODE.EQ.12) THEN
-
- IDEC=59
-
- ENDIF
-
-C---STATUS, IDs AND POINTERS
-
- ISTHEP(IHIG)=195
-
- DO 600 I=1,2
-
- ISTHEP(NHEP+I)=193
-
- IDHW(NHEP+I)=IDEC
-
- IDHEP(NHEP+I)=IDPDG(IDEC)
-
- JDAHEP(I,IHIG)=NHEP+I
-
- JMOHEP(1,NHEP+I)=IHIG
-
- JMOHEP(2,NHEP+I)=NHEP+(3-I)
-
- JDAHEP(2,NHEP+I)=NHEP+(3-I)
-
- PHEP(5,NHEP+I)=RMASS(IDEC)
-
- IDEC=IDEC+6
-
- IF (IDEC.EQ.204) IDEC=199
-
- IF (IDEC.EQ.206) IDEC=200
-
- IF (IDEC.EQ. 65) IDEC= 59
-
- 600 CONTINUE
-
-C---ALLOW W/Z TO BE OFF-SHELL
-
- IF (IMODE.EQ.10.OR.IMODE.EQ.11) THEN
-
- IF (IMODE.EQ.10) THEN
-
- EMB=EMW
-
- GAMB=GAMW
-
- ELSE
-
- EMB=EMZ
-
- GAMB=GAMZ
-
- ENDIF
-
-C---STANDARD MASS DISTRIBUTION
-
- 700 TMIN=ATAN(-EMB/GAMB)
-
- TMAX1=ATAN((EMH**2/EMB-EMB)/GAMB)
-
- EM1=HWUSQR(EMB*(GAMB*TAN(HWRUNI(0,TMIN,TMAX1))+EMB))
-
- TMAX2=ATAN(((EMH-EM1)**2/EMB-EMB)/GAMB)
-
- EM2=HWUSQR(EMB*(GAMB*TAN(HWRUNI(0,TMIN,TMAX2))+EMB))
-
- X1=(EM1/EMH)**2
-
- X2=(EM2/EMH)**2
-
-C---CORRECT MASS DISTRIBUTION
-
- PROB=HWUSQR(1+X1**2+X2**2-2*X1-2*X2-2*X1*X2)
-
- & * ((X1+X2-1)**2 + 8*X1*X2)
-
- IF (.NOT.HWRLOG(PROB)) GOTO 700
-
-C---CALCULATE SPIN DENSITY MATRIX
-
- RHOHEP(1,NHEP+1)=4*X1*X2 / (8*X1*X2 + (X1+X2-1)**2)
-
- RHOHEP(2,NHEP+1)=(X1+X2-1)**2 / (8*X1*X2 + (X1+X2-1)**2)
-
- RHOHEP(3,NHEP+1)=RHOHEP(1,NHEP+1)
-
-C---SYMMETRIZE DISTRIBUTIONS IN PARTICLES 1,2
-
- IF (HWRLOG(HALF)) THEN
-
- PHEP(5,NHEP+1)=EM1
-
- PHEP(5,NHEP+2)=EM2
-
- ELSE
-
- PHEP(5,NHEP+1)=EM2
-
- PHEP(5,NHEP+2)=EM1
-
- ENDIF
-
- ENDIF
-
-C---DO DECAY
-
- PCM=HWUPCM(EMH,PHEP(5,NHEP+1),PHEP(5,NHEP+2))
-
- IF (PCM.LT.ZERO) CALL HWWARN('HWDHIG',105,*999)
-
- CALL HWDTWO(PHEP(1,IHIG),PHEP(1,NHEP+1),PHEP(1,NHEP+2),
-
- & PCM,TWO,.TRUE.)
-
- NHEP=NHEP+2
-
-C---IF QUARK DECAY, HADRONIZE
-
- IF (IMODE.LE.6) THEN
-
- ISTHEP(NHEP-1)=113
-
- ISTHEP(NHEP)=114
-
- CALL HWBGEN
-
- CALL HWDHOB
-
- CALL HWCFOR
-
- CALL HWCDEC
-
- ENDIF
-
- 999 END
-
-CDECK ID>, HWDHOB.
-
-*CMZ :- -20/10/99 09:46:43 by Peter Richardson
-
-*-- Author : Ian Knowles & Bryan Webber
-
-C-----------------------------------------------------------------------
-
- SUBROUTINE HWDHOB
-
-C-----------------------------------------------------------------------
-
-C Performs decays of heavy objects (heavy quarks & SUSY particles)
-
-C MODIFIED TO INCLUDE R-PARITY VIOLATING SUSY PR 9/4/99
-
-C-----------------------------------------------------------------------
-
- INCLUDE 'HERWIG61.INC'
-
- DOUBLE PRECISION HWUMBW,HWUPCM,HWR,SDKM,RN,BF,PCM,
-
- & EMMX,EMWSQ,GMWSQ,EMLIM,PW(5),EMTST,HWDPWT,HWDWWT,HWULDO,PDW(5,3)
-
- INTEGER IST(3),IHEP,IS,ID,IM,I,JHEP,KHEP,LHEP,MHEP,NPR,ISM,JCM,
-
- & MTRY,NTRY,IDM,IDM2,THEP,CLSAVE(2),WHEP,RHEP
-
- LOGICAL FOUND
-
- EXTERNAL HWR,HWDPWT,HWDWWT
-
- DATA IST/113,114,114/
-
- IF (IERROR.NE.0) RETURN
-
- 10 FOUND=.FALSE.
-
- CLSAVE(1) = 0
-
- CLSAVE(2) = 0
-
- DO 60 IHEP=1,NMXHEP
-
- IS=ISTHEP(IHEP)
-
- ID=IDHW(IHEP)
-
- IF (.NOT.RSTAB(ID).AND.(ID.EQ.6.OR.ID.EQ.12.OR.
-
- & (ID.GE.203.AND.ID.LE.218).OR.ABS(IDPDG(ID)).GT.1000000).AND.
-
- & (IS.EQ.190.OR.(IS.GE.147.AND.IS.LE.151))) THEN
-
- FOUND=.TRUE.
-
- IF(.NOT.RPARTY) THEN
-
- NHEP = NHEP+1
-
- ISTHEP(NHEP) = 3
-
- IDHW(NHEP) = 20
-
- IDHEP(NHEP) = 0
-
- CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP))
-
- CALL HWVEQU(4,VHEP(1,IHEP),VHEP(1,NHEP))
-
- JMOHEP(1,NHEP)=JMOHEP(1,IHEP)
-
- JMOHEP(2,NHEP)=JMOHEP(2,IHEP)
-
- JDAHEP(1,NHEP)=JDAHEP(1,IHEP)
-
- JDAHEP(2,NHEP)=JDAHEP(2,IHEP)
-
- ENDIF
-
-C Make a copy of decaying object
-
- NHEP=NHEP+1
-
- ISTHEP(NHEP)=155
-
- IDHW(NHEP)=IDHW(IHEP)
-
- IDHEP(NHEP)=IDHEP(IHEP)
-
- CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP))
-
- CALL HWVEQU(4,VHEP(1,IHEP),VHEP(1,NHEP))
-
- JMOHEP(1,NHEP)=JMOHEP(1,IHEP)
-
- JMOHEP(2,NHEP)=JMOHEP(2,IHEP)
-
- MTRY=0
-
- 15 MTRY=MTRY+1
-
-C Select decay mode
-
- RN=HWR()
-
- BF=0.
-
- IM=LSTRT(ID)
-
- DO 20 I=1,NMODES(ID)
-
- BF=BF+BRFRAC(IM)
-
- IF (BF.GE.RN) GOTO 30
-
- 20 IM=LNEXT(IM)
-
- CALL HWWARN('HWDHOB',50,*30)
-
- 30 IF (NHEP+5.GT.NMXHEP) CALL HWWARN('HWDHOB',100,*999)
-
- NPR=NPRODS(IM)
-
- JDAHEP(1,NHEP)=NHEP+1
-
- JDAHEP(2,NHEP)=NHEP+NPR
-
-C Reset colour pointers (if set)
-
- JHEP=JMOHEP(2,IHEP)
-
- IF (JHEP.GT.0) THEN
-
- IF (JDAHEP(2,JHEP).EQ.IHEP) JDAHEP(2,JHEP)=NHEP
-
- IF(.NOT.RPARTY.AND.ISTHEP(JHEP).EQ.155
-
- & .AND.ABS(IDHEP(JHEP)).GT.1000000
-
- & .AND.JDAHEP(2,JHEP-1).EQ.IHEP) JDAHEP(2,JHEP-1) = NHEP
-
- ENDIF
-
- JHEP=JDAHEP(2,IHEP)
-
- IF (JHEP.GT.0) THEN
-
- IF (JMOHEP(2,JHEP).EQ.IHEP) JMOHEP(2,JHEP)=NHEP
-
- IF(.NOT.RPARTY.AND.ISTHEP(JHEP).EQ.155
-
- & .AND.ABS(IDHEP(JHEP)).GT.1000000
-
- & .AND.JMOHEP(2,JHEP-1).EQ.IHEP) JMOHEP(2,JHEP-1) = NHEP
-
- ENDIF
-
-C--Reset colour pointers if baryon number violated
-
- IF(.NOT.RPARTY) THEN
-
- DO JHEP=1,NHEP
-
- IF(ISTHEP(JHEP).EQ.155
-
- & .AND.ABS(IDHEP(JHEP)).GT.1000000.AND.
-
- & JDAHEP(2,JHEP-1).EQ.IHEP) JDAHEP(2,JHEP-1)= NHEP
-
- IF(JDAHEP(2,JHEP).EQ.IHEP) JDAHEP(2,JHEP)=NHEP
-
- IF(JMOHEP(2,JHEP).EQ.IHEP) JMOHEP(2,JHEP)=NHEP
-
- ENDDO
-
- IF(HRDCOL(1,1).EQ.IHEP) HRDCOL(1,1)=NHEP
-
- ENDIF
-
-C Relabel original track
-
- ISTHEP(IHEP)=3
-
- JMOHEP(2,IHEP)=JMOHEP(1,IHEP)
-
- JDAHEP(1,IHEP)=NHEP
-
- JDAHEP(2,IHEP)=NHEP
-
-C Label decay products and choose masses
-
- LHEP=NHEP
-
- MHEP=LHEP+1
-
- NTRY=0
-
- 35 NTRY=NTRY+1
-
- SDKM=PHEP(5,NHEP)
-
- DO 40 I=1,NPR
-
- NHEP=NHEP+1
-
- IDHW(NHEP)=IDKPRD(I,IM)
-
- IDHEP(NHEP)=IDPDG(IDKPRD(I,IM))
-
- ISTHEP(NHEP)=IST(I)
-
- JMOHEP(1,NHEP)=LHEP
-
- JDAHEP(1,NHEP)=0
-
- PHEP(5,NHEP)=HWUMBW(IDKPRD(I,IM))
-
- 40 SDKM=SDKM-PHEP(5,NHEP)
-
- IF (SDKM.LT.ZERO) THEN
-
- NHEP=NHEP-NPR
-
- IF (NTRY.LE.NETRY) GO TO 35
-
- CALL HWWARN('HWDHOB',1,*45)
-
- 45 IF (MTRY.LE.NETRY) GO TO 15
-
- CALL HWWARN('HWDHOB',101,*999)
-
- ENDIF
-
-C Assign production vertices to decay products
-
- CALL HWUDKL(ID,PHEP(1,IHEP),VHEP(1,MHEP))
-
- CALL HWVSUM(4,VHEP(1,IHEP),VHEP(1,MHEP),VHEP(1,MHEP))
-
- CALL HWVEQU(4,VHEP(1,MHEP),VHEP(1,NHEP))
-
- IF (NPR.EQ.2) THEN
-
-C Two body decay: LHEP -> MHEP + NHEP
-
- PCM=HWUPCM(PHEP(5,IHEP),PHEP(5,MHEP),PHEP(5,NHEP))
-
- CALL HWDTWO(PHEP(1,IHEP),PHEP(1,MHEP),
-
- & PHEP(1,NHEP),PCM,TWO,.FALSE.)
-
- ELSEIF (NPR.EQ.3) THEN
-
-C Three body decay: LHEP -> KHEP + MHEP + NHEP
-
- KHEP=MHEP
-
- MHEP=MHEP+1
-
-C Provisional colour self-connection of KHEP
-
- JMOHEP(2,KHEP)=KHEP
-
- JDAHEP(2,KHEP)=KHEP
-
- IF (NME(IM).EQ.100) THEN
-
-C Generate decay momenta using full (V-A)*(V-A) matrix element
-
- EMMX=PHEP(5,IHEP)-PHEP(5,NHEP)
-
- EMWSQ=RMASS(198)**2
-
- GMWSQ=(RMASS(198)*GAMW)**2
-
- EMLIM=GMWSQ
-
- IF (EMMX.LT.RMASS(198)) EMLIM=EMLIM+(EMWSQ-EMMX**2)**2
-
- 50 CALL HWDTHR(PHEP(1,IHEP),PHEP(1,MHEP),
-
- & PHEP(1,KHEP),PHEP(1,NHEP),HWDWWT)
-
- CALL HWVSUM(4,PHEP(1,KHEP),PHEP(1,MHEP),PW)
-
- PW(5)=HWULDO(PW,PW)
-
- EMTST=(EMWSQ-PW(5))**2
-
- IF ((EMTST+GMWSQ)*HWR().GT.EMLIM) GOTO 50
-
- PW(5)=SQRT(PW(5))
-
-C Assign production vertices to 1 and 2
-
- CALL HWUDKL(198,PW,VHEP(1,KHEP))
-
- CALL HWVSUM(4,VHEP(1,NHEP),VHEP(1,KHEP),VHEP(1,KHEP))
-
- ELSEIF(NME(IM).EQ.300) THEN
-
-C Generate momenta using 3-body RPV matrix element
-
- CALL HWDRME(LHEP,KHEP)
-
- ELSE
-
-C Three body phase space decay
-
- CALL HWDTHR(PHEP(1,IHEP),PHEP(1,MHEP),
-
- & PHEP(1,KHEP),PHEP(1,NHEP),HWDPWT)
-
- ENDIF
-
- CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,MHEP))
-
- ELSEIF(NPR.EQ.4) THEN
-
-C Four body decay: LHEP -> KHEP + RHEP + MHEP + NHEP
-
- KHEP = MHEP
-
- RHEP = MHEP+1
-
- MHEP = MHEP+2
-
-C Provisional colour connections of KHEP and RHEP
-
- JMOHEP(2,KHEP)=RHEP
-
- JDAHEP(2,KHEP)=RHEP
-
- JMOHEP(2,RHEP)=KHEP
-
- JDAHEP(2,RHEP)=KHEP
-
-C Four body phase space decay
-
- CALL HWDFOR(PHEP(1,IHEP),PHEP(1,KHEP),PHEP(1,RHEP),
-
- & PHEP(1,MHEP),PHEP(1,NHEP))
-
- CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,RHEP))
-
- CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,MHEP))
-
- ELSE
-
- CALL HWWARN('HWDHOB',102,*999)
-
- ENDIF
-
-C Colour connections
-
- IF (ID.EQ.6.OR.ID.EQ.12.OR.(ID.GE.209.AND.ID.LE.212)
-
- & .OR.(ID.GE.215.AND.ID.LE.218)) THEN
-
- IF (NPR.EQ.3.AND.NME(IM).EQ.100) THEN
-
-C usual heavy quark decay
-
- JMOHEP(2,KHEP)=MHEP
-
- JDAHEP(2,KHEP)=MHEP
-
- JMOHEP(2,MHEP)=KHEP
-
- JDAHEP(2,MHEP)=KHEP
-
- JMOHEP(2,NHEP)=LHEP
-
- JDAHEP(2,NHEP)=LHEP
-
- ELSEIF (ABS(IDHEP(MHEP)).EQ.37) THEN
-
-C heavy quark to charged Higgs
-
- JMOHEP(2,MHEP)=MHEP
-
- JDAHEP(2,MHEP)=MHEP
-
- JMOHEP(2,NHEP)=LHEP
-
- JDAHEP(2,NHEP)=LHEP
-
- ELSEIF (ABS(IDHEP(NHEP)).EQ.37) THEN
-
- JMOHEP(2,MHEP)=LHEP
-
- JDAHEP(2,MHEP)=LHEP
-
- JMOHEP(2,NHEP)=NHEP
-
- JDAHEP(2,NHEP)=NHEP
-
- ELSE
-
- CALL HWWARN('HWDHOB',103,*999)
-
- ENDIF
-
- ELSE
-
- IF(.NOT.RPARTY.AND.
-
- & ((NPR.EQ.2.AND.ID.GE.401.AND.ID.LT.448.AND.
-
- & IDHW(MHEP).LE.132.AND.IDHW(NHEP).LE.132)
-
- & .OR.(NPR.EQ.3.AND.ID.GE.449.AND.ID.LE.457.AND.
-
- & IDHW(MHEP).LE.132.AND.IDHW(NHEP).LE.132.AND.
-
- & IDHW(MHEP-1).LE.132))) THEN
-
-C R-parity violating SUSY decays
-
- IF(NPR.EQ.2) THEN
-
-C--Rparity slepton colour connections
-
- IF(ID.GE.425.AND.ID.LE.448) THEN
-
- IF(IDHW(MHEP).GT.12) THEN
-
- JMOHEP(2,MHEP) = MHEP
-
- JDAHEP(2,MHEP) = MHEP
-
- JMOHEP(2,NHEP) = NHEP
-
- JDAHEP(2,NHEP) = NHEP
-
- ELSE
-
- JMOHEP(2,MHEP) = NHEP
-
- JDAHEP(2,MHEP) = NHEP
-
- JMOHEP(2,NHEP) = MHEP
-
- JDAHEP(2,NHEP) = MHEP
-
- ENDIF
-
-C--Rparity squark colour connections
-
- ELSE
-
- IF(IDHEP(LHEP).GT.0) THEN
-
-C--LQD decay colour connections
-
- IF(IDHW(MHEP).GT.12) THEN
-
- JMOHEP(2,MHEP) = MHEP
-
- JDAHEP(2,MHEP) = MHEP
-
- JMOHEP(2,NHEP) = LHEP
-
- JDAHEP(2,NHEP) = LHEP
-
- ELSE
-
-C--UDD decay colour connections
-
- HVFCEN = .TRUE.
-
- CALL HWDRCL(LHEP,MHEP,CLSAVE)
-
- ENDIF
-
- ELSE
-
-C--Antisquark connections
-
- IF(IDHW(MHEP).GT.12) THEN
-
- JMOHEP(2,MHEP) = MHEP
-
- JDAHEP(2,MHEP) = MHEP
-
- JMOHEP(2,NHEP) = LHEP
-
- JDAHEP(2,NHEP) = LHEP
-
- ELSE
-
- HVFCEN = .TRUE.
-
- CALL HWDRCL(LHEP,MHEP,CLSAVE)
-
- ENDIF
-
- ENDIF
-
- ENDIF
-
- ELSE
-
- IF(ID.GE.450.AND.ID.LE.457) THEN
-
-C--Rparity Neutralino/Chargino colour connection
-
- IF(IDHW(MHEP-1).LE.12.AND.IDHW(MHEP).LE.12.
-
- & AND.IDHW(NHEP).LE.12) THEN
-
- HVFCEN = .TRUE.
-
- CALL HWDRCL(LHEP,MHEP,CLSAVE)
-
- ELSE
-
- JMOHEP(2,MHEP) = NHEP
-
- JDAHEP(2,MHEP) = NHEP
-
- JMOHEP(2,NHEP) = MHEP
-
- JDAHEP(2,NHEP) = MHEP
-
- ENDIF
-
-C--Rparity gluino colour connections
-
- ELSEIF(ID.EQ.449) THEN
-
- IF(IDHW(MHEP-1).LE.12.AND.IDHW(MHEP).LE.12.
-
- & AND.IDHW(NHEP).LE.12) THEN
-
- HVFCEN = .TRUE.
-
- CALL HWDRCL(LHEP,MHEP,CLSAVE)
-
-C--Now the lepton number violating decay
-
- ELSE
-
- IF(IDHW(MHEP).LE.6) THEN
-
- JMOHEP(2,MHEP) = LHEP
-
- JDAHEP(2,MHEP) = NHEP
-
- JMOHEP(2,NHEP) = MHEP
-
- JDAHEP(2,NHEP) = LHEP
-
- ELSE
-
- JMOHEP(2,MHEP) = NHEP
-
- JDAHEP(2,MHEP) = LHEP
-
- JMOHEP(2,NHEP) = LHEP
-
- JDAHEP(2,NHEP) = MHEP
-
- ENDIF
-
- ENDIF
-
- ELSE
-
- CALL HWWARN('HWDHOB',104,*999)
-
- ENDIF
-
- ENDIF
-
- ELSE
-
-C Normal SUSY decays
-
- IF (ID.LE.448.AND.ID.GT.207) THEN
-
-C Squark (or slepton)
-
- IF (IDHW(MHEP).EQ.449) THEN
-
- IF (IDHEP(LHEP).GT.0) THEN
-
- JMOHEP(2,MHEP)=LHEP
-
- JDAHEP(2,MHEP)=NHEP
-
- JMOHEP(2,NHEP)=MHEP
-
- JDAHEP(2,NHEP)=LHEP
-
- ELSE
-
- JMOHEP(2,MHEP)=NHEP
-
- JDAHEP(2,MHEP)=LHEP
-
- JMOHEP(2,NHEP)=LHEP
-
- JDAHEP(2,NHEP)=MHEP
-
- ENDIF
-
- ELSE
-
- IF(NPR.EQ.3.AND.IDHW(MHEP).LE.12) THEN
-
- JMOHEP(2,MHEP)=NHEP
-
- JDAHEP(2,MHEP)=NHEP
-
- JMOHEP(2,NHEP)=MHEP
-
- JDAHEP(2,NHEP)=MHEP
-
- ELSE
-
- JMOHEP(2,MHEP)=MHEP
-
- JDAHEP(2,MHEP)=MHEP
-
- JMOHEP(2,NHEP)=LHEP
-
- JDAHEP(2,NHEP)=LHEP
-
- ENDIF
-
- ENDIF
-
- ELSEIF (ID.EQ.449) THEN
-
-C Gluino
-
- IF (IDHW(NHEP).EQ.13) THEN
-
- JMOHEP(2,MHEP)=MHEP
-
- JDAHEP(2,MHEP)=MHEP
-
- JMOHEP(2,NHEP)=LHEP
-
- JDAHEP(2,NHEP)=LHEP
-
- ELSEIF (IDHEP(MHEP).GT.0) THEN
-
- JMOHEP(2,MHEP)=LHEP
-
- JDAHEP(2,MHEP)=NHEP
-
- JMOHEP(2,NHEP)=MHEP
-
- JDAHEP(2,NHEP)=LHEP
-
- ELSE
-
- JMOHEP(2,MHEP)=NHEP
-
- JDAHEP(2,MHEP)=LHEP
-
- JMOHEP(2,NHEP)=LHEP
-
- JDAHEP(2,NHEP)=MHEP
-
- ENDIF
-
- ELSE
-
-C Gaugino or Higgs
-
- JMOHEP(2,MHEP)=NHEP
-
- JDAHEP(2,MHEP)=NHEP
-
- JMOHEP(2,NHEP)=MHEP
-
- JDAHEP(2,NHEP)=MHEP
-
- ENDIF
-
- ENDIF
-
- ENDIF
-
-C---SPECIAL CASE FOR THREE-BODY TOP DECAYS:
-
-C RELABEL THEM AS TWO TWO-BODY DECAYS FOR PARTON SHOWERING
-
- IF ((ID.EQ.6.OR.ID.EQ.12).AND.NPR.EQ.3.AND.NME(IM).EQ.100) THEN
-
-C---STORE W DECAY PRODUCTS
-
- CALL HWVEQU(10,PHEP(1,KHEP),PDW)
-
-C---BOOST THEM INTO W REST FRAME
-
- CALL HWULOF(PW,PDW(1,1),PDW(1,3))
-
-C---REPLACE THEM BY W
-
- CALL HWVEQU(5,PW,PHEP(1,KHEP))
-
- WHEP=KHEP
-
- IDHW(KHEP)=198
-
- IF (ID.EQ.12) IDHW(KHEP)=199
-
- IDHEP(KHEP)=IDPDG(IDHW(KHEP))
-
- JMOHEP(2,KHEP)=KHEP
-
- JDAHEP(2,KHEP)=KHEP
-
- CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,KHEP))
-
-C---AND MOVE B UP
-
- CALL HWVEQU(5,PHEP(1,NHEP),PHEP(1,MHEP))
-
- IDHW(MHEP)=IDHW(NHEP)
-
- IDHEP(MHEP)=IDHEP(NHEP)
-
- JDAHEP(2,LHEP)=MHEP
-
- JMOHEP(2,MHEP)=JMOHEP(2,NHEP)
-
- JDAHEP(2,MHEP)=JDAHEP(2,NHEP)
-
- CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,MHEP))
-
- NHEP=MHEP
-
-C---DO PARTON SHOWER
-
- EMSCA=PHEP(5,IHEP)
-
- CALL HWBGEN
-
- IF (IERROR.NE.0) RETURN
-
-C---FIND BOOSTED W MOMENTUM
-
- NTRY=0
-
- 41 NTRY=NTRY+1
-
- IF (NTRY.GT.NHEP.OR.WHEP.LE.0.OR.WHEP.GT.NHEP)
-
- $ CALL HWWARN('HWDHOB',101,*999)
-
- WHEP=JDAHEP(1,WHEP)
-
- IF (ISTHEP(WHEP).NE.190) GOTO 41
-
-C---AND HENCE ITS CHILDRENS MOMENTA
-
- CALL HWULOB(PHEP(1,WHEP),PDW(1,3),PHEP(1,NHEP+1))
-
- CALL HWVDIF(4,PHEP(1,WHEP),PHEP(1,NHEP+1),PHEP(1,NHEP+2))
-
- PHEP(5,NHEP+2)=PDW(5,2)
-
-C---LABEL THEM
-
- ISTHEP(WHEP)=195
-
- DO 51 I=1,2
-
- IDHW(NHEP+I)=IDKPRD(I,IM)
-
- IDHEP(NHEP+I)=IDPDG(IDHW(NHEP+I))
-
- ISTHEP(NHEP+I)=112+I
-
- JDAHEP(I,WHEP)=NHEP+I
-
- JMOHEP(1,NHEP+I)=WHEP
-
- JMOHEP(2,NHEP+I)=NHEP+3-I
-
- JDAHEP(2,NHEP+I)=NHEP+3-I
-
- 51 CONTINUE
-
- NHEP=NHEP+2
-
-C---ASSIGN PRODUCTION VERTICES TO 1 AND 2
-
- CALL HWUDKL(198,PW,VHEP(1,NHEP))
-
- CALL HWVSUM(4,VHEP(1,WHEP),VHEP(1,NHEP),VHEP(1,NHEP))
-
- CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-1))
-
-C---DO PARTON SHOWERS
-
- EMSCA=PW(5)
-
- CALL HWBGEN
-
- IF (IERROR.NE.0) RETURN
-
- ELSE
-
-C Do parton showers
-
- EMSCA=PHEP(5,IHEP)
-
- CALL HWBGEN
-
- IF (IERROR.NE.0) RETURN
-
- ENDIF
-
- ENDIF
-
-C--New to correct colour connections in Rslash
-
- IF(CLSAVE(1).NE.0) THEN
-
- THEP = MHEP+1
-
- ID = IDHW(CLSAVE(1))
-
- IDM = IDHW(JMOHEP(1,CLSAVE(1)))
-
- IDM2 = IDHW(LHEP)
-
- IF(IDM.EQ.15) ID=IDHW(JMOHEP(1,JMOHEP(1,CLSAVE(1))))
-
- IF((ID.LE.6.AND.((IDM.GE.419.AND.IDM.LE.424).OR.IDM.EQ.411.OR.
-
- & IDM.EQ.412).
-
- & AND.((IDM2.GE.413.AND.IDM2.LE.418)
-
- & .OR.IDM2.EQ.449).OR.IDM2.EQ.405.OR.IDM2.EQ.406)
-
- & .OR.(ID.LE.6.AND.IDM.EQ.449.AND.
-
- & (((IDM2.GE.413.AND.IDM2.LE.418).OR.IDM2.EQ.405.OR.IDM2.EQ.406)
-
- & .OR.IDM2.EQ.449)).OR.
-
- & (IDM.EQ.15.AND.ID.LE.12.AND.ID.GE.7.AND.((IDM2.GE.413.AND.
-
- & IDM2.LE.418).OR.IDM2.EQ.449.OR.IDM2.
-
- & EQ.405.OR.IDM2.EQ.406))) THEN
-
- IF(JMOHEP(2,CLSAVE(1)).EQ.MHEP) THEN
-
- IF(IDHW(CLSAVE(1)).NE.13.AND.IDHW(CLSAVE(1)).NE.449)
-
- & JMOHEP(2,CLSAVE(2)) = THEP
-
- JDAHEP(2,MHEP) = CLSAVE(1)
-
- JDAHEP(2,THEP) = CLSAVE(2)
-
- ELSE
-
- IF(IDHW(CLSAVE(1)).NE.13.AND.IDHW(CLSAVE(1)).NE.449)
-
- & JMOHEP(2,CLSAVE(2)) = MHEP
-
- JDAHEP(2,MHEP) = CLSAVE(2)
-
- JDAHEP(2,THEP) = CLSAVE(1)
-
- ENDIF
-
- ELSEIF((ID.GT.6.AND.ID.LE.12.
-
- & AND.((IDM.GE.413.AND.IDM.LE.418).OR.IDM.EQ.405.OR.
-
- & IDM.EQ.406).AND.
-
- & ((IDM2.GE.419.AND.IDM2.LE.424).OR.IDM2.EQ.449.OR.
-
- & IDM2.EQ.411.OR.IDM2.EQ.412)).OR.
-
- & (ID.GT.6.AND.ID.LE.12.AND.IDM.EQ.449.
-
- & AND.((IDM2.GE.419.AND.IDM2.LE.424).OR.IDM2.EQ.449.OR.
-
- & IDM2.EQ.411.OR.IDM2.EQ.412)).OR.
-
- & (IDM.EQ.15.AND.ID.LE.6.AND.((IDM2.GE.419.AND.
-
- & IDM2.LE.424).OR.IDM2.EQ.449.OR.IDM2.EQ.411.OR.
-
- & IDM2.EQ.412))) THEN
-
- IF(JDAHEP(2,CLSAVE(1)).EQ.MHEP) THEN
-
- JDAHEP(2,CLSAVE(2))=THEP
-
- JMOHEP(2,MHEP)=CLSAVE(1)
-
- JMOHEP(2,THEP)=CLSAVE(2)
-
- ELSE
-
- JDAHEP(2,CLSAVE(2))=MHEP
-
- JMOHEP(2,MHEP)=CLSAVE(2)
-
- JMOHEP(2,THEP)=CLSAVE(1)
-
- ENDIF
-
- ENDIF
-
- COLUPD = .FALSE.
-
- CALL HWBCON
-
- ENDIF
-
- IF (IHEP.EQ.NHEP) GOTO 70
-
- 60 CONTINUE
-
- 70 IF (FOUND) THEN
-
-C Fix any SUSY colour disconnections
-
- DO 80 IHEP=1,NHEP
-
- IF (ISTHEP(IHEP).GE.147.AND.ISTHEP(IHEP).LE.151
-
- & .AND.JDAHEP(2,IHEP).EQ.0) THEN
-
- IM=JMOHEP(1,IHEP)
-
-C Chase connection back through SUSY decays
-
- 75 IM=JMOHEP(1,IM)
-
- ISM=ISTHEP(IM)
-
- IF (ISM.EQ.120) GOTO 80
-
- IF (ISM.NE.123.AND.ISM.NE.124.AND.ISM.NE.155) GOTO 75
-
-C Look for unclustered parton to connect
-
- DO JHEP=1,NHEP
-
- IF (ISTHEP(JHEP).GE.147.AND.ISTHEP(JHEP).LE.151) THEN
-
- JCM=JMOHEP(2,JHEP)
-
- IF (JCM.EQ.IM) THEN
-
-C Found it: connect
-
- JMOHEP(2,JHEP)=IHEP
-
- JDAHEP(2,IHEP)=JHEP
-
- GOTO 80
-
- ENDIF
-
- ENDIF
-
- ENDDO
-
-C Not found: need to go further back
-
- GOTO 75
-
- ENDIF
-
- 80 CONTINUE
-
-C Go back to check for further heavy decay products
-
- GOTO 10
-
- ENDIF
-
- 999 END
-
-CDECK ID>, HWDHVY.
-
-*CMZ :- -26/04/91 12.19.24 by Federico Carminati
-
-*-- Author : Ian Knowles & Bryan Webber
-
-C-----------------------------------------------------------------------
-
- SUBROUTINE HWDHVY
-
-C-----------------------------------------------------------------------
-
-C Performs partonic decays of hadrons containing heavy quark(s):
-
-C either, meson/baryon spectator model weak decays;
-
-C or, quarkonia -> 2-gluons, q-qbar, 3-gluons, or 2-gluons + photon.
-
-C-----------------------------------------------------------------------
-
- INCLUDE 'HERWIG61.INC'
-
- DOUBLE PRECISION HWULDO,HWR,XS,XB,EMWSQ,GMWSQ,EMLIM,PW(4),
-
- & EMTST,X1,X2,X3,TEST,HWDWWT,HWDPWT
-
- INTEGER IST(3),I,IHEP,IM,ID,IDQ,IQ,IS,J
-
- EXTERNAL HWR,HWDWWT,HWDPWT,HWULDO
-
- DATA IST/113,114,114/
-
- IF (IERROR.NE.0) RETURN
-
- DO 100 I=1,NMXQDK
-
- IF (I.GT.NQDK) THEN
-
- NQDK=0
-
- RETURN
-
- ENDIF
-
- IHEP=LOCQ(I)
-
- IF (ISTHEP(IHEP).EQ.199) GOTO 100
-
- IM=IMQDK(I)
-
- IF (NHEP+NPRODS(IM).GT.NMXHEP) CALL HWWARN('HWDHVY',100,*999)
-
- IF (IDKPRD(4,IM).NE.0) THEN
-
-C Weak decay of meson or baryon
-
-C Idenitify decaying heavy quark and spectator
-
- ID=IDHW(IHEP)
-
- IF (ID.EQ.136.OR.ID.EQ.140.OR.ID.EQ.144.OR.
-
- & ID.EQ.150.OR.ID.EQ.155.OR.ID.EQ.158.OR.ID.EQ.161.OR.
-
- & (ID.EQ.254.AND.IDKPRD(4,IM).EQ.11)) THEN
-
-C c hadron or c decay of B_c+
-
- IDQ=4
-
- IQ=NHEP+1
-
- IS=NHEP+2
-
- ELSEIF (ID.EQ.171.OR.ID.EQ.175.OR.ID.EQ.179.OR.
-
- & ID.EQ.185.OR.ID.EQ.190.OR.ID.EQ.194.OR.ID.EQ.196.OR.
-
- & (ID.EQ.230.AND.IDKPRD(4,IM).EQ.5)) THEN
-
-C cbar hadron or cbar decay of B_c-
-
- IDQ=10
-
- IS=NHEP+1
-
- IQ=NHEP+2
-
- ELSEIF ((ID.GE.221.AND.ID.LE.229).OR.
-
- & (ID.EQ.230.AND.IDKPRD(4,IM).EQ.10)) THEN
-
-C b hadron or b decay of B_c-
-
- IDQ=5
-
- IQ=NHEP+1
-
- IS=NHEP+2
-
- ELSEIF ((ID.GE.245.AND.ID.LE.253).OR.
-
- & (ID.EQ.254.AND.IDKPRD(4,IM).EQ.4)) THEN
-
-C bbar hadron or bbar decay of B_c+
-
- IDQ=11
-
- IS=NHEP+1
-
- IQ=NHEP+2
-
- ELSE
-
-C Decay not recognized
-
- CALL HWWARN('HWDHVY',101,*999)
-
- ENDIF
-
-C Label constituents
-
- IF (NHEP+5.GT.NMXHEP) CALL HWWARN('HWDHVY',102,*999)
-
- ISTHEP(IHEP)=199
-
- JDAHEP(1,IHEP)=NHEP+1
-
- JDAHEP(2,IHEP)=NHEP+2
-
- IDHW(IQ)=IDQ
-
- IDHW(IS)=IDKPRD(4,IM)
-
- IDHEP(IQ)=IDPDG(IDQ)
-
- IDHEP(IS)=IDPDG(IDKPRD(4,IM))
-
- ISTHEP(IQ)=155
-
- ISTHEP(IS)=115
-
- JMOHEP(1,IQ)=IHEP
-
- JMOHEP(2,IQ)=IS
-
- JDAHEP(1,IQ)=NHEP+3
-
- JDAHEP(2,IQ)=NHEP+5
-
- JMOHEP(1,IS)=IHEP
-
- JMOHEP(2,IS)=NHEP+5
-
- JDAHEP(1,IS)=0
-
- JDAHEP(2,IS)=NHEP+5
-
- NHEP=NHEP+2
-
-C and weak decay product jets
-
- DO 10 J=1,3
-
- NHEP=NHEP+1
-
- IDHW(NHEP)=IDKPRD(J,IM)
-
- IDHEP(NHEP)=IDPDG(IDKPRD(J,IM))
-
- ISTHEP(NHEP)=IST(J)
-
- JMOHEP(1,NHEP)=IQ
-
- JDAHEP(1,NHEP)=0
-
- 10 PHEP(5,NHEP)=RMASS(IDKPRD(J,IM))
-
- JMOHEP(2,NHEP-2)=NHEP-1
-
- JDAHEP(2,NHEP-2)=NHEP-1
-
- JMOHEP(2,NHEP-1)=NHEP-2
-
- JDAHEP(2,NHEP-1)=NHEP-2
-
- JMOHEP(2,NHEP )=IQ
-
- JDAHEP(2,NHEP )=IQ
-
-C Share momenta in ratio of masses, preserving specator mass
-
- XS=RMASS(IDHW(IS))/PHEP(5,IHEP)
-
- XB=ONE-XS
-
- CALL HWVSCA(5,XB,PHEP(1,IHEP),PHEP(1,IQ))
-
- CALL HWVSCA(5,XS,PHEP(1,IHEP),PHEP(1,IS))
-
- IF (NME(IM).EQ.100) THEN
-
-C Generate decay momenta using full (V-A)*(V-A) matrix element
-
- EMWSQ=RMASS(198)**2
-
- GMWSQ=(RMASS(198)*GAMW)**2
-
- EMLIM=GMWSQ+(EMWSQ-(PHEP(5,IQ)-PHEP(5,NHEP))**2)**2
-
- 20 CALL HWDTHR(PHEP(1,IQ ),PHEP(1,NHEP-1),
-
- & PHEP(1,NHEP-2),PHEP(1,NHEP),HWDWWT)
-
- CALL HWVSUM(4,PHEP(1,NHEP-2),PHEP(1,NHEP-1),PW)
-
- EMTST=(HWULDO(PW,PW)-EMWSQ)**2
-
- IF ((EMTST+GMWSQ)*HWR().GT.EMLIM) GOTO 20
-
- ELSE
-
-C Use phase space
-
- CALL HWDTHR(PHEP(1,IQ ),PHEP(1,NHEP-2),
-
- & PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
-
- CALL HWVSUM(4,PHEP(1,NHEP-2),PHEP(1,NHEP-1),PW)
-
- ENDIF
-
-C Set up production vertices
-
- CALL HWVZRO(4,VHEP(1,IQ))
-
- CALL HWVEQU(4,VHEP(1,IQ),VHEP(1,IS))
-
- CALL HWVEQU(4,VHEP(1,IQ),VHEP(1,NHEP))
-
- CALL HWUDKL(198,PW,VHEP(1,NHEP-2))
-
- CALL HWVSUM(4,VHEP(1,IQ),VHEP(1,NHEP-2),VHEP(1,NHEP-2))
-
- CALL HWVEQU(4,VHEP(1,NHEP-2),VHEP(1,NHEP-1))
-
- EMSCA=PHEP(5,IQ)
-
- ELSE
-
-C Quarkonium decay
-
-C Label products
-
- ISTHEP(IHEP)=199
-
- JDAHEP(1,IHEP)=NHEP+1
-
- DO 30 J=1,NPRODS(IM)
-
- NHEP=NHEP+1
-
- IDHW(NHEP)=IDKPRD(J,IM)
-
- IDHEP(NHEP)=IDPDG(IDKPRD(J,IM))
-
- ISTHEP(NHEP)=IST(J)
-
- JMOHEP(1,NHEP)=IHEP
-
- JDAHEP(1,NHEP)=0
-
- PHEP(5,NHEP)=RMASS(IDKPRD(J,IM))
-
- 30 CALL HWVZRO(4,VHEP(1,NHEP))
-
- JDAHEP(2,IHEP)=NHEP
-
-C Establish colour connections and select momentum configuration
-
- IF (NPRODS(IM).EQ.3) THEN
-
- IF (IDKPRD(3,IM).EQ.13) THEN
-
-C 3-gluon decay
-
- JMOHEP(2,NHEP-2)=NHEP
-
- JMOHEP(2,NHEP-1)=NHEP-2
-
- JMOHEP(2,NHEP )=NHEP-1
-
- JDAHEP(2,NHEP-2)=NHEP-1
-
- JDAHEP(2,NHEP-1)=NHEP
-
- JDAHEP(2,NHEP )=NHEP-2
-
- ELSE
-
-C or 2-gluon + photon decay
-
- JMOHEP(2,NHEP-2)=NHEP-1
-
- JMOHEP(2,NHEP-1)=NHEP-2
-
- JMOHEP(2,NHEP )=NHEP
-
- JDAHEP(2,NHEP-2)=NHEP-1
-
- JDAHEP(2,NHEP-1)=NHEP-2
-
- JDAHEP(2,NHEP )=NHEP
-
- ENDIF
-
- IF (NME(IM).EQ.130) THEN
-
-C Use Ore & Powell orthopositronium matrix element
-
- 40 CALL HWDTHR(PHEP(1,IHEP),PHEP(1,NHEP-2),
-
- & PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
-
- X1=TWO*HWULDO(PHEP(1,IHEP),PHEP(1,NHEP-2))/PHEP(5,IHEP)**2
-
- X2=TWO*HWULDO(PHEP(1,IHEP),PHEP(1,NHEP-1))/PHEP(5,IHEP)**2
-
- X3=TWO-X1-X2
-
- TEST=((X1*(ONE-X1))**2+(X2*(ONE-X2))**2+(X3*(ONE-X3))**2)
-
- & /(X1*X2*X3)**2
-
- IF (TEST.LT.TWO*HWR()) GOTO 40
-
- ELSE
-
-C Use phase space
-
- CALL HWDTHR(PHEP(1,IHEP),PHEP(1,NHEP-2),
-
- & PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
-
- ENDIF
-
- ELSE
-
-C Parapositronium 2-gluon or q-qbar decay
-
- JMOHEP(2,NHEP-1)=NHEP
-
- JMOHEP(2,NHEP )=NHEP-1
-
- JDAHEP(2,NHEP-1)=NHEP
-
- JDAHEP(2,NHEP )=NHEP-1
-
- CALL HWDTWO(PHEP(1,IHEP),PHEP(1,NHEP-1),
-
- & PHEP(1,NHEP),CMMOM(IM),TWO,.FALSE.)
-
- ENDIF
-
- EMSCA=PHEP(5,IHEP)
-
- ENDIF
-
-C Process this new hard scatter
-
- CALL HWVEQU(4,VTXQDK(1,I),VTXPIP)
-
- CALL HWBGEN
-
- CALL HWCFOR
-
- CALL HWCDEC
-
- CALL HWDHAD
-
- 100 CONTINUE
-
- NQDK=0
-
- 999 END
-
-CDECK ID>, HWDRCL.
-
-*CMZ :- -20/07/99 10:56:12 by Peter Richardson
-
-*-- Author : Peter Richardson
-
-C-----------------------------------------------------------------------
-
- SUBROUTINE HWDRCL(IHEP,MHEP,CLSAVE)
-
-C-----------------------------------------------------------------------
-
-C Sets the colour connections in Baryon number violating decays
-
-C-----------------------------------------------------------------------
-
- INCLUDE 'HERWIG61.INC'
-
- INTEGER IHEP,MHEP,ID,ID2,IDM2,IDM3,COLCON(2,2,3),FLACON(2,3),JHEP,
-
- & DECAY,COLANT,KHEP,IDM,IDMB,IDMB2,IDMB3,IDMB4,QHEP,IDM4,
-
- & CLSAVE(2),XHEP,I,HWRINT,THEP
-
- LOGICAL CONBV
-
-C--Colour connections for the decays
-
- DATA COLCON/-1,1,-1,-2,-2,1,-3,-1,-1,1,-2,-1/
-
- DATA FLACON/1,-1,1,-1,-1,0/
-
-C--identify the decay
-
- IF(IERROR.NE.0) RETURN
-
- ID = IDHW(IHEP)
-
- ID2 = IDHW(MHEP)
-
- IF(ID.GE.450.AND.ID.LE.457) THEN
-
- DECAY = 1
-
- ELSEIF(ID.EQ.449) THEN
-
- DECAY = 2
-
- ELSEIF((ID.GE.411.AND.ID.LE.424).OR.ID.EQ.405.OR.ID.EQ.406) THEN
-
- DECAY = 3
-
- ELSE
-
-C--UNKNOWN DECAY
-
- CALL HWWARN('HWDRCL',100,*999)
-
- ENDIF
-
- COLANT = 1
-
-C--identify the colour partner
-
- IF(DECAY.GT.1.AND.ID2.LE.6) THEN
-
-C--colour partner
-
- COLANT = 2
-
- KHEP = JDAHEP(2,IHEP-1)
-
- ELSEIF(DECAY.GT.1.AND.ID2.GE.7) THEN
-
-C--anticolour partner
-
- COLANT = 3
-
- KHEP = JMOHEP(2,IHEP)
-
- ELSE
-
- KHEP=IHEP
-
- ENDIF
-
- IDM = IDHW(JMOHEP(1,KHEP))
-
- IF(ABS(IDPDG(IDM)).GT.1000000.OR.IDM.EQ.15) THEN
-
- IDM2 = IDHW(JDAHEP(1,JMOHEP(1,KHEP)))
-
- IDM3 = IDHW(JDAHEP(2,JMOHEP(1,KHEP)))
-
- IDM4 = IDHW(JDAHEP(2,JMOHEP(1,KHEP))-1)
-
- QHEP = JMOHEP(1,KHEP)
-
- IDMB = IDHW(JMOHEP(1,QHEP))
-
- IDMB2 = IDHW(JMOHEP(2,QHEP))
-
- IDMB3 = IDHW(JDAHEP(1,QHEP))
-
- IDMB4 = IDHW(JDAHEP(2,QHEP))
-
- ENDIF
-
-C--Now decide if the colour partner decayed via BV
-
- IF(COLANT.EQ.2.AND.((((IDM.GE.413.AND.IDM.LE.418).OR.
-
- & IDM.EQ.449.OR.IDM.EQ.405.OR.IDM.EQ.406).AND.
-
- & (IDM2.GE.7.AND.IDM2.LE.12.AND.
-
- & IDM3.GE.7.AND.IDM3.LE.12.AND.
-
- & IDM4.GE.7.AND.IDM4.LE.12)).OR.
-
- & (IDM.EQ.15.AND.IDMB.LE.6.AND.IDMB2.LE.6.AND.
-
- & ((IDMB3.GE.7.AND.IDMB4.GE.12.AND.IDMB4.EQ.449).OR.
-
- & (IDMB3.GE.198.AND.IDMB3.LE.207.AND.
-
- & ABS(IDPDG(IDMB4)).GT.1000000))))) THEN
-
- CONBV = .TRUE.
-
- COLUPD = .TRUE.
-
- HVFCEN = .FALSE.
-
- XHEP = JMOHEP(2,JDAHEP(2,JMOHEP(1,KHEP)))
-
- ELSEIF(COLANT.EQ.3.AND.((((IDM.GE.419.AND.IDM.LE.424).OR.
-
- & IDM.EQ.449.OR.IDM.EQ.411.OR.IDM.EQ.412).AND.
-
- & (IDM2.LE.6.AND.IDM3.LE.6.AND.IDM4.LE.6)).OR.
-
- & (IDM.EQ.15.AND.IDMB.GE.7.AND.IDMB.LE.12.AND.
-
- & IDMB2.GE.7.AND.IDMB2.LE.12.AND.((IDMB3.LE.6.AND.
-
- & IDMB4.EQ.449).OR.(ABS(IDPDG(IDMB4)).GT.1000000
-
- & .AND.IDMB3.GE.198.AND.IDMB3.LE.207))))) THEN
-
- CONBV = .TRUE.
-
- COLUPD = .TRUE.
-
- HVFCEN = .FALSE.
-
- XHEP = JDAHEP(2,JDAHEP(2,JMOHEP(1,KHEP)))
-
- ELSE
-
- CONBV = .FALSE.
-
- COLUPD = .FALSE.
-
- XHEP = 0
-
- ENDIF
-
- IF(CONBV) THEN
-
- IF(IDM.NE.15) THEN
-
- CLSAVE(1) = JDAHEP(2,JMOHEP(1,KHEP))-1
-
- CLSAVE(2) = CLSAVE(1)+1
-
- ELSE
-
- IF(IDMB4.EQ.449) THEN
-
- DO I=1,2
-
- CLSAVE(I) = JMOHEP(I,JMOHEP(1,KHEP))
-
- IF(CLSAVE(I).EQ.XHEP) CLSAVE(I)=JDAHEP(1,JMOHEP(1,KHEP))
-
- ENDDO
-
- ELSE
-
- CLSAVE(1) = JMOHEP(1,JMOHEP(1,KHEP))
-
- CLSAVE(2) = JMOHEP(2,JMOHEP(1,KHEP))
-
- ENDIF
-
- ENDIF
-
- ELSE
-
- CLSAVE(1)=0
-
- CLSAVE(2)=0
-
- ENDIF
-
-C--Now set the colours for angular ordering
-
- THEP = MHEP-1
-
- IF(DECAY.EQ.1) THEN
-
- IF(ID2.LE.6) THEN
-
- JMOHEP(2,THEP) = THEP+HWRINT(1,2)
-
- JDAHEP(2,THEP) = THEP
-
- ELSE
-
- JMOHEP(2,THEP) = THEP
-
- JDAHEP(2,THEP) = THEP+HWRINT(1,2)
-
- ENDIF
-
- ELSEIF(DECAY.EQ.2) THEN
-
- IF(ID2.LE.6) THEN
-
- JMOHEP(2,THEP) = IHEP
-
- JDAHEP(2,THEP) = THEP
-
- ELSE
-
- JMOHEP(2,THEP) = THEP
-
- JDAHEP(2,THEP) = IHEP
-
- ENDIF
-
- ENDIF
-
-C--Colour of the second two
-
- DO JHEP=1,2
-
- IF(ID2.LE.6) THEN
-
- JMOHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+
-
- & COLCON(HWRINT(1,2),JHEP,DECAY)
-
- JDAHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+FLACON(JHEP,DECAY)
-
- ELSE
-
- JDAHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+
-
- & COLCON(HWRINT(1,2),JHEP,DECAY)
-
- JMOHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+FLACON(JHEP,DECAY)
-
- ENDIF
-
- ENDDO
-
-C--Now set the colours of the colour partner
-
- IF(DECAY.GT.1.AND..NOT.CONBV) THEN
-
- IF(ID2.LE.6) JMOHEP(2,KHEP) = MHEP+HWRINT(0,1)
-
- IF(ID2.GE.7) JDAHEP(2,KHEP) = MHEP+HWRINT(0,1)
-
- ELSEIF(CONBV) THEN
-
- IF(ID2.GT.6) THEN
-
- JMOHEP(2,CLSAVE(1)) = MHEP+HWRINT(0,1)
-
- IF(JMOHEP(2,CLSAVE(1)).EQ.MHEP) THEN
-
- JMOHEP(2,CLSAVE(2)) = MHEP+1
-
- ELSE
-
- JMOHEP(2,CLSAVE(2)) = MHEP
-
- ENDIF
-
- ELSE
-
- JDAHEP(2,CLSAVE(1)) = MHEP+HWRINT(0,1)
-
- IF(JDAHEP(2,CLSAVE(1)).EQ.MHEP) THEN
-
- JDAHEP(2,CLSAVE(2)) = MHEP+1
-
- ELSE
-
- JDAHEP(2,CLSAVE(2)) = MHEP
-
- ENDIF
-
- ENDIF
-
- ENDIF
-
- 999 END
-
-CDECK ID>, HWDRME.
-
-*CMZ :- -20/07/99 10:56:12 by Peter Richardson
-
-*-- Author : Peter Richardson
-
-C-----------------------------------------------------------------------
-
- SUBROUTINE HWDRME(LHEP,MHEP)
-
-C-----------------------------------------------------------------------
-
-C SUBROUTINE TO IMPLEMENT ALL RPARITY DECAY MATRIX ELEMENTS
-
-C-----------------------------------------------------------------------
-
- INCLUDE 'HERWIG61.INC'
-
- DOUBLE PRECISION SM(6),SW(6),HWULDO,INFCOL,AM, M12SQ,M23SQ,MSGN,
-
- & M13SQ,A(6),B(6),SWEAK,MW,DECMOM(5),TEST(4),EPS,
-
- & M12SQT(6),M23SQT(6),M13SQT(6),LIMIT,M(4),RAND,
-
- & MC(2),MX2(6),MX(6),HWDPWT,HWR,HWDRM1,LAMD(3)
-
- EXTERNAL HWDRM1,HWULDO,HWDPWT,HWR
-
- INTEGER K,SN(3),LHEP,CSP,I,SB(3),J,ND,LTRY,MHEP,NSP,ID(3),IG,
-
- & IDHWTP,IDHPTP,MTRY
-
- PARAMETER(EPS=1D-20)
-
- IF(IERROR.NE.0) RETURN
-
-C--Electroweak parameters, etc
-
- SWEAK = SQRT(SWEIN)
-
- MW = RMASS(198)
-
- M(4) = PHEP(5,LHEP)
-
- IG = IDHW(LHEP)
-
-C--Find the masses of the final state and zero parameters
-
- DO K=1,3
-
- ID(K) = IDHW(MHEP+K-1)
-
- IF(ID(K).LE.12) THEN
-
- SN(K)=ID(K)
-
- ELSE
-
- SN(K)=ID(K)-120
-
- ENDIF
-
- IF(SN(K).GT.6) SN(K)=SN(K)-6
-
- M(K) = PHEP(5,LHEP+K)
-
- SB(K)=SN(K)
-
- LAMD(K) = ZERO
-
- ENDDO
-
- DO J=1,6
-
- MX2(J) = ZERO
-
- MX(J) = ZERO
-
- M13SQT(J) = ZERO
-
- M23SQT(J) = ZERO
-
- M12SQT(J) = ZERO
-
- ENDDO
-
-C--Evaluate the coefficents for the mode we want
-
- IF(IG.GE.450.AND.IG.LE.453) THEN
-
-C--NEUTRALINO
-
- NSP = IG-449
-
- AM = RMASS(IG)
-
- MSGN = ZSGNSS(NSP)
-
- MC(1) = ZMIXSS(NSP,3)/(2*MW*COSB*SWEAK)
-
- MC(2) = ZMIXSS(NSP,4)/(2*MW*SINB*SWEAK)
-
-C--Calculate the combinations of couplings needed
-
- IF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
-
-C--first for the UDD modes
-
- DO J=1,2
-
- A(J) = M(1)*MC(2)*QMIXSS(SN(1),2,J)
-
- & +SLFCH(SN(1),NSP)*QMIXSS(SN(1),1,J)
-
- B(J) = MSGN*(M(1)*MC(2)*QMIXSS(SN(1),1,J)
-
- & +SRFCH(SN(1),NSP)*QMIXSS(SN(1),2,J))
-
- MX2(J) = QMIXSS(SN(1),2,J)
-
- A(J+2) = M(2)*MC(1)*QMIXSS(SN(2),2,J)
-
- & +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
-
- B(J+2) = MSGN*(M(2)*MC(1)*QMIXSS(SN(2),1,J)
-
- & +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
-
- MX2(J+2) = QMIXSS(SN(2),2,J)
-
- A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
-
- & +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
-
- B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
-
- & +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
-
- MX2(J+2) = QMIXSS(SN(3),2,J)
-
- ENDDO
-
- DO K=1,3
-
- SN(K) = SN(K)+400
-
- SB(K) = SB(K)+412
-
- ENDDO
-
- ELSEIF(ID(1).GE.121.AND.ID(2).GE.121.AND.ID(3).GE.121) THEN
-
-C--Now for the LLE modes
-
- DO J=1,2
-
- A(J) = MSGN*(M(1)*MC(1)*LMIXSS(SN(1),1,J)
-
- & +SRFCH(10+SN(1),NSP)*LMIXSS(SN(1),2,J))
-
- B(J) = M(1)*MC(1)*LMIXSS(SN(1),2,J)
-
- & +SLFCH(10+SN(1),NSP)*LMIXSS(SN(2),1,J)
-
- MX2(J)= LMIXSS(SN(1),1,J)
-
- A(J+2) = ZERO
-
- B(J+2) = SLFCH(10+SN(2),NSP)*LMIXSS(SN(2),1,J)
-
- MX2(J+2) = LMIXSS(SN(2),1,J)
-
- A(J+4) = M(3)*MC(1)*LMIXSS(SN(3),2,J)
-
- & +SLFCH(10+SN(3),NSP)*LMIXSS(SN(3),1,J)
-
- B(J+4) = MSGN*(M(3)*MC(1)*LMIXSS(SN(3),1,J)
-
- & +SRFCH(10+SN(3),NSP)*LMIXSS(SN(3),2,J))
-
- MX2(4+J) = LMIXSS(SN(3),2,J)
-
- ENDDO
-
- DO J=1,3
-
- SN(J) = SN(J) + 424
-
- SB(J) = SB(J) + 436
-
- ENDDO
-
- ELSE
-
-C--Now for both types of LQD modes
-
- IF(MOD(SN(1),2).EQ.0) THEN
-
-C--First the neutrino,down,antidown mode
-
- DO J=1,2
-
- A(J) = ZERO
-
- B(J) = SLFCH(10+SN(1),NSP)*
-
- & LMIXSS(SN(1),1,J)
-
- MX2(J) = LMIXSS(SN(1),1,J)
-
- A(J+2) = MSGN*(M(2)*MC(1)*QMIXSS(SN(2),1,J)
-
- & +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
-
- B(J+2) = M(2)*MC(1)*QMIXSS(SN(2),2,J)
-
- & +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
-
- MX2(2+J) = QMIXSS(SN(2),1,J)
-
- A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
-
- & +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
-
- B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
-
- & +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
-
- MX2(J+4) = QMIXSS(SN(3),2,J)
-
- ENDDO
-
- ELSE
-
-C--Now the charged lepton, antiup,down modes
-
- DO J=1,2
-
- A(J) = MSGN*(M(1)*MC(1)*LMIXSS(SN(1),1,J)
-
- & +SRFCH(10+SN(1),NSP)*LMIXSS(SN(1),2,J))
-
- B(J) = M(1)*MC(1)*LMIXSS(SN(1),2,J)
-
- & +SLFCH(10+SN(1),NSP)*LMIXSS(SN(1),1,J)
-
- MX2(J) = LMIXSS(SN(1),1,J)
-
- A(J+2) =MSGN*(M(2)*MC(2)*QMIXSS(SN(2),1,J)
-
- & +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
-
- B(J+2) = M(2)*MC(2)*QMIXSS(SN(2),2,J)
-
- & +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
-
- MX2(2+J) = QMIXSS(SN(2),1,J)
-
- A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
-
- & +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
-
- B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
-
- & +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
-
- MX2(J+4) = QMIXSS(SN(3),2,J)
-
- ENDDO
-
- ENDIF
-
- SN(1) = SN(1) + 424
-
- SB(1) = SB(1) + 436
-
- DO J=2,3
-
- SN(J) = SN(J) + 400
-
- SB(J) = SB(J) + 412
-
- ENDDO
-
- ENDIF
-
- DO K=1,3
-
- SM(2*K-1) = RMASS(SN(K))
-
- SM(2*K) = RMASS(SB(K))
-
- SW(2*K-1) = HBAR/RLTIM(SN(K))
-
- SW(2*K) = HBAR/RLTIM(SB(K))
-
- ENDDO
-
- ND = 3
-
- DO K=1,3
-
- LAMD(K) = ONE
-
- ENDDO
-
- INFCOL = ONE
-
- ELSEIF(IG.EQ.449) THEN
-
-C--GLUINO
-
-C--First obtian the masses and widths needed
-
- AM = RMASS(IG)
-
- ND = 3
-
-C--Calculate the combinations of couplings needed
-
- IF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
-
-C--first for the UDD modes
-
- INFCOL = -0.5D0
-
-C--Couplings
-
- DO I=1,3
-
- DO J=1,2
-
- A(2*I-2+J) = -QMIXSS(SN(I),1,J)
-
- B(2*I-2+J) = QMIXSS(SN(I),2,J)
-
- MX2(2*I-2+J) = QMIXSS(SN(I),2,J)
-
- ENDDO
-
- SN(I) = SN(I)+400
-
- SB(I) = SB(I)+412
-
- ENDDO
-
- ELSE
-
- INFCOL = ONE
-
-C--Now for both types of LQD modes
-
- IF(MOD(SN(1),2).EQ.0) THEN
-
-C--First the neutrino,down,antidown mode
-
- DO J=1,2
-
- A(J) = ZERO
-
- B(J) = ZERO
-
- MX2(J) = ZERO
-
- A(J+2) = QMIXSS(SN(2),2,J)
-
- B(J+2) = -QMIXSS(SN(2),1,J)
-
- MX2(J+2) = QMIXSS(SN(2),1,J)
-
- A(J+4) = -QMIXSS(SN(3),1,J)
-
- B(J+4) = QMIXSS(SN(3),2,J)
-
- MX2(4+J) = QMIXSS(SN(3),2,J)
-
- ENDDO
-
- ELSEIF(MOD(SN(1),2).EQ.1) THEN
-
-C--Now the charged lepton, antiup,down modes
-
- DO J=1,2
-
- A(J) = ZERO
-
- B(J) = ZERO
-
- MX2(J) = ZERO
-
- A(J+2) = QMIXSS(SN(2),2,J)
-
- B(J+2) = -QMIXSS(SN(2),1,J)
-
- MX2(J+2) = QMIXSS(SN(2),1,J)
-
- A(J+4) = -QMIXSS(SN(3),1,J)
-
- B(J+4) = QMIXSS(SN(3),2,J)
-
- MX2(J+4) = QMIXSS(SN(3),2,J)
-
- ENDDO
-
- ENDIF
-
- SN(1) = SN(1) + 424
-
- SB(1) = SB(1) + 436
-
- DO K=2,3
-
- SN(K) = SN(K) + 400
-
- SB(K) = SB(K) + 412
-
- ENDDO
-
- ENDIF
-
- DO K=1,3
-
- SM(2*K-1) = RMASS(SN(K))
-
- SM(2*K) = RMASS(SB(K))
-
- SW(2*K-1) = HBAR/RLTIM(SN(K))
-
- SW(2*K) = HBAR/RLTIM(SB(K))
-
- ENDDO
-
- DO K=1,3
-
- LAMD(K) = ONE
-
- ENDDO
-
- ELSEIF(IG.GE.454.AND.IG.LE.457) THEN
-
-C--CHARGINO
-
- CSP = IG-453
-
- IF(CSP.GT.2) CSP = CSP-2
-
- AM = RMASS(IG)
-
- INFCOL = -ONE
-
- MSGN = WSGNSS(CSP)
-
- MC(1) = ONE/(SQRT(2.0D0)*MW*COSB)
-
- MC(2) = ONE/(SQRT(2.0D0)*MW*SINB)
-
-C--Calculate the combinations of the couplings needed
-
- IF(ID(1).GT.120.AND.ID(2).GT.120.AND.ID(3).GT.120) THEN
-
-C--first for the LLE modes, three modes
-
- IF(MOD(SN(1),2).EQ.0.AND.MOD(SN(3),2).EQ.0) THEN
-
-C--the one diagram mode nubar,positron, nu
-
- DO J=1,2
-
- A(J+4) = LMIXSS(SN(3)-1,1,J)*WMXUSS(CSP,1)
-
- & -RMASS(SN(3)+119)*MC(1)*LMIXSS(SN(3)-1,2,J)*WMXUSS(CSP,2)
-
- B(J+4) = ZERO
-
- MX2(J+4) = LMIXSS(SN(3)-1,2,J)
-
- ENDDO
-
- ND = 1
-
- SN(3) = SN(3)+423
-
- SB(3) = SB(3)+435
-
- ELSEIF(MOD(SN(1),2).EQ.0.AND.MOD(SN(2),2).EQ.0) THEN
-
-C--the first two diagram mode nu, nu, positron
-
- DO J=1,2
-
- A(J) = ZERO
-
- B(J) = LMIXSS(SN(1)-1,1,J)*WMXUSS(CSP,1)
-
- & -RMASS(SN(1)+119)*MC(1)*LMIXSS(SN(1)-1,2,J)*WMXUSS(CSP,2)
-
- A(J+2) = ZERO
-
- B(J+2) = LMIXSS(SN(2)-1,1,J)*WMXUSS(CSP,1)
-
- & -RMASS(SN(2)+119)*MC(1)*LMIXSS(SN(2)-1,2,J)*WMXUSS(CSP,2)
-
- MX2(J) = LMIXSS(SN(1)-1,1,J)
-
- MX2(J+2) = LMIXSS(SN(2)-1,1,J)
-
- ENDDO
-
- ND = 2
-
- DO J=1,2
-
- SN(J) = SN(J)+423
-
- SB(J) = SB(J)+435
-
- ENDDO
-
- ELSE
-
-C--the second two diagram mode positron, positron, electron
-
- DO J=1,2
-
- A(J) = -M(1)*WMXUSS(CSP,2)*MC(1)*LMIXSS(SN(1)+1,1,J)
-
- B(J) = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(1)+1,1,J)
-
- A(J+2) = -M(2)*WMXUSS(CSP,2)*MC(1)*LMIXSS(SN(2)+1,1,J)
-
- B(J+2) = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(2)+1,1,J)
-
- MX2(J) = LMIXSS(SN(1)+1,1,J)
-
- MX2(J+2) = LMIXSS(SN(2)+1,1,J)
-
- ENDDO
-
- DO J=1,2
-
- SN(J) = SN(J)+425
-
- SB(J) = SB(J)+437
-
- ENDDO
-
- ND = 2
-
- ENDIF
-
- DO K=1,3
-
- LAMD(K) = ONE
-
- ENDDO
-
- ELSEIF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
-
-C--now for the UDD
-
- IF(MOD(SN(1),2).EQ.0) THEN
-
-C--two diagram mode
-
- LAMD(1) = LAMDA3(SN(2)/2,SN(1)/2,(SN(3)+1)/2)
-
- LAMD(2) = LAMDA3(SN(1)/2,SN(2)/2,(SN(3)+1)/2)
-
- DO J=1,2
-
- A(J) = WMXUSS(CSP,1)*QMIXSS(SN(1)-1,1,J)
-
- & -RMASS(SN(1)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(1)-1,2,J)
-
- B(J) = -MSGN*M(2)*WMXVSS(CSP,2)*QMIXSS(SN(1)-1,1,J)
-
- A(J+2) = WMXUSS(CSP,1)*QMIXSS(SN(2)-1,1,J)
-
- & -RMASS(SN(2)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(2)-1,2,J)
-
- B(J+2) = -MSGN*M(2)*WMXVSS(CSP,2)*QMIXSS(SN(2)-1,1,J)
-
- MX2(J) = QMIXSS(SN(1)-1,2,J)
-
- MX2(J+2) = QMIXSS(SN(2)-1,2,J)
-
- ENDDO
-
- DO J=1,2
-
- SN(J) = SN(J) + 399
-
- SB(J) = SB(J) + 411
-
- ENDDO
-
- ND = 2
-
- ELSE
-
-C--three diagram mode
-
- LAMD(1) = LAMDA3((SN(1)+1)/2,(SN(2)+1)/2,(SN(3)+1)/2)
-
- LAMD(2) = LAMDA3((SN(2)+1)/2,(SN(1)+1)/2,(SN(3)+1)/2)
-
- LAMD(3) = LAMDA3((SN(3)+1)/2,(SN(2)+1)/2,(SN(1)+1)/2)
-
- DO I=1,3
-
- DO J=1,2
-
- A(J+2*I-2) = MSGN*(WMXVSS(CSP,1)*QMIXSS(SN(I)+1,1,J)
-
- & -RMASS(SN(I)+1)*MC(2)*WMXVSS(CSP,2)*QMIXSS(SN(I)+1,2,J))
-
- B(J+2*I-2) = -M(I)*MC(1)*WMXUSS(CSP,2)
-
- & *QMIXSS(SN(I)+1,1,J)
-
- MX2(J+2*I-2) = QMIXSS(SN(I)+1,2,J)
-
- ENDDO
-
- SN(I) = SN(I) + 401
-
- SB(I) = SB(I) + 413
-
- ENDDO
-
- ND = 3
-
- ENDIF
-
- ELSE
-
-C--now for the LQD modes
-
- IF(MOD(SN(2),2).EQ.1.AND.MOD(SN(3),2).EQ.0) THEN
-
-C--first one diagram mode nubar, dbar, up
-
- DO J=1,2
-
- A(J+4) = -MSGN*M(3)*WMXVSS(CSP,2)*MC(2)*
-
- & QMIXSS(SN(3)-1,1,J)
-
- B(J+4) = WMXUSS(CSP,1)*QMIXSS(SN(3)-1,1,J)
-
- & -RMASS(SN(3)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(3)-1,2,1)
-
- MX2(J+4) = QMIXSS(SN(3)-1,2,J)
-
- ENDDO
-
- SN(3) = SN(3) + 399
-
- SB(3) = SB(3) + 411
-
- ND = 1
-
- ELSEIF(MOD(SN(2),2).EQ.0.AND.MOD(SN(3),2).EQ.0) THEN
-
-C--second one diagram mode positron, ubar, up
-
- DO J=1,2
-
- A(J+4) = -MSGN*M(3)*WMXVSS(CSP,2)*MC(2)*
-
- & QMIXSS(SN(3)-1,1,J)
-
- B(J+4) = WMXUSS(CSP,1)*QMIXSS(SN(3)-1,1,J)
-
- & -RMASS(SN(3)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(3)-1,2,1)
-
- MX2(J+4) = QMIXSS(SN(3)-1,2,J)
-
- ENDDO
-
- SN(3) = SN(3) + 399
-
- SB(3) = SB(3) + 411
-
- ND = 1
-
- ELSEIF(MOD(SN(2),2).EQ.1.AND.MOD(SN(3),2).EQ.1) THEN
-
-C--first two diagram mode positron, dbar, down
-
- DO J=1,2
-
- A(J) = -M(1)*MC(1)*WMXUSS(CSP,2)*LMIXSS(SN(1)+1,1,J)
-
- B(J) = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(2)+1,1,J)
-
- A(J+2) = -M(2)*WMXUSS(CSP,2)*MC(1)*QMIXSS(SN(2)+1,1,J)
-
- B(J+2) = MSGN*(WMXVSS(CSP,1)*QMIXSS(SN(2)+1,1,J)
-
- & -RMASS(SN(2)+1)*MC(2)*WMXVSS(CSP,2)*QMIXSS(SN(2)+1,2,J))
-
- MX2(J) = LMIXSS(SN(1)+1,1,J)
-
- MX2(J+2) = QMIXSS(SN(2)+1,1,J)
-
- ENDDO
-
- SN(1) = SN(1) + 425
-
- SB(1) = SB(1) + 437
-
- SN(2) = SN(2) + 401
-
- SB(2) = SB(2) + 413
-
- ND = 2
-
- ELSE
-
-C--second two diagram mode nu, up, dbar
-
- DO J=1,2
-
- A(J) = ZERO
-
- B(J) = WMXUSS(CSP,1)*LMIXSS(SN(1)-1,1,J)
-
- & -RMASS(119+SN(1))*MC(1)*WMXUSS(CSP,2)*LMIXSS(SN(1)-1,2,J)
-
- A(J+2) = -MSGN*M(2)*MC(2)*WMXVSS(CSP,2)*
-
- & QMIXSS(SN(2)-1,1,J)
-
- B(J+2) = WMXUSS(CSP,1)*QMIXSS(SN(2)-1,1,J)
-
- & -RMASS(SN(2)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(2)-1,2,J)
-
- MX2(J) = LMIXSS(SN(1)-1,1,J)
-
- MX2(J+2) = QMIXSS(SN(2)-1,1,J)
-
- ENDDO
-
- SN(1) = SN(1) + 423
-
- SB(1) = SB(1) + 435
-
- SN(2) = SN(2) + 399
-
- SB(2) = SB(2) + 411
-
- ND = 2
-
- ENDIF
-
- DO K=1,3
-
- LAMD(K) = ONE
-
- ENDDO
-
- ENDIF
-
- IF(ND.EQ.1) THEN
-
- DO K=1,2
-
- SM(2*K-1) = 0.0D0
-
- SM(2*K) = 0.0D0
-
- SW(2*K-1) = 0.0D0
-
- SW(2*K) = 0.0D0
-
- ENDDO
-
- SM(5) = RMASS(SN(3))
-
- SM(6) = RMASS(SB(3))
-
- SW(5) = HBAR/RLTIM(SN(3))
-
- SW(6) = HBAR/RLTIM(SB(3))
-
- ELSE
-
- DO K=1,2
-
- SM(2*K-1) = RMASS(SN(K))
-
- SM(2*K) = RMASS(SB(K))
-
- SW(2*K-1) = HBAR/RLTIM(SN(K))
-
- SW(2*K) = HBAR/RLTIM(SB(K))
-
- SM(4+K) = ZERO
-
- SW(4+K) = ZERO
-
- ENDDO
-
- ENDIF
-
- ELSE
-
-C--UNKNOWN
-
- CALL HWWARN('HWDRME',500,*999)
-
- ENDIF
-
-C--Set mixing to zero if diagram not available
-
- IF((AM.LT.(ABS(SM(1))+M(1)).OR.ABS(SM(1)).LT.(M(2)+M(3)))
-
- & .AND.ABS(MX2(1)).GT.ZERO.AND.ND.NE.1) MX(1) = MX2(1)*LAMD(1)
-
- IF((AM.LT.(ABS(SM(2))+M(1)).OR.ABS(SM(2)).LT.(M(2)+M(3)))
-
- & .AND.ABS(MX2(2)).GT.ZERO.AND.ND.NE.1) MX(2) = MX2(2)*LAMD(1)
-
- IF((AM.LT.(ABS(SM(3))+M(2)).OR.ABS(SM(3)).LT.(M(1)+M(3)))
-
- & .AND.ABS(MX2(3)).GT.ZERO.AND.ND.NE.1) MX(3) = MX2(3)*LAMD(2)
-
- IF((AM.LT.(ABS(SM(4))+M(2)).OR.ABS(SM(4)).LT.(M(1)+M(3)))
-
- & .AND.ABS(MX2(4)).GT.ZERO.AND.ND.NE.1) MX(4) = MX2(4)*LAMD(2)
-
- IF((AM.LT.(ABS(SM(5))+M(3)).OR.ABS(SM(5)).LT.(M(1)+M(2)))
-
- & .AND.ABS(MX2(5)).GT.ZERO.AND.ND.NE.2) MX(5) = MX2(5)*LAMD(3)
-
- IF((AM.LT.(ABS(SM(6))+M(3)).OR.ABS(SM(6)).LT.(M(1)+M(2)))
-
- & .AND.ABS(MX2(6)).GT.ZERO.AND.ND.NE.2) MX(6) = MX2(6)*LAMD(3)
-
-C--Calculate the limiting points
-
- DO J=1,2
-
- IF(ND.NE.1) THEN
-
- IF(ABS(MX(J)).GT.EPS) CALL HWDRM5(M23SQT(J),M13SQT(J),
-
- & M12SQT(J),A(J),B(J),M(2),M(3),M(1),M(4),SM(J),SW(J))
-
- IF(ABS(MX(J+2)).GT.EPS) CALL HWDRM5(M13SQT(2+J),M23SQT(2+J),
-
- & M12SQT(2+J),A(2+J),B(2+J),M(1),M(3),M(2),M(4),SM(2+J),SW(2+J))
-
- ENDIF
-
- IF(ND.NE.2) THEN
-
- IF(ABS(MX(J+4)).GT.EPS) CALL HWDRM5(M12SQT(4+J),M23SQT(4+J),
-
- & M13SQT(4+J),A(4+J),B(4+J),M(1),M(2),M(3),M(4),SM(4+J),SW(4+J))
-
- ENDIF
-
- ENDDO
-
-C--Now evaluate the limit using these points
-
- LIMIT = ZERO
-
- DO 100 I=1,6
-
- IF(ABS(MX(I)).LT.EPS) GOTO 100
-
- LIMIT = LIMIT+HWDRM1(TEST,M12SQT(I),M13SQT(I),M23SQT(I),A,B,MX,
-
- & M,SM,SW,INFCOL,AM,0,ND)
-
- 100 CONTINUE
-
- LIMIT = TWO*LIMIT
-
-C--Now evaluate at a random point
-
- MTRY = 0
-
- 25 MTRY = MTRY+1
-
- LTRY = 0
-
- 35 LTRY = LTRY+1
-
- CALL HWDTHR(PHEP(1,LHEP),PHEP(1,MHEP),
-
- & PHEP(1,MHEP+1),PHEP(1,MHEP+2),HWDPWT)
-
-C--Now calculate the m12sq etc for the actual point
-
- M12SQ = M(1)**2+M(2)**2+2*HWULDO(PHEP(1,MHEP),PHEP(1,MHEP+1))
-
- M13SQ = M(1)**2+M(3)**2+2*HWULDO(PHEP(1,MHEP),PHEP(1,MHEP+2))
-
- M23SQ = M(2)**2+M(3)**2+2*HWULDO(PHEP(1,MHEP+1),PHEP(1,MHEP+2))
-
-C--Now calulate the matrix element
-
- TEST(4) = HWDRM1(TEST,M12SQ,M13SQ,M23SQ,A,B,MX,
-
- & M,SM,SW,INFCOL,AM,1,ND)
-
-C--Now test the value againest the limit
-
- RAND = HWR()*LIMIT
-
- IF(TEST(4).GT.LIMIT) THEN
-
- LIMIT = 1.1D0*TEST(4)
-
- CALL HWWARN('HWDRME',51,*150)
-
- ENDIF
-
- 150 IF(TEST(4).LT.RAND.AND.LTRY.LT.NETRY) THEN
-
- GOTO 35
-
- ELSEIF(LTRY.GE.NETRY) THEN
-
- IF(MTRY.LE.NETRY) THEN
-
- LIMIT = LIMIT*0.9D0
-
- CALL HWWARN('HWDRME',52,*25)
-
- ELSE
-
- CALL HWWARN('HWDRME',100,*999)
-
- ENDIF
-
- ENDIF
-
-C--Reorder the particles in gluino decay to get angular ordering right
-
- IF(IG.EQ.449.AND.ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
-
- DO LTRY=1,3
-
- IF(TEST(LTRY).GT.RAND) THEN
-
- IF(LTRY.EQ.2) THEN
-
- IDHWTP = IDHW(MHEP)
-
- IDHW(MHEP) = IDHW(MHEP+1)
-
- IDHW(MHEP+1) = IDHWTP
-
- IDHPTP = IDHEP(MHEP)
-
- IDHEP(MHEP) = IDHEP(MHEP+1)
-
- IDHEP(MHEP+1) = IDHPTP
-
- CALL HWVEQU(5,PHEP(1,MHEP),DECMOM)
-
- CALL HWVEQU(5,PHEP(1,MHEP+1),PHEP(1,MHEP))
-
- CALL HWVEQU(5,DECMOM,PHEP(1,MHEP+1))
-
- ELSEIF(LTRY.EQ.3) THEN
-
- IDHWTP = IDHW(MHEP)
-
- IDHW(MHEP) = IDHW(MHEP+2)
-
- IDHW(MHEP+2) = IDHWTP
-
- IDHPTP = IDHEP(MHEP)
-
- IDHEP(MHEP) = IDHEP(MHEP+2)
-
- IDHEP(MHEP+2) = IDHPTP
-
- DO I=1,5
-
- CALL HWVEQU(5,PHEP(1,MHEP),DECMOM)
-
- CALL HWVEQU(5,PHEP(1,MHEP+2),PHEP(1,MHEP))
-
- CALL HWVEQU(5,DECMOM,PHEP(1,MHEP+2))
-
- ENDDO
-
- ENDIF
-
- GOTO 52
-
- ENDIF
-
- RAND=RAND-TEST(LTRY)
-
- ENDDO
-
- ENDIF
-
- 52 CONTINUE
-
- 999 END
-
-CDECK ID>, HWDRM1.
-
-*CMZ :- -20/07/99 10:56:12 by Peter Richardson
-
-*-- Author : Peter Richardson
-
-C-----------------------------------------------------------------------
-
- FUNCTION HWDRM1(TEST,M12SQ,M13SQ,M23SQ,A,B,MX,M,SM,SW
-
- & ,INFCOL,AM,LM,ND)
-
-C-----------------------------------------------------------------------
-
-C FUNCTION TO GIVE THE R-PARITY VIOLATING MATRIX ELEMENT AT A GIVEN
-
-C PHASE-SPACE POINT
-
-C-----------------------------------------------------------------------
-
- IMPLICIT NONE
-
- DOUBLE PRECISION M12SQ,M13SQ,M23SQ,MX(6),A(6),B(6),SM(6),SW(6),
-
- & INFCOL,AM,TERM(21),TEST(4),PLN,NPLN,ZERO,
-
- & M(4),HWDRM1,HWDRM2,HWDRM3,HWDRM4
-
- PARAMETER (ZERO=0)
-
- EXTERNAL HWDRM2,HWDRM3,HWDRM4
-
- INTEGER LM,K,ND
-
-C--Zero the array
-
- DO K=1,21
-
- TERM(K) = 0.0D0
-
- ENDDO
-
- HWDRM1 = 0.0D0
-
-C--The amplitude
-
- IF(ABS(MX(1)).GT.ZERO.AND.ND.NE.1) THEN
-
- TERM(1) = MX(1)**2*HWDRM2(M23SQ,M(2),M(3),M(1),M(4),SM(1),
-
- & SW(1),A(1),B(1))
-
- IF(ABS(MX(2)).GT.ZERO) TERM(7)= MX(1)*MX(2)*HWDRM3(M23SQ,M(2),
-
- & M(3),M(1),M(4),SM(1),SM(2),SW(1),SW(2),A(1),A(2),B(1),B(2))
-
- IF(ABS(MX(3)).GT.ZERO) TERM(10)=-MX(1)*MX(3)*HWDRM4(M13SQ,M23SQ,
-
- & M(1),M(3),M(2),M(4),SM(3),SM(1),SW(3),SW(1),A(1),A(3),B(1),B(3))
-
- IF(ABS(MX(4)).GT.ZERO) TERM(11)=-MX(1)*MX(4)*HWDRM4(M13SQ,M23SQ,
-
- & M(1),M(3),M(2),M(4),SM(4),SM(1),SW(4),SW(1),A(1),A(4),B(1),B(4))
-
- IF(ABS(MX(5)).GT.ZERO) TERM(12)=-MX(1)*MX(5)*HWDRM4(M23SQ,M12SQ,
-
- & M(3),M(2),M(1),M(4),SM(1),SM(5),SW(1),SW(5),A(5),A(1),B(5),B(1))
-
- IF(ABS(MX(6)).GT.ZERO) TERM(13)=-MX(1)*MX(6)*HWDRM4(M23SQ,M12SQ,
-
- & M(3),M(2),M(1),M(4),SM(1),SM(6),SW(1),SW(6),A(6),A(1),B(6),B(1))
-
- ENDIF
-
- IF(ABS(MX(2)).GT.ZERO.AND.ND.NE.1) THEN
-
- TERM(2) = MX(2)**2*HWDRM2(M23SQ,M(2),M(3),M(1),M(4),SM(2),
-
- & SW(2),A(2),B(2))
-
- IF(ABS(MX(3)).GT.ZERO) TERM(14)=-MX(2)*MX(3)*HWDRM4(M13SQ,M23SQ,
-
- & M(1),M(3),M(2),M(4),SM(3),SM(2),SW(3),SW(2),A(2),A(3),B(2),B(3))
-
- IF(ABS(MX(4)).GT.ZERO) TERM(15)=-MX(2)*MX(4)*HWDRM4(M13SQ,M23SQ,
-
- & M(1),M(3),M(2),M(4),SM(4),SM(2),SW(4),SW(2),A(2),A(4),B(2),B(4))
-
- IF(ABS(MX(5)).GT.ZERO) TERM(16)=-MX(2)*MX(5)*HWDRM4(M23SQ,M12SQ,
-
- & M(3),M(2),M(1),M(4),SM(2),SM(5),SW(2),SW(5),A(5),A(2),B(5),B(2))
-
- IF(ABS(MX(6)).GT.ZERO) TERM(17)=-MX(2)*MX(6)*HWDRM4(M23SQ,M12SQ,
-
- & M(3),M(2),M(1),M(4),SM(2),SM(6),SW(2),SW(6),A(6),A(2),B(6),B(2))
-
- ENDIF
-
- IF(ABS(MX(3)).GT.ZERO.AND.ND.NE.1) THEN
-
- TERM(3) = MX(3)**2*HWDRM2(M13SQ,M(1),M(3),M(2),M(4),SM(3),
-
- & SW(3),A(3),B(3))
-
- IF(ABS(MX(4)).GT.ZERO) TERM(8)= MX(3)*MX(4)*HWDRM3(M13SQ,M(1),
-
- & M(3),M(2),M(4),SM(3),SM(4),SW(3),SW(4),A(3),A(4),B(3),B(4))
-
- IF(ABS(MX(5)).GT.ZERO) TERM(18)=-MX(3)*MX(5)*HWDRM4(M12SQ,M13SQ,
-
- & M(2),M(1),M(3),M(4),SM(5),SM(3),SW(5),SW(3),A(3),A(5),B(3),B(5))
-
- IF(ABS(MX(6)).GT.ZERO) TERM(19)=-MX(3)*MX(6)*HWDRM4(M12SQ,M13SQ,
-
- & M(2),M(1),M(3),M(4),SM(6),SM(3),SW(6),SW(3),A(3),A(6),B(3),B(6))
-
- ENDIF
-
- IF(ABS(MX(4)).GT.ZERO.AND.ND.NE.1) THEN
-
- TERM(4) = MX(4)**2*HWDRM2(M13SQ,M(1),M(3),M(2),M(4),SM(4),
-
- & SW(4),A(4),B(4))
-
- IF(ABS(MX(5)).GT.ZERO) TERM(20)=-MX(4)*MX(5)*HWDRM4(M12SQ,M13SQ,
-
- & M(2),M(1),M(3),M(4),SM(5),SM(4),SW(5),SW(4),A(4),A(5),B(4),B(5))
-
- IF(ABS(MX(6)).GT.ZERO) TERM(21)=-MX(4)*MX(6)*HWDRM4(M12SQ,M13SQ,
-
- & M(2),M(1),M(3),M(4),SM(6),SM(4),SW(6),SW(4),A(4),A(6),B(4),B(6))
-
- ENDIF
-
- IF(ABS(MX(5)).GT.ZERO.AND.ND.NE.2) THEN
-
- TERM(5) = MX(5)**2*HWDRM2(M12SQ,M(1),M(2),M(3),M(4),SM(5),
-
- & SW(5),A(5),B(5))
-
- IF(ABS(MX(6)).GT.ZERO) TERM(9)= MX(5)*MX(6)*HWDRM3(M12SQ,M(1),
-
- & M(2),M(3),M(4),SM(5),SM(6),SW(5),SW(6),A(5),A(6),B(5),B(6))
-
- ENDIF
-
- IF(ABS(MX(6)).GT.ZERO.AND.ND.NE.2) TERM(6) = MX(6)**2*
-
- & HWDRM2(M12SQ,M(1),M(2),M(3),M(4),SM(6),SW(6),A(6),B(6))
-
- DO K=10,21
-
- TERM(K)=TERM(K)*INFCOL
-
- ENDDO
-
-C--Add them up
-
- DO K=1,21
-
- HWDRM1 = HWDRM1+TERM(K)
-
- ENDDO
-
-C--Different colour flows for the gluino
-
- IF(LM.NE.0) THEN
-
- NPLN = 0.0D0
-
- PLN = 0.0D0
-
- DO K=1,9
-
- PLN = PLN+TERM(K)
-
- ENDDO
-
- DO K=10,21
-
- NPLN= NPLN+TERM(K)
-
- ENDDO
-
- DO K=1,3
-
- TEST(K) = (TERM(2*K-1)+TERM(2*K)+TERM(6+K))*(1+NPLN/PLN)
-
- ENDDO
-
- ELSE
-
- DO K=1,3
-
- TEST(K) = 0.0D0
-
- ENDDO
-
- ENDIF
-
- IF(TEST(4).LT.ZERO) CALL HWWARN('HWDRM1',50,*999)
-
- 999 END
-
-CDECK ID>, HWDRM2.
-
-*CMZ :- -20/07/99 10:56:12 by Peter Richardson
-
-*-- Author : Peter Richardson
-
-C-----------------------------------------------------------------------
-
- FUNCTION HWDRM2(X,MA,MB,MC,MD,MR1,GAM1,A,B)
-
-C-----------------------------------------------------------------------
-
-C Function to compute the matrix element squared part of a 3-body
-
-C R-parity decay
-
-C-----------------------------------------------------------------------
-
- IMPLICIT NONE
-
- DOUBLE PRECISION X,MA,MB,MC,MD,A,B,HWDRM2,MR1,GAM1
-
- HWDRM2 = (X - MA**2 - MB**2)*(4*A*B*MC*MD +
-
- & (A**2 + B**2)*(-X + MC**2 + MD**2))/
-
- & ((X-MR1**2)**2+GAM1**2*MR1**2)
-
- END
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff