Merge lp:~maddevelopers/mg5amcnlo/release_candidate into lp:~madteam/mg5amcnlo/trunk

[Johan Alwall]

Too many changes to detail here...

[Johan Alwall (johan-alwall)]

Please check out this version and the new web interface on http://cp3wks05.fynu.ucl.ac.be/.
We are hoping to do the merge by tomorrow afternoon (April 12), as soon as I get up :-)

[Olivier Mattelaer (olivier-mattelaer)]

> Please check out this version and the new web interface on
> http://cp3wks05.fynu.ucl.ac.be/.
> We are hoping to do the merge by tomorrow afternoon (April 12), as soon as I
> get up :-)

Hi Johan,

Concerning the '+' problem with pythia8. I think this sould be resolve in the pythia 8 interface.
i.e made a string replacement for that. This type of quite stupid rule (you cann't name things with...) should be avoid as much as possible.

Otherwise, I agree to pass in golden master :-)

[Olivier Mattelaer (olivier-mattelaer)]

another (very) minor point:
one of the logger is not beautiful:
root: Finding symmetric diagrams for process g g > g g QED=0
could you try to pass in the normal format?

[lp:~maddevelopers/mg5amcnlo/release_candidate]

138. By Olivier Mattelaer on 2011-04-12

add nmssm and correct logging format in write_symmetry

139. By Johan Alwall on 2011-04-12

Allow for + in model names with C++ output by replacing + with _plus_ in class and file names

140. By Johan Alwall on 2011-04-12

Updated VERSION to 1.0.0

141. By Johan Alwall on 2011-04-12

Ensure that general param_cards can be read by model, in particular MG4 sm param_card (except for yukawa couplings, which are in different blocks)

=== added file '.bzrignore'
--- .bzrignore	1970-01-01 00:00:00 +0000
+++ .bzrignore	2011-04-12 19:56:30 +0000
@@ -0,0 +1,1 @@
+models/RS/.DS_Store

=== added directory 'DECAY'
=== added file 'DECAY/DECAYVersion.txt'
--- DECAY/DECAYVersion.txt	1970-01-01 00:00:00 +0000
+++ DECAY/DECAYVersion.txt	2011-04-12 19:56:30 +0000
@@ -0,0 +1,1 @@
+1.0.5

=== added file 'DECAY/README'
--- DECAY/README	1970-01-01 00:00:00 +0000
+++ DECAY/README	2011-04-12 19:56:30 +0000
@@ -0,0 +1,13 @@
+IMPORTANT: this is strictly a SM Decay routine!!!
+AUTHOR: Fabio Maltoni
+NAME  : decay
+DESCRIPTION:
+c***********************************************************************
+c     1. reads events from file                                        *
+c     2. asks  which particle to decay                                 *
+c     3. asks  the decay mode                                          *
+c     4. decays particle                                               *
+c     5. write out decayed events.                                     *
+c***********************************************************************
+
+

=== added file 'DECAY/UpdateNotes.txt'
--- DECAY/UpdateNotes.txt	1970-01-01 00:00:00 +0000
+++ DECAY/UpdateNotes.txt	2011-04-12 19:56:30 +0000
@@ -0,0 +1,48 @@
+Update notes for DECAY package
+
+c----------------------------------------------------------------------*
+c     FIRST VERSION 16-May-2003                                        *
+c----------------------------------------------------------------------*
+c     LAST UPDATE  26-Sep-2003:                                        *
+c     - W->anything and t->b anything decays added.                    *
+c     - rnd number generator seed changes in sequential runs.          *
+c----------------------------------------------------------------------*
+c     LAST UPDATE  07-Nov-2003:                                        *
+c     - exactly turns unweighted evts in unweighted evts.              *
+c     - decaying identical particles in one event fixed.               *
+c     - bug in the Z->jets quark fractions corrected.                  *
+c     - error trap routines added                                      *
+c----------------------------------------------------------------------*
+c     LAST UPDATE 25-Feb-2004:                                         *
+c     - H->ZZ>leptons : wrong extra factor of two deleted              *
+c----------------------------------------------------------------------*
+c     LAST UPDATE 12-Dec-2004:                                         *
+c     - Merging it to the official web template                        *
+c----------------------------------------------------------------------*
+c     LAST UPDATE  June-10-2005:                                       *
+c     - Updated for use with new rw_events.f format                    *
+c     LAST UPDATE  August-18-2005:                                     *
+c     - Fixed a bug in  the previous update                            *
+c----------------------------------------------------------------------*
+c     LAST UPDATE 25-Jan-2006:                                         *
+c     - Template change structure, moved to directory Decay            *
+c----------------------------------------------------------------------*
+c     LAST UPDATE 27-July-2006:                                        *
+c     - Compliant with MadGraph_v4.0                                   *
+c----------------------------------------------------------------------*
+c     LAST UPDATE Jan 2007                                             *
+c     - Compliant with MadGraph_v4.1                                   *
+c----------------------------------------------------------------------*
+c     LAST UPDATE Mar 2008                                             *
+c     - RF: Fixed bug in writing rnd number seed in event file banner  *
+c----------------------------------------------------------------------*
+c     LAST UPDATE May 2008                                             *
+c     - FM: Fixed bug in Higgs 4 body decay                            *
+c----------------------------------------------------------------------*
+c     LAST UPDATE June 2008                                            *
+c     - RF: Fixed bug if there is more than one event number in the    *
+c           init block of event file.                                  *
+c----------------------------------------------------------------------*
+c     LAST UPDATE Sep 2009                                             *
+c     - OM: Modification for gfortran compilation compatibility        *
+c***********************************************************************

=== added file 'DECAY/alfas.inc'
--- DECAY/alfas.inc	1970-01-01 00:00:00 +0000
+++ DECAY/alfas.inc	2011-04-12 19:56:30 +0000
@@ -0,0 +1,11 @@
+c***********************************************************************
+c     this files contains the common blocks for the 
+c     the alpha_s settings
+c
+c     asmz   = alpha_s(Mz) is set based on the pdf chosen in setcuts.f
+c     nloop  = order of the running of alpha_s based on the pdf chosen  
+c***********************************************************************
+      integer nloop
+      double precision asmz
+      common/a_block/asmz,nloop
+

=== added file 'DECAY/alfas_functions.f'
--- DECAY/alfas_functions.f	1970-01-01 00:00:00 +0000
+++ DECAY/alfas_functions.f	2011-04-12 19:56:30 +0000
@@ -0,0 +1,280 @@
+C
+C-----------------------------------------------------------------------------
+C
+      double precision function alfa(alfa0,qsq )
+C
+C-----------------------------------------------------------------------------
+C
+C	This function returns the 1-loop value of alpha.
+C
+C	INPUT: 
+C		qsq   = Q^2
+C
+C-----------------------------------------------------------------------------
+C
+      implicit none
+      double precision  qsq,alfa0
+c
+c constants
+c
+      double precision  One, Three, Pi,zmass
+      parameter( One = 1.0d0, Three = 3.0d0 )
+      parameter( Pi = 3.14159265358979323846d0 )
+      parameter( zmass = 91.188d0 )
+cc
+      alfa = alfa0 / ( 1.0d0 - alfa0*dlog( qsq/zmass**2 ) /Three /Pi )
+ccc
+      return
+      end
+
+C
+C-----------------------------------------------------------------------------
+C
+      double precision function alfaw(alfaw0,qsq,nh )
+C
+C-----------------------------------------------------------------------------
+C
+C	This function returns the 1-loop value of alpha_w.
+C
+C	INPUT: 
+C		qsq = Q^2
+C               nh  = # of Higgs doublets
+C
+C-----------------------------------------------------------------------------
+C
+      implicit none
+      double precision  qsq, alphaw, dum,alfaw0
+      integer  nh, nq
+c
+c	  include
+c
+	  
+c
+c constants
+c
+      double precision  Two, Four, Pi, Twpi, zmass,tmass
+      parameter( Two = 2.0d0, Four = 4.0d0 )
+      parameter( Pi = 3.14159265358979323846d0 )
+      parameter( Twpi = 3.0d0*Four*Pi )
+      parameter( zmass = 91.188d0,tmass=174d0 )
+cc
+      if ( qsq.ge.tmass**2 ) then
+         nq = 6
+      else
+         nq = 5
+      end if
+      dum = ( 22.0d0 - Four*nq - nh/Two ) / Twpi
+      alfaw = alfaw0 / ( 1.0d0 + dum*alfaw0*dlog( qsq/zmass**2 ) )
+ccc
+      return
+      end
+
+C-----------------------------------------------------------------------------
+C
+      DOUBLE PRECISION FUNCTION ALPHAS(Q)
+c
+c     Evaluation of strong coupling constant alpha_S
+c     Author: R.K. Ellis
+c
+c     q -- scale at which alpha_s is to be evaluated
+c
+c-- common block alfas.inc
+c     asmz -- value of alpha_s at the mass of the Z-boson
+c     nloop -- the number of loops (1,2, or 3) at which beta 
+c
+c     function is evaluated to determine running.
+c     the values of the cmass and the bmass should be set
+c     in common block qmass.
+C-----------------------------------------------------------------------------
+
+      IMPLICIT NONE
+c
+      include 'alfas.inc'
+      DOUBLE PRECISION Q,T,AMZ0,AMB,AMC
+      DOUBLE PRECISION AS_OUT
+      INTEGER NLOOP0,NF3,NF4,NF5
+      PARAMETER(NF5=5,NF4=4,NF3=3)
+C
+      REAL*8       CMASS,BMASS
+      COMMON/QMASS/CMASS,BMASS
+      DATA CMASS,BMASS/1.42D0,4.7D0/  ! HEAVY QUARK MASSES FOR THRESHOLDS
+C
+      REAL*8 ZMASS
+      DATA ZMASS/91.188D0/
+C
+      SAVE AMZ0,NLOOP0,AMB,AMC
+      DATA AMZ0,NLOOP0/0D0,0/
+      IF (Q .LE. 0D0) THEN 
+         WRITE(6,*) 'q .le. 0 in alphas'
+         WRITE(6,*) 'q= ',Q
+         STOP
+      ENDIF
+      IF (asmz .LE. 0D0) THEN 
+         WRITE(6,*) 'asmz .le. 0 in alphas',asmz
+c         WRITE(6,*) 'continue with asmz=0.1185'
+         STOP
+         asmz=0.1185D0
+      ENDIF
+      IF (CMASS .LE. 0.3D0) THEN 
+         WRITE(6,*) 'cmass .le. 0.3GeV in alphas',CMASS
+c         WRITE(6,*) 'continue with cmass=1.5GeV?'
+         STOP
+         CMASS=1.42D0
+      ENDIF
+      IF (BMASS .LE. 0D0) THEN 
+         WRITE(6,*) 'bmass .le. 0 in alphas',BMASS
+         WRITE(6,*) 'COMMON/QMASS/CMASS,BMASS'
+         STOP
+         BMASS=4.7D0
+      ENDIF
+c--- establish value of coupling at b- and c-mass and save
+      IF ((asmz .NE. AMZ0) .OR. (NLOOP .NE. NLOOP0)) THEN
+         AMZ0=asmz
+         NLOOP0=NLOOP
+         T=2D0*DLOG(BMASS/ZMASS)
+         CALL NEWTON1(T,asmz,AMB,NLOOP,NF5)
+         T=2D0*DLOG(CMASS/BMASS)
+         CALL NEWTON1(T,AMB,AMC,NLOOP,NF4)
+      ENDIF
+
+c--- evaluate strong coupling at scale q
+      IF (Q  .LT. BMASS) THEN
+           IF (Q  .LT. CMASS) THEN
+             T=2D0*DLOG(Q/CMASS)
+             CALL NEWTON1(T,AMC,AS_OUT,NLOOP,NF3)
+           ELSE
+             T=2D0*DLOG(Q/BMASS)
+             CALL NEWTON1(T,AMB,AS_OUT,NLOOP,NF4)
+           ENDIF
+      ELSE
+      T=2D0*DLOG(Q/ZMASS)
+      CALL NEWTON1(T,asmz,AS_OUT,NLOOP,NF5)
+      ENDIF
+      ALPHAS=AS_OUT
+      RETURN
+      END
+
+
+      SUBROUTINE NEWTON1(T,A_IN,A_OUT,NLOOP,NF)
+C     Author: R.K. Ellis
+
+c---  calculate a_out using nloop beta-function evolution 
+c---  with nf flavours, given starting value as-in
+c---  given as_in and logarithmic separation between 
+c---  input scale and output scale t.
+c---  Evolution is performed using Newton's method,
+c---  with a precision given by tol.
+
+      IMPLICIT NONE
+      INTEGER NLOOP,NF
+      REAL*8 T,A_IN,A_OUT,AS,TOL,F2,F3,F,FP,DELTA
+      REAL*8 B0(3:5),C1(3:5),C2(3:5),DEL(3:5)
+      PARAMETER(TOL=5.D-4)
+C---     B0=(11.-2.*NF/3.)/4./PI
+      DATA B0/0.716197243913527D0,0.66314559621623D0,0.61009394851893D0/
+C---     C1=(102.D0-38.D0/3.D0*NF)/4.D0/PI/(11.D0-2.D0/3.D0*NF)
+      DATA C1/.565884242104515D0,0.49019722472304D0,0.40134724779695D0/
+C---     C2=(2857.D0/2.D0-5033*NF/18.D0+325*NF**2/54)
+C---     /16.D0/PI**2/(11.D0-2.D0/3.D0*NF)
+      DATA C2/0.453013579178645D0,0.30879037953664D0,0.14942733137107D0/
+C---     DEL=SQRT(4*C2-C1**2)
+      DATA DEL/1.22140465909230D0,0.99743079911360D0,0.66077962451190D0/
+      F2(AS)=1D0/AS+C1(NF)*LOG((C1(NF)*AS)/(1D0+C1(NF)*AS))
+      F3(AS)=1D0/AS+0.5D0*C1(NF)
+     & *LOG((C2(NF)*AS**2)/(1D0+C1(NF)*AS+C2(NF)*AS**2))
+     & -(C1(NF)**2-2D0*C2(NF))/DEL(NF)
+     & *ATAN((2D0*C2(NF)*AS+C1(NF))/DEL(NF))
+
+           
+      A_OUT=A_IN/(1D0+A_IN*B0(NF)*T)
+      IF (NLOOP .EQ. 1) RETURN
+      A_OUT=A_IN/(1D0+B0(NF)*A_IN*T+C1(NF)*A_IN*LOG(1D0+A_IN*B0(NF)*T))
+      IF (A_OUT .LT. 0D0) AS=0.3D0
+ 30   AS=A_OUT
+
+      IF (NLOOP .EQ. 2) THEN
+      F=B0(NF)*T+F2(A_IN)-F2(AS)
+      FP=1D0/(AS**2*(1D0+C1(NF)*AS))
+      ENDIF
+      IF (NLOOP .EQ. 3) THEN
+      F=B0(NF)*T+F3(A_IN)-F3(AS)
+      FP=1D0/(AS**2*(1D0+C1(NF)*AS+C2(NF)*AS**2))
+      ENDIF
+      A_OUT=AS-F/FP
+      DELTA=ABS(F/FP/AS)
+      IF (DELTA .GT. TOL) GO TO 30
+      RETURN
+      END
+
+
+C-----------------------------------------------------------------------------
+C
+      double precision function mfrun(mf,scale,asmz,nloop)
+C
+C-----------------------------------------------------------------------------
+C
+C	This function returns the 2-loop value of a MSbar fermion mass
+C       at a given scale.
+C
+C	INPUT: mf    = MSbar mass of fermion at MSbar fermion mass scale 
+C	       scale = scale at which the running mass is evaluated
+C	       asmz  = AS(MZ) : this is passed to alphas(scale,asmz,nloop)
+C              nloop = # of loops in the evolution
+C       
+C
+C
+C	EXTERNAL:      double precision alphas(scale,asmz,nloop)
+C                      
+C-----------------------------------------------------------------------------
+C
+      implicit none
+C
+C     ARGUMENTS
+C
+      double precision  mf,scale,asmz
+      integer           nloop
+C
+C     LOCAL
+C
+      double precision  beta0, beta1,gamma0,gamma1
+      double precision  A1,as,asmf,l2
+      integer  nf
+C
+C     EXTERNAL
+C
+      double precision  alphas
+      external          alphas
+c
+c     CONSTANTS
+c
+      double precision  One, Two, Three, Pi
+      parameter( One = 1.0d0, Two = 2.0d0, Three = 3.0d0 )
+      parameter( Pi = 3.14159265358979323846d0) 
+      double precision tmass
+      parameter(tmass=174d0)
+cc
+C
+C
+      if ( mf.gt.tmass ) then
+         nf = 6
+      else
+         nf = 5
+      end if
+
+      beta0 = ( 11.0d0 - Two/Three *nf )/4d0
+      beta1 = ( 102d0  - 38d0/Three*nf )/16d0
+      gamma0= 1d0
+      gamma1= ( 202d0/3d0  - 20d0/9d0*nf )/16d0
+      A1    = -beta1*gamma0/beta0**2+gamma1/beta0
+      as    = alphas(scale)
+      asmf  = alphas(mf)
+      l2    = (1+ A1*as/Pi)/(1+ A1*asmf/Pi)
+      
+      
+      mfrun = mf * (as/asmf)**(gamma0/beta0)
+
+      if(nloop.eq.2) mfrun =mfrun*l2
+ccc
+      return
+      end
+

=== added file 'DECAY/calc_values.inc'
--- DECAY/calc_values.inc	1970-01-01 00:00:00 +0000
+++ DECAY/calc_values.inc	2011-04-12 19:56:30 +0000
@@ -0,0 +1,17 @@
+c
+c     Calculated intermediate values that are printed out
+c
+      double precision  decw, w_w_nl, w_w_tau, w_w_ud, w_w_cs
+      double precision  decz, w_z_nn, w_z_ll, w_z_tau
+      double precision  w_z_uu, w_z_dd, w_z_cc, w_z_bb
+      double precision  mt_h,mb_h,mc_h
+      double precision  w_h_tt,w_h_bb,w_h_tau,w_h_cc
+      double precision  w_h_ww,w_h_zz,w_h_WLWL,w_h_ZLZL
+      common/calc_values/
+     &                   decw, w_w_nl, w_w_tau, w_w_ud, w_w_cs,
+     &                   decz, w_z_nn, w_z_ll, w_z_tau,
+     &                   w_z_uu, w_z_dd, w_z_cc, w_z_bb,
+     &                   mt_h, mb_h, mc_h,
+     &                   w_h_tt,w_h_bb,w_h_tau,w_h_cc,
+     &                   w_h_ww,w_h_zz,w_h_WLWL,w_h_ZLZL
+

=== added file 'DECAY/coupl.inc'
--- DECAY/coupl.inc	1970-01-01 00:00:00 +0000
+++ DECAY/coupl.inc	2011-04-12 19:56:30 +0000
@@ -0,0 +1,50 @@
+c====================================================================
+c
+c  Define common block containing all coupling constants and masses
+c  which are used in the HELAS routines.
+c
+c  These include masses, widths and real/complex couplings.
+c
+c  This file can be built automatically from particles.dat and
+c  interactions.dat
+c
+c====================================================================
+c
+c     QCD 
+c
+      double complex       gg(2)
+      double precision     g
+      common /COUPL_QCD/   g,gg   
+c
+c    kinematical masses
+c
+      double precision     hmass, wmass, zmass, amass,
+     &                     tmass, bmass, lmass, cmass
+      common /COUPL_MASS/  hmass, wmass, zmass, amass,
+     &                     tmass, bmass, lmass, cmass
+c
+c   widths
+c
+      double precision     hwidth, wwidth, zwidth, 
+     &                     twidth, lwidth, awidth
+      common /COUPL_WIDTH/ hwidth, wwidth, zwidth, 
+     &                     twidth, lwidth, awidth
+c
+c  couplings in the feynman rules
+c
+      double complex       gal(2), gad(2), gau(2), gwf(2),
+     &                     gzn(2), gzl(2), gzd(2), gzu(2)
+      double precision     gw, gwwa, gwwz
+      common /COUPL_GAUGE/ gal   , gad   , gau   , gwf   ,
+     &                     gzn   , gzl   , gzd   , gzu   ,
+     &                     gw, gwwa, gwwz
+c
+      double complex       gwfc(2),  gwfs(2), gwfm(2)
+      common /coupl_ckm/   gwfc,     gwfs   , gwfm	
+c
+      double complex       gwwh, gzzh, gwwhh, gzzhh, ghhh, ghhhh
+      common /COUPL_SCAL/  gwwh, gzzh, gwwhh, gzzhh, ghhh, ghhhh
+c
+      double complex       ghtop(2), ghbot(2), ghtau(2), ghcha(2)
+      common /COUPL_YUK/   ghtop   , ghbot   , ghtau   , ghcha
+

