diff --git a/Decay/General/SSVDecayer.cc b/Decay/General/SSVDecayer.cc
--- a/Decay/General/SSVDecayer.cc
+++ b/Decay/General/SSVDecayer.cc
@@ -1,355 +1,343 @@
 // -*- C++ -*-
 //
 // SSVDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2019 The Herwig Collaboration
 //
 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
 //
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the SSVDecayer class.
 //
 
 #include "SSVDecayer.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "ThePEG/PDT/DecayMode.h"
 #include "Herwig/Utilities/Kinematics.h"
 #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
 #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
 #include "Herwig/Decay/GeneralDecayMatrixElement.h"
 using namespace Herwig;
 using namespace ThePEG::Helicity;
 
 IBPtr SSVDecayer::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr SSVDecayer::fullclone() const {
   return new_ptr(*this);
 }
 
 void SSVDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
 			      vector<VertexBasePtr> vertex,
 			      map<ShowerInteraction,VertexBasePtr> & inV,
 			      const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
 			      map<ShowerInteraction,VertexBasePtr> fourV) {
   decayInfo(incoming,outgoing);
   for(auto vert : vertex)
     vertex_            .push_back(dynamic_ptr_cast<AbstractVSSVertexPtr>(vert));
   vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
   for(auto & inter : itemp) {
     incomingVertex_[inter]  = dynamic_ptr_cast<AbstractVSSVertexPtr>(inV.at(inter));
     fourPointVertex_[inter] = dynamic_ptr_cast<AbstractVVSSVertexPtr>(fourV.at(inter));
     outgoingVertexS_[inter] = AbstractVSSVertexPtr();
-    outgoingVertexV_[inter] = AbstractVVVVertexPtr();
+    outgoingVertexV_[inter] = AbstractVVVVertexPtr();  
     if(outV[0].at(inter)) {
       if (outV[0].at(inter)->getName()==VertexType::VSS)
 	outgoingVertexS_[inter]   = dynamic_ptr_cast<AbstractVSSVertexPtr>(outV[0].at(inter));
       else
 	outgoingVertexV_[inter]   = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[0].at(inter));
     }
     if(outV[1].at(inter)) {
       if (outV[1].at(inter)->getName()==VertexType::VSS)
 	outgoingVertexS_[inter]   = dynamic_ptr_cast<AbstractVSSVertexPtr>(outV[1].at(inter));
-      else
+      else 
 	outgoingVertexV_[inter]   = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[1].at(inter));
     }
   }
 }
 
 void SSVDecayer::persistentOutput(PersistentOStream & os) const {
   os << vertex_
      << incomingVertex_   << outgoingVertexS_
      << outgoingVertexV_  << fourPointVertex_;
 }
 
 void SSVDecayer::persistentInput(PersistentIStream & is, int) {
   is >> vertex_
      >> incomingVertex_   >> outgoingVertexS_
      >> outgoingVertexV_  >> fourPointVertex_;
 }
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<SSVDecayer,GeneralTwoBodyDecayer>
 describeHerwigSSVDecayer("Herwig::SSVDecayer", "Herwig.so");
 
 void SSVDecayer::Init() {
 
   static ClassDocumentation<SSVDecayer> documentation
     ("This implements the decay of a scalar to a vector and a scalar");
 
 }
 
 void SSVDecayer::
 constructSpinInfo(const Particle & part, ParticleVector decay) const {
   unsigned int isc(0),ivec(1);
   if(decay[0]->dataPtr()->iSpin() != PDT::Spin0) swap(isc,ivec);
   ScalarWaveFunction::
     constructSpinInfo(const_ptr_cast<tPPtr>(&part),incoming,true);
   ScalarWaveFunction::
     constructSpinInfo(decay[isc],outgoing,true);
   VectorWaveFunction::
     constructSpinInfo(vector_,decay[ivec],outgoing,true,false);
 }
 
 double SSVDecayer::me2(const int,const Particle & part,
 		       const tPDVector & outgoing,
 		       const vector<Lorentz5Momentum> & momenta,
 		       MEOption meopt) const {
   unsigned int isc(0),ivec(1);
   if(outgoing[0]->iSpin() != PDT::Spin0) swap(isc,ivec);
   if(!ME()) {
     if(ivec==1)
       ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin1)));
     else
       ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1,PDT::Spin0)));
   }
   if(meopt==Initialize) {
     ScalarWaveFunction::
       calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&part),incoming);
     swave_ = ScalarWaveFunction(part.momentum(),part.dataPtr(),incoming);
     // fix rho if no correlations
     fixRho(rho_);
   }
-  bool massless = decay[1]->id()==ParticleID::gamma || decay[1]->id()==ParticleID::g;
-  if(meopt==Terminate) {
-    ScalarWaveFunction::
-      constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),incoming,true);
-    ScalarWaveFunction::
-      constructSpinInfo(decay[isc],outgoing,true);
-    VectorWaveFunction::
-      constructSpinInfo(vector_,decay[ivec],outgoing,true,massless);
-  }
   VectorWaveFunction::
-    calculateWaveFunctions(vector_,decay[ivec],outgoing,massless);
-
-  ScalarWaveFunction sca(decay[isc]->momentum(),decay[isc]->dataPtr(),outgoing);
-  Energy2 scale(sqr(inpart.mass()));
+    calculateWaveFunctions(vector_,momenta[ivec],outgoing[ivec],Helicity::outgoing,false);
+  ScalarWaveFunction sca(momenta[isc],outgoing[isc],Helicity::outgoing);
+  Energy2 scale(sqr(part.mass()));
   //make sure decay matrix element is in the correct order
   double output(0.);
   if(ivec == 0) {
     for(unsigned int ix = 0; ix < 3; ++ix) {
       (*ME())(0, ix, 0) = 0.;
-      if(massless && ix==1) continue;
       for(auto vert : vertex_)
 	(*ME())(0, ix, 0) += vert->evaluate(scale,vector_[ix],sca, swave_);
     }
   }
   else {
     for(unsigned int ix = 0; ix < 3; ++ix) {
       (*ME())(0, 0, ix) = 0.;
-      if(massless && ix==1) continue;
       for(auto vert : vertex_)
 	(*ME())(0, 0, ix) += vert->evaluate(scale,vector_[ix],sca,swave_);
     }
   }
   output = (ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
   // colour and identical particle factors
   output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]);
   // return the answer
   return output;
 }
 
-Energy SSVDecayer:: partialWidth(PMPair inpart, PMPair outa,
+Energy SSVDecayer:: partialWidth(PMPair inpart, PMPair outa, 
 				 PMPair outb) const {
   if( inpart.second < outa.second + outb.second  ) return ZERO;
   return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
 }
 
 
 double  SSVDecayer::threeBodyME(const int , const Particle & inpart,
 				const ParticleVector & decay,
 				ShowerInteraction inter, MEOption meopt) {
   int iscal (0), ivect (1), iglu (2);
   // get location of outgoing scalar/vector
   if(decay[1]->dataPtr()->iSpin()==PDT::Spin0) swap(iscal,ivect);
 
   if(meopt==Initialize) {
     // create scalar wavefunction for decaying particle
     ScalarWaveFunction::calculateWaveFunctions(rho3_,const_ptr_cast<tPPtr>(&inpart),incoming);
     swave3_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),incoming);
   }
   // setup spin information when needed
   if(meopt==Terminate) {
     ScalarWaveFunction::
       constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),incoming,true);
     ScalarWaveFunction::
       constructSpinInfo(decay[iscal],outgoing,true);
      VectorWaveFunction::
       constructSpinInfo(vector3_,decay[ivect],outgoing,true,false);
     VectorWaveFunction::
       constructSpinInfo(gluon_,  decay[iglu ],outgoing,true,false);
     return 0.;
   }
   // calculate colour factors and number of colour flows
   unsigned int nflow;
   vector<DVector> cfactors = getColourFactors(inpart, decay, nflow);
   vector<GeneralDecayMEPtr> ME(nflow,new_ptr(GeneralDecayMatrixElement(PDT::Spin0, PDT::Spin0,
 								       PDT::Spin1, PDT::Spin1)));
 
   // create wavefunctions
   ScalarWaveFunction scal(decay[iscal]->momentum(),  decay[iscal]->dataPtr(),outgoing);
   VectorWaveFunction::calculateWaveFunctions(vector3_,decay[ivect],outgoing,false);
   VectorWaveFunction::calculateWaveFunctions(gluon_,  decay[iglu ],outgoing,true );
 
   // gauge invariance test
 #ifdef GAUGE_CHECK
   gluon_.clear();
   for(unsigned int ix=0;ix<3;++ix) {
     if(ix==1) gluon_.push_back(VectorWaveFunction());
     else {
       gluon_.push_back(VectorWaveFunction(decay[iglu ]->momentum(),
   					  decay[iglu ]->dataPtr(),10,
   					  outgoing));
     }
   }
 #endif
 
   if (! ((incomingVertex_[inter]  && (outgoingVertexS_[inter] || outgoingVertexV_[inter])) ||
 	 (outgoingVertexS_[inter] &&  outgoingVertexV_[inter])))
     throw Exception()
       << "Invalid vertices for radiation in SSV decay in SSVDecayer::threeBodyME"
       << Exception::runerror;
 
 
   // sort out colour flows
   int S(1), V(2);
-  if (decay[iscal]->dataPtr()->iColour()==PDT::Colour3bar &&
+  if (decay[iscal]->dataPtr()->iColour()==PDT::Colour3bar && 
       decay[ivect]->dataPtr()->iColour()==PDT::Colour8)
     swap(S,V);
-  else if (decay[ivect]->dataPtr()->iColour()==PDT::Colour3 &&
+  else if (decay[ivect]->dataPtr()->iColour()==PDT::Colour3 && 
 	   decay[iscal]->dataPtr()->iColour()==PDT::Colour8)
     swap(S,V);
 
   Energy2 scale(sqr(inpart.mass()));
 
   const GeneralTwoBodyDecayer::CFlow & colourFlow
         = colourFlows(inpart, decay);
   double gs(0.);
   bool couplingSet(false);
 #ifdef GAUGE_CHECK
   double total=0.;
 #endif
   for(unsigned int iv = 0; iv < 3; ++iv) {
     for(unsigned int ig = 0; ig < 2; ++ig) {
       // radiation from the incoming scalar
       if((inpart.dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
 	 (inpart.dataPtr()->charged()  && inter==ShowerInteraction::QED) ) {
 	assert(incomingVertex_[inter]);
-	ScalarWaveFunction scalarInter =
+	ScalarWaveFunction scalarInter = 
 	  incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),
 					   gluon_[2*ig],swave3_,inpart.mass());
-
+	
 	assert(swave3_.particle()->id()==scalarInter.particle()->id());
 	Complex diag = 0.;
 	for(auto vertex : vertex_)
 	  diag += vertex->evaluate(scale,vector3_[iv],scal,scalarInter);
 	if(!couplingSet) {
 	  gs = abs(incomingVertex_[inter]->norm());
 	  couplingSet = true;
 	}
 	for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
-	  (*ME[colourFlow[0][ix].first])(0, 0, iv, ig) +=
-	    colourFlow[0][ix].second*diag;
+	  (*ME[colourFlow[0][ix].first])(0, 0, iv, ig) += 
+	    colourFlow[0][ix].second*diag; 
 	}
 #ifdef GAUGE_CHECK
 	total+=norm(diag);
 #endif
       }
       // radiation from the outgoing scalar
       if((decay[iscal]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
 	 (decay[iscal]->dataPtr()->charged()  && inter==ShowerInteraction::QED) ) {
 	assert(outgoingVertexS_[inter]);
 	// ensure you get correct outgoing particle from first vertex
 	tcPDPtr off = decay[iscal]->dataPtr();
 	if(off->CC()) off = off->CC();
-	ScalarWaveFunction scalarInter =
+	ScalarWaveFunction scalarInter = 
 	  outgoingVertexS_[inter]->evaluate(scale,3,off,gluon_[2*ig],scal,decay[iscal]->mass());
-
+	
 	assert(scal.particle()->id()==scalarInter.particle()->id());
 
 	if(!couplingSet) {
 	  gs = abs(outgoingVertexS_[inter]->norm());
 	  couplingSet = true;
 	}
 	Complex diag = 0.;
 	for(auto vertex : vertex_)
 	  diag += vertex->evaluate(scale,vector3_[iv],scalarInter,swave3_);
 	for(unsigned int ix=0;ix<colourFlow[S].size();++ix) {
-	  (*ME[colourFlow[S][ix].first])(0, 0, iv, ig) +=
+	  (*ME[colourFlow[S][ix].first])(0, 0, iv, ig) += 
 	    colourFlow[S][ix].second*diag;
 	}
 #ifdef GAUGE_CHECK
 	total+=norm(diag);
 #endif
       }
 
       // radiation from outgoing vector
       if((decay[ivect]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
 	 (decay[ivect]->dataPtr()->charged()  && inter==ShowerInteraction::QED) ) {
 	assert(outgoingVertexV_[inter]);
 	// ensure you get correct outgoing particle from first vertex
 	tcPDPtr off = decay[ivect]->dataPtr();
 	if(off->CC()) off = off->CC();
-	VectorWaveFunction  vectorInter =
+	VectorWaveFunction  vectorInter = 
 	  outgoingVertexV_[inter]->evaluate(scale,3,off,gluon_[2*ig],
 					    vector3_[iv],decay[ivect]->mass());
-
+	
 	assert(vector3_[iv].particle()->id()==vectorInter.particle()->id());
-
+	
 	if(!couplingSet) {
 	  gs = abs(outgoingVertexV_[inter]->norm());
 	  couplingSet = true;
-	}
+	}	
 	Complex diag = 0.;
 	for(auto vertex : vertex_)
 	  diag += vertex->evaluate(scale,vectorInter,scal,swave3_);
 	for(unsigned int ix=0;ix<colourFlow[V].size();++ix) {
-	  (*ME[colourFlow[V][ix].first])(0, 0, iv, ig) +=
+	  (*ME[colourFlow[V][ix].first])(0, 0, iv, ig) += 
 	    colourFlow[V][ix].second*diag;
 	}
 #ifdef GAUGE_CHECK
 	total+=norm(diag);
 #endif
       }
       // radiation from 4 point vertex
       if (fourPointVertex_[inter]) {
 	Complex diag =  fourPointVertex_[inter]->evaluate(scale, gluon_[2*ig], vector3_[iv],
 							  scal, swave3_);
 	for(unsigned int ix=0;ix<colourFlow[3].size();++ix) {
-	  (*ME[colourFlow[3][ix].first])(0, 0, iv, ig) +=
+	  (*ME[colourFlow[3][ix].first])(0, 0, iv, ig) += 
 	     colourFlow[3][ix].second*diag;
 	}
 #ifdef GAUGE_CHECK
 	total+=norm(diag);
 #endif
       }
     }
   }
-
+  
   // contract matrices
   double output=0.;
   for(unsigned int ix=0; ix<nflow; ++ix){
     for(unsigned int iy=0; iy<nflow; ++iy){
       output+=cfactors[ix][iy]*(ME[ix]->contract(*ME[iy],rho3_)).real();
     }
   }
   // divide by alpha_(S,EM)
   output*=(4.*Constants::pi)/sqr(gs);
 #ifdef GAUGE_CHECK
   double ratio = output/total;
   if(abs(ratio)>1e-20) {
     generator()->log() << "Test of gauge invariance in decay\n" << inpart << "\n";
     for(unsigned int ix=0;ix<decay.size();++ix)
       generator()->log() << *decay[ix] << "\n";
     generator()->log() << "Test of gauge invariance " << ratio << "\n";
   }
 #endif
   // return the answer
   return output;
 }