diff --git a/Models/Feynrules/python/ufo2peg/particles.py b/Models/Feynrules/python/ufo2peg/particles.py
--- a/Models/Feynrules/python/ufo2peg/particles.py
+++ b/Models/Feynrules/python/ufo2peg/particles.py
@@ -1,611 +1,611 @@
 from __future__ import print_function
 from string import Template
 import os
 import numpy as np
 from FR_Parameters import *
 
 # ignore these, they're in Hw++ already # TODO reset Hw++ settings instead
 SMPARTICLES = {
 
 1:'d',
 2:'u',
 3:'s',
 4:'c',
 5:'b',
 6:'t', # think about this one later
 
 11:'e-',
 12:'nu_e',
 13:'mu-',
 14:'nu_mu',
 15:'tau-',
 16:'nu_tau',
 
 21:'g',
 22:'gamma',
 23:'Z0',
 24:'W+',
 
 -1:'dbar',
 -2:'ubar',
 -3:'sbar',
 -4:'cbar',
 -5:'bbar',
 -6:'tbar',
 
 -11:'e+',
 -12:'nu_ebar',
 -13:'mu+',
 -14:'nu_mubar',
 -15:'tau+',
 -16:'nu_taubar',
 
 -24:'W-'
 }
 
 
 
 particleT = Template(
 """
 create ThePEG::ParticleData $name
 # values set to 999999 are recalculated later from other model parameters
 setup $name $pdg_code $name $mass $width $wcut $ctau $charge $color $spin 0
 """
 )
 
 
 class ParticleConverter:
     'Convert a FR particle to extract the information ThePEG needs.'
     def __init__(self,p,parmsubs,modelparameters):
         self.name = p.name
         self.pdg_code = p.pdg_code
         self.spin = p.spin
         self.color = p.color
         if self.color == 1:
             self.color = 0
         self.selfconjugate = 0
 
         self.mass = parmsubs[str(p.mass)]
         if type(self.mass) == str:
             value = modelparameters[self.mass]
             try:
                 value = value.real
             except:
                 pass
             newname = '%s_ABS' % self.mass
             self.mass = '${%s}' % newname
             modelparameters[newname] = abs(value)
         else:
             try:
                 self.mass = self.mass.real
             except:
                 pass
             self.mass = 999999. # abs(self.mass)
 
         hbarc = 197.3269631e-15 # GeV mm (I hope ;-) )
         self.width = parmsubs[str(p.width)]
         if type(self.width) == str:
             width = modelparameters[self.width]
             ctau = (hbarc / width) if width != 0 else 0
             newname = '%s_CTAU' % self.width
             self.ctau = '${%s}' % newname
             modelparameters[newname] = ctau
 
             wcut = 10 * width
             newname = '%s_WCUT' % self.width
             self.wcut = '${%s}' % newname
             modelparameters[newname] = wcut
 
             self.width = '${%s}' % self.width
         else:
             self.ctau = 999999. # (hbarc / self.width) if self.width != 0 else 0
             self.wcut = 999999. #10.0 * self.width
             self.width = 999999. # was blank line before
 
         self.charge = int(3 * p.charge)
 
     def subs(self):
         return self.__dict__
 
 def check_effective_vertex(FR,p,ig) :
     for vertex in FR.all_vertices:
         if(len(vertex.particles) != 3) : continue
         if(p not in vertex.particles ) : continue
         ng=0
         for part in vertex.particles :
             if(part.pdg_code==ig) : ng+=1
         if(ng==2) :
             return False
     return True
 
 # finds all dim-Dimention vectices in FR that involve pIn
 def sort_vertices(FR,pIn,dim):
     possibleVertices = []
     for V in FR.all_vertices:
         # only keep 1 -> dim-1 vetrices
         if len(V.particles) != dim :
             continue
         # keep vertices with pIn
         if pIn == str(V.particles[0]) or pIn == str(V.particles[1]) or pIn == str(V.particles[2]):
             possibleVertices.append(V)
     return possibleVertices
 
 # extracts all possible splittings for incoming particle p
 def sort_splittings(FR,Vertices,p):
     pSplittings = []
     for V in Vertices:
         lorz = V.lorentz
         # calculate the total coupling value for this splitting
         coup = V.couplings
         keys = coup.keys()
         coupling_value = [0.,0.]
         for k in keys :
             color_idx, lorentz_idx = k
             # https://arxiv.org/pdf/2304.09883.pdf
             if 'FFS' in str(lorz[lorentz_idx]):
                 # distinguish CP-even/-odd couplings, each couplings in the square bracket correspond to [1,Gamma5] respectively
                 if str(lorz[lorentz_idx]) == 'FFS1': # ProjM
                     coupling_value[0] += eval(coup[k].value)/2
                     coupling_value[1] -= eval(coup[k].value)/2
                 elif str(lorz[lorentz_idx]) == 'FFS2': # ProjM - ProjP
                     coupling_value[1] -= eval(coup[k].value)
                 elif str(lorz[lorentz_idx]) == 'FFS3': # ProjP
                     coupling_value[0] += eval(coup[k].value)/2
                     coupling_value[1] += eval(coup[k].value)/2
                 elif str(lorz[lorentz_idx]) == 'FFS4': # ProjP + ProjM
                     coupling_value[0] += eval(coup[k].value)
             elif 'FFV' in str(lorz[lorentz_idx]):
                 # distinguish left-/right-couplings, each couplings in the square bracket correspond to [P_L,P_R] respectively
                 if str(lorz[lorentz_idx]) == 'FFV1': # Gamma
                     coupling_value[0] += eval(coup[k].value)
                     coupling_value[1] += eval(coup[k].value)
                 elif str(lorz[lorentz_idx]) == 'FFV2': # Gamma*ProjM
                     coupling_value[0] += eval(coup[k].value)
                 elif str(lorz[lorentz_idx]) == 'FFV3': # Gamma*(ProjM-2*ProjP)
                     coupling_value[0] += eval(coup[k].value)
                     coupling_value[1] -= 2*eval(coup[k].value)
                 elif str(lorz[lorentz_idx]) == 'FFV4': # Gamma*(ProjM+2*ProjP)
                     coupling_value[0] += eval(coup[k].value)
                     coupling_value[1] += 2*eval(coup[k].value)
                 elif str(lorz[lorentz_idx]) == 'FFV5': # Gamma*(ProjM+4*ProjP)
                     coupling_value[0] += eval(coup[k].value)
                     coupling_value[1] += 4*eval(coup[k].value)
             else:
                 coupling_value[0] += eval(coup[k].value)
 
         # extract splitting format as p -> p1, p2
         p0set = False
         p1set = False
         p2set = False
         for particle in V.particles:
             if particle == p and not p0set:
                 p0set = True
             elif not p1set :
                 p1 = particle
                 p1set = True
             elif not p2set :
                 p2 = particle
                 p2set = True
         if not p0set or not p1set or not p2set :
             continue
 
         def isGoldstone(p) :
             def gstest(name):
                 try:
                     return getattr(p,name)
                 except AttributeError:
                     return False
             gsnames = ['goldstone','goldstoneboson','GoldstoneBoson']
             if any(map(gstest, gsnames)):
                 return True
             return False
 
         def isGhost(p) :
             try:
                 getattr(p,'GhostNumber')
             except AttributeError:
                 return False
             return p.GhostNumber != 0
 
         if isGoldstone(p1) or isGoldstone(p2) :
             continue
         if isGhost(p1) or isGhost(p2) :
             continue
         id1 = abs(p1.pdg_code)
         id2 = abs(p2.pdg_code)
         # put the bigger spin last
         if p1.spin > p2.spin :
             p1, p2 = p2, p1
         pp1p2 = [p,p1,p2] + coupling_value
         pp2p1 = [p,p2,p1] + coupling_value
         if pp1p2 not in pSplittings and pp2p1 not in pSplittings:
             pSplittings.append(pp1p2)
     return pSplittings
 
 def extract_mass(FR,Vertex) :
     m = [0.,0.,0.]
     m[0] = Vertex[0].mass.value
     m[1] = Vertex[1].mass.value
     m[2] = Vertex[2].mass.value
     for i in range(len(m)) :
         if isinstance(m[i],str) :
             m[i] = eval(m[i])
     return m
 
 def split_name(Vertex,split=False) :
     p = [Vertex[0].name, Vertex[1].name, Vertex[2].name]
     for i in range(0,3) :
         if Vertex[i].pdg_code in SMPARTICLES :
             p[i] = SMPARTICLES[Vertex[i].pdg_code]
         else :
             p[i] = Vertex[i].name
     if not split :
         splitname = p[0] + p[1] + p[2]
         splitname = splitname.replace("+", "p")
         splitname = splitname.replace("-", "m")
         return splitname
     else :
         return p[0], p[1], p[2]
 
 def isQuark(particle) :
     if abs(particle.pdg_code) >= 1 and abs(particle.pdg_code) <= 6 :
         return True
     else :
         return False
 
 def isLepton(particle) :
     if abs(particle.pdg_code) >= 11 and abs(particle.pdg_code) <= 16 :
         return True
     else :
         return False
 
 def isScalar(particle) :
     if particle.spin==1 :
         return True
     else :
         return False
 
 def isVector(particle) :
     if particle.spin==3 :
         return True
     else :
         return False
 
 def isGVB(particle) :
     if abs(particle.pdg_code)==24 :
         return True
     elif particle.pdg_code==23 :
         return True
     elif particle.pdg_code==22 :
         return True
     elif particle.pdg_code==21 :
         return True
     else :
         return False
 
 def isBSMVB(particle) :
     if particle.spin==3 :
         pdgid=abs(particle.pdg_code)
         if 20 < pdgid and pdgid < 26:
             return False
         else :
             return True
     else :
         return False
 
 def antiparticle(FR,particle) :
     for anti in FR.all_particles :
         if particle.pdg_code == -anti.pdg_code :
             return anti
     return particle
 
 def thepeg_particles(FR,parameters,modelname,modelparameters,forbidden_names,hw_higgs,\
                     enable_bsm_shower,allow_fcnc) :
     plist = []
     antis = {}
     names = []
     splittings = []
     done_splitting_QCD = []
     done_splitting_QED = []
     done_splitting_EW  = []
 
     for p in FR.all_particles:
         if p.spin == -1:
             continue
 
         gsnames = ['goldstone',
                    'goldstoneboson',
                    'GoldstoneBoson']
 
         def gstest(name):
             try:
                 return getattr(p,name)
             except AttributeError:
                 return False
 
         if any(map(gstest, gsnames)):
             continue
 
         if p.pdg_code in SMPARTICLES:
             continue
 
         if p.pdg_code == 25 and not hw_higgs:
             plist.append(
 """
 set /Herwig/Particles/h0:Mass_generator NULL
 set /Herwig/Particles/h0:Width_generator NULL
 rm /Herwig/Masses/HiggsMass
 rm /Herwig/Widths/hWidth
 """
 )
         if p.name in forbidden_names:
             print('RENAMING PARTICLE',p.name,'as ',p.name+'_UFO')
             p.name +="_UFO"
         subs = ParticleConverter(p,parameters,modelparameters).subs()
         if not (p.pdg_code == 25 and hw_higgs) :
             plist.append( particleT.substitute(subs) )
 
         pdg, name = subs['pdg_code'],  subs['name']
         names.append(name)
         if -pdg in antis:
             plist.append( 'makeanti %s %s\n' % (antis[-pdg], name) )
 
         elif not (p.pdg_code == 25 and hw_higgs) :
             plist.append( 'insert /Herwig/NewPhysics/NewModel:DecayParticles 0 %s\n' % name )
             plist.append( 'insert /Herwig/Shower/ShowerHandler:DecayInShower 0 %s #  %s' % (abs(pdg), name) )
             antis[pdg] = name
             selfconjugate = 1
 
         class SkipMe(Exception):
             pass
 
         def spin_name(s):
             spins = { 1 : 'Zero',
                       2 : 'Half',
                       3 : 'One' }
             if s not in spins:
                 raise SkipMe()
             else:
                 return spins[s]
 
         def col_name(c):
             cols = { 3 : 'Triplet',
                      6 : 'Sextet',
                      8 : 'Octet' }
             return cols[c]
 
         try:
             # QCD splitting functions
-            if p.color in [3,6,8] and abs(pdg) not in done_splitting_QCD: # which colors?
+            if p.color in [3,6,8] and abs(pdg) not in done_splitting_QCD and (enable_bsm_shower or p.pdg_code in SMPARTICLES): # which colors?
                 done_splitting_QCD.append(abs(pdg))
                 splitname = '{name}SplitFnQCD'.format(name=p.name)
                 sudname = '{name}SudakovQCD'.format(name=p.name)
                 splittings.append(
 """
 create Herwig::{s}{s}OneSplitFn {name}
 set {name}:InteractionType QCD
 set {name}:ColourStructure {c}{c}Octet
 cp /Herwig/Shower/SudakovCommon {sudname}
 set {sudname}:SplittingFunction {name}
 do /Herwig/Shower/SplittingGenerator:AddFinalSplitting {pname}->{pname},g; {sudname}
 """.format(s=spin_name(p.spin), name=splitname,
            c=col_name(p.color), pname=p.name, sudname=sudname)
                 )
         except SkipMe:
             pass
         # QED splitting functions
         try:
-            if p.charge != 0 and abs(pdg) not in done_splitting_QED:
+            if p.charge != 0 and abs(pdg) not in done_splitting_QED and (enable_bsm_shower or p.pdg_code in SMPARTICLES):
                 done_splitting_QED.append(abs(pdg))
                 splitname = '{name}SplitFnQED'.format(name=p.name)
                 sudname = '{name}SudakovQED'.format(name=p.name)
                 splittings.append(
 """
 create Herwig::{s}{s}OneSplitFn {name}
 set {name}:InteractionType QED
 set {name}:ColourStructure ChargedChargedNeutral
 cp /Herwig/Shower/SudakovCommon {sudname}
 set {sudname}:SplittingFunction {name}
 set {sudname}:Alpha /Herwig/Shower/AlphaQED
 do /Herwig/Shower/SplittingGenerator:AddFinalSplitting {pname}->{pname},gamma; {sudname}
 """.format(s=spin_name(p.spin), name=splitname, pname=p.name, sudname=sudname)
                 )
         except SkipMe:
             pass
 
 
         # EW and BSM splitting functions
         if enable_bsm_shower :
             try:
                 Vertices = sort_vertices(FR,str(p.name),3)
                 pSplittings = sort_splittings(FR,Vertices,p)
                 for Vertex in pSplittings :
                     # do not do anything for CouplingValue < 1e-6
                     if abs(Vertex[3].real) < 1e-6 and abs(Vertex[3].imag) < 1e-6 and abs(Vertex[4].real) < 1e-6 and abs(Vertex[4].imag) < 1e-6 :
                         continue
                     # do not include QCD splittings
                     if Vertex[1].pdg_code == 21 or Vertex[2].pdg_code == 21 :
                         continue
                     # do not include QED splittings
                     if Vertex[2].pdg_code == 22 and abs(Vertex[0].pdg_code) == abs(Vertex[1].pdg_code) or\
                        Vertex[1].pdg_code == 22 and abs(Vertex[0].pdg_code) == abs(Vertex[2].pdg_code):
                         continue
                     # skip lepton vertices
                     if isLepton(Vertex[1]) or isLepton(Vertex[2]) : ###HERE###
                         continue
                     # loop over all possible configurations in the splitting
                     for pos in range(0,3) :
                         # rearrange to all possible cases
                         V=[Vertex[0], Vertex[1], Vertex[2], Vertex[3], Vertex[4]]
                         if pos==0:
                             V[0], V[1], V[2] = Vertex[0], Vertex[1], Vertex[2]
                         elif pos==1:
                             V[0], V[1], V[2] = Vertex[1], Vertex[2], Vertex[0]
                         else :
                             V[0], V[1], V[2] = Vertex[2], Vertex[0], Vertex[1]
                         # don't allow photon as progenitor
                         if V[0].pdg_code == 22 : ###HERE###
                             continue
                         # for a generic splitting m0 < m1+m2, otherwise it's a decay
                         m = extract_mass(FR,V)
                         # filter out m+i*zero
                         for ix in range(len(m)) :
                             if isinstance(m[ix], complex) :
                                 m[ix] = m[ix].real
                         if m[0] > m[1] + m[2] :
                             continue
 
                         """
                         possible EW splittings:
                             - V > V  H   with SM GVBs and BSM (charged and nuetral) Scalars
                             - q > q' H   with SM quarks and BSM (charged and nuetral) Scalars, may include FCNC-inducing splittings
                             - H > H' V   with BSM (charged and nuetral) Scalars and SM GVBs (W,Z,photon)
                             - V > H  H'  with SM GVBs and BSM (charged and nuetral) Scalars
                             - H > H' H'' Higgs splittings
                         """
 
                         if isVector(V[0]) : ###HERE###
                             # allow V > V' H
                             if isVector(V[1]) and isScalar(V[2]) :
                                 pass
                             elif isVector(V[2]) and isScalar(V[1]) :
                                 V[0], V[1], V[2] = V[0], V[2], V[1]
                             # allow V > H H'
                             elif isScalar(V[1]) and isScalar(V[2]):
                                 pass
                             # allow V > V' V''
                             elif isVector(V[1]) and isVector(V[2]):
                                 pass
                             # nothing else with a GVB progenitor
                             else :
                                 continue
                         elif isQuark(V[0]) :
                             # allow q > q' H (including FCNC splittings)
                             if isQuark(V[1]) and isScalar(V[2]) :
                                 pass
                             elif isQuark(V[2]) and isScalar(V[1]) :
                                 V[0], V[1], V[2] = V[0], V[2], V[1]
                             elif isQuark(V[1]) and isVector(V[2]) :
                                 pass
                             elif isQuark(V[2]) and isVector(V[1]) :
                                 V[0], V[1], V[2] = V[0], V[2], V[1]
                             # nothing else with a quark progenitor
                             else :
                                 continue
                         elif isScalar(V[0]) :
                             # allow H > H' H''
                             if isScalar(V[1]) and isScalar(V[2]) :
                                 pass
                             # allow H > H' V
                             elif isScalar(V[1]) and isVector(V[2]) :
                                 pass
                             elif isScalar(V[2]) and isVector(V[1]) :
                                 V[0], V[1], V[2] = V[0], V[2], V[1]
                             # nothing else with a scalar progenitor
                             else :
                                 continue
                         # nothing else
                         else :
                             continue
 
                         # getting the electric charge right
                         if V[0].charge != V[1].charge+V[2].charge :
                             if V[0].pdg_code*V[1].pdg_code < 0 :
                                 if V[0].pdg_code > 0 :
                                     V[1] = antiparticle(FR,V[1])
                                 else :
                                     V[0] = antiparticle(FR,V[0])
                             elif V[0].pdg_code*V[2].pdg_code < 0 :
                                 if V[0].pdg_code > 0 :
                                     V[2] = antiparticle(FR,V[2])
                                 else :
                                     V[0] = antiparticle(FR,V[0])
                         if V[0].charge != V[1].charge+V[2].charge :
                             if abs(V[0].charge-V[1].charge-V[2].charge) < 1e-10 :
                                 pass
                             else :
                                 continue
 
                         # deal with FCNC inducing splittings
                         if not allow_fcnc :
                             if isQuark(V[0]) and isQuark(V[1]) and isScalar(V[2]) and\
                                V[2].charge==0. and V[0].pdg_code!=V[1].pdg_code :
                                 print("Omitting",V[0],"->",V[1],",",V[2], "FCNC-inducing splitting.",\
                                      "Use --allow-fcnc flag if you wish to keep this.")
                                 continue
 
                         # set up the splitting
                         SPname = split_name(V,False)
                         # no scalar > quark, antiquark splitting
                         s = [spin_name(V[0].spin),spin_name(V[1].spin),spin_name(V[2].spin)]
 
                         # If the real the coupling value is small, then ignore
                         if abs(V[3].real) < 1e-6:
                             V[3] = complex(0.,V[3].imag)
                         if abs(V[3].imag) < 1e-6:
                             V[3] = complex(V[3].real,0.)
                         if abs(V[4].real) < 1e-6:
                             V[4] = complex(0.,V[4].imag)
                         if abs(V[4].imag) < 1e-6:
                             V[4] = complex(V[4].real,0.)
 
                         if SPname not in done_splitting_EW and (V[3]!=0. or V[4]!=0.):
                             done_splitting_EW.append(SPname)
                             splitname = '{name}SplitFnEW'.format(name=SPname)
                             sudname = '{name}SudakovEW'.format(name=SPname)
                             p0name, p1name, p2name = split_name(V,True)
                             splittings.append(
 """
 create Herwig::{s0}{s1}{s2}EWSplitFn {name}
 set {name}:InteractionType EW
 set {name}:ColourStructure EW
 """.format(s0=s[0],s1=s[1],s2=s[2],name=splitname)
                             )
                             if s[0]=='Half' and s[1]=='Half' and s[2]=='Zero':
                                 splittings.append(
 """set {name}:CouplingValue.CP0.Im {i}
 set {name}:CouplingValue.CP0.Re {j}
 set {name}:CouplingValue.CP1.Im {k}
 set {name}:CouplingValue.CP1.Re {l}
 """.format(name=splitname,i=V[3].imag,j=V[3].real,k=V[4].imag,l=V[4].real)
                                 )
                             elif s[0]=='Half' and s[1]=='Half' and s[2]=='One':
                                 splittings.append(
 """set {name}:CouplingValue.Left.Im {i}
 set {name}:CouplingValue.Left.Re {j}
 set {name}:CouplingValue.Right.Im {k}
 set {name}:CouplingValue.Right.Re {l}
 """.format(name=splitname,i=V[3].imag,j=V[3].real,k=V[4].imag,l=V[4].real)
                                 )
                             else:
                                 splittings.append(
 """set {name}:CouplingValue.Im {i}
 set {name}:CouplingValue.Re {j}
 """.format(name=splitname,i=V[3].imag,j=V[3].real)
                                 )
                             splittings.append(
 """cp /Herwig/Shower/SudakovCommon {sudname}
 set {sudname}:SplittingFunction {name}
 set {sudname}:Alpha /Herwig/Shower/AlphaEW
 do /Herwig/Shower/SplittingGenerator:AddFinalSplitting {p0}->{p1},{p2}; {sudname}
 """.format(name=splitname,p0=p0name,p1=p1name,p2=p2name,sudname=sudname)
                             )
             except SkipMe:
                 pass
 
 
         if p.charge == 0 and p.color == 1 and p.spin == 1 and not (p.pdg_code == 25 and hw_higgs) :
             if(check_effective_vertex(FR,p,21)) :
                 plist.append(
 """
 insert /Herwig/{ModelName}/V_GenericHGG:Bosons 0 {pname}
 """.format(pname=p.name, ModelName=modelname)
                 )
             if(check_effective_vertex(FR,p,22)) :
                 plist.append(
 """
 insert /Herwig/{ModelName}/V_GenericHPP:Bosons 0 {pname}
 """.format(pname=p.name, ModelName=modelname)
                 )
 
     return ''.join(plist)+''.join(splittings), names
diff --git a/Shower/QTilde/SplittingFunctions/OneZeroZeroEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/OneZeroZeroEWSplitFn.cc
--- a/Shower/QTilde/SplittingFunctions/OneZeroZeroEWSplitFn.cc
+++ b/Shower/QTilde/SplittingFunctions/OneZeroZeroEWSplitFn.cc
@@ -1,213 +1,212 @@
 // -*- C++ -*-
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the OneZeroZeroEWSplitFn class.
 //
 
 #include "OneZeroZeroEWSplitFn.h"
 #include "ThePEG/StandardModel/StandardModelBase.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/PDT/ParticleData.h"
 #include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
 #include "Herwig/Models/StandardModel/SMFFHVertex.h"
 #include "ThePEG/Interface/Parameter.h"
 
 using namespace Herwig;
 
 IBPtr OneZeroZeroEWSplitFn::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr OneZeroZeroEWSplitFn::fullclone() const {
   return new_ptr(*this);
 }
 
 void OneZeroZeroEWSplitFn::persistentOutput(PersistentOStream & os) const {
   os << _couplingValueIm << _couplingValueRe;
 }
 
 void OneZeroZeroEWSplitFn::persistentInput(PersistentIStream & is, int) {
   is >> _couplingValueIm >> _couplingValueRe;
 }
 
 // The following static variable is needed for the type description system in ThePEG.
 DescribeClass<OneZeroZeroEWSplitFn,SplittingFunction>
 describeHerwigOneZeroZeroEWSplitFn("Herwig::OneZeroZeroEWSplitFn", "HwShower.so");
 
 
 void OneZeroZeroEWSplitFn::Init() {
 
   static ClassDocumentation<OneZeroZeroEWSplitFn> documentation
     ("The OneZeroZeroEWSplitFn class implements purly beyond SM electroweak splittings V->HH'");
 
   static Parameter<OneZeroZeroEWSplitFn, double> interfaceCouplingValueIm
     ("CouplingValue.Im",
      "The numerical value (imaginary part) of the splitting coupling to be imported for BSM splittings",
      &OneZeroZeroEWSplitFn::_couplingValueIm, 0.0, -1.0E6, +1.0E6,
      false, false, Interface::limited);
 
   static Parameter<OneZeroZeroEWSplitFn, double> interfaceCouplingValueRe
     ("CouplingValue.Re",
      "The numerical value (real part) of the splitting coupling to be imported for BSM splittings",
      &OneZeroZeroEWSplitFn::_couplingValueRe, 0.0, -1.0E6, +1.0E6,
      false, false, Interface::limited);
 
 }
 
 
 void OneZeroZeroEWSplitFn::doinit() {
   SplittingFunction::doinit();
 }
 
 
 void OneZeroZeroEWSplitFn::getCouplings(Complex & g, const IdList &) const {
   g = Complex(_couplingValueRe,_couplingValueIm);
 }
 
 
 double OneZeroZeroEWSplitFn::P(const double z, const Energy2 t,
 			       const IdList &ids, const bool mass, const RhoDMatrix & rho) const {
   Complex gvhh(0.,0.);
   getCouplings(gvhh,ids);
   double rho00 = abs(rho(0,0));
   double rho11 = abs(rho(1,1));
   double rho22 = abs(rho(2,2));
   // the splitting in the massless limit
   double val = (1.-z)*z*(rho00+rho22);
   // the massive limit
   if(mass){
     // get the running mass
     double m0t2 = sqr(getParticleData(ids[0]->id())->mass())/t;
     double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
     double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
-    val += (m0t2*sqr(-1.+2.*z)*rho11)/2.+(-(m1t2*(1.-z))-m2t2*z+m0t2*(1.-z)*z)
-        *(rho00+rho22);
+    val += (rho00+rho22)*(z*(1.-z)*m0t2-(1.-z)*m1t2-z*m2t2)+rho11*sqr(1.-2.*z)/2.*m0t2;
   }
   return norm(gvhh)*val;
 }
 
 
 double OneZeroZeroEWSplitFn::overestimateP(const double z,
 					   const IdList & ids) const {
   Complex gvhh(0.,0.);
   getCouplings(gvhh,ids);
   return norm(gvhh)*(1.-z)*z;
 }
 
 
 double OneZeroZeroEWSplitFn::ratioP(const double z, const Energy2 t,
 				    const IdList & ids, const bool mass,
 				    const RhoDMatrix & rho) const {
   double rho00 = abs(rho(0,0));
   double rho11 = abs(rho(1,1));
   double rho22 = abs(rho(2,2));
   // ratio in the massless limit
   double val = rho00+rho22;
   // the massive limit
   if(mass){
     // get the running mass
     double m0t2 = sqr(ids[0]->mass())/t;
     double m1t2 = sqr(ids[1]->mass())/t;
     double m2t2 = sqr(ids[2]->mass())/t;
     val += ((m0t2*sqr(-1.+2.*z)*rho11)/2.+(-(m1t2*(1.-z))-m2t2*z+m0t2*(1.-z)*z)
         *(rho00+rho22))/((1.-z)*z);
   }
   return val;
 }
 
 
 double OneZeroZeroEWSplitFn::integOverP(const double z,
 				      const IdList & ids,
 				      unsigned int PDFfactor) const {
   Complex gvhh(0.,0.);
   getCouplings(gvhh,ids);
   double pre = norm(gvhh);
   switch (PDFfactor) {
   case 0:
     return pre/6.*(3.-2.*z)*sqr(z);
   case 1:
   case 2:
   case 3:
   default:
     throw Exception() << "OneZeroZeroEWSplitFn::integOverP() invalid PDFfactor = "
 		      << PDFfactor << Exception::runerror;
   }
 }
 
 
 double OneZeroZeroEWSplitFn::invIntegOverP(const double r, const IdList & ids,
 					   unsigned int PDFfactor) const {
   Complex gvhh(0.,0.);
   getCouplings(gvhh,ids);
   double pre = norm(gvhh);
   switch (PDFfactor) {
   case 0:
     return (1.-pre/pow(-pow(pre,3)+12.*pow(pre,2)*r+2.*sqrt(6.)
             *sqrt(-(pow(pre,4)*(pre-6.*r)*r)),1./3.)- pow(-pow(pre,3)+12.
             *pow(pre,2)*r+2.*sqrt(6)*sqrt(-(pow(pre,4)*(pre-6.*r)*r)),1./3.)/pre)/2.;
   case 1:
   case 2:
   case 3:
   default:
     throw Exception() << "OneZeroZeroEWSplitFn::invIntegOverP() invalid PDFfactor = "
 		      << PDFfactor << Exception::runerror;
   }
 }
 
 
 bool OneZeroZeroEWSplitFn::accept(const IdList &ids) const {
   if(ids.size()!=3)
     return false;
   if(_couplingValueIm==0.&&_couplingValueRe==0.)
     return false;
   if(ids[0]->iCharge()!=ids[1]->iCharge()+ids[2]->iCharge())
     return false;
   if(ids[0]->iSpin()==PDT::Spin1 && ids[1]->iSpin()==PDT::Spin0 && ids[2]->iSpin()==PDT::Spin0)
     return true;
   else
     return false;
 }
 
 
 vector<pair<int, Complex> >
 OneZeroZeroEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & ,
 				       const RhoDMatrix &) {
   // no dependence on the spin density matrix, dependence on off-diagonal terms cancels
   // and rest = splitting function for Tr(rho)=1 as required by defn
   return vector<pair<int, Complex> >(1,make_pair(0,1.));
 }
 
 
 vector<pair<int, Complex> >
 OneZeroZeroEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & ,
 					const RhoDMatrix &) {
   // no dependence on the spin density matrix, dependence on off-diagonal terms cancels
   // and rest = splitting function for Tr(rho)=1 as required by defn
   return vector<pair<int, Complex> >(1,make_pair(0,1.));
 }
 
 
 DecayMEPtr OneZeroZeroEWSplitFn::matrixElement(const double z, const Energy2 t,
                                              const IdList & ids, const double phi,
                                              bool) {
   Complex gvhh(0.,0.);
   getCouplings(gvhh,ids);
   // calculate the kernal
   DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,PDT::Spin0,PDT::Spin0)));
   double m0t2 = sqr(getParticleData(ids[0]->id())->mass())/t;
   double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
   double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
   Complex phase  = exp(Complex(0.,1.)*phi);
   Complex cphase = conj(phase);
   double sqrtmass = sqrt(m0t2*z*(1.-z)-m1t2*(1.-z)-m2t2*z+z*(1.-z));
   // assign kernel
   (*kernal)(0,0,0) = -cphase*sqrtmass;                        // 111
   (*kernal)(1,0,0) =  sqrt(m0t2)*(1.-2.*z)/sqrt(2.);         // 211 -> 411
   (*kernal)(2,0,0) =  phase*sqrtmass;                         // 311
   // return the answer
   return kernal;
 }