=== added file 'DECAY/decay.f'
--- DECAY/decay.f	1970-01-01 00:00:00 +0000
+++ DECAY/decay.f	2011-04-12 19:56:30 +0000
@@ -0,0 +1,1613 @@
+      program decay
+c***********************************************************************
+c     1. reads events from file                                        *
+c     2. asks  which particle to decay                                 *
+c     3. asks  the decay mode                                          *
+c     4. decays particle                                               *
+c     5. write out decayed events.                                     *    
+c----------------------------------------------------------------------*
+c     FIRST VERSION 16-May-2003                                        *
+c----------------------------------------------------------------------*
+c     LAST UPDATE  26-Sep-2003:                                        *
+c     - W->anything and t->b anything decays added.                    *
+c     - rnd number generator seed changes in sequential runs.          *
+c----------------------------------------------------------------------*
+c     LAST UPDATE  07-Nov-2003:                                        *
+c     - exactly turns unweighted evts in unweighted evts.              * 
+c     - decaying identical particles in one event fixed.               * 
+c     - bug in the Z->jets quark fractions corrected.                  *
+c     - error trap routines added                                      *
+c----------------------------------------------------------------------*
+c     LAST UPDATE 25-Feb-2004:                                         *
+c     - H->ZZ>leptons : wrong extra factor of two deleted              * 
+c----------------------------------------------------------------------*
+c     LAST UPDATE 12-Dec-2004:                                         *
+c     - Merging it to the official web template                        *
+c----------------------------------------------------------------------*
+c     LAST UPDATE  June-10-2005:                                       *
+c     - Updated for use with new rw_events.f format                    *
+c     LAST UPDATE  August-18-2005:                                     *
+c     - Fixed a bug in  the previous update                            * 
+c----------------------------------------------------------------------*
+c     LAST UPDATE  July-27-2006:                                        *
+c     - Make it compliant with MadGraph v 4.0                          *
+c***********************************************************************
+      implicit none
+      include 'decay.inc'
+c
+c     parameters
+c
+      real*8     smallmass
+      parameter (smallmass=1d0) 
+c     
+c     Local
+c
+c
+      integer nexternal, ic(7,MaxParticles),idecay(MaxParticles)
+      double precision P(0:4,MaxParticles),wgt
+      integer i,k,iten
+      integer nevent,nfound(0:MaxParticles)
+      data nevent/0/
+c     
+      real*8 sum(0:maxievent),maxwgt(0:maxievent)
+      character*50 dec_name
+      character*3 name
+      character*132 buff
+      character*140 buff2
+      integer iseed
+      real*8 aa
+      real*8 scale,aqcd,aqed
+      integer ievent,eventnr,eventnumber
+      integer maxpup
+      parameter(maxpup=100)
+      integer idbmup(2),pdfgup(2),pdfsup(2),idwtup,nprup,lprup
+      double precision ebmup(2),xsecup,xerrup,xmaxup
+C      
+      logical done,firsttime,fopenend,lwrite,newbanner
+      LOGICAL DEBUG
+      data firsttime/.true./
+      data done/.false./
+      DATA DEBUG/.TRUE./
+c
+c     external
+c     
+      real*8 dot
+      real xran1
+c
+c     Global
+c
+      include 'coupl.inc'
+      data id /201*' '/
+c
+      character*70 infile,outfile
+      integer run_mode
+      common/mode/run_mode,infile
+c
+      integer iunit
+      common/to_lh/iunit
+
+      integer sizeievent
+      common/ievents/sizeievent
+c
+c     Date
+c
+c      CHARACTER*20 DATE
+c      EXTERNAL DATE_Y2KBUG
+      
+c-----
+c   Begin Code
+c-----
+c
+      write(*,*) '*****************************************************'
+      write(*,*) '*                    DECAY                          *'
+      write(*,*) '*              a MadEvent program                   *'
+      write(*,*) '*        for decaying unstable particles            *'
+      write(*,*) '*             in the Standar Model                  *'  
+      write(*,*) '*        ---------------------------------          *'
+      write(*,*) '*        version compliant with MG_ME_V4.0          *'
+      write(*,*) '*                                                   *'
+      write(*,*) '*                27-July-2006                       *'
+      write(*,*) '*****************************************************'
+c
+c     run mode: 0 calculates partial widths
+c               1 decays event in file
+c
+      write(*,*)
+      write(*,*) '  Input run mode:'
+      write(*,*) '  ---------------'
+      write(*,*)
+      write(*,*) '  0 = calculates decay widths'
+      write(*,*) '  1 = decay events in file'
+      write(*,*)
+            
+      read(*,'(i1)')  run_mode
+
+c
+c     initialize rnd number for decay 
+c
+      call rnd_init(lunrnd,iseed)  
+c
+c
+c     open scratch file which will contain the param_card.dat
+c
+      open(lunp,status='scratch')
+c
+      if(run_mode.eq.0) then
+c
+c     calculates the decay width at LO
+c
+       write(*,*) '*****************************************************'
+       write(*,*) '*  run_mode=0 => calculting decay widths            *'
+       write(*,*) '*****************************************************' 
+       write(*,*)
+       write(*,*) '  Input parameter file (e.g. ../Cards/param_card.dat)'
+       write(*,*) '  ---------------------------------------------------'
+       write(*,*)
+       read(*,'(a)')  infile
+       open(lunr,file=infile,status='old',err=102)
+       write(*,*)
+       write(*,*) '*****************************************************'
+       write(*,*) '* Using file:                                       *'
+       write(*,'(1x,a1,a25,a1)') ' ',infile,' '
+       write(*,*) '* for the input params.                             *'
+       write(*,*) '*                                                   *'
+       write(*,*) '* >>>>>>Total widths are recalculated here<<<<<     *'
+       write(*,*) '*****************************************************'
+       write(*,*)
+
+       goto 103 
+ 102   write(*,*) 'file', infile ,' not found : stopping here'
+       stop
+ 103   continue
+
+c
+c      write the param information into the scratch file
+c
+      done=.false.
+      do while(.not.done)
+         read(lunr,'(a132)',err=101,end=101) buff
+         write(lunp,'(a132)') buff
+      enddo
+c
+ 101  continue
+       
+      elseif(run_mode.eq.1) then
+
+       write(*,*) '*****************************************************'
+       write(*,*) '*     run_mode=1 => decaying events                 *'
+       write(*,*) '*                                                   *'
+       write(*,*) '* Using the param_card.dat in the banner for        *'
+       write(*,*) '* the input params.                                 *'
+       write(*,*) '*                                                   *'
+       write(*,*) '* >>>>>>Total widths are recalculated here<<<<<     *'
+       write(*,*) '*****************************************************'
+c 
+      write(*,*) 'input event file: (e.g. events.lhe)'
+      read(*,'(a)')  infile
+c
+      write(*,*) 'name for output file: (e.g. dec-events.lhe)'
+      read(*,'(a)')  outfile
+c
+      open(lunr,file=infile,status='old')
+      open(lunw,status='scratch')
+      open(luni,file=outfile,status='unknown')
+c
+c     copy old banner into new banner and the param_card.dat into the scratch file
+c
+      done=.false.
+      lwrite=.false.
+      newbanner=.false.
+      read(lunr,'(a132)',err=99) buff
+      do while(index(buff,'<init>') .eq. 0)
+         if(index(buff,"<header>") .ne. 0) newbanner=.true.
+         if (index(buff,"<slha>") .ne. 0 .or.
+     $        index(buff,"Begin param_card.dat") .ne. 0) lwrite=.true.
+         if (index(buff,"</slha>") .ne. 0 .or.
+     $        index(buff,"End param_card.dat") .ne. 0) lwrite=.false.
+c         
+         if(index(buff,'</header>') .eq. 0 .and.
+     $        (newbanner .or. index(buff,'-->') .eq. 0 ))
+     $      write(luni,'(a)') buff(1:len_trim(buff))
+         if(lwrite) write(lunp,'(a)') buff(1:len_trim(buff))
+c         if(lwrite) write(*,'(a50)') 'found in param_card: ',buff
+c         if(.not.lwrite) write(*,'(a50)') 'found in banner: ',buff
+         read(lunr,'(a132)',err=99) buff
+      enddo
+c
+      rewind(lunr)
+c
+      endif
+
+      call setpara
+      call printout
+      call set_mass
+      call set_id
+c
+c
+c
+      write(*,*) '*****************************************************'
+         aa=xran1(iseed)
+         write(*,'(a22,f8.6,a24)') ' *  first rnd number= ',aa    ,
+     &        '                 *'
+      write(*,*) '*****************************************************'
+c
+c     Input the particle to be decayed
+c
+      write(*,*)
+      write(*,*) ' Implemented decays are for:'
+      write(*,*) ' ---------------------------'
+      write(*,*)
+      write(*,*) ' Leptons: ta- ta+' 
+      write(*,*) ' Quarks : t t~' 
+      write(*,*) ' Bosons : z w+ w- h'        
+      write(*,*) ' Input particle to be decayed (e.g. t~):'
+      read(*,'(a3)',err=99)  name
+      call case_trap(name,3)
+c
+c     identify the particle ip
+c
+      ip=0
+      i=-25
+      do while(i.le.25.and.ip.eq.0)
+         if(name.eq.id(i)) ip=i
+         i=i+1
+      enddo
+c
+      if(ip.eq.0) goto 99  ! particle not existing
+c
+c     find out whether the decay requested is implemented.
+c     In case it is, printout the possibilities and
+c     ask the user to choose.
+c
+      call decay_modes(dec_name)
+      if(imode.eq.0) goto 99  ! decay not implemented
+c
+c     setting up branching ratio, widths and so on..
+c
+      call init
+c
+c     write out decay information
+c
+      write(*,*) 
+      write(*,*) '----------------------------------------------'
+      write(*,*) '              Decay information               '
+      write(*,*) '----------------------------------------------'
+      write(*,*) ' '
+      write(*,*) '  particle name   : ',id(ip)
+      write(*,*) '  decay mode      : ',dec_name      
+c
+c     calculate the MC partial width for normalization purposes
+c
+      call tot_decay
+c
+      if(run_mode.eq.0) goto 98  !If only decay width is needed, I'm done 
+c
+c     first check that the particle to be decayed
+c     is present in the events and gather some info
+c     on the events
+c
+
+      write(*,*) ' '
+      write(*,*) '  Wait....Reading In Event File '
+      write(*,*) ' '      
+      done=.false.
+      nevent=0
+      sizeievent=0
+      do i=0,maxievent
+         sum(i)=0d0
+         maxwgt(i)=-1d0
+      enddo
+      do i=0,MaxParticles
+         nfound(i)=0
+      enddo
+      do while (.not. done)
+	 call read_event(lunr,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
+	 if(.not.done) then 
+            call find_particle(ip,nexternal,ic,k,idecay)         
+            nfound(k)=nfound(k)+1
+            eventnumber=eventnr(ievent)
+            sum(eventnumber)=sum(eventnumber)+wgt
+            maxwgt(eventnumber) = max(wgt,maxwgt(eventnumber))
+            nevent=nevent+1
+         endif
+      enddo
+      do i=1, sizeievent
+         sum(0)=sum(0)+sum(i)
+         maxwgt(0)=maxwgt(0)+maxwgt(i)
+      enddo
+
+
+      rewind(lunr)
+
+      write(*,*) '----------------------------------------------'
+      write(*,*) '        Input events information              '
+      write(*,*) '----------------------------------------------'
+      write(*,*) ' '
+      write(*,*) '  Input Event file             :  ',infile
+      write(*,*) '  Number of Events             :  ',nevent
+      write(*,*) '  Integrated weight (pb)       :  ',sum(0)
+      write(*,*) '  Max wgt                      :  ',maxwgt(0)
+      write(*,*) '  Average wgt                  :  ',sum(0)/nevent
+      write(*,*)
+
+      if(nfound(0).eq.nevent) then
+         write(*,*) '  There are no events with particle ',
+     &        id(ip),' to decay !!'
+         stop
+      else
+         do i=0,MaxParticles
+            if(nfound(i).gt.0) then
+      write(*,'(1x,a14,i2,a1,a2,a13,i10)') 
+     & '  Events with ',i,' ',id(ip),'           :',nfound(i)
+            endif   
+         enddo                       
+      endif
+
+      
+c
+c     read each event, decay it, and write it out
+c     note that it is NOT assumed that all the events
+c     have the same number of particles in the final
+c     state or the same type. Each event is treated
+c     independently from the others.
+c
+      write(*,*) ' '
+      write(*,*) '----------------------------------------------'
+      write(*,*) '   Decay events in file : in progress         '      
+      write(*,*) '----------------------------------------------'
+      write(*,*) ' '
+      done=.false.
+      i   = 0
+      do while (.not. done)
+         call read_event(lunr,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
+         if (.not. done) then
+            i=i+1
+            call decay_event(P,wgt,nexternal,ic)
+            call write_event(lunw,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2)
+         endif
+         iten=int(real(nevent)/10e0)+1
+         if(real(i)/real(iten).eq.real(i/iten)) then
+         write(*,'(10x,i3,a)')int(real(i)/real(nevent)*100),' % done'
+         endif
+      enddo
+      write(*,*) '           100% done'
+c
+c     decayed event information
+c
+c-- statistics of the decay_event routine
+c
+      write(*,*) 
+      write(*,*) '----------------------------------------------'
+      write(*,*) '   Statistics of the unweighting              '      
+      write(*,*) '----------------------------------------------'
+      write(*,*) ' '
+      write(*,*) '  Events weighted   = ' ,n_wei
+      write(*,*) '  Events unweighted = ' ,nevent
+      write(*,*) '  Width  (wgt)      = ' ,s_wei/real(n_wei), ' GeV'
+      write(*,*) '  Width  (unwgt)    = ' ,s_unw/real(n_wei), ' GeV'
+      IF(DEBUG) then
+       write(*,*) '  Events   over     = ' ,n_ove
+       write(*,*) '  Integral over     = ' ,s_ove/real(n_wei), ' GeV'
+      ENDIF
+      write(*,*) ' '
+c
+      write(*,*) ' '
+      write(*,*) '  Wait....Writing Out Decayed Event file '
+      write(*,*) ' '      
+      rewind(lunw)
+      done=.false.
+      nevent=0
+      do i=0,sizeievent
+         sum(i)=0d0
+         maxwgt(i)=-1d0
+      enddo
+      do while (.not. done)
+         call read_event(lunw,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
+         if(.not.done) then
+            eventnumber=eventnr(ievent)
+            sum(eventnumber)=sum(eventnumber)+wgt
+            maxwgt(eventnumber) = max(wgt,maxwgt(eventnumber))
+            nevent=nevent+1
+         endif
+      enddo
+      sum(0)=0d0
+      maxwgt(0)=-1d0
+      do i=1, sizeievent
+         sum(0)=sum(0)+sum(i)
+         maxwgt(0)=maxwgt(0)+maxwgt(i)
+      enddo
+
+c RF: To write the correct iseed in the banner, read it again form file iseed.dat
+      open(lunrnd,file='./iseed.dat',status="old",err=300)
+      read(lunrnd,err=300,end=300,fmt='(i10)') iseed
+      iseed=iseed+1
+      close(lunrnd)
+ 300  continue
+cend RF
+
+
+      write(*,*) '----------------------------------------------'
+      write(*,*) '        Output events information             '
+      write(*,*) '----------------------------------------------'
+      write(*,*) ' '
+      write(*,*) '  Output Event file            :  ',outfile
+      write(*,*) '  Number of Events             :  ',nevent
+      write(*,*) '  Integrated weight (pb)       :  ',sum(0)
+      write(*,*) '  Max wgt                      :  ',maxwgt(0)
+      write(*,*) '  Average wgt                  :  ',sum(0)/nevent
+      write(*,*) ' '
+c
+      if(newbanner) then
+         write(luni,'(a)') "<MGDecayInfo>"
+      else
+         write(luni,'(a)') "#********************************************************************"
+         write(luni,'(a)') '#    particle decayed               '
+         write(luni,'(a)') "#********************************************************************"
+      endif
+      write(luni,'(a,a3)')  '#  particle name   : ',id(ip)
+      write(luni,'(a,a)')   '#  decay mode      : ',dec_name      
+      write(luni,'(a,e10.5)')   '#  MC partial width: ',MC_width
+      write(luni,'(a,i10  )') '#  Rnd seed              : ',iseed
+      write(luni,'(a,i10  )') '#  Number of Events      : ',nevent
+      write(luni,'(a,e10.5)') '#  Integrated weight (pb): ',sum(0)
+      write(luni,'(a,e10.5)') '#  Max wgt               : ',maxwgt(0)
+      write(luni,'(a,e10.5)') '#  Average wgt           : ',sum(0)/nevent
+      if(newbanner) then
+         write(luni,'(a)') '</MGDecayInfo>'
+         write(luni,'(a)') '</header>'
+      else
+         write(luni,'(a)') '#********************************************************************' 
+         write(luni,'(a)') '-->'
+      endif
+
+
+C   Write out compulsory init info
+      write(luni,'(a)') '<init>'
+
+      rewind(lunr)
+      do while(index(buff,'</init') .eq. 0)
+         read(lunr,'(a132)',err=99) buff
+         if (index(buff,'<init') .eq. 0) cycle
+         read(lunr,*) (idbmup(i),i=1,2),(ebmup(i),i=1,2),(pdfgup(i),i=1,2),
+     $        (pdfsup(i),i=1,2),idwtup,nprup
+         write(luni,90) (idbmup(i),i=1,2),(ebmup(i),i=1,2),(pdfgup(i),i=1,2),
+     $        (pdfsup(i),i=1,2),idwtup,nprup
+         do i=1,sizeievent
+            read(lunr,*) xsecup,xerrup,xmaxup,lprup 
+            write(luni,91) sum(i),xerrup*sum(i)/xsecup,maxwgt(i),lprup
+         enddo
+      enddo
+      write(luni,'(a)') '</init>'
+
+ 90   FORMAT(2i6,2e19.11,2i2,2i6,i2,i3)
+ 91   FORMAT(3e19.11,i4)
+
+      rewind(lunw)
+      i=0
+      done = .false.
+      do while (.not. done)
+         call read_event(lunw,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
+
+      if (.not. done) then
+         call write_event(luni,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2)
+         
+      endif
+         i=i+1
+      enddo
+
+
+ 98   continue
+
+      write(luni,'(a)')'</LesHouchesEvents>'
+
+c
+c     close datafiles
+c
+      close(luni)
+      close(lunr)
+      close(lunw)
+
+      stop
+
+ 99   write(*,*) 'error'
+      stop
+
+      end 
+
+
+      
+      integer function eventnr(number)
+c*************************************************************************
+c     function that should be given a number and returns
+c     the number:
+c     1   for the first number (i.e. the first time this
+c         function gets executed)
+c     2   for the second number if the second number is
+c         diffent from the first, otherwise it returns 1
+c     3   for the third number, if the third number is
+c         different from the first or second, otherwise
+c         it returns 1, or 2 depending if the third number
+c         is equal to the first or the second.
+c     4   etc.
+c     It updates also a common block with the total number
+c     of different values used as the argument of this function
+c
+c     before the first call to this function sizeievent should
+c     should be set to 0.
+c*************************************************************************
+      include "decay.inc"
+      integer i,number,test(maxievent)
+      logical exist
+      
+      integer sizeievent
+      common/ievents/sizeievent
+
+      save test
+
+      exist = .false.
+      if (sizeievent.lt.1)goto 188
+      do i=1, sizeievent
+         if (number.eq.test(i))then
+            exist=.true.
+            eventnr=i
+         endif
+      enddo
+
+ 188  continue
+      if (.not.exist)then
+         sizeievent=sizeievent+1
+         eventnr=sizeievent
+         test(sizeievent)=number
+      endif
+
+      if (sizeievent.ge.maxievent)then
+         write (*,*)
+     &   'Error, too many different event numbers in event file'
+         return
+      endif
+
+      return
+      end
+
+
+
+
+
+
+      subroutine set_id
+c****************************************
+c     to set the particles ids
+*****************************************
+      implicit none
+      include 'decay.inc'
+c
+c
+c-Quarks: d u s c b t d~ u~ s~ c~ b~ t~ 
+      id(1) ='d'
+      id(2) ='u'
+      id(3) ='s'
+      id(4) ='c'
+      id(5) ='b'
+      id(6) ='t'
+      id(-1)='d~'
+      id(-2)='u~'
+      id(-3)='s~'
+      id(-4)='c~'
+      id(-5)='b~'
+      id(-6)='t~'
+c-Leptons: e- mu- ta- ve vm vt e+ mu+ ta+ ve~ vm~ vt~ 
+      id(11) ='e-'
+      id(12) ='ve'
+      id(13) ='mu-'
+      id(14) ='vm'
+      id(15) ='ta-'
+      id(16) ='vt'
+      id(-11)='e+'
+      id(-12)='ve~'
+      id(-13)='mu+'
+      id(-14)='vm~'
+      id(-15)='ta+'
+      id(-16)='vt~'
+c-Bosons: g a z w+ w- h        
+      id(21)  ='g'
+      id(22)  ='a'
+      id(23)  ='z'
+      id(24)  ='w+'
+      id(-24) ='w-'
+      id(25)  ='h'
+      return
+      end
+ 
+
+      subroutine set_mass
+c****************************************
+c     to set the particles ids
+*****************************************
+      implicit none
+      include 'decay.inc'
+      include 'coupl.inc'
+c
+c
+c-Quarks: d u s c b t d~ u~ s~ c~ b~ t~ 
+      pdgmass(1) =0d0
+      pdgmass(2) =0d0
+      pdgmass(3) =0d0
+      pdgmass(4) =cmass
+      pdgmass(5) =bmass
+      pdgmass(6) =tmass
+      pdgmass(-1)=0d0
+      pdgmass(-2)=0d0
+      pdgmass(-3)=0d0
+      pdgmass(-4)=cmass
+      pdgmass(-5)=bmass
+      pdgmass(-6)=tmass
+c-Leptons: e- mu- ta- ve vm vt e+ mu+ ta+ ve~ vm~ vt~ 
+      pdgmass(11) =0d0
+      pdgmass(12) =0d0
+      pdgmass(13) =0d0
+      pdgmass(14) =0d0
+      pdgmass(15) =lmass
+      pdgmass(16) =0d0
+      pdgmass(-11)=0d0
+      pdgmass(-12)=0d0
+      pdgmass(-13)=0d0
+      pdgmass(-14)=0d0
+      pdgmass(-15)=lmass
+      pdgmass(-16)=0d0
+c-Bosons: g a z w+ w- h        
+      pdgmass(21)  =0d0
+      pdgmass(22)  =0d0
+      pdgmass(23)  =zmass
+      pdgmass(24)  =wmass
+      pdgmass(-24) =wmass
+      pdgmass(25)  =hmass
+      return
+      end
+      
+      subroutine decay_modes(dec_name)
+c*********************************************
+c     to set the particles ids
+c*********************************************
+      implicit none
+      include 'decay.inc'
+c
+c     argument
+c
+      character*50 dec_name
+c
+c     common
+c
+      include 'coupl.inc'
+c
+c     local
+c
+      character*50 string(30)
+      integer i
+c----------
+c     START
+c----------
+
+      write(*,*) 'particle to decay :',id(ip)
+      write(*,*) 
+      write(*,*) 'Implemented decay modes:'
+      write(*,*) '------------------------'
+c
+      if(ip.eq.6) then  !top
+         string(1) = ' t -> b  w+ '
+         string(2) = ' t -> b  ve e+'
+         string(3) = ' t -> b  vm mu+'
+         string(4) = ' t -> b  vt ta+'
+         string(5) = ' t -> b  vl l+    (e+mu)'
+         string(6) = ' t -> b  vl l+    (e+mu+ta)'
+         string(7) = ' t -> b  j   j    (ud+cs)'
+         string(8) = ' t -> b  anything (e+mu+ta+ud+cs)'
+
+         write(*,'(i2,2x,a40)') (i,string(i),i=1,8)
+         write(*,*) ' '
+         write(*,*) 'your choice is:'   
+         read (*,*) imode
+         if(.not.(imode.ge.1.and.imode.le.8)) then
+            write(*,*) 'choice not implemented'
+            imode=0
+            return
+         endif
+         if(imode.eq.1) then !number of particles in the decay
+            nd=2
+         else
+            nd=3
+         endif
+         dec_name=string(imode)
+         m1=tmass
+         return
+      endif
+c
+      if(ip.eq.-6) then  !anti-top
+         string(1) =  ' t~ -> b~  w- '
+         string(2) =  ' t~ -> b~  e- ve~'
+         string(3) =  ' t~ -> b~ mu- vm~'
+         string(4) =  ' t~ -> b~ ta- vt~'
+         string(5) =  ' t~ -> b~  l- vl~  (e+mu)'
+         string(6) =  ' t~ -> b~  l- vl~  (e+mu+ta)'
+         string(7) =  ' t~ -> b~  j   j   (ud+cs)'
+         string(8) =  ' t~ -> b~ anything (e+mu+ta+ud+cs)'
+         write(*,'(i2,2x,a40)') (i,string(i),i=1,8)
+         write(*,*) ' '
+         write(*,*) 'your choice is:'   
+         read (*,*) imode
+         if(.not.(imode.ge.1.and.imode.le.8)) then
+            write(*,*) 'choice not implemented'
+            imode=0
+            return
+         endif
+         if(imode.eq.1) then    !number of particles in the decay 
+            nd=2
+         else
+            nd=3
+         endif
+         m1=tmass
+         dec_name=string(imode)
+         return
+      endif
+c
+      if(ip.eq.15) then  ! tau-
+         string(1) =  ' ta- -> vt   e- ve~  '
+         string(2) =  ' ta- -> vt  mu- vm~  '
+         string(3) =  ' ta- -> vt   l- vl~ (e+mu) '
+         string(4) =  ' ta- -> vt   pi-     '
+         string(5) =  ' ta- -> vt  rho(770)-'
+         write(*,'(i2,2x,a30)') (i,string(i),i=1,5)
+         write(*,*) ' '
+         write(*,*) 'your choice is:'   
+         read (*,*) imode
+         if(.not.(imode.ge.1.and.imode.le.5)) then
+            write(*,*) 'choice not implemented'
+            imode=0
+            return
+         endif
+         if(imode.le.3) then    !number of particles in the decay
+            nd=3
+         else
+            nd=2
+         endif
+         m1=lmass
+         dec_name=string(imode)
+         return
+      endif
+c
+      if(ip.eq.-15) then  ! tau+
+         string(1) =  ' ta+ -> vt~  ve e+   '
+         string(2) =  ' ta+ -> vt~  vm mu+   '
+         string(3) =  ' ta+ -> vt~  vl  l+  (e+mu)'
+         string(4) =  ' ta+ -> vt~   pi+     '
+         string(5) =  ' ta+ -> vt~  rho(770)+'
+         write(*,'(i2,2x,a30)') (i,string(i),i=1,5)
+         write(*,*) ' '
+         write(*,*) 'your choice is:'   
+         read (*,*) imode
+         if(.not.(imode.ge.1.and.imode.le.5)) then
+            write(*,*) 'choice not implemented'
+            imode=0
+         endif
+         if(imode.le.3) then !number of particles in the decay
+            nd=3
+         else
+            nd=2
+         endif
+         m1=lmass
+         dec_name=string(imode)
+         return
+      endif
+c
+      if(ip.eq.23) then  ! z
+         string(1) =  ' z -> e-   e+ '
+         string(2) =  ' z -> mu- mu+ '
+         string(3) =  ' z -> ta- ta+ '
+         string(4) =  ' z -> l-   l+  (e+mu)'
+         string(5) =  ' z -> l-   l+  (e+mu+ta )'
+         string(6) =  ' z -> vl  vl~  (ve+vm+vt)'
+         string(7) =  ' z -> b    b~            '
+         string(8) =  ' z -> c    c~            '
+         string(9) =  ' z -> j    j~  (u+d+c+s) '
+         string(10)=  ' z -> j    j~  (u+d+c+s+b)'
+         string(11)=  ' z -> visible  (e+mu+ta+u+d+c+s+b)'
+         string(12)=  ' z -> anything'
+
+         write(*,'(i2,2x,a40)') (i,string(i),i=1,10)
+         write(*,*) ' '
+         write(*,*) 'your choice is:'   
+         read (*,*) imode
+         if(.not.(imode.ge.1.and.imode.le.12)) then
+            write(*,*) 'choice not implemented'
+            imode=0
+            return
+         endif
+         m1=zmass
+         nd=2                   !number of particles in the decay
+         dec_name=string(imode)
+         return
+      endif
+c
+      if(ip.eq.24) then  ! w+
+         string(1) =  ' w+ -> e+  ve'
+         string(2) =  ' w+ -> mu+ vm'
+         string(3) =  ' w+ -> ta+ vt '
+         string(4) =  ' w+ -> l+  vl (e+mu)'
+         string(5) =  ' w+ -> l+  vl (e+mu+ta)'
+         string(6) =  ' w+ -> j   j  (ud+cs)'
+         string(7) =  ' w+ -> anything (e+mu+ta+ud+cs)'
+
+         write(*,'(i2,2x,a40)') (i,string(i),i=1,7)
+         write(*,*) 'your choice is:'   
+         write(*,*) ' '
+         read (*,*) imode
+         if(.not.(imode.ge.1.and.imode.le.7)) then
+            write(*,*) 'choice not implemented'
+            imode=0
+         endif
+         nd=2                   !number of particles in the decay
+         m1=wmass
+         dec_name=string(imode)
+         return
+      endif
+c
+      if(ip.eq.-24) then  ! w-
+         string(1) =  ' w- -> e-  ve~'
+         string(2) =  ' w- -> mu- vm~'
+         string(3) =  ' w- -> ta- vt~ '
+         string(4) =  ' w- -> l-  vl~ (e+mu)'
+         string(5) =  ' w- -> l-  vl~ (e+mu+ta)'
+         string(6) =  ' w- -> j   j   (ud+cs)'
+         string(7) =  ' w- -> anything (e+mu+ta+ud+cs)'
+         write(*,'(i2,2x,a40)') (i,string(i),i=1,7)
+         write(*,*) 'your choice is:'   
+         write(*,*) ' '
+         read (*,*) imode
+         if(.not.(imode.ge.1.and.imode.le.7)) then
+            write(*,*) 'choice not implemented'
+            imode=0
+            return
+         endif
+         m1=wmass
+         nd=2                   !number of particles in the decay
+         dec_name=string(imode)
+         return
+      endif
+c
+      if(ip.eq.25) then  ! h
+        string(1) =' h -> b   b~ '
+        string(2) =' h -> ta- ta+'
+        string(3) =' h -> mu- mu+'
+        string(4) =' h -> c   c~ '
+        string(5) =' h -> t   t~ (when m_h>2*m_t)'
+        string(6) =' h -> g   g  '
+        string(7) =' h -> a   a  '
+        string(8) =' h -> z   a  (when m_h>  m_z)'
+        string(9) =' h -> w+  w- (when m_h>2*m_w)'
+        string(10)=" h -> w*  w -> l  vl  l' vl' (l,l'=e,mu)" 
+        string(11)=" h -> w*  w -> l  vl  l' vl' (l,l'=e,mu,ta)" 
+        string(12)=" h -> w*  w -> j  j   l  vl  (jj=ud,cs;l=e,mu)"
+        string(13)=" h -> w*  w -> j  j   l  vl  (jj=ud,cs;l=e,mu,ta)"
+        string(14)=' h -> z   z  (when m_h>2*m_z)'
+        string(15)=" h -> z*  z -> l- l+  l-' l+'(l,l'=e,mu)"
+        string(16)=" h -> z*  z -> l- l+  l-' l+'(l,l'=e,mu,ta)"
+        string(17)=" h -> z*  z -> j  j~  l-  l+ (j=u,d,c,s;l=e,mu )"
+        string(18)=" h -> z*  z -> j  j~  l-  l+ (j=u,d,c,s;l=e,mu,ta)"
+        string(19)=" h -> z*  z -> b  b~  l-  l+ (l=e,mu )"
+        string(20)=" h -> z*  z -> b  b~  l-  l+ (l=e,mu,ta )"
+        string(21)=" h -> z*  z -> vl vl~ l-' l+'(l=e,mu,ta;l'=e,mu) "
+        string(22)=" h -> z*  z -> vl vl~ l-' l+'(l,l'=e,mu,ta) "
+
+         write(*,'(i2,2x,a50)') (i,string(i),i=1,22)
+         write(*,*) 'your choice is:'   
+         write(*,*) ' '
+         read (*,*) imode
+         if(.not.(imode.ge.1.and.imode.le.22)) then
+            write(*,*) 'choice not implemented'
+            imode=0
+            return
+         endif
+c--   number of decay products
+         ND=4
+         m1=hmass
+         if(imode.le.9.or.imode.eq.14) ND=2
+         dec_name=string(imode)
+         return
+      endif
+      
+      write(*,*) 'no decay mode implemented for particle ',id(ip)
+      
+
+      return
+      end
+      
+
+      double precision function dot(p1,p2)
+C****************************************************************************
+C     4-Vector Dot product
+C****************************************************************************
+      implicit none
+      double precision p1(0:3),p2(0:3)
+      dot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
+      end
+
+
+      subroutine init
+c********************************************************
+c     sets up : the branching ratios (best=HDECAY or PDG)
+c               LO partial widths 
+c               ND: the number of decay products
+c********************************************************
+      implicit none
+c
+c     Include
+c
+      include 'decay.inc'
+      include 'coupl.inc'
+      include 'calc_values.inc'
+c
+c-----
+cstart
+c-----
+c      
+
+c
+c     PDG2002 values for the Branching ratios:
+c
+
+c--   tau
+      br_ta_lv =0.175d0
+      br_ta_pi =0.111d0
+      br_ta_ro =0.254d0
+c--   w
+      br_w_lv =0.1068d0
+      br_w_jj =1d0-3d0*br_w_lv
+c--   z
+      br_z_ll =0.0336d0
+      br_z_vv =0.2000d0
+      br_z_cc =0.1176d0
+      br_z_bb =0.1514d0
+      br_z_jj =0.6991d0-br_z_cc-br_z_bb
+      
+      write(*,*)
+      write(*,*) '----------------------------------------------'
+      write(*,*) '         Relevant branching ratios            '
+      write(*,*) '----------------------------------------------'
+      write(*,*)
+
+c-------------------------------------------------------------------
+      if(abs(ip).eq.6) then     !top
+
+
+         write(*,16) 'BR (t -> b w    ) =  ', 1d0
+         bratio=1d0
+         if(imode.eq.2) then
+         write(*,16) 'BR (w -> e ve   ) =  ', br_w_lv
+         bratio=br_w_lv
+         elseif(imode.eq.3) then
+         write(*,16) 'BR (w -> mu vm  ) =  ', br_w_lv
+         bratio=br_w_lv
+         elseif(imode.eq.4) then
+         write(*,16) 'BR (w -> tau vt ) =  ', br_w_lv
+         bratio=br_w_lv
+         elseif(imode.eq.5) then
+         write(*,15) 'BR (w -> l vl   ) =  ', 2d0*br_w_lv,' (l=e,mu)'
+         bratio=2d0*br_w_lv
+         elseif(imode.eq.6) then
+         write(*,15) 'BR (w -> l vl   ) =  ', 3d0*br_w_lv,' (l=e,mu,ta)'
+         bratio=3d0*br_w_lv
+         elseif(imode.eq.7) then
+         write(*,15) 'BR (w -> j j    ) =  ', br_w_jj   ,' (jj=ud+cs)'
+         bratio=br_w_jj
+         elseif(imode.eq.8) then
+         write(*,15) 'BR (w ->anything) =  ',1d0
+         bratio=1d0
+         endif
+         
+         If(imode.eq.1) then    !t -> b w+
+            calc_width=twidth
+         elseif(imode.eq.2.or.imode.eq.3.or.imode.eq.5) then !  t -> b vl l+
+            calc_width=twidth*(w_w_nl/wwidth)
+            if(imode.eq.5) calc_width=calc_width*2d0
+         elseif(imode.eq.4) then !     t -> b vt tau+
+            calc_width=twidth*(w_w_tau/wwidth)
+         elseif(imode.eq.6) then !     t -> b vl l+
+            calc_width=twidth*((2d0*w_w_nl+w_w_tau)/wwidth)
+         elseif(imode.eq.7) then !    t -> b j j 
+            calc_width=twidth*(w_w_ud+w_w_cs)/wwidth
+         elseif(imode.eq.8) then !    t -> b anything 
+            calc_width=twidth
+         endif
+                  
+      endif
+c-------------------------------------------------------------------
+      if(abs(ip).eq.15) then    !tau
+         if(imode.eq.1)     then
+            write(*,16) 'BR (ta -> vt  e ve) =  ', br_ta_lv
+         bratio=br_ta_lv
+         elseif(imode.eq.2)     then
+            write(*,16) 'BR (ta -> vt mu vm) =  ', br_ta_lv
+         bratio=br_ta_lv
+         elseif(imode.eq.3)     then
+        write(*,15) 'BR (ta -> vt  l vl) =  ', 2d0*br_ta_lv,' (l=e,mu)'
+         bratio=2d0*br_ta_lv
+         elseif(imode.eq.4) then
+            write(*,16) 'BR (ta -> vt  pi  ) =  ', br_ta_pi
+         bratio=br_ta_pi
+         elseif(imode.eq.5) then
+            write(*,16) 'BR (ta -> vt  rho ) =  ', br_ta_ro
+         bratio=br_ta_ro
+         endif
+         
+         if(imode.eq.1.or.imode.eq.2) then !ta -> vt vl l
+            calc_width=lwidth*br_ta_lv
+         elseIf(imode.eq.3) then
+               calc_width=2d0*lwidth*br_ta_lv 
+         elseIf(imode.eq.5) then !     tau -> vt rho
+            calc_width=lwidth*br_ta_ro
+         elseif(imode.eq.4) then !    tau -> vt pi
+            calc_width=lwidth*br_ta_pi
+         endif
+         
+      endif
+
+c-------------------------------------------------------------------
+      if(ip.eq.23) then         ! z
+         if    (imode.eq.1) then
+            write(*,16) 'BR (z -> e- e+  ) =  ', br_z_ll
+         bratio=br_z_ll
+         elseif(imode.eq.2) then
+            write(*,16) 'BR (z -> mu- mu+) =  ', br_z_ll
+         bratio=br_z_ll
+         elseif(imode.eq.3) then
+            write(*,16) 'BR (z -> ta- ta+) =  ', br_z_ll
+         bratio=br_z_ll
+         elseif(imode.eq.4) then
+        write(*,15) 'BR (z -> l- l+  ) =  ', 2d0*br_z_ll,' (l=e,mu)'
+         bratio=2d0*br_z_ll
+         elseif(imode.eq.5) then
+        write(*,15) 'BR (z -> l- l+  ) =  ', 3d0*br_z_ll,' (l=e,mu,ta)'
+         bratio=3d0*br_z_ll
+         elseif(imode.eq.6) then
+            write(*,15) 'BR (z -> vl vl~ ) =  ', br_z_vv,' (l=e,mu,ta)'
+         bratio=br_z_vv
+         elseif(imode.eq.7) then
+             write(*,16) 'BR (z -> b  b~  ) =  ', br_z_bb
+         bratio=br_z_bb
+         elseif(imode.eq.8)then
+            write(*,16) 'BR (z -> c  c~  ) =  ', br_z_cc
+         bratio=br_z_cc
+         elseif(imode.eq.9)then
+            write(*,15) 'BR (z -> j  j~  ) =  ', br_z_jj+br_z_cc,
+     &           ' (j=u+d+c+s)'
+         bratio=br_z_jj+br_z_cc
+         elseif(imode.eq.10)then
+            write(*,15) 'BR (z -> j  j~  ) =  ', 
+     &           br_z_jj+br_z_bb+br_z_cc,' (j=u+d+c+s+b)'
+         bratio=br_z_jj+br_z_bb+br_z_cc
+         elseif(imode.eq.11)then
+            write(*,15) 'BR (z -> visible) =  ', 
+     &           br_z_jj+br_z_bb+br_z_cc+3d0*br_z_ll,
+     &           ' (j=l+j)'
+         bratio=br_z_jj+br_z_bb+br_z_cc+3d0*br_z_ll
+         elseif(imode.eq.12)then
+            write(*,15) 'BR (z -> anyth. ) =  ', 
+     &           1d0,' (j=l+v+j)'
+         bratio=1d0
+         endif
+         
+         
+         If(imode.eq.1.or.imode.eq.2) then !     z->e- e+
+            calc_width=w_z_ll
+         elseif(imode.eq.4) then !     z->e- e+,mu-mu+,ta-ta+
+            calc_width=2d0*w_z_ll
+         elseif(imode.eq.3) then !     z->ta- ta+
+            calc_width=w_z_tau
+         elseif(imode.eq.5) then !     z->e- e+,mu-mu+,ta-ta+
+            calc_width=2d0*w_z_ll+w_z_tau
+         elseif(imode.eq.6) then !    z-> vl  vl~
+            calc_width=w_z_nn*3d0
+         elseif(imode.eq.7) then !    z->b b~
+            calc_width=w_z_bb
+         elseif(imode.eq.8) then !     z->c c~
+            calc_width=w_z_cc
+         elseif(imode.eq.9) then !     z->jj (u+d+c+s)
+            calc_width=Two*w_z_dd + w_z_uu + w_z_cc 
+         elseif(imode.eq.10) then !*     z->jj (u+d+c+s+b)
+            calc_width=Two*w_z_dd + w_z_uu + w_z_cc + w_z_bb 
+         elseif(imode.eq.11) then !*     z->any (l+u+d+c+s+b)
+            calc_width=zwidth-w_z_nn*3d0
+         elseif(imode.eq.12) then !*     total width
+            calc_width=zwidth
+         endif
+      endif
+      
+c-------------------------------------------------------------------
+      if(abs(ip).eq.24) then    ! w
+         if    (imode.eq.1) then
+         write(*,16) 'BR (w -> e ve   ) =  ', br_w_lv
+         bratio=br_w_lv
+         elseif(imode.eq.2) then
+         write(*,16) 'BR (w -> mu vm  ) =  ', br_w_lv
+         bratio=br_w_lv
+         elseif(imode.eq.3) then
+         write(*,16) 'BR (w -> tau vt ) =  ', br_w_lv
+         bratio=br_w_lv
+         elseif(imode.eq.4) then
+         write(*,15) 'BR (w -> l vl   ) =  ', 2d0*br_w_lv,' (l=e,mu)'
+         bratio=2d0*br_w_lv
+         elseif(imode.eq.5) then
+         write(*,15) 'BR (w -> l vl   ) =  ', 3d0*br_w_lv,' (l=e,mu,ta)'
+         bratio=3d0*br_w_lv
+         elseif(imode.eq.6) then
+         write(*,15) 'BR (w -> j j    ) =  ', br_w_jj   ,' (jj=ud+cs)'
+         bratio=br_w_jj
+         elseif(imode.eq.7) then
+         write(*,15) 'BR (w ->anything) =  1  (=e+mu+ta+ud+cs)'
+         bratio=1d0
+         endif
+         
+         If(imode.eq.1.or.imode.eq.2) then !w-> l ve
+            calc_width=w_w_nl
+         elseif(imode.eq.3) then !     w -> ta vt
+            calc_width=w_w_tau
+         elseif(imode.eq.4) then !     w -> l vl (e,mu)
+            calc_width=2d0*w_w_nl
+         elseif(imode.eq.5) then !    w -> l vl (e,mu,tau)
+            calc_width=2d0*w_w_nl+w_w_tau
+         elseif(imode.eq.6) then !     w -> ud+cs
+            calc_width=w_w_ud+w_w_cs
+         elseif(imode.eq.7) then !     w -> anything
+            calc_width=wwidth
+         endif
+
+      endif
+
+c-------------------------------------------------------------------
+      if(ip.eq.25) then         ! higgs
+         
+         write(*,15) 'Higgs mass        =  ', hmass,  ' GeV'
+         write(*,15) 'Higgs tot width   =  ', SMWDTH, ' GeV'
+         if(imode.eq.1) then
+            write(*,16) 'BR (h -> b b~   ) =  ', SMBRB
+            bratio=SMBRB
+         elseif(imode.eq.2) then
+            write(*,16) 'BR (h -> ta+ ta-) =  ', SMBRL
+            bratio=SMBRL
+         elseif(imode.eq.3) then
+            write(*,16) 'BR (h -> mu+ mu-) =  ', SMBRM
+            bratio=SMBRM
+         elseif(imode.eq.4) then 
+            write(*,16) 'BR (h -> c  c~  ) =  ', SMBRC
+            bratio=SMBRC
+         elseif(imode.eq.5) then 
+            write(*,16) 'BR (h -> t  t~  ) =  ', SMBRT
+            bratio=SMBRT
+         elseif(imode.eq.6) then 
+            write(*,16) 'BR (h -> g  g   ) =  ', SMBRG
+            bratio=SMBRG
+         elseif(imode.eq.7) then 
+            write(*,16) 'BR (h -> a  a   ) =  ', SMBRGA
+            bratio=SMBRGA
+         elseif(imode.eq.8) then 
+            write(*,16) 'BR (h -> z  a   ) =  ', SMBRZGA
+            bratio=SMBRZGA
+         elseif(imode.ge.9.and.imode.le.13) then 
+            write(*,16) 'BR (h -> w+ w-  ) =  ', SMBRW
+            bratio=SMBRW
+         elseif(imode.ge.14.and.imode.le.22) then 
+            write(*,16) 'BR (h -> z  z   ) =  ', SMBRZ    
+            bratio=SMBRZ    
+         endif
+
+         if(imode.eq.10) then
+         write(*,15) 'BR (w -> l vl   ) =  ', 2d0*br_w_lv,' (l=e,mu)'
+            bratio=bratio*(2d0*br_w_lv)**2
+         elseif(imode.eq.11) then
+         write(*,15) 'BR (w -> l vl   ) =  ', 3d0*br_w_lv,' (l=e,mu,ta)'
+            bratio=bratio*(3d0*br_w_lv)**2
+         elseif(imode.eq.12) then
+         write(*,15) 'BR (w -> j j    ) =  ', br_w_jj   ,' (jj=ud+cs)'
+         write(*,15) 'BR (w -> l vl   ) =  ', 2d0*br_w_lv,' (l=e,mu)'
+            bratio=bratio*2d0*br_w_jj*2d0*br_w_lv
+         elseif(imode.eq.13) then
+         write(*,15) 'BR (w -> j j    ) =  ', br_w_jj   ,' (jj=ud+cs)'
+         write(*,15) 'BR (w -> l vl   ) =  ', 3d0*br_w_lv,' (l=e,mu,ta)'
+            bratio=bratio*2d0*br_w_jj*3d0*br_w_lv
+         endif
+
+         if(imode.eq.15) then
+        write(*,15) 'BR (z -> l- l+  ) =  ', 2d0*br_z_ll,' (l=e,mu)'
+            bratio=bratio*(2d0*br_z_ll)**2
+         elseif(imode.eq.16) then
+        write(*,15) 'BR (z -> l- l+  ) =  ', 3d0*br_z_ll,' (l=e,mu,ta)'
+            bratio=bratio*(3d0*br_z_ll)**2
+         elseif(imode.eq.17) then
+        write(*,15) 'BR (z -> j  j~  ) =  ', br_z_jj+br_z_cc,
+     &           ' (j=u+d+c+s)'
+        write(*,15) 'BR (z -> l- l+  ) =  ', 2d0*br_z_ll,' (l=e,mu)'
+            bratio=bratio*(br_z_jj+br_z_cc)*2d0*br_z_ll*2d0
+         elseif(imode.eq.18) then
+        write(*,15) 'BR (z -> j  j~  ) =  ', br_z_jj+br_z_cc,
+     &           ' (j=u+d+c+s)'
+        write(*,15) 'BR (z -> l- l+  ) =  ', 3d0*br_z_ll,' (l=e,mu,ta)'
+            bratio=bratio*(br_z_jj+br_z_cc)*3d0*br_z_ll*2d0
+         elseif(imode.eq.19) then
+        write(*,16) 'BR (z -> b  b~  ) =  ', br_z_bb
+        write(*,15) 'BR (z -> l- l+  ) =  ', 2d0*br_z_ll,' (l=e,mu)'
+            bratio=bratio*br_z_bb*2d0*br_z_ll*2d0
+         elseif(imode.eq.20) then
+        write(*,16) 'BR (z -> b  b~  ) =  ', br_z_bb
+        write(*,15) 'BR (z -> l- l+  ) =  ', 3d0*br_z_ll,' (l=e,mu,ta)'
+            bratio=bratio*br_z_bb*3d0*br_z_ll*2d0
+         elseif(imode.eq.21) then
+        write(*,15) 'BR (z -> vl vl~ ) =  ', br_z_vv,' (l=e,mu,ta)'
+        write(*,15) 'BR (z -> l- l+  ) =  ', 2d0*br_z_ll,' (l=e,mu)'
+            bratio=bratio*br_z_vv*2d0*br_z_ll*2d0
+         elseif(imode.eq.22) then
+        write(*,15) 'BR (z -> vl vl~ ) =  ', br_z_vv,' (l=e,mu,ta)'
+        write(*,15) 'BR (z -> l- l+  ) =  ', 3d0*br_z_ll,' (l=e,mu,ta)'
+            bratio=bratio*br_z_vv*3d0*br_z_ll*2d0
+         endif
+
+
+        write(*,16) 'Best tot Br Ratio =  ', bratio
+
+         If(imode.eq.1) then    !    h->b b~
+            calc_width=w_h_bb
+         elseif(imode.eq.2) then ! h->ta- ta+
+            calc_width=w_h_tau
+         elseif(imode.eq.3) then ! h->mu- mu+
+            calc_width=w_h_tau/(lmass/0.105658389d0)**2
+         elseif(imode.eq.4) then !   h-> c c~
+            calc_width=w_h_cc
+         elseif(imode.eq.5) then !    h-> t t~
+            if(hmass.le.2d0*tmass) then
+               write(*,*) 
+               write(*,*) 'Error: this decay mode assumes m_h>2 m_t'
+               write(*,*) 
+               stop
+            endif
+             calc_width=w_h_tt
+         elseif(imode.eq.6) then ! h-> g g 
+            calc_width=SMBRG*SMWDTH 
+         elseif(imode.eq.7) then !   h-> a a 
+            calc_width= SMBRGA*SMWDTH 
+         elseif(imode.eq.8) then !     h-> z a   
+            if(hmass.le.zmass) then
+               write(*,*) 
+               write(*,*) 'Error: this decay mode assumes m_h>m_Z'
+               write(*,*) 
+               stop
+            endif
+            calc_width= SMBRZGA*SMWDTH 
+         elseif(imode.eq.9) then !    h-> w w
+            if(hmass.le.2d0*wmass) then
+               write(*,*) 
+               write(*,*) 'Error: this decay mode assumes m_h>2*m_W'
+               write(*,*) 
+               stop
+            endif
+            calc_width=w_h_ww
+         elseif(imode.eq.10) then !  h -> w*  w -> e  ve mu vmu 
+            calc_width= SMBRW*SMWDTH*(2d0*w_w_nl/wwidth)**2
+         elseif(imode.eq.11) then ! h -> w*  w -> l  vl  l' vl' (l,l'=e,mu,ta)
+            calc_width= SMBRW*SMWDTH*((2d0*w_w_nl+w_w_tau)/wwidth)**2
+         elseif(imode.eq.12) then !    h -> w*  w -> j  j   l  vl  (jj=ud,cs;l=e,mu)
+        calc_width= SMBRW*SMWDTH*2d0*(2d0*w_w_nl*(w_w_ud+w_w_cs))/wwidth**2
+         elseif(imode.eq.13) then !h -> w*  w -> j  j   l  vl  (jj=ud,cs;l=e,mu,ta)
+            calc_width= SMBRW*SMWDTH*2d0*
+     &           ((2d0*w_w_nl+w_w_tau)*(w_w_ud+w_w_cs))/wwidth**2
+         elseif(imode.eq.14) then ! h-> z z 
+            if(hmass.le.2d0*zmass) then
+               write(*,*) 
+               write(*,*) 'Error: this decay mode assumes m_h>2*m_Z'
+               write(*,*) 
+               stop
+            endif
+            calc_width=w_h_zz
+         elseif(imode.eq.15) then ! h -> z*  z -> l- l+  l-' l+'(l,l'=e,mu)
+            calc_width= SMBRZ*SMWDTH*(2d0*w_z_ll/zwidth)**2
+         elseif(imode.eq.16) then ! h -> z*  z -> l- l+  l-' l+'(l,l'=e,mu,ta)
+            calc_width= SMBRZ*SMWDTH*
+     &                  ((2d0*w_z_ll+w_z_tau)/zwidth)**2
+         elseif(imode.eq.17) then ! h -> z*  z -> j  j~  l-  l+ (j=u,d,c,s;l=e,mu )
+            calc_width= SMBRZ*SMWDTH*2d0*
+     &      (2d0*w_z_ll*(Two*w_z_dd + w_z_uu + w_z_cc))/zwidth**2
+         elseif(imode.eq.18) then ! h -> z*  z -> j  j~  l-  l+ (j=u,d,c,s;l=e,mu,ta)
+            calc_width= SMBRZ*SMWDTH*2d0*
+     &  ((2d0*w_z_ll+w_z_tau)*
+     &  (Two*w_z_dd + w_z_uu + w_z_cc))/zwidth**2
+         elseif(imode.eq.19) then ! h -> z*  z -> b  b~  l-  l+ (l=e,mu )
+            calc_width=SMBRZ*SMWDTH*2d0*
+     &                 w_z_ll*w_z_bb*2d0/zwidth**2
+         elseif(imode.eq.20) then ! h -> z*  z -> b  b~  l-  l+ (l=e,mu,ta )
+            calc_width=SMBRZ*SMWDTH*2d0*
+     &           (2d0*w_z_ll+w_z_tau)*w_z_bb/zwidth**2
+         elseif(imode.eq.21) then ! h -> z*  z -> vl vl~ l-' l+'(l=e,mu,ta;l'=e,mu)
+            calc_width=SMBRZ*SMWDTH*2d0*
+     &                 3d0*w_z_nn*2d0*w_z_ll/zwidth**2
+         elseif(imode.eq.22) then ! h -> z*  z -> vl vl~ l-' l+'(l,l'=e,mu,ta)
+            calc_width=SMBRZ*SMWDTH*2d0*
+     &                 3d0*w_z_nn*(2d0*w_z_ll+w_z_tau)/zwidth**2
+         endif   
+ 
+         calc_br=calc_width/SMWDTH
+         write(*,16) 'Calc Br Ratio     =  ', calc_br
+         
+      endif
+
+
+
+
+
+
+ 15   format( 3x,a,f9.5,a )
+ 16   format( 3x,a,f9.5 )
+
+   
+
+      return
+      end
+
+
+      subroutine tot_decay
+c**********************************************************    
+c     vegas evaluation of the width
+c**********************************************************
+      implicit none
+      include 'decay.inc'
+c
+c     LOCAL
+c
+      INTEGER init,itmx,ncall,nprn,ndim,i
+      REAL*4 region(2*mxdim),fxn,tgral,chi2a,sd
+      real*8 diff,eff,ave_wei
+      LOGICAL WRITEOUT
+      DATA WRITEOUT /.TRUE./
+c
+c     EXTERNAL
+c
+      EXTERNAL FXN 
+c
+c-----
+c     Begin 
+c-----     
+c
+c
+c     write out that calculation is going on
+c
+      write(*,*) 
+      write(*,*) '----------------------------------------------'
+      write(*,*) '      MC Evaluation of the Partial Width      '
+      write(*,*) '----------------------------------------------'
+c
+c-MC
+c
+      ndim  =3*ND-4  !number of dimensions of the phase space
+c
+c-warm up
+c
+      init  =0
+      ncall =100000
+      itmx  =5
+      nprn  =0
+      do i=1,ndim
+         region(i)=0E0
+         region(ndim+i)=1E0
+      enddo
+      call vegas(region,ndim,fxn,init,ncall,itmx,nprn,tgral,sd,chi2a)
+c
+c-final run
+c
+      mxwgt =0d0 !reset the maximum weight to zero
+      ncall =maxpoints
+      itmx  =1   !only one large iteration
+      nprn  =0
+      init  =1
+      call vegas(region,ndim,fxn,init,ncall,itmx,nprn,tgral,sd,chi2a)
+c
+c-result
+c
+      MC_width=tgral
+      mxwgt=mxwgt*real(ncall)
+      eff=MC_width/mxwgt*100e0
+C
+C     write results on the screen
+C      
+      write(*,*) ' '
+      write(*,*) '  LO partial width = ' ,calc_width, ' GeV'
+      write(*,*) '  MC partial width = ' ,MC_width,' +-',sd,' GeV'
+      write(*,*) '  Max. wgt         = ' ,mxwgt
+      write(*,1) '  Unwgt. Eff.(%)   = ' ,eff
+ 1    format(1x,a21,1x,f4.1)
+c
+c     write the result on a log file when WRITEOUT=.true.
+c      
+      IF(WRITEOUT) THEN
+         open(lunout,file='all_decay.dat',status='unknown')
+         call toend(lunout)
+         diff=0d0
+         if(calc_width.gt.0d0) then
+            diff=(calc_width-MC_width)/(calc_width+MC_width)*100d0
+         endif
+         write(lunout,'(a3,2x,i2,2x,e12.5,2x,e10.5,
+     &              1x,f5.2,3x,e10.5,3x,f5.1)') 
+     &        id(ip),imode,calc_width,MC_width,diff,mxwgt,eff
+         close(lunout)
+      ENDIF   
+      
+
+      return
+      end
+
+
+      real*4 function fxn(xx,wgt)
+c**********************************************************    
+c     vegas evaluation of the width
+c**********************************************************
+      implicit none
+      include 'decay.inc'
+c
+c     arguments
+c
+      real*4 xx(mxdim),wgt
+c                  
+c     Local
+c
+      integer IDS(MaxDecPart),ICOL(2,MaxDecPart)
+      integer jhel(MaxDecPart)
+      double precision pd(0:3,MaxDecPart),px(0:3),weight
+      integer i
+c
+c     External
+c
+      real xran1
+      integer iseed
+      data iseed/1/   !this value is irrelevant
+      integer get_hel
+c      
+c     Global
+c
+      include 'coupl.inc'
+      include 'calc_values.inc'
+c
+c-----
+c     Begin 
+c-----
+c--   momentum-random if you want to check Lorentz Invariance
+      PD(0,1)=M1
+      PD(1,1)=0d0
+      PD(2,1)=0d0
+      PD(3,1)=0d0
+      PD(1,1)=xran1(iseed)*10d0
+      PD(2,1)=xran1(iseed)*10d0
+      PD(3,1)=xran1(iseed)*10d0 
+      PD(0,1)=dsqrt(M1**2+PD(1,1)**2+PD(2,1)**2+PD(3,1)**2)
+
+c--   helicity
+      if(abs(ip).eq.24.or.ip.eq.23) then
+         jhel(1) = get_hel(xran1(iseed),3)
+      else
+         jhel(1) = get_hel(xran1(iseed),2)
+      endif
+c--   pass information to the common block
+      do i=1,mxdim
+         x(i)=xx(i)
+      enddo
+C      
+      CALL EVENT(PD,JHEL,IDS,ICOL,WEIGHT)
+C      
+c      write(*,*) 'from fxn: weight=',weight
+      mxwgt   = max(weight*wgt,mxwgt) !max weight = wgt * fxn
+      fxn     = REAL(weight)
+C     
+      return
+      end
+
+
+      subroutine toend(iunit)
+c**********************************************************    
+c read a file until the end
+c**********************************************************
+      ios = 0
+      dowhile(ios.eq.0)
+         read(unit=iunit,fmt='(1x)',iostat=ios)
+      enddo
+      end
+
+      subroutine find_particle(ipp,nexternal,ic,n,idecay)
+c***********************************************************    
+c    looks for particle ip in the event ic and
+c    returns the n instances found which
+c    are not decayed yet. idecay(i) gives the position
+c    of the i-th particle to be decayed.
+c***********************************************************
+      implicit none
+      include 'decay.inc'
+c
+c     Arguments
+c
+      integer n,ipp,nexternal,ic(7,MaxParticles),idecay(MaxParticles)
+c
+c     Local
+c
+      integer i
+
+c-----
+c  Begin Code
+c-----
+      
+      n=0
+      do i=1,nexternal
+         idecay(i)=0
+         if (ic(1,i) .eq. ipp .and. ic(6,i) .eq. 1 ) then
+            n=n+1
+            idecay(n)=i
+         endif
+      enddo
+      
+      return
+      end
+      
+
+
+      subroutine rnd_init(iunit,iseed)
+c***********************************************************    
+c initialize rnd number generator xran1 by reading iseed.dat
+c***********************************************************
+      implicit none
+c
+c     Arguments
+c
+      integer iunit,iseed
+c
+c     External
+c
+      real xran1
+c
+c     Local
+c
+      real*8 aa
+      character*50 adummy
+      logical done
+      character*(132) buff
+c-----
+c  Begin Code
+c-----
+
+c
+c     calculating decay rates only
+c
+c- try to open the iseed.dat file 
+      open(iunit,file='./iseed.dat',status="old",err=200)
+      read(iunit,err=200,end=200,fmt='(i10)') iseed
+      write(*,*) '*****************************************************'
+      write(*,*) '*      reading seed from iseed.dat                  *'
+      write(*,'(a22,i6,a26)')   
+     &           ' *  rnd number seed = ',iseed ,'               *'
+
+      rewind(iunit)
+      write(iunit,fmt='(i10)') iseed-1
+      close(iunit)
+      write(*,*) '*****************************************************'
+      goto 201
+c- if iseed.dat does not exist then set it to -1.
+ 200  iseed=-1
+      open(iunit,file='./iseed.dat',status="new")
+      write(iunit,fmt='(i10)') iseed-1
+      write(*,*) '*****************************************************'
+      write(*,*) '*      no iseed.dat  => iseed=-1                    *'
+      write(*,*) '*      iseed.dat  now written                       *'
+      write(*,*) '*****************************************************'
+ 201  continue
+      close(iunit)
+      return
+
+      return
+ 99   write (*,*) 'error'
+
+      return
+      end
+
+
+

=== added file 'DECAY/decay.inc'
--- DECAY/decay.inc	1970-01-01 00:00:00 +0000
+++ DECAY/decay.inc	2011-04-12 19:56:30 +0000
@@ -0,0 +1,105 @@
+c*****************************************************
+c     PARAMETERS
+c*****************************************************     	
+c
+c     VEGAS MAX DIMENSIONS
+c
+      integer 	 mxdim  	!to be the same as in vegas.f 
+      parameter (mxdim=10)      !rebin.f
+c                               
+      integer    MaxParticles	 ! max number of particles in 
+      parameter (MaxParticles=25)! the event
+      integer    MaxDecPart	 ! max particles in the decay
+      parameter (MaxDecPart=5)   ! 1->2,..,MaxDecPart  
+      integer    MaxiEvent	 ! max nr of processes
+      parameter (MaxiEvent=100)  ! in the proc_card
+c
+      integer MaxPoints  	! maxpoints used in the integration
+      parameter(MaxPoints=500000)
+c
+c     Numbers
+c
+      double precision  Zero, One, Two, Three, Four, Half,Rt2
+      parameter( Zero = 0.0d0, One = 1.0d0, Two = 2.0d0 )
+      parameter( Three = 3.0d0, Four = 4.0d0, Half = 0.5d0 )
+      parameter( Rt2   = 1.414213562d0 )
+      double precision  Pi, Fourpi
+      parameter( Pi = 3.14159265358979323846d0 )
+      parameter( Fourpi = Four * Pi )
+c
+c     Units	
+c
+      integer luni,lunr,lunw,lunp,lunrnd,lunout
+      parameter (luni=15,lunr=16,lunw=17,lunp=18,lunrnd=19,lunout=20) 	
+c
+c*****************************************************
+c     GLOBAL
+c*****************************************************     	
+c
+c     IP    = PDG CODE FOR THE PARTICLE TO DECAY
+c     IMODE = DECAY MODE
+c     ND    = NUMBER OF DECAY PRODUCTS
+c
+      INTEGER       IP,IMODE,ND	
+      COMMON/DEC_ID/IP,IMODE,ND
+c
+c     masses and couplings for the decay routines
+c
+      REAL*8          M1,M2,M3,M4,M5,MV,GV
+      COMMON/m_decays/M1,M2,M3,M4,M5,MV,GV
+      DOUBLE COMPLEX       GX,GXX(2),GXX1(2)  
+      COMMON/WEAK_COUPL/   GX,GXX   ,GXX1
+c
+c    MC width, analytic width and best branching ratio
+c     	
+      real*8          MC_width,calc_width,bratio,calc_br
+      common/to_width/MC_width,calc_width,bratio,calc_br
+c
+c    Particle's ID
+c     	
+      character*3 id(-100:100)
+      real*8 pdgmass(-100:100)
+      common/to_id/pdgmass,id
+c
+c    Color label
+c	
+      integer cindex 
+      common/c_ind/cindex
+c
+c    Vegas 	
+c
+      real*4          x(MXDIM)
+      COMMON/TO_VEGAS/x
+      integer                                 n_wei,n_ove
+      real*8          mxwgt,s_wei,s_unw,s_ove
+      common/to_event/mxwgt,s_wei,s_unw,s_ove,n_wei,n_ove
+c
+c     HDECAY branching ratios
+c
+      integer           ihiggs, nnlo, ipole
+      common /flag/     ihiggs, nnlo, ipole  
+      integer           ionsh, ionwz, iofsusy
+      common /onshell/  ionsh, ionwz, iofsusy 
+      double precision  gf, alph, almass, ammuon, amz, amw
+      common /param/    gf, alph, almass, ammuon, amz, amw 
+      double precision  ams, amc, amb, amt
+      common /masses/   ams, amc, amb, amt 
+      double precision  gamc0, gamt0, gamt1, gamw, gamz
+      common /wzwdth/   gamc0, gamt0, gamt1, gamw, gamz 
+      double precision  SMBRB,SMBRL,SMBRM,SMBRS,SMBRC,SMBRT,
+     &                  SMBRG,SMBRGA,SMBRZGA,SMBRW,SMBRZ,SMWDTH
+      common /widthsm/  SMBRB,SMBRL,SMBRM,SMBRS,SMBRC,SMBRT,
+     &                  SMBRG,SMBRGA,SMBRZGA,SMBRW,SMBRZ,SMWDTH
+c
+      double precision  hdals, hdmhbeg, hdmhend, tgbet
+      common /hdparms/  hdals, hdmhbeg, hdmhend, tgbet
+c
+c     PDG branching ratios
+c
+      double precision br_ta_lv,br_ta_pi,br_ta_ro
+      double precision br_w_lv,br_w_jj
+      double precision br_z_ll,br_z_vv,br_z_cc,br_z_bb,br_z_jj 
+      common/brs/      br_ta_lv,br_ta_pi,br_ta_ro,
+     &                 br_w_lv,br_w_jj,
+     &                 br_z_ll,br_z_vv,br_z_cc,br_z_bb,br_z_jj 
+

=== added file 'DECAY/decay_couplings.f'
--- DECAY/decay_couplings.f	1970-01-01 00:00:00 +0000
+++ DECAY/decay_couplings.f	2011-04-12 19:56:30 +0000
@@ -0,0 +1,725 @@
+c----------------------------------------------------------------------
+C  couplings.f 
+c----------------------------------------------------------------------
+c  This files takes the inputs of the standard model from a Les Houches 
+c  file (param_card.dat) and calculates all the couplings that end up
+c  in the Feynman rules, i.e., in the HELAS routines.
+c   
+c  With the convention of the New Web Generation MadEvent in this
+c  part no non-trivial computation is made. The SM model secondary
+c  parameters have been all calculated by the SM calculator, SMCalc
+c  which produced param_card.dat.
+c
+c  The only exception to the above rule is for alpha_S. In the case
+c  where a pdf is used in the initial state, then the values as(MZ)
+c  set in the les houches card is superseeded.
+c  Scales and pdf's are set in the run_card.dat.
+c
+c This file contains the following routines:
+c 
+c- madgraph original call to set the parameters
+c- lh_readin is called from here.
+c  It talks to the lh_readin through the common block values.
+c      subroutine setpara
+c
+c-routine to read the LH parameters in
+c      subroutine lh_readin
+c
+c-to set
+c      subroutine set_it(npara,ivalue,value,name,id,
+c      subroutine case_trap(string,length)
+c      subroutine no_spaces(buff,nchars)
+c---------------------------------------------------------------------- 
+
+
+      subroutine setpara
+c***********************************************************************
+c This subroutine sets up the HELAS couplings of the STANDARD MODEL.
+c***********************************************************************
+      implicit none
+c
+c local
+c
+      integer i
+      real*8 dum
+c
+c     common file with the couplings
+c
+      include 'coupl.inc'
+      include 'decay.inc'
+      include 'calc_values.inc'
+c
+c     local
+c
+      double precision  v
+      double precision  ee, ee2, ez, ey, sw, cw, sc2
+      double precision  gwne, gwud, lambda, lam4, xt, rew, rqcd
+      double precision  alphas, alfa, alfaw, mfrun
+      external          alphas, alfa, alfaw, mfrun
+c
+c     Common to lh_readin and printout
+c
+      double precision  alpha, sin2w, gfermi, alfas
+      double precision  mtMS,mbMS,mcMS,mtaMS!MSbar masses
+      double precision  Vud,Vus             !CKM matrix elements
+      common/values/    alpha,sin2w,gfermi,alfas,   
+     &                  mtMS,mbMS,mcMS,mtaMS,
+     &                  Vud
+c
+c constants
+c
+      double complex  ci
+      parameter( ci = ( 0.0d0, 1.0d0 ) )
+c
+c     alfas and run
+c
+      include 'alfas.inc'
+      include 'run.inc'
+c
+c     Auxiliary functions
+c
+      REAL*8 XX,Y,Z
+      REAL*8 HFF,BETA,HEAVY
+
+      BETA(XX) = DSQRT(One-Four*XX)
+      HFF(XX,Y)  = One/(Two*Four*PI)*XX*(BETA(Y))**3
+      HEAVY(XX,Y,Z)= ( One - Half*(y**2+z**2) - Half*(y**2-z**2)**2
+     &                 + Three*y*z*(XX**2 - One)/(XX**2 + One)        )
+     &               * sqrt( (One-y**2-z**2)**2 - Four * y**2 * z**2 )
+c
+c------------------------------------------
+c Start calculating the couplings for HELAS
+c------------------------------------------
+c
+      call lh_readin
+c     
+c     Strong coupling
+c
+c     As a rule we first check if a pdf has been chosen in the    
+c     run_card.dat (which has been already read at this stage).
+c     If there pdfs in the initial state, then the alpha_s(MZ) used
+c     is set to the corresponding value.  
+
+      if(scale.le.1d0)    scale=zmass   
+   
+      if(lpp(1).ne.0.or.lpp(2).ne.0) then    
+          if(asmz .le.0.01d0) asmz =0.118d0   
+      else 
+	      asmz=alfas   !value read from the param_card.dat
+c          nloop=2 
+      endif      
+      if(nloop.eq.0)      nloop=1   
+
+      G = DSQRT(4d0*PI*ALPHAS(SCALE)) ! use setting of the param_card.dat @ NLO
+      GG(1) = -G
+      GG(2) = -G     
+c
+c auxiliary local values
+c
+      cw  = sqrt( One - sin2w )
+      ee2 = alpha * Fourpi
+      sw  = sqrt( sin2w )
+      ee  = sqrt( ee2 )
+      ez  = ee/(sw*cw)
+      ey  = ee*(sw/cw)
+      sc2 = sin2w*( One - sin2w )
+      v   = Two*wmass*sw/ee   ! the wmass is used to calculate v
+      lambda = hmass**2 / (Two * v**2)
+c
+c vector boson couplings
+c
+      gw   = ee/sw
+      gwwa = ee
+      gwwz = ee*cw/sw
+c
+c gauge & higgs boson coupling constants
+c
+      gwwh  = dcmplx( ee2/sin2w*Half*v, Zero )
+      gzzh  = dcmplx( ee2/sc2*Half*v, Zero )
+      ghhh  = dcmplx( -hmass**2/v*Three, Zero )
+      gwwhh = dcmplx( ee2/sin2w*Half, Zero )
+      gzzhh = dcmplx( ee2/sc2*Half, Zero)
+      ghhhh = ghhh/v
+c
+c fermion-fermion-vector couplings
+c
+      gal(1) = dcmplx(  ee          , Zero )
+      gal(2) = dcmplx(  ee          , Zero )
+      gau(1) = dcmplx( -ee*Two/Three, Zero )
+      gau(2) = dcmplx( -ee*Two/Three, Zero )
+      gad(1) = dcmplx(  ee/Three    , Zero )
+      gad(2) = dcmplx(  ee/Three    , Zero )
+
+      gwf(1) = dcmplx( -ee/sqrt(Two*sin2w), Zero )
+      gwf(2) = dcmplx(  Zero              , Zero )
+
+      gzn(1) = dcmplx( -ez*Half                     , Zero )
+      gzn(2) = dcmplx(  Zero                        , Zero )
+      gzl(1) = dcmplx( -ez*(-Half + sin2w)          , Zero )
+      gzl(2) = dcmplx( -ey                          , Zero )
+      gzu(1) = dcmplx( -ez*( Half - sin2w*Two/Three), Zero )
+      gzu(2) = dcmplx(  ey*Two/Three                , Zero )
+      gzd(1) = dcmplx( -ez*(-Half + sin2w/Three)    , Zero )
+      gzd(2) = dcmplx( -ey/Three                    , Zero )
+c
+c fermion-fermion-Higgs couplings (complex) hff(2)
+c
+c NOTE: the running mass is evaluated @ the same order 
+c nloop of alpha_s set by the PDF choice
+c 
+
+      if(mtMS.gt.1d0) then
+         ghtop(1) = dcmplx( -mtMS/v, Zero )
+      else
+         ghtop(1) = dcmplx( Zero,Zero)
+      endif
+      ghtop(2) = ghtop(1)
+
+      if(mbMS.gt.1d0) then
+         ghbot(1) = dcmplx( -mbMS/v, Zero )
+      else
+         ghbot(1) = dcmplx( Zero, Zero )
+      endif
+      ghbot(2) = ghbot(1)
+      
+      if(mcMS.gt.1d0) then
+         ghcha(1) = dcmplx( -mcMS/v, Zero )
+      else
+         ghcha(1) = dcmplx( Zero, Zero )
+      endif
+      ghcha(2) = ghcha(1)
+
+      ghtau(1) = dcmplx( -mtaMS/v, Zero )
+      ghtau(2) = ghtau(1)
+c
+c     CKM matrix: 
+c     symmetric 3x3 matrix, Vud=Vcs, Vus=Vcd Vcb=Vub=0
+c
+c     >>>>>>>>>>>>>>>***** NOTE****<<<<<<<<<<<<<<<<<<<<<<<<<
+c     these couplings matter only when interaction_CKM.dat
+c     is used to generate all the diagrams with off-diagonal
+c     couplings. The default of MadEvent is a diagonal
+c     CKM matrix.
+
+	  Vus=DSQRT(1d0-Vud**2)
+      do i=1,2
+         gwfc(i) = gwf(i)*Vud
+         gwfs(i) = gwf(i)*Vus
+         gwfm(i) =-gwf(i)*Vus
+      enddo
+
+c---------------------------------------------------------
+c Set Photon Width to Zero, used by symmetry optimization
+c---------------------------------------------------------
+
+      awidth = 0d0
+c     
+c Z boson partial widths
+c     
+      decz = zmass / ( 24.0d0 * Pi )
+      w_z_nn = decz * ( gzn(1)**2 + gzn(2)**2 )
+      w_z_ll = decz * ( gzl(1)**2 + gzl(2)**2 )
+      decz = decz * Three 
+      w_z_uu = decz * ( gzu(1)**2 + gzu(2)**2 )
+      w_z_dd = decz * ( gzd(1)**2 + gzd(2)**2 )
+      dum = dble( (gzl(2)+gzl(1))/(gzl(2)-gzl(1)) )
+      w_z_tau = w_z_ll * heavy( dum, lmass/zmass, lmass/zmass )
+      dum = dble( (gzu(2)+gzu(1))/(gzu(2)-gzu(1)) )
+      w_z_cc = w_z_uu *  heavy( dum, cmass/zmass, cmass/zmass )
+      dum = dble( (gzd(2)+gzd(1))/(gzd(2)-gzd(1)) )
+      w_z_bb = w_z_dd *  heavy( dum, bmass/zmass, bmass/zmass )
+      
+         zwidth =   Three*w_z_nn + Two*w_z_ll + w_z_tau
+     &        + Two*w_z_dd + w_z_uu + w_z_cc + w_z_bb
+
+c
+c W boson partial widths
+c     
+      decw = wmass / ( 24.0d0 * Pi )
+      w_w_nl = decw * ( gwf(1)**2 + gwf(2)**2 )
+      dum = dble( (gwf(2)+gwf(1))/(gwf(2)-gwf(1)) )
+      w_w_tau = w_w_nl * heavy( dum, lmass/wmass, Zero )
+      w_w_ud = w_w_nl * Three 
+      w_w_cs = w_w_ud * heavy( dum, cmass/wmass, Zero )
+      
+      wwidth = Two*w_w_nl + w_w_tau + w_w_ud + w_w_cs
+
+c
+c top quark width
+c
+      call topwid(tmass,wmass,bmass,wwidth,gw,twidth) 
+
+c
+c tau width
+c
+      lwidth  = 2.36d-12 !tau width, PDG value
+c
+c     LO withds of the Higgs into tt~,bb~,tau tau~.
+c
+      if(hmass.gt.2d0*tmass) then
+         w_h_tt =3d0*cdabs(ghtop(1)**2)*hff(hmass,(tmass/hmass)**2)
+      else
+         w_h_tt =0d0
+      endif
+
+      w_h_bb =3d0*cdabs(ghbot(1)**2)*hff(hmass,(bmass/hmass)**2)
+      w_h_tau=cdabs(ghtau(1)**2)*hff(hmass,(lmass/hmass)**2)
+      w_h_bb =3d0*cdabs(ghbot(1)**2)*hff(hmass,(bmass/hmass)**2)
+      w_h_cc =3d0*cdabs(ghcha(1)**2)*hff(hmass,(cmass/hmass)**2)
+      w_h_tau=cdabs(ghtau(1)**2)*hff(hmass,(lmass/hmass)**2)
+
+      
+      if(hmass.gt.2d0*wmass) then         
+         w_h_ww=gfermi/8d0/pi/rt2*hmass**3*
+     &        dsqrt(1-4d0*(wmass/hmass)**2)*
+     &        (12d0*(wmass/hmass)**4 -4d0*(wmass/hmass)**2+1d0)         
+         w_h_WLWL=gw**2/64d0/pi*hmass**3/wmass**2* !longitudinal W's
+     &        dsqrt(1-4d0*(wmass/hmass)**2)*
+     &        (1-2d0*(wmass/hmass)**2)**2 
+      else
+         w_h_ww=0d0
+         w_h_WLWL=0d0
+      endif
+      
+      if(hmass.gt.2d0*zmass) then
+         w_h_zz=gfermi/16d0/pi/rt2*hmass**3*
+     &        dsqrt(1-4d0*(zmass/hmass)**2)*
+     &        (12d0*(zmass/hmass)**4 -4d0*(zmass/hmass)**2+1d0)
+         w_h_ZLZL=gw**2/128d0/pi*hmass**3/wmass**2* !longitudinal Z's
+     &        dsqrt(1-4d0*(zmass/hmass)**2)*
+     &        (1-2d0*(zmass/hmass)**2)**2         
+      else
+         w_h_zz=0d0
+         w_h_zLzL=0d0
+      endif
+            
+c--------------------------
+c start interface to HDECAY
+c--------------------------
+c do not change things here unless you exactly know what you are doing
+c
+      ihiggs = 0
+      nnlo   = 0
+      ipole  = 0
+
+      ionsh   = 0
+      ionwz   = 0
+      iofsusy = 1
+
+      hdals  = asmz
+c-- do not set masses to zero here
+      ams    = 0.190d0    !strange pole mass
+
+      if(cmass.gt.0d0) then 
+         amc    = cmass
+      else
+         amc    = 1.42d0        !charm   pole mass      
+      endif
+
+      if(bmass.gt.0d0) then 
+         amb    = bmass
+      else
+         amb    = 4.7d0          !bottom   pole mass      
+      endif
+
+      if(tmass.gt.0d0) then 
+         amt    = tmass
+      else
+         amt    = 174.3d0        !top pole mass      
+      endif
+
+      if(lmass.gt.0d0) then 
+         almass    = lmass
+      else
+         almass    = 1.777d0     !tau  mass      
+      endif
+
+      ammuon = 0.105658389d0 !muon mass
+
+      alph   = 137.0359895D0 !alpha_EM
+      gf     = gfermi    
+
+      if(wwidth.gt.0d0) then
+         gamw   = wwidth
+      else
+         gamw=2.12d0
+      endif
+      
+      if(zwidth.gt.0d0) then
+         gamz   = zwidth
+      else
+         gamz=2.495d0
+      endif
+      
+      amz    = zmass
+      amw    = wmass
+
+      hdmhbeg = hmass
+      hdmhend = hmass
+
+c 
+c     this calculates the branching ratios of the Higgs
+c     at the best of our knowledge 
+c
+      call hdecay     
+      hwidth = SMWDTH
+
+c------------------------
+c end interface to HDECAY
+c------------------------
+
+
+c----------------------------
+c end subroutine coupsm
+c----------------------------
+
+
+      return
+      end
+ 
+
+      subroutine lh_readin
+c----------------------------------------------------------------------
+c Read the parameters from the lh file
+c
+c 1. Input values for the EW sector 
+c 2. Higgs mass and width
+c 3. Fermion masses (pole and MSbar) and widths
+c---------------------------------------------------------------------- 
+      implicit none
+c
+c     parameters
+c
+      integer maxpara
+      parameter (maxpara=1000)
+c
+c     local
+c
+      integer npara,l1,l2
+      character*132 buff
+      real*8 real_value 
+      real*8 value(maxpara)
+      integer ivalue(maxpara),n
+      character*20 name(maxpara),bn,dumstring
+      logical block_found,done,fopened	  
+      integer i,name_length,idum
+
+c
+c	block info
+c      
+      character*20 block_name
+c
+c   Common
+c
+      include 'coupl.inc'
+      include 'decay.inc'
+c
+c     Common to lh_readin and printout
+c
+      double precision  alpha, sin2w, gfermi, alfas
+      double precision  mtMS,mbMS,mcMS,mtaMS!MSbar masses
+      double precision  Vud,Vus             !CKM matrix elements
+      common/values/    alpha,sin2w,gfermi,alfas,   
+     &                  mtMS,mbMS,mcMS,mtaMS,
+     &                  Vud
+c
+c----------
+c     start
+c----------
+c
+        n=0
+       rewind(lunp)
+       done=.false.
+
+       do while(.not.done)
+          block_found=.false.
+c
+c looks for the blocks or for decay
+c
+      do while(.not.block_found)
+       read(lunp,'(a132)',end=99,err=99) buff
+c--	change to lower case
+	   call case_trap(buff,20)
+       if(buff(1:5).eq.'block') then
+         if(index(buff,"#").ne.0) l1=index(buff,"#")-1 
+         block_name=buff(6:min(l1,26)) 
+         call no_spaces(block_name,name_length)
+c        write(*,*) block_name(1:name_length)
+         block_found=.true.        
+	   elseif(buff(1:5).eq.'decay') then
+               n=n+1
+               l1=30
+               if(index(buff,"#").ne.0) l1=index(buff,"#")-1 ! ignore comments       
+               read(buff(6:l1),*) ivalue(n),value(n)
+               name(n)="decay" 
+       endif		   
+      end do
+c
+c     
+ 
+      if(block_found) then
+ 	  do while(.true.)
+          read(lunp,'(a132)',end=99,err=99) buff
+     	  call case_trap(buff,20)
+          if(buff(1:1).eq.'b'.or.buff(1:1).eq.'d') then
+          	backspace lunp
+          	exit
+          endif	         	
+            if(buff(1:1).ne.'#') then  !if it not a comment
+              n=n+1	       
+              l1=30
+              if(index(buff,"#").ne.0) l1=index(buff,"#")-1 ! ignore comments       
+c
+c  WARNING:... not all blocks have the same sintax!! You need to change it
+c  depending on the block you are reading
+c
+          
+             if(block_name(1:5).eq."mgckm") then
+                read(buff(1:l1),*) ivalue(n),idum,value(n)
+              elseif (block_name(1:6).eq."spinfo".or.
+     &block_name(1:6).eq."dcinfo") then
+                read(buff(1:l1),*) ivalue(n),dumstring
+              else
+                read(buff(1:l1),*) ivalue(n),value(n)
+              endif  
+              name(n)=block_name(1:name_length)
+c             write(*,"(1x,i2,2x,e16.8,1x,a)") 
+c     &        ivalue(n),value(n),name(n)
+           	  endif
+      end do ! do while in the block
+      else
+	  done=.true.
+	  endif
+	  end do ! do while the entire file
+	  
+
+99	continue      
+	
+	   bn="sminputs"
+       call set_it(n,ivalue,value,name,1,bn,alpha,128.9d0)
+       alpha=1d0/alpha
+       call set_it(n,ivalue,value,name,2,bn,gfermi,0.1166d-4)
+       call set_it(n,ivalue,value,name,3,bn,alfas,0.119d0)
+       call set_it(n,ivalue,value,name,4,bn,zmass,91.188d0)
+       call set_it(n,ivalue,value,name,6,bn,tmass,174.3d0)
+       call set_it(n,ivalue,value,name,7,bn,lmass,1.777d0)
+	   bn="mgsmparam"
+       call set_it(n,ivalue,value,name,1,bn,sin2w,0.2312d0)
+       call set_it(n,ivalue,value,name,2,bn,wmass,80.419d0)
+       bn="mgyukawa"
+       call set_it(n,ivalue,value,name,4,bn,mcMS,1.25d0)
+       call set_it(n,ivalue,value,name,5,bn,mbMS,4.2d0)
+       call set_it(n,ivalue,value,name,6,bn,mtMS,174d0)
+       call set_it(n,ivalue,value,name,15,bn,mtaMS,1.777d0)
+       bn="mgckm"
+       call set_it(n,ivalue,value,name,1,bn,vud,1d0)
+       bn="mass"
+       call set_it(n,ivalue,value,name,4,bn,cmass,1.4d0)
+       call set_it(n,ivalue,value,name,5,bn,bmass,4.7d0)
+       call set_it(n,ivalue,value,name,6,bn,tmass,tmass*1d0)
+       call set_it(n,ivalue,value,name,15,bn,lmass,lmass*1d0)
+       call set_it(n,ivalue,value,name,25,bn,hmass,120d0)
+       call set_it(n,ivalue,value,name,23,bn,zmass,zmass*1d0)
+       call set_it(n,ivalue,value,name,24,bn,wmass,wmass*1d0)
+
+c       bn="decay"
+c       call set_it(n,ivalue,value,name,6,bn,twidth,1.5083d0)
+c       call set_it(n,ivalue,value,name,25,bn,hwidth,0.0037d0)
+c       call set_it(n,ivalue,value,name,23,bn,zwidth,2.441d0)
+c       call set_it(n,ivalue,value,name,24,bn,wwidth,2.0476d0)
+
+       write(*,*)
+       write(*,*) ' >>> Widths in param_card.dat are ignored <<<'
+       write(*,*)
+
+      return
+      end
+
+      
+      subroutine set_it(npara,ivalue,value,name,id,
+     &                  block_name,var,def_value)
+c----------------------------------------------------------------------------------
+c     finds the parameter value  in block_name and associate var to it.
+c     If it is not found a default is given.
+c----------------------------------------------------------------------------------
+      implicit none
+
+c
+c     parameters
+c
+      integer maxpara
+      parameter (maxpara=100)
+c
+c     arguments
+c
+      integer npara,ivalue(maxpara),id
+      character*20  block_name,name(maxpara)
+      real*8 var,def_value,value(maxpara)
+c
+c     local
+c
+      logical found
+      integer i
+c
+c     start
+c
+	  found=.false.
+      do i=1,npara
+         found = (id.eq.ivalue(i)).and.(name(i).eq.block_name)
+ 	               if(found) then
+         	var=value(i)
+            exit
+          endif	
+      enddo
+      
+      if (.not.found) then
+c         write (*,*) "Warning: parameter not found"
+c         write (*,*) "         setting it to default value ",def_value
+         var=def_value
+      endif
+      return
+
+      end
+      
+      
+      subroutine case_trap(string,length)
+c**********************************************************    
+c change string to lowercase if the input is not
+c**********************************************************
+      implicit none
+c
+c     ARGUMENT
+c      
+      character*(*) string
+      integer length
+c
+c     LOCAL
+c
+      integer i,k
+
+      do i=1,length
+         k=ichar(string(i:i))
+         if(k.ge.65.and.k.le.90) then  !upper case A-Z
+            k=ichar(string(i:i))+32   
+            string(i:i)=char(k)        
+         endif
+      enddo
+
+      return
+      end
+
+
+      subroutine no_spaces(buff,nchars)
+c**********************************************************************
+c     Given buff a buffer of words separated by spaces
+c     returns it where all space are moved to the right
+c     returns also the length of the single word.
+c     maxlength is the length of the buffer
+c**********************************************************************
+      implicit none
+c
+c     Constants
+c
+      integer    maxline
+      parameter (maxline=20)
+      character*1 null
+      parameter  (null=' ')
+c
+c     Arguments
+c
+      character*(maxline) buff
+      integer nchars,maxlength
+c
+c     Local
+c
+      integer i,j
+      character*(maxline) temp
+c-----
+c  Begin Code
+c-----
+      nchars=0
+c      write (*,*) "buff=",buff(1:maxlength)
+      do i=1,maxline
+         if(buff(i:i).ne.null) then
+            nchars=nchars+1
+            temp(nchars:nchars)=buff(i:i)
+         endif
+c         write(*,*) i,":",buff(1:maxlength),":",temp(1:nchars),":"
+      enddo
+      buff=temp      
+      end
+
+      
+      SUBROUTINE TOPWID(RMT,RMW,RMB,RGW,GW,RGT)
+c*************************************************************************
+c     THE TOTAL WEAK DECAY WIDTH OF THE TOP QUARK, INCLUDING
+c     THE EFFECTS OF BOTTOM MASS AND, IF IGW=1,  A FINITE W WIDTH.
+c     From James Stirling 6-10-94
+c
+c     RMT=TOP MASS
+c     RMW=W   MASS
+c     RMB=B   MASS
+c     RGW=W   WIDTH
+c     GW =WEAK COUPLING
+c
+c     RGT=TOP WIDTH
+c
+c*************************************************************************
+      IMPLICIT COMPLEX*16(A-H,O-Z)
+      REAL*8 RMT,RMB,RMW,XW,XB,RGW,RGT,GW
+*
+      PI=4.*DATAN(1.D0)
+      XGW=dcmplx(GW/2d0/dsqrt(2d0))
+*
+      XB=RMB/RMT
+      XW=RMW/RMT
+*
+      RM=XB**2
+      OM=1.+RM-DCMPLX(RMW**2,RMW*RGW)/RMT**2
+      Y1=OM+CDSQRT(OM*OM-4.*RM)
+      Y0=OM-CDSQRT(OM*OM-4.*RM)
+      Z1=2.
+      Z0=2.*CDSQRT(RM)
+*
+      D0=(-Y0**8+3.*Y0**7*RM+3.*Y0**7-8.*Y0**6*RM-12.*Y0**5*RM**
+     . 2-12.*Y0**5*RM+96.*Y0**4*RM**2-48.*Y0**3*RM**3-48.*Y0**3*
+     . RM**2-128.*Y0**2*RM**3+192.*Y0*RM**4+192.*Y0*RM**3-256.*
+     . RM**4)/(24.*Y0**4*(Y1-Y0))
+      D1=(-Y1**8+3.*Y1**7*RM+3.*Y1**7-8.*Y1**6*RM-12.*Y1**5*RM**
+     . 2-12.*Y1**5*RM+96.*Y1**4*RM**2-48.*Y1**3*RM**3-48.*Y1**3*
+     . RM**2-128.*Y1**2*RM**3+192.*Y1*RM**4+192.*Y1*RM**3-256.*
+     . RM**4)/(24.*Y1**4*(Y1-Y0))
+      A4=(32.*RM**4*(Y1-Y0))/(3.*Y1*Y0*(Y1-Y0))
+      A3=(8.*RM**3*(-3.*Y1**2*Y0*RM-3.*Y1**2*Y0+4.*Y1**2*RM+3.*
+     . Y1*Y0**2*RM+3.*Y1*Y0**2-4.*Y0**2*RM))/(3.*Y1**2*Y0**2*(Y1
+     . -Y0))
+      A2=(8.*RM**3*(2.*Y1**3*Y0**2-3.*Y1**3*Y0*RM-3.*Y1**3*Y0+4.
+     . *Y1**3*RM-2.*Y1**2*Y0**3+3.*Y1*Y0**3*RM+3.*Y1*Y0**3-4.*Y0
+     . **3*RM))/(3.*Y1**3*Y0**3*(Y1-Y0))
+      A1=(2.*RM**2*(3.*Y1**4*Y0**3*RM+3.*Y1**4*Y0**3+8.*Y1**4*Y0
+     . **2*RM-12.*Y1**4*Y0*RM**2-12.*Y1**4*Y0*RM+16.*Y1**4*RM**2
+     . -3.*Y1**3*Y0**4*RM-3.*Y1**3*Y0**4-8.*Y1**2*Y0**4*RM+12.*
+     . Y1*Y0**4*RM**2+12.*Y1*Y0**4*RM-16.*Y0**4*RM**2))/(3.*Y1**
+     . 4*Y0**4*(Y1-Y0))
+      B0=(Y1**3-3.*Y1**2*RM-3.*Y1**2+8.*Y1*RM-Y0**3+3.*Y0**2*RM+
+     . 3.*Y0**2-8.*Y0*RM)/(24.*(Y1-Y0))
+      B1=(Y1+Y0-3.*RM-3.)/24.
+      B2=1./24.
+*
+      RINT=D0*CDLOG((Z1-Y0)/(Z0-Y0))
+     .    -D1*CDLOG((Y1-Z1)/(Y1-Z0))
+     .    -A4/3.*(1./Z1**3-1./Z0**3)
+     .    -A3/2.*(1./Z1**2-1./Z0**2)
+     .    -A2   *(1./Z1   -1./Z0   )
+     .    +A1*CDLOG(Z1/Z0)
+     .    +B0   *(Z1   -Z0   )
+     .    +B1/2.*(Z1**2-Z0**2)
+     .    +B2/3.*(Z1**3-Z0**3)
+*
+      XGW4=XGW**4
+*
+* TOTAL WIDTH INCLUDES FLAVOUR & COLOUR FACTORS
+      RGT=RMT**3/(RMW*RGW)*XGW4/(8.*PI**3)*DIMAG(RINT)
+      RGT=9.*RGT
+      RETURN
+      END

=== added file 'DECAY/decay_event.f'
--- DECAY/decay_event.f	1970-01-01 00:00:00 +0000
+++ DECAY/decay_event.f	2011-04-12 19:56:30 +0000
@@ -0,0 +1,2537 @@
+      subroutine decay_event(P5,wgt,ni,ic)
+c********************************************************************
+c     Decay particle(s) in the event 
+c
+c     input: P(0,MaxParticles) :four momenta
+c            wgt               :event weight
+c            ni                :number of particle in the event 
+c                               before the decay
+c            ic(6,MaxParticles):particle labels       
+c
+c     output:P(0,MaxParticles) :four momenta
+c            wgt               :event weight
+c            ic(6,MaxParticles):particle labels       
+c
+c     which particle has to be decayed is passed through
+c     a common block DEC_ID in decay.inc
+c
+c********************************************************************
+      implicit none
+      include 'decay.inc'
+c
+c     Arguments
+c      
+      integer ni, ic(7,MaxParticles)
+      double precision P5(0:4,MaxParticles),wgt
+c
+c     Local
+c
+      integer i,k,jj, jhel(MaxDecPart), idecay, idec(MaxParticles),j, im
+      integer JD(MaxDecPart),ICOL(2,MaxDecPart)
+      integer maxcol
+      double precision pcms(0:3), pd(0:3,MaxDecPart),pswgt
+      double precision p(0:3,MaxParticles)
+      double precision totwidth, dwgt, goal_wgt
+      double precision weight,r
+      logical done
+      real wgt_vegas
+c
+c     External
+c
+      real xran1
+      integer iseed
+      data iseed/1/  !no effect if already called
+c      
+c     Global
+c
+      data s_unw,s_wei,s_ove/3*0d0/
+      data n_wei,n_ove/2*0/
+c
+      include 'coupl.inc'
+      include 'calc_values.inc'
+c
+c     first map the momenta into objects 0:3
+c
+      do i=1,MaxParticles
+         do j=0,3
+            p(j,i)=p5(j,i)
+         enddo
+      enddo
+
+c
+c--   First find out how many particles there are in the 
+c     event to decay
+c
+      call find_particle(ip,ni,ic,k,idec)         
+      if (k.eq.0) return !no particle to decay
+c
+c--   Loop over particles to decay
+c
+      do jj=1,k
+         idecay=idec(jj)
+c--   find the largest label present for the color flow
+         maxcol=max(ic(4,1),ic(5,1))
+         do i=2,ni
+            maxcol=max(maxcol,ic(4,i),ic(5,i))
+         enddo   
+         cindex=maxcol+1        !available color index
+c--   setup color of the decay particle
+         icol(1,1) =  ic(4,idecay) ! color      of particle
+         icol(2,1) =  ic(5,idecay) ! anti-color of particle
+c--   Boost to the momentum of particle to the CMS frame 
+         do i=0,3
+            pcms(i) = -p(i,1) - p(i,2)
+         enddo
+         pcms(0)=-pcms(0)
+         call boostx(p(0,idecay),pcms,pd(0,1))
+         jhel(1) =  ic(7,idecay) ! helicity   of particle 
+c--   Unweigthed decay
+         done=.false.
+         do while (.not. done)
+            CALL GET_X(X,WGT_VEGAS)
+            CALL EVENT(PD(0,1),JHEL(1),JD(1),ICOL,WEIGHT)
+            weight=weight*real(wgt_vegas)
+            n_wei=n_wei+1
+            s_wei=s_wei+weight
+            r = xran1(iseed)
+            if(weight / r .gt. mxwgt) then
+               done=.true.         
+               s_unw=s_unw+mxwgt
+               if (weight .gt. mxwgt) then
+                  s_ove = s_ove-mxwgt+weight
+                  n_ove = n_ove + 1
+               endif
+            endif
+         enddo
+c--   Multiply by the branching ratio
+         wgt = wgt * bratio     
+c--   Boost back to lab frame and put momentum into P(0,x)
+         do i=1,3
+            pcms(i)=-pcms(i)
+         enddo
+         do i=2,MaxDecPart      !loop over decay products
+            call boostx(pd(0,i),pcms,p(0,i-1+ni))
+         enddo
+c--   Set the new particle information
+c--   The info on the decayed particles are kept intact but 
+c--   status changed to decayed
+         do i=1,ND
+            ic(1,ni+i)=jd(i+1)  !particle's IDs
+            ic(2,ni+i)=idecay   !Mother info
+            ic(3,ni+i)=idecay   !Mother info
+            ic(4,ni+i)=icol(1,i+1) !     color info
+            ic(5,ni+i)=icol(2,i+1) !anti-color info
+            ic(6,ni+i)=+1       !Final Particle
+            ic(7,ni+i)=jhel(1+i) !Helicity info
+         enddo
+         ni = ni + nd           !Keep intermediate state decaying particle
+         ic(6,idecay)=+2        !This is decayed particle      
+c--   end loop over particles to decay
+      enddo                     
+c
+c     finally map the P back to 0:4
+c
+      do i=1,ni
+         do j=0,3
+            p5(j,i)=p(j,i)
+         enddo
+c         p5(4,i)=pdgmass(ic(1,i))
+      enddo
+
+      do i=ni-nd,ni
+          p5(4,i)=pdgmass(ic(1,i))
+      enddo
+
+      return
+      end
+
+
+      SUBROUTINE EVENT(PD,JHEL,JD,ICOL,DWGT)
+c**********************************************************    
+c     Calculation of the decay event quantities
+c
+c     input: pd(0:3,1) decay particle momentum
+c
+c     ouput: pd(0:3,MaxDecPart) momenta
+c            jhel(MaxDecPart)   helicities
+c            jd(MaxDecPart)     particle id's
+c            icol(2,MaxDecPart) color labels
+c            DWGT               pswgt*emmesq*factor
+c
+c**********************************************************
+      implicit none
+      include 'decay.inc'
+c
+c     Arguments
+c
+      real*8 pd(0:3,MaxDecPart)
+      integer jhel(MaxDecPart),JD(MaxDecPart),ICOL(2,MaxDecPart)
+      real*8 dwgt
+c
+c     Global
+c
+      include 'coupl.inc'
+      include 'calc_values.inc'
+c
+c     local
+c
+      real*8 pswgt,emmesq,fudge,factor
+      real*8 aa,aa1,aa2,aa3
+      integer isign,jl,jvl,i,j
+      logical hadro_dec
+      integer multi_decay
+      real*8 xprob1,xprob2
+c
+c     EXTERNAL
+c
+      integer get_hel
+      real xran1
+      integer iseed
+      data iseed/1/ !no effect if already called
+
+C------
+C START
+C------
+
+      DWGT   =0d0
+      emmesq =0d0
+      isign  =ip/abs(ip)
+      aa     =0d0
+      aa1    =0d0
+      aa2    =0d0
+      do i=2,MaxDecPart
+         icol(1,i)=0
+         icol(2,i)=0
+      enddo
+      
+c      write(*,*) 'from event: imode,ip',imode,ip
+
+      If(abs(ip).eq.6) then
+*-----------------------------------------------------
+*     top decays
+*-----------------------------------------------------
+
+         
+      If(imode.eq.1) then 
+*------------------------
+*     t  -> b  w+
+*     t~ -> b~ w-
+*------------------------
+
+c--   masses
+      M1=tmass
+      M2=bmass
+      M3=wmass
+c--   id's
+      jd(2)=isign*5  !b  or b~
+      jd(3)=isign*24 !w+ or w-
+c--   color
+      icol(1,2)=icol(1,1)  !     color of the top
+      icol(2,2)=icol(2,1)  !anti-color of the top
+c--   couplings
+      GXX(1)=gwf(1)
+      GXX(2)=gwf(2)
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*6d0    
+c--   helicities 
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  get_hel(xran1(iseed),3)
+c--   phase space
+      call phasespace(pd,pswgt)     
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      if(isign.eq. 1) call emme_ffv(pd,jhel,emmesq)      
+      if(isign.eq.-1) call emme_fxfv(pd,jhel,emmesq)      
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+
+      return
+      endif !end imode=1
+
+      if(imode.ge.2.and.imode.le.8) then 
+*------------------------
+*     t  -> b  vl l+  ,jj
+*     t~ -> b~ l  vl~ ,jj
+*------------------------
+
+c--   masses
+      M1=tmass
+      M2=bmass
+      M3=0d0
+      M4=0d0
+      MV=WMASS
+      GV=WWIDTH
+c--   color
+      icol(1,2)=icol(1,1)  !     color of the top
+      icol(2,2)=icol(2,1)  !anti-color of the top
+c--   id's
+      if(isign.eq.1) then       !t
+         jl=4
+         jvl=3
+      else                      !t~
+         jl=3
+         jvl=4
+      endif   
+      jd(2)=isign*5             !b  or b~
+c--   rnd number 
+      if(imode.ge.5) aa =dble(xran1(iseed))
+      if(imode.eq.8) aa1=dble(xran1(iseed))
+c
+c     various cases imode=2,3,4,5,6,7,8
+c
+      if(imode.eq.2) then
+*-------------------
+*     t  -> b  ve e+
+*-------------------
+         jd(jvl)=isign*12       !ve or ve~
+         jd(jl) =isign*(-11)    !e+ or e-
+
+      elseif(imode.eq.3) then   
+*--------------------
+*     t  -> b  vm mu+
+*--------------------
+         jd(jvl)=isign*14       !vm or vm~
+         jd(jl) =isign*(-13)    !mu+ or mu-
+
+      elseif(imode.eq.4) then
+*-------------------
+*     t  -> b  vt ta+
+*-------------------
+         if(isign.eq.1) then    !t
+            M3=0D0
+            M4=lmass
+         else                   !t~
+            M3=lmass
+            M4=0d0
+         endif   
+         jd(jvl)=isign*16       !vt or vt~
+         jd(jl) =isign*(-15)    !ta+ or ta-
+
+      elseif(imode.eq.5) then
+*-------------------
+*     t -> b  vl l+  (e+mu) 
+*-------------------
+        if(aa.lt.Half) then
+            jd(jvl)=isign*12    !ve or ve~
+            jd(jl) =isign*(-11) !e+ or e-
+         else
+            jd(jvl)=isign*14    !vm or vm~
+            jd(jl) =isign*(-13) !mu+ or mu-
+         endif
+
+      elseif(imode.eq.6) then
+*-------------------
+*     t -> b  vl l+  (e+mu+ta) 
+*-------------------
+         if(aa.lt.one/Three) then
+            jd(jvl)=isign*12    !ve or ve~
+            jd(jl) =isign*(-11) !e+ or e-
+         elseif(aa.lt.two/three) then
+            jd(jvl)=isign*14    !vm or vm~
+            jd(jl) =isign*(-13) !mu+ or mu-
+         else
+            if(isign.eq.1) then !t
+               M3=0D0
+               M4=lmass
+            else                !t~
+               M3=lmass
+               M4=0d0
+            endif   
+            jd(jvl)=isign*16    !vt or vt~
+            jd(jl) =isign*(-15) !ta+ or ta-
+         endif
+
+      elseif(imode.eq.7) then
+*-------------------
+*     t -> b  j j   (ud,cs) 
+*-------------------
+c--   color
+         icol(1,3)=cindex   !position 3 is always a particle  
+         icol(2,3)=0
+         icol(1,4)=0        !position 4 is always a anti-particle
+         icol(2,4)=cindex
+
+         if(aa.lt..5d0) then
+            jd(jvl)=isign*2    !u  or u~
+            jd(jl) =isign*(-1) !d~ or d
+         else
+            if(isign.eq.1) then !t
+               M3=cmass
+               M4=0d0
+            else                !t~
+               M3=0d0
+               M4=cmass
+            endif   
+            jd(jvl)=isign*4     !c  or c~
+            jd(jl) =isign*(-3)  !s~ or s
+         endif
+
+      elseif(imode.eq.8) then
+*-------------------
+*     t -> b  anything
+*-------------------
+c
+c     first decide if W decays leptonically or hadronically
+c
+         if(aa1.lt.3d0*br_w_lv) then  !leptonic decay
+            hadro_dec=.false.
+            if(aa.lt.one/Three) then
+               jd(jvl)=isign*12 !ve or ve~
+               jd(jl) =isign*(-11) !e+ or e-
+            elseif(aa.lt.two/three) then
+               jd(jvl)=isign*14 !vm or vm~
+               jd(jl) =isign*(-13) !mu+ or mu-
+            else
+               if(isign.eq.1) then !t
+                  M3=0D0
+                  M4=lmass
+               else             !t~
+                  M3=lmass
+                  M4=0d0
+               endif   
+               jd(jvl)=isign*16 !vt or vt~
+               jd(jl) =isign*(-15) !ta+ or ta-
+            endif
+            
+         else                   ! hadronic decay
+            hadro_dec=.true.
+c--   color
+            icol(1,3)=cindex    !position 3 is always a particle  
+            icol(2,3)=0
+            icol(1,4)=0         !position 4 is always a anti-particle
+            icol(2,4)=cindex
+            
+            if(aa.lt..5d0) then
+               jd(jvl)=isign*2  !u  or u~
+               jd(jl) =isign*(-1) !d~ or d
+            else
+               if(isign.eq.1) then !t
+                  M3=cmass
+                  M4=0d0
+               else             !t~
+                  M3=0d0
+                  M4=cmass
+               endif   
+               jd(jvl)=isign*4  !c  or c~
+               jd(jl) =isign*(-3) !s~ or s
+            endif
+         endif
+
+      endif !imode(from 2 to 8)
+      endif !imode
+ 
+c--   couplings
+      GXX(1)=gwf(1)
+      GXX(2)=gwf(2)
+c--   color*(bose factor)*number of helicities
+      factor=(3d0/3d0)*1d0*4d0    
+      if(imode.eq.7) factor=factor*3d0 !quark color in jj
+      if(imode.eq.8.and.hadro_dec) factor=factor*3d0 !quark color in jj
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  get_hel(xran1(iseed),2)
+      jhel(4) =  -jhel(3)
+c--   phase space
+      if(GV.eq.0d0) call error_trap("GV")
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      if(isign.eq. 1) call emme_f3f(pd,jhel,emmesq) !t
+      if(isign.eq.-1) call emme_fx3f(pd,jhel,emmesq) !t~
+c--   weight
+      dwgt=pswgt*emmesq*factor
+      if(imode.eq.5) dwgt=dwgt*2d0 !l=e,mu
+      if(imode.eq.6) dwgt=dwgt*3d0 !l=e,mu,tau
+      if(imode.eq.7) dwgt=dwgt*2d0 !jj=ud,cs      
+      if(imode.eq.8) then
+         if(hadro_dec) then
+            dwgt=dwgt*2d0/br_w_jj  !jj=ud,cs      
+         else
+            dwgt=dwgt*3d0/(3d0*br_w_lv)  !l=e,mu,tau
+         endif
+      endif
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+      write(*,*) 'non-existing imode for top decays'
+      endif ! top decays
+
+
+
+      If(abs(ip).eq.15) then
+*-----------------------------------------------------
+*     tau decays
+*-----------------------------------------------------
+
+
+      if(imode.le.3) then 
+*------------------------
+*     ta -> vt vl l
+*------------------------
+
+c--   masses
+      M1=lmass
+      M2=0d0
+      M3=0d0
+      M4=0d0
+      MV=WMASS
+      GV=WWIDTH
+c--   id's
+      if(isign.eq.1) then       !tau-
+         jl=3
+         jvl=4
+      else                      !tau+
+         jl=4
+         jvl=3
+      endif   
+      jd(2)=isign*16   !vt or vt~
+c--   couplings
+      GXX(1)=gwf(1)
+      GXX(2)=gwf(2)
+c--   rnd number 
+      if(imode.eq.3) aa=dble(xran1(iseed))
+
+      if    (imode.eq.1) then
+*-------------------
+*     ta  -> vt  e  ve
+*-------------------
+         jd(jvl)=isign*(-12)   !ve~ or ve
+         jd(jl) =isign*(11)    !e-  or e+
+
+      elseif(imode.eq.2) then   
+*-----------------------
+*     tau  -> vt  mu vm
+*-----------------------
+         jd(jvl)=isign*(-14)    !vm~ or vm
+         jd(jl) =isign*( 13)    !mu- or mu+
+
+      elseif(imode.eq.3) then
+*----------------------
+*     tau -> vt  l vl  (e+mu) 
+*----------------------
+        if(aa.lt.Half) then
+           jd(jvl)=isign*(-12)  !ve~ or ve
+           jd(jl) =isign*(11)   !e-  or e+
+         else
+            jd(jvl)=isign*(-14) !vm~ or vm
+            jd(jl) =isign*( 13) !mu- or mu+
+         endif
+
+      endif   
+
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*1d0    
+c--   helicities: by definition of 3rd and 4th particle 
+      jhel(2) =  -isign
+      jhel(3) =  -1
+      jhel(4) =   1
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      if(isign.eq. 1) call emme_f3f(pd,jhel,emmesq)  !tau-
+      if(isign.eq.-1) call emme_fx3f(pd,jhel,emmesq) !tau+
+c--   weight
+      dwgt=pswgt*emmesq*factor
+      if(imode.eq.3) dwgt=dwgt*2d0 !l=e,mu
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+
+      elseif(imode.eq.4) then 
+*------------------------
+*     tau -> vt pi
+*------------------------
+
+c--   masses
+      M1=lmass
+      M2=0d0
+      M3=0.1395d0
+c--   id's
+      jd(2)=isign*16            !vt or vt~
+      jd(3)=isign*(-211)        !pi-  pi+
+c--   couplings
+      GXX(1)=gwf(1)
+      GXX(2)=gwf(2)
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*2d0    
+c--   fudge to normalize
+      fudge=lwidth*br_ta_pi/0.00373
+      factor=factor*fudge
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  0              !not used
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      if(isign.eq. 1) call emme_ffs (pd,jhel,emmesq)
+      if(isign.eq.-1) call emme_fxfs(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseIf(imode.eq.5) then 
+*------------------------
+*     tau -> vt rho
+*------------------------
+
+c--   masses
+      M1=lmass
+      M2=0d0
+      M3=0.770d0
+c--   id's
+      jd(2)=isign*16            !vt or vt~
+      jd(3)=isign*(-213)        !rho- rho +
+c--   couplings
+      GXX(1)=gwf(1)
+      GXX(2)=gwf(2)
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*6d0    
+c--   fudge to normalize
+      fudge=lwidth*br_ta_ro/0.0182
+      factor=factor*fudge
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  get_hel(xran1(iseed),3)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      if(isign.eq. 1) call emme_ffv (pd,jhel,emmesq)
+      if(isign.eq.-1) call emme_fxfv(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+      endif  !imodes
+      write(*,*) 'non-existing imode for tau decays'
+      endif !tau
+
+
+
+      If(ip.eq.23) then
+*-----------------------------------------------------
+*     Z decays
+*-----------------------------------------------------
+
+c--   masses
+      M1=zmass
+      M2=0D0
+      M3=0D0
+c--   couplings
+      GXX(1)=gzl(1)
+      GXX(2)=gzl(2)
+c--   rnd number for flavour sums
+      aa =dble(xran1(iseed))
+      aa1=dble(xran1(iseed))
+
+
+      if    (imode.eq.1) then
+*-------------------
+*     z->e- e+
+*-------------------
+         jd(2) = 11   
+         jd(3) =-11
+
+      elseif(imode.eq.2) then   
+*--------------------
+*     z->mu- mu+
+*--------------------
+         jd(2) = 13   
+         jd(3) =-13
+
+      elseif(imode.eq.3) then
+*-------------------
+*     z->ta- ta+
+*-------------------
+         M2=lmass
+         jd(2) = 15   
+         jd(3) =-15
+
+      elseif(imode.eq.4) then
+*-------------------
+*     z->e- e+,mu-mu+
+*-------------------
+         if(aa.lt.Half) then
+            jd(2) = 11   
+            jd(3) =-11
+         else
+            jd(2) = 13   
+            jd(3) =-13
+         endif
+
+      elseif(imode.eq.5) then
+*-------------------
+*     z->e- e+,mu-mu+,ta-ta+
+*-------------------
+         if(aa.lt.one/Three) then
+            jd(2) = 11   
+            jd(3) =-11
+         elseif(aa.lt.two/three) then
+            jd(2) = 13   
+            jd(3) =-13
+         else
+            M2=lmass
+            jd(2) = 15   
+            jd(3) =-15
+         endif
+
+      elseif(imode.eq.6) then
+*------------------------
+*     z->vl vl~
+*------------------------
+
+c--   couplings
+         GXX(1)=gzn(1)
+         GXX(2)=gzn(2)
+
+         if(aa.lt.one/Three) then
+            jd(2) = 12   
+            jd(3) =-12
+         elseif(aa.lt.two/three) then
+            jd(2) = 14   
+            jd(3) =-14
+         else
+            jd(2) = 16   
+            jd(3) =-16
+         endif
+
+
+      elseif(imode.eq.7) then
+*------------------------
+*     z-> b b~
+*------------------------
+c--   color
+         icol(1,2)=cindex   !position 2 is always a particle  
+         icol(2,2)=0
+         icol(1,3)=0        !position 3 is always a anti-particle
+         icol(2,3)=cindex
+
+c--   couplings
+         GXX(1)=gzd(1)
+         GXX(2)=gzd(2)
+         M2=bmass
+         jd(2) = 5  
+         jd(3) =-5
+
+      elseif(imode.eq.8) then
+*------------------------
+*     z-> c c~
+*------------------------
+c--   color
+         icol(1,2)=cindex   !position 2 is always a particle  
+         icol(2,2)=0
+         icol(1,3)=0        !position 3 is always a anti-particle
+         icol(2,3)=cindex
+
+c--   couplings
+         GXX(1)=gzu(1)
+         GXX(2)=gzu(2)
+         M2=cmass
+         jd(2) = 4  
+         jd(3) =-4
+
+
+      elseif(imode.eq.9) then 
+*------------------------
+*     z-> u,d,c,s
+*------------------------
+c--   color
+         icol(1,2)=cindex   !position 2 is always a particle  
+         icol(2,2)=0
+         icol(1,3)=0        !position 3 is always a anti-particle
+         icol(2,3)=cindex
+
+c-- probability of a uc or ds decay
+         
+         xprob1=w_z_uu/(w_z_dd + w_z_uu)
+
+         if(aa1.lt.xprob1) then ! decay into ups
+            multi_decay=1
+            
+            if(aa.lt. .5d0) then !u
+               M2=0d0
+               jd(2) = 2  
+               jd(3) =-2
+               GXX(1)=gzu(1)
+               GXX(2)=gzu(2)
+            else                !c
+               jd(2) = 4  
+               jd(3) =-4
+               M2=cmass
+               GXX(1)=gzu(1)
+               GXX(2)=gzu(2)
+            endif
+
+         else                   !decay into downs
+            multi_decay=2
+            if(aa.lt..5d0) then !d
+               M2=0d0
+               jd(2) = 1  
+               jd(3) =-1
+               GXX(1)=gzd(1)
+               GXX(2)=gzd(2)
+            else                !s
+               jd(2) = 3  
+               jd(3) =-3
+               M2=0d0
+               GXX(1)=gzd(1)
+               GXX(2)=gzd(2)
+            endif
+         
+         endif                  !which decay: in ups or downs
+
+      elseif(imode.eq.10) then 
+*------------------------
+*     z-> u,d,c,s,b
+*------------------------
+c--   color
+         icol(1,2)=cindex   !position 2 is always a particle  
+         icol(2,2)=0
+         icol(1,3)=0        !position 3 is always a anti-particle
+         icol(2,3)=cindex
+
+c-- probability of a uc or ds decay
+         
+         xprob1=2d0*w_z_uu/(2d0*w_z_dd + 2d0*w_z_uu+w_z_bb)
+         xprob2=xprob1+2d0*w_z_dd/(2d0*w_z_dd + 2d0*w_z_uu+w_z_bb)
+         
+         if(aa1.lt.xprob1) then ! decay into ups
+            multi_decay=1
+            if(aa.lt.0.5d0) then !u
+               jd(2) = 2  
+               jd(3) =-2
+               M2=0d0
+               GXX(1)=gzu(1)
+               GXX(2)=gzu(2)
+            else                !c 
+               jd(2) = 4  
+               jd(3) =-4
+               M2=cmass
+               GXX(1)=gzu(1)
+               GXX(2)=gzu(2)
+            endif
+         
+         elseif(aa1.lt.xprob2) then ! decay into downs
+            multi_decay=2
+            if(aa.lt. .5d0) then !d
+               M2=0d0
+               jd(2) = 1  
+               jd(3) =-1
+               GXX(1)=gzd(1)
+               GXX(2)=gzd(2)
+            else                !s      
+               jd(2) = 3  
+               jd(3) =-3
+               M2=0d0
+               GXX(1)=gzd(1)
+               GXX(2)=gzd(2)
+            endif
+         else                   ! decay into b's
+            multi_decay=3
+            jd(2) = 5  
+            jd(3) =-5
+            M2=bmass
+            GXX(1)=gzd(1)
+            GXX(2)=gzd(2)
+         endif   
+         
+      else
+         write(*,*) 'non-existing imode for Z decays'
+      endif                     ! imode
+      
+      M3=M2
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*4d0    
+      if(imode.ge.7) factor=factor*3d0 !quark color in jj
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  get_hel(xran1(iseed),2)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_vff(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+      if(imode.eq.4)  dwgt=dwgt*2d0 !l=e,mu
+      if(imode.eq.5)  dwgt=dwgt*3d0 !l=e,mu,tau
+      if(imode.eq.6)  dwgt=dwgt*3d0 !vl=ve,vm,vt
+c
+      if(imode.eq.9) then       !j=u,d,c,s
+         if(multi_decay.eq.1) then
+            dwgt=dwgt*2d0/xprob1 !j=u,c 
+         else
+            dwgt=dwgt*2d0/(1.-xprob1) !j=d,s
+         endif
+      endif   
+c     
+      if(imode.eq.10) then      !j=u,d,c,s,b
+         if(multi_decay.eq.1) then
+            dwgt=dwgt*2d0/xprob1 !j=u,c 
+         elseif(multi_decay.eq.2) then 
+            dwgt=dwgt*2d0/(xprob2-xprob1) !j=d,s
+         elseif(multi_decay.eq.3) then 
+            dwgt=dwgt/(1.-xprob2) !j=b
+         endif 
+      endif   
+
+      return
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      endif ! Z
+
+
+      If(abs(ip).eq.24) then
+*-----------------------------------------------------
+*     W decays
+*-----------------------------------------------------
+
+c--   masses
+      M1=wmass
+      M2=0d0
+      M3=0d0
+c--   couplings
+      GXX(1)=gwf(1)
+      GXX(2)=gwf(2)
+c--   id's
+      if(isign.eq.1) then       !W+
+         jl=3
+         jvl=2
+      else                      !W-
+         jl=2
+         jvl=3
+      endif   
+c--   rnd number 
+      if(imode.ge.4) aa =dble(xran1(iseed))
+      if(imode.eq.7) aa1=dble(xran1(iseed))
+
+
+      if(imode.eq.1) then 
+*------------------------
+*     w-> e ve
+*------------------------
+         jd(jvl)=isign*12       !ve or ve~
+         jd(jl) =isign*(-11)    !e+ or e-
+
+      elseif(imode.eq.2) then   
+*--------------------
+*     w -> mu vm
+*--------------------
+         jd(jvl)=isign*14       !vm or vm~
+         jd(jl) =isign*(-13)    !mu+ or mu-
+
+      elseif(imode.eq.3) then
+*-------------------
+*     w -> ta vt
+*-------------------
+         if(isign.eq.1) then    !t
+            M2=0D0
+            M3=lmass
+         else                   !t~
+            M2=lmass
+            M3=0d0
+         endif   
+         jd(jvl)=isign*16       !vt or vt~
+         jd(jl) =isign*(-15)    !ta+ or ta-
+
+      elseif(imode.eq.4) then
+*-------------------
+*     w -> vl l+  (e+mu) 
+*-------------------
+        if(aa.lt.Half) then
+            jd(jvl)=isign*12    !ve or ve~
+            jd(jl) =isign*(-11) !e+ or e-
+         else
+            jd(jvl)=isign*14    !vm or vm~
+            jd(jl) =isign*(-13) !mu+ or mu-
+         endif
+
+      elseif(imode.eq.5) then
+*-------------------
+*     w -> vl l+  (e+mu+ta) 
+*-------------------
+         if(aa.lt.one/Three) then
+            jd(jvl)=isign*12    !ve or ve~
+            jd(jl) =isign*(-11) !e+ or e-
+         elseif(aa.lt.two/three) then
+            jd(jvl)=isign*14    !vm or vm~
+            jd(jl) =isign*(-13) !mu+ or mu-
+         else
+            if(isign.eq.1) then !t
+               M2=0D0
+               M3=lmass
+            else                !t~
+               M2=lmass
+               M3=0d0
+            endif   
+            jd(jvl)=isign*16    !vt or vt~
+            jd(jl) =isign*(-15) !ta+ or ta-
+         endif
+
+      elseif(imode.eq.6) then
+*-------------------
+*     w -> j j   (ud,cs) 
+*-------------------
+c--   color
+         icol(1,2)=cindex   !position 2 is always a particle  
+         icol(2,2)=0
+         icol(1,3)=0        !position 3 is always a anti-particle
+         icol(2,3)=cindex
+
+         if(aa.lt..5d0) then
+            jd(jvl)=isign*2    !u  or u~
+            jd(jl) =isign*(-1) !d~ or d
+         else
+            if(isign.eq.1) then !t
+               M2=cmass
+               M3=0d0
+            else                !t~
+               M2=0d0
+               M3=cmass
+            endif   
+            jd(jvl)=isign*4     !c  or c~
+            jd(jl) =isign*(-3)  !s~ or s
+         endif
+
+      elseif(imode.eq.7) then
+*-------------------
+*     w -> anything
+*-------------------
+c     
+c     first decide if W decays leptonically or hadronically
+c     
+         if(aa1.lt.3d0*br_w_lv) then !leptonic decay
+            hadro_dec=.false.
+            if(aa.lt.one/Three) then
+               jd(jvl)=isign*12 !ve or ve~
+               jd(jl) =isign*(-11) !e+ or e-
+            elseif(aa.lt.two/three) then
+               jd(jvl)=isign*14 !vm or vm~
+               jd(jl) =isign*(-13) !mu+ or mu-
+            else
+               if(isign.eq.1) then !t
+                  M2=0D0
+                  M3=lmass
+               else             !t~
+                  M2=lmass
+                  M3=0d0
+               endif   
+               jd(jvl)=isign*16 !vt or vt~
+               jd(jl) =isign*(-15) !ta+ or ta-
+            endif
+         else
+            hadro_dec=.true.    !hadronic decay
+c--   color
+            icol(1,2)=cindex    !position 2 is always a particle  
+            icol(2,2)=0
+            icol(1,3)=0         !position 3 is always a anti-particle
+            icol(2,3)=cindex
+            
+            if(aa.lt..5d0) then
+               jd(jvl)=isign*2  !u  or u~
+               jd(jl) =isign*(-1) !d~ or d
+            else
+               if(isign.eq.1) then !t
+                  M2=cmass
+                  M3=0d0
+               else             !t~
+                  M2=0d0
+                  M3=cmass
+               endif   
+               jd(jvl)=isign*4  !c  or c~
+               jd(jl) =isign*(-3) !s~ or s
+            endif            
+         endif
+      endif
+c
+c     done with the modes: now I start the calculation
+c      
+     
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*2d0    
+      if(imode.eq.6) factor=factor*3d0 !quark color in jj
+      if(imode.eq.7.and.hadro_dec) factor=factor*3d0 !quark color in jj
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  -jhel(2)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_vff(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+      if(imode.eq.4)  dwgt=dwgt*2d0 !l=e,mu
+      if(imode.eq.5)  dwgt=dwgt*3d0 !l=e,mu,tau
+      if(imode.eq.6)  dwgt=dwgt*2d0 !j=ud+cs
+      if(imode.eq.7) then
+         if(hadro_dec) then
+            dwgt=dwgt*2d0/br_w_jj        !jj=ud,cs      
+         else
+            dwgt=dwgt*3d0/(3d0*br_w_lv)  !l=e,mu,tau
+         endif
+      endif
+
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+
+      return
+      endif ! w
+
+
+      If(ip.eq.25) then
+*-----------------------------------------------------
+*     Higgs decays
+*-----------------------------------------------------
+
+      If(imode.eq.1) then 
+*------------------------
+*     h->b b~
+*------------------------
+c--   masses
+      M1=hmass
+      M2=bmass
+      M3=bmass
+c--   id's
+      jd(2) = 5   
+      jd(3) =-5
+c--   color
+      icol(1,2)=cindex          !position 2 is always a particle  
+      icol(2,2)=0
+      icol(1,3)=0               !position 3 is always a anti-particle
+      icol(2,3)=cindex
+c--   couplings
+      GXX(1)=ghbot(1)
+      GXX(2)=ghbot(2)
+c--   color*(bose factor)*number of helicities
+      factor=3d0*1d0*2d0    
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  jhel(2)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_hff(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.2) then 
+*------------------------
+*     h->ta- ta+
+*------------------------
+c--   masses
+      M1=hmass
+      M2=lmass
+      M3=lmass
+c--   id's
+      jd(2) = 15   
+      jd(3) =-15
+c--   couplings
+      GXX(1)=ghtau(1)
+      GXX(2)=ghtau(2)
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*2d0    
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  jhel(2)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_hff(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+
+      return
+
+      elseif(imode.eq.3) then 
+*------------------------
+*     h->mu- mu+
+*------------------------
+c--   masses
+      M1=hmass
+      M2=0d0
+      M3=0d0
+c--   id's
+      jd(2) = 13   
+      jd(3) =-13
+c--   couplings
+      GXX(1)=ghtau(1)/lmass*0.105658389d0
+      GXX(2)=ghtau(2)/lmass*0.105658389d0
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*2d0    
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  jhel(2)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_hff(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.4) then 
+*------------------------
+*     h-> c c~
+*------------------------
+c--   masses
+      M1=hmass
+      M2=cmass
+      M3=cmass
+c--   id's
+      jd(2) = 4   
+      jd(3) =-4
+c--   color
+      icol(1,2)=cindex          !position 2 is always a particle  
+      icol(2,2)=0
+      icol(1,3)=0               !position 3 is always a anti-particle
+      icol(2,3)=cindex
+c--   couplings
+      GXX(1)=ghcha(1)
+      GXX(2)=ghcha(2)
+c--   color*(bose factor)*number of helicities
+      factor=3d0*1d0*4d0    
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  get_hel(xran1(iseed),2)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_hff(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.5) then 
+*------------------------
+*     h-> t t~
+*------------------------
+c--   masses
+      M1=hmass
+      M2=tmass
+      M3=tmass
+c--   id's
+      jd(2) = 6   
+      jd(3) =-6
+c--   couplings
+      GXX(1)=ghtop(1)
+      GXX(2)=ghtop(2)
+c--   color*(bose factor)*number of helicities
+      factor=3d0*1d0*4d0    
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  get_hel(xran1(iseed),2)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_hff(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.6) then 
+*------------------------
+*     h-> g g 
+*------------------------
+c--   masses
+      M1=hmass
+      M2=0d0
+      M3=0d0
+c--   id's
+      jd(2) = 21   
+      jd(3) = 21
+c--   color
+      icol(1,2)=cindex 
+      icol(2,2)=cindex+1
+      icol(1,3)=cindex+1      
+      icol(2,3)=cindex
+c--   couplings
+      GX=cmplx(1d0)
+c--   color*(bose factor)*number of helicities
+      factor=8d0*.5d0*4d0 
+c--   fudge factor for normalization
+      factor=factor*hmass*2d0*pi*(SMBRG*SMWDTH)
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  get_hel(xran1(iseed),2)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_hvv(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.7) then 
+*------------------------
+*     h-> a a 
+*------------------------
+
+c--   masses
+      M1=hmass
+      M2=0d0
+      M3=0d0
+c--   id's
+      jd(2) = 22   
+      jd(3) = 22
+c--   couplings
+      GX=cmplx(1d0)
+c--   color*(bose factor)*number of helicities
+      factor=1d0*.5d0*4d0    
+c--   fudge factor for normalization
+      factor=factor*hmass*16d0*pi*(SMBRGA*SMWDTH)
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),2)
+      jhel(3) =  get_hel(xran1(iseed),2)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_hvv(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.8) then 
+*------------------------
+*     h-> z a   
+*------------------------
+c--   masses
+      M1=hmass
+      M2=zmass
+      M3=0d0
+c--   id's
+      jd(2) = 23   
+      jd(3) = 22
+c--   couplings
+      GX=cmplx(1d0)
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*6d0    
+c--   fudge factor for normalization
+      factor=factor*SMBRZGA*SMWDTH*
+     &8d0*pi*hmass/(1d0-(zmass/hmass)**2)
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),3)        !-1,0,1
+      jhel(3) =  get_hel(xran1(iseed),2)
+c--   phase space
+      call twomom  (pd(0,1),pd(0,2),pd(0,3),pswgt)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_hvv(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.9) then 
+*------------------------
+*     h-> w w
+*------------------------
+c--   masses
+      M1=hmass
+      M2=wmass
+      M3=wmass
+c--   id's
+      jd(2) = 24   
+      jd(3) =-24
+c--   couplings
+      GX=gwwh
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*9d0    
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),3)
+      jhel(3) =  get_hel(xran1(iseed),3)
+      jhel(2) =  INT(3e0*xran1(iseed)) - 1
+      jhel(3) =  INT(3e0*xran1(iseed)) - 1
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_hvv(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+
+      elseif(imode.eq.10) then 
+*--------------------------------
+*     h -> w*  w -> l  vl  l' vl' (l,l'=e,mu)
+*--------------------------------
+
+c--   masses
+      M1=hmass
+      M2=0d0
+      M3=0d0
+      M4=0d0
+      M5=0d0
+      MV=WMASS
+      GV=WWIDTH
+c--   id's
+      aa1=dble(xran1(iseed))
+      aa2=dble(xran1(iseed))
+c--   W*-
+      if(aa1.lt.half) then
+         jd(2) =  11            !e-  
+         jd(3) = -12            !ve~ 
+      else
+         jd(2) =  13            !mu- 
+         jd(3) = -14            !vm~ 
+      endif
+c--   W*+
+      if(aa2.lt.half) then
+         jd(4) =  12            !ve 
+         jd(5) = -11            !e+ 
+      else
+         jd(4) =  14            !vm  
+         jd(5) = -13            !mu+ 
+      endif
+c--   couplings
+      GX     =gwwh
+      GXX(1) =gwf(1)
+      GXX(2) =gwf(2)
+      GXX1(1)=gwf(1)
+      GXX1(2)=gwf(2)
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0 
+c--   helicities
+      jhel(2)=-1
+      jhel(3)=1
+      jhel(4)=-1
+      jhel(5)=1
+c--   phase space
+      if(GV.eq.0d0) call error_trap("GV")
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_h4f(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   sum of flavours
+      dwgt=dwgt*4d0
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.11) then 
+*--------------------------------
+*     h -> w*  w -> l  vl  l' vl' (l,l'=e,mu,ta)
+*--------------------------------
+
+c--   masses
+      M1=hmass
+      M2=0d0
+      M3=0d0
+      M4=0d0
+      M5=0d0
+      MV=WMASS
+      GV=WWIDTH
+c--   id's
+      aa1=dble(xran1(iseed))
+      aa2=dble(xran1(iseed))
+c--   W*-
+      if(aa1.lt.one/three) then
+         jd(2) =  11            !e-  
+         jd(3) = -12            !ve~ 
+      elseif(aa1.lt.two/three) then
+         jd(2) =  13            !mu- 
+         jd(3) = -14            !vm~ 
+      else
+         M2    =  lmass
+         jd(2) =  15            !ta- 
+         jd(3) = -16            !vt~ 
+      endif
+c--   W*+
+      if(aa2.lt.one/three) then
+         jd(4) =  12            !ve 
+         jd(5) = -11            !e+ 
+      elseif(aa2.lt.two/three) then
+         jd(4) =  14            !vm  
+         jd(5) = -13            !mu+ 
+      else
+         M5    =  lmass
+         jd(4) =  16            !vt 
+         jd(5) = -15            !ta+ 
+      endif
+
+c--   couplings
+      GX     =gwwh
+      GXX(1) =gwf(1)
+      GXX(2) =gwf(2)
+      GXX1(1)=gwf(1)
+      GXX1(2)=gwf(2)
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0    
+c--   helicities
+      jhel(2)=-1
+      jhel(3)=1
+      jhel(4)=-1
+      jhel(5)=1
+c--   phase space
+      if(GV.eq.0d0) call error_trap("GV")
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_h4f(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   sum of flavours
+      dwgt=dwgt*9d0
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.12) then 
+*--------------------------------
+*     h -> w*  w -> j  j   l  vl  (jj=ud,cs;l=e,mu)
+*--------------------------------
+
+c--   masses
+      M1=hmass
+      M2=0d0
+      M3=0d0
+      M4=0d0
+      M5=0d0
+      MV=WMASS
+      GV=WWIDTH
+c--   color
+      icol(1,2)=cindex          !position 2 is always a particle  
+      icol(2,2)=0
+      icol(1,3)=0               !position 3 is always a anti-particle
+      icol(2,3)=cindex
+c--   id's
+      aa1=dble(xran1(iseed))
+      aa2=dble(xran1(iseed))
+      aa3=dble(xran1(iseed))
+
+      if(aa1.lt.half) then
+c--   W*-
+         if(aa2.lt.half) then
+            jd(2) =  1          !d 
+            jd(3) = -2          !u~
+         else
+            M3    =  cmass
+            jd(2) =  3          !s 
+            jd(3) = -4          !c~ 
+         endif         
+c--   W*+
+         if(aa3.lt.half) then
+            jd(4) =  12         !ve 
+            jd(5) = -11         !e+ 
+         else
+            jd(4) =  14         !vm  
+            jd(5) = -13         !mu+ 
+         endif
+
+      else
+c--   W*+
+         if(aa2.lt.half) then
+            jd(2) =  2          !u
+            jd(3) = -1          !d~
+         else
+            M2    =  cmass
+            jd(2) =  4          !c 
+            jd(3) = -3          !s~ 
+         endif         
+c--   W*-
+         if(aa3.lt.half) then
+            jd(4) =  11         !e- 
+            jd(5) = -12         !ve~ 
+         else
+            jd(4) =  13         !mu- 
+            jd(5) = -14         !vm~ 
+         endif
+      endif
+
+c--   couplings
+      GX     =gwwh
+      GXX(1) =gwf(1)
+      GXX(2) =gwf(2)
+      GXX1(1)=gwf(1)
+      GXX1(2)=gwf(2)
+c--   color*(bose factor)*number of helicities
+      factor=3d0*1d0  
+c--   helicities
+      jhel(2)=-1
+      jhel(3)=1
+      jhel(4)=-1
+      jhel(5)=1
+c--   phase space
+      if(GV.eq.0d0) call error_trap("GV")
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_h4f(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   sum of flavours
+      dwgt=dwgt*4d0
+c--   sum of W+ and W- possibilities
+      dwgt=dwgt*2d0
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.13) then 
+*--------------------------------
+*     h -> w*  w -> j  j   l  vl  (jj=ud,cs;l=e,mu,ta)
+*--------------------------------
+
+c--   masses
+      M1=hmass
+      M2=0d0
+      M3=0d0
+      M4=0d0
+      M5=0d0
+      MV=WMASS
+      GV=WWIDTH
+c--   color
+      icol(1,2)=cindex          !position 2 is always a particle  
+      icol(2,2)=0
+      icol(1,3)=0               !position 3 is always a anti-particle
+      icol(2,3)=cindex
+c--   id's
+      aa1=dble(xran1(iseed))
+      aa2=dble(xran1(iseed))
+      aa3=dble(xran1(iseed))
+
+      if(aa1.lt.half) then
+c--   W*-
+         if(aa2.lt.half) then
+            jd(2) =  1          !d 
+            jd(3) = -2          !u~
+         else
+            M3    =  cmass
+            jd(2) =  3          !s 
+            jd(3) = -4          !c~ 
+         endif         
+c--   W*+
+         if(aa3.lt.one/three) then
+            jd(4) =  12         !ve 
+            jd(5) = -11         !e+ 
+         elseif(aa2.lt.two/three) then
+            jd(4) =  14         !vm  
+            jd(5) = -13         !mu+ 
+         else
+            M5    =  lmass
+            jd(4) =  16         !vt 
+            jd(5) = -15         !ta+ 
+         endif
+
+      else
+c--   W*+
+         if(aa2.lt.half) then
+            jd(2) =  2          !u
+            jd(3) = -1          !d~
+         else
+            M2    =  cmass
+            jd(2) =  4          !c 
+            jd(3) = -3          !s~ 
+         endif         
+c--   W*-
+         if(aa3.lt.one/three) then
+            jd(4) =  11         !e-  
+            jd(5) = -12         !ve~ 
+         elseif(aa3.lt.two/three) then
+            jd(4) =  13         !mu- 
+            jd(5) = -14         !vm~ 
+         else
+            M4    =  lmass
+            jd(4) =  15         !ta- 
+            jd(5) = -16         !vt~ 
+         endif
+      endif
+
+c--   couplings
+      GX     =gwwh
+      GXX(1) =gwf(1)
+      GXX(2) =gwf(2)
+      GXX1(1)=gwf(1)
+      GXX1(2)=gwf(2)
+c--   color*(bose factor)*number of helicities
+      factor=3d0*1d0   
+c--   helicities
+      jhel(2)=-1
+      jhel(3)=1
+      jhel(4)=-1
+      jhel(5)=1
+c--   phase space
+      if(GV.eq.0d0) call error_trap("GV")
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_h4f(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   sum of flavours
+      dwgt=dwgt*6d0
+c--   sum of W+ and W- possibilities
+      dwgt=dwgt*2d0
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.14) then 
+*------------------------
+*     h-> z z 
+*------------------------
+c--   masses
+      M1=hmass
+      M2=zmass
+      M3=zmass
+c--   id's
+      jd(2) = 23  
+      jd(3) = 23
+c--   couplings
+      GX=gzzh
+c--   color*(bose factor)*number of helicities
+      factor=1d0*.5d0*9d0    
+c--   helicities
+      jhel(2) =  get_hel(xran1(iseed),3)
+      jhel(3) =  get_hel(xran1(iseed),3)
+c--   phase space
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_hvv(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.15) then 
+*---------------------------------
+*     h -> z*  z -> l- l+  l-' l+'(l,l'=e,mu)
+*---------------------------------
+
+c--   masses
+      M1=hmass
+      M2=0d0
+      M3=0d0
+      M4=0d0
+      M5=0d0
+      MV=ZMASS
+      GV=ZWIDTH
+c--   id's
+      aa1=dble(xran1(iseed))
+      aa2=dble(xran1(iseed))
+c     Z1
+      if(aa1.lt.half) then
+         jd(2) =  11            !e-  
+         jd(3) = -11            !e+
+      else
+         jd(2) =  13            !mu- 
+         jd(3) = -13            !mu+ 
+      endif
+c     Z2
+      if(aa2.lt.half) then
+         jd(4) =  11            !e- 
+         jd(5) = -11            !e+ 
+      else
+         jd(4) =  13            !mu-  
+         jd(5) = -13            !mu+ 
+      endif
+c--   couplings
+      GX    =gzzh
+      GXX(1)=gzl(1)
+      GXX(2)=gzl(2)
+      GXX1(1)=gzl(1)
+      GXX1(2)=gzl(2)
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*4d0/2d0   
+c--   helicities
+      jhel(2)=+1
+      if(xran1(iseed) .gt. .5e0) jhel(2)=-1
+      jhel(3)=-jhel(2)
+      jhel(4)=+1
+      if(xran1(iseed) .gt. .5e0) jhel(4)=-1
+      jhel(5)=-jhel(4)
+c--   phase space
+      if(GV.eq.0d0) call error_trap("GV")
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_h4f(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   sum of flavours
+      dwgt=dwgt*4d0
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+      
+
+      elseif(imode.eq.16) then 
+*---------------------------------
+*     h -> z*  z -> l- l+  l-' l+'(l,l'=e,mu,ta)
+*---------------------------------
+
+c--   masses
+      M1=hmass
+      M2=0d0
+      M3=0d0
+      M4=0d0
+      M5=0d0
+      MV=ZMASS
+      GV=ZWIDTH
+c--   id's
+      aa1=dble(xran1(iseed))
+      aa2=dble(xran1(iseed))
+c     Z
+      if(aa1.lt.one/three) then
+         jd(2) =  11            !e-  
+         jd(3) = -11            !e+
+      elseif(aa1.lt.two/three) then
+         jd(2) =  13            !mu- 
+         jd(3) = -13            !mu+ 
+      else
+         M2    =  lmass
+         M3    =  lmass
+         jd(2) =  15            !ta- 
+         jd(3) = -15            !ta+ 
+      endif
+c     Z
+      if(aa2.lt.one/three) then
+         jd(4) =  11            !e-
+         jd(5) = -11            !e+ 
+      elseif(aa2.lt.two/three) then
+         jd(4) =  13            !mu-  
+         jd(5) = -13            !mu+ 
+      else
+         M4    =  lmass
+         M5    =  lmass
+         jd(4) =  15            !ta-
+         jd(5) = -15            !ta+ 
+      endif
+c--   couplings
+      GX    =gzzh
+      GXX(1)=gzl(1)
+      GXX(2)=gzl(2)
+      GXX1(1)=gzl(1)
+      GXX1(2)=gzl(2)
+c--   color*(bose factor)*number of helicities
+      factor=1d0*1d0*4d0/2d0    
+c--   helicities
+      jhel(2)=+1
+      if(xran1(iseed) .gt. .5e0) jhel(2)=-1
+      jhel(3)=-jhel(2)
+      jhel(4)=+1
+      if(xran1(iseed) .gt. .5e0) jhel(4)=-1
+      jhel(5)=-jhel(4)
+c--   phase space
+      if(GV.eq.0d0) call error_trap("GV")
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_h4f(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   sum of flavours
+      dwgt=dwgt*9d0
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.17) then 
+*---------------------------------
+*     h -> z*  z -> j  j~  l-  l+ (j=u,d,c,s;l=e,mu )
+*---------------------------------
+
+c--   masses
+      M1=hmass
+      M2=0d0
+      M3=0d0
+      M4=0d0
+      M5=0d0
+      MV=ZMASS
+      GV=ZWIDTH
+c--   color
+      icol(1,2)=cindex          !position 2 is always a particle  
+      icol(2,2)=0
+      icol(1,3)=0               !position 3 is always a anti-particle
+      icol(2,3)=cindex
+c--   couplings
+      GX    =gzzh
+      GXX(1)=gzl(1)
+      GXX(2)=gzl(2)
+      GXX1(1)=gzl(1)
+      GXX1(2)=gzl(2)
+c--   id's
+      aa =dble(xran1(iseed)) 
+      aa1=dble(xran1(iseed))
+      aa2=dble(xran1(iseed))
+c     Z1
+c-- probability of a uc or ds decay
+         
+         xprob1=w_z_uu/(w_z_dd + w_z_uu)
+
+         if(aa1.lt.xprob1) then ! decay into ups
+            multi_decay=1
+            
+            if(aa.lt. .5d0) then !u
+               M2=0d0
+               jd(2) = 2  
+               jd(3) =-2
+               GXX(1)=gzu(1)
+               GXX(2)=gzu(2)
+            else                !c
+               jd(2) = 4  
+               jd(3) =-4
+               M2=cmass
+               GXX(1)=gzu(1)
+               GXX(2)=gzu(2)
+            endif
+
+         else                   !decay into downs
+            multi_decay=2
+            if(aa.lt..5d0) then !d
+               M2=0d0
+               jd(2) = 1  
+               jd(3) =-1
+               GXX(1)=gzd(1)
+               GXX(2)=gzd(2)
+            else                !s
+               jd(2) = 3  
+               jd(3) =-3
+               M2=0d0
+               GXX(1)=gzd(1)
+               GXX(2)=gzd(2)
+            endif
+         
+         endif                  !which decay: in ups or downs
+      M3=M2
+c     Z2
+      if(aa2.lt.half) then
+         jd(4) =  11            !e- 
+         jd(5) = -11            !e+ 
+      else
+         jd(4) =  13            !mu-  
+         jd(5) = -13            !mu+ 
+      endif
+
+c--   color*(bose factor)*number of helicities
+      factor=3d0*1d0*4d0/2d0  
+c--   helicities
+      jhel(2)=+1
+      if(xran1(iseed) .gt. .5e0) jhel(2)=-1
+      jhel(3)=-jhel(2)
+      jhel(4)=+1
+      if(xran1(iseed) .gt. .5e0) jhel(4)=-1
+      jhel(5)=-jhel(4)
+c--   phase space
+      if(GV.eq.0d0) call error_trap("GV")
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_h4f(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   sum of flavours
+      dwgt=dwgt*2d0
+      if(multi_decay.eq.1) then
+         dwgt=dwgt*2d0/xprob1   !j=u,c 
+      else
+         dwgt=dwgt*2d0/(1.-xprob1) !j=d,s
+      endif
+c--   sum of two possibilities for the decay of a Z
+      dwgt=dwgt*2d0
+c--   check that dwgt is a reasonable number
+      call check_nan(dwgt)
+      return
+
+      elseif(imode.eq.18) then 
+*---------------------------------
+*     h -> z*  z -> j  j~  l-  l+ (j=u,d,c,s;l=e,mu,ta)
+*---------------------------------
+
+c--   masses
+      M1=hmass
+      M2=0d0
+      M3=0d0
+      M4=0d0
+      M5=0d0
+      MV=ZMASS
+      GV=ZWIDTH
+c--   color
+      icol(1,2)=cindex          !position 2 is always a particle  
+      icol(2,2)=0
+      icol(1,3)=0               !position 3 is always a anti-particle
+      icol(2,3)=cindex
+c--   couplings
+      GX    =gzzh
+      GXX(1)=gzl(1)
+      GXX(2)=gzl(2)
+      GXX1(1)=gzl(1)
+      GXX1(2)=gzl(2)
+c--   id's
+      aa =dble(xran1(iseed))
+      aa1=dble(xran1(iseed))
+      aa2=dble(xran1(iseed))
+c     Z1
+      xprob1=w_z_uu/(w_z_dd + w_z_uu)
+      if(aa1.lt.xprob1) then    ! decay into ups
+         multi_decay=1
+         
+         if(aa.lt. .5d0) then   !u
+            M2=0d0
+            jd(2) = 2  
+            jd(3) =-2
+            GXX(1)=gzu(1)
+            GXX(2)=gzu(2)
+         else                   !c
+            jd(2) = 4  
+            jd(3) =-4
+            M2=cmass
+            GXX(1)=gzu(1)
+            GXX(2)=gzu(2)
+         endif
+         
+      else                      !decay into downs
+         multi_decay=2
+         if(aa.lt..5d0) then    !d
+            M2=0d0
+            jd(2) = 1  
+            jd(3) =-1
+            GXX(1)=gzd(1)
+            GXX(2)=gzd(2)
+         else                   !s
+            jd(2) = 3  
+            jd(3) =-3
+            M2=0d0
+            GXX(1)=gzd(1)
+            GXX(2)=gzd(2)
+         endif
+         
+      endif                     !which decay: in ups or downs
+      M3=M2
+c     Z2
+      if(aa2.lt.one/three) then
+         jd(4) =  11            !e-
+         jd(5) = -11            !e+ 
+      elseif(aa2.lt.two/three) then
+         jd(4) =  13            !mu-  
+         jd(5) = -13            !mu+ 
+      else
+         M4    =  lmass
+         M5    =  lmass
+         jd(4) =  15            !ta-
+         jd(5) = -15            !ta+ 
+      endif
+      
+
+c--   color*(bose factor)*number of helicities
+      factor=3d0*1d0*4d0/2d0    
+c--   helicities
+      jhel(2)=+1
+      if(xran1(iseed) .gt. .5e0) jhel(2)=-1
+      jhel(3)=-jhel(2)
+      jhel(4)=+1
+      if(xran1(iseed) .gt. .5e0) jhel(4)=-1
+      jhel(5)=-jhel(4)
+c--   phase space
+      if(GV.eq.0d0) call error_trap("GV")
+      call phasespace(pd,pswgt)
+      if(pswgt.eq.0d0) return
+c--   matrix element
+      call emme_h4f(pd,jhel,emmesq)
+c--   weight
+      dwgt=pswgt*emmesq*factor
+c--   sum of flavours
+      dwgt=dwgt*3d0
+      if(multi_decay.eq.1) then
+         dwgt=dwgt*2d0/xprob1   !j=u,c