No OneTemporary
Actions

Size

95 KB

Subscribers

None

View Options

	diff --git a/Decay/General/SFFDecayer.cc b/Decay/General/SFFDecayer.cc
	--- a/Decay/General/SFFDecayer.cc
	+++ b/Decay/General/SFFDecayer.cc
	@@ -1,493 +1,500 @@
	// -- C++ --
	//
	// SFFDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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 SFFDecayer class.
	//

	#include "SFFDecayer.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 "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
	#include "Herwig/Utilities/Kinematics.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"

	using namespace Herwig;
	using namespace ThePEG::Helicity;

	IBPtr SFFDecayer::clone() const {
	return new_ptr(*this);
	}

	IBPtr SFFDecayer::fullclone() const {
	return new_ptr(*this);
	}

	void SFFDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
	vector<VertexBasePtr> vertex,
	map<ShowerInteraction,VertexBasePtr> & inV,
	const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
	map<ShowerInteraction,VertexBasePtr> ) {
	decayInfo(incoming,outgoing);
	for(auto vert : vertex) {
	vertex_ .push_back(dynamic_ptr_cast<AbstractFFSVertexPtr>(vert));
	perturbativeVertex_ .push_back(dynamic_ptr_cast<FFSVertexPtr> (vert));
	}
	vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
	for(auto & inter : itemp) {
	incomingVertex_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(inV.at(inter));
	outgoingVertex1_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[0].at(inter));
	outgoingVertex2_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[1].at(inter));
	}
	}

	void SFFDecayer::persistentOutput(PersistentOStream & os) const {
	os << vertex_ << perturbativeVertex_
	<< incomingVertex_ << outgoingVertex1_
	<< outgoingVertex2_;
	}

	void SFFDecayer::persistentInput(PersistentIStream & is, int) {
	is >> vertex_ >> perturbativeVertex_
	>> incomingVertex_ >> outgoingVertex1_
	>> outgoingVertex2_;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	-DescribeClass<SFFDecayer,GeneralTwoBodyDecayer>
	+DescribeClass<SFFDecayer,GeneralTwoBodyDecayer2>
	describeHerwigSFFDecayer("Herwig::SFFDecayer", "Herwig.so");

	void SFFDecayer::Init() {

	static ClassDocumentation<SFFDecayer> documentation
	("This class implements to decay of a scalar to 2 fermions");

	}
	-
	-double SFFDecayer::me2(const int , const Particle & inpart,
	- const ParticleVector & decay,MEOption meopt) const {
	- if(!ME())
	- ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1Half,PDT::Spin1Half)));
	- // work out which is the fermion and antifermion
	+void SFFDecayer::
	+constructSpinInfo(const Particle & part, ParticleVector decay) const {
	int iferm(1),ianti(0);
	int itype[2];
	for(unsigned int ix=0;ix<2;++ix) {
	if(decay[ix]->dataPtr()->CC()) itype[ix] = decay[ix]->id()>0 ? 0:1;
	else itype[ix] = 2;
	}
	if(itype[0]==0\|\|itype[1]==1\|\|(itype[0]==2&&itype[1]==2)) swap(iferm,ianti);
	+ ScalarWaveFunction::
	+ constructSpinInfo(const_ptr_cast<tPPtr>(&part),incoming,true);
	+ SpinorBarWaveFunction::
	+ constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
	+ SpinorWaveFunction::
	+ constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
	+}
	+double SFFDecayer::me2(const int,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	+ MEOption meopt) const {
	+ if(!ME())
	+ ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1Half,PDT::Spin1Half)));
	+ // work out which is the fermion and antifermion
	+ int iferm(1),ianti(0);
	+ int itype[2];
	+ for(unsigned int ix=0;ix<2;++ix) {
	+ if(outgoing[ix]->CC()) itype[ix] = outgoing[ix]->id()>0 ? 0:1;
	+ else itype[ix] = 2;
	+ }
	+ if(itype[0]==0\|\|itype[1]==1\|\|(itype[0]==2&&itype[1]==2)) swap(iferm,ianti);

	if(meopt==Initialize) {
	ScalarWaveFunction::
	- calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&inpart),incoming);
	- swave_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),incoming);
	+ calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&part),incoming);
	+ swave_ = ScalarWaveFunction(part.momentum(),part.dataPtr(),incoming);
	// fix rho if no correlations
	fixRho(rho_);
	}
	- if(meopt==Terminate) {
	- ScalarWaveFunction::
	- constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),incoming,true);
	- SpinorBarWaveFunction::
	- constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
	- SpinorWaveFunction::
	- constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
	- return 0.;
	- }
	SpinorBarWaveFunction::
	- calculateWaveFunctions(wavebar_,decay[iferm],outgoing);
	+ calculateWaveFunctions(wavebar_,momenta[iferm],outgoing[iferm],Helicity::outgoing);
	SpinorWaveFunction::
	- calculateWaveFunctions(wave_ ,decay[ianti],outgoing);
	- Energy2 scale(sqr(inpart.mass()));
	+ calculateWaveFunctions(wave_ ,momenta[ianti],outgoing[ianti],Helicity::outgoing);
	+ Energy2 scale(sqr(part.mass()));
	for(unsigned int ifm = 0; ifm < 2; ++ifm){
	for(unsigned int ia = 0; ia < 2; ++ia) {
	if(iferm > ianti) (*ME())(0, ia, ifm) = 0.;
	else (*ME())(0, ifm, ia) = 0.;
	for(auto vert : vertex_) {
	if(iferm > ianti){
	(*ME())(0, ia, ifm) += vert->evaluate(scale,wave_[ia],
	wavebar_[ifm],swave_);
	}
	else {
	(*ME())(0, ifm, ia) += vert->evaluate(scale,wave_[ia],
	wavebar_[ifm],swave_);
	}
	}
	}
	}

	double output = (ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
	// colour and identical particle factors
	- output *= colourFactor(inpart.dataPtr(),decay[0]->dataPtr(),
	- decay[1]->dataPtr());
	+ output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]);
	// return the answer
	return output;
	}

	Energy SFFDecayer::partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const {
	if( inpart.second < outa.second + outb.second ) return ZERO;
	if(perturbativeVertex_.size()==1 &&
	perturbativeVertex_[0]) {
	tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
	perturbativeVertex_[0]->setCoupling(sqr(inpart.second), outb.first, outa.first,
	in);
	double mu1(outa.second/inpart.second),mu2(outb.second/inpart.second);
	double c2 = norm(perturbativeVertex_[0]->norm());
	Complex al(perturbativeVertex_[0]->left()), ar(perturbativeVertex_[0]->right());
	double me2 = -c2( (norm(al) + norm(ar))( sqr(mu1) + sqr(mu2) - 1.)
	+ 2.(arconj(al) + alconj(ar)).real()mu1*mu2 );
	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second, outa.second,
	outb.second);
	Energy output = me2pcm/(8Constants::pi);
	// colour factor
	output *= colourFactor(inpart.first,outa.first,outb.first);
	return output;
	}
	else {
	- return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
	+ return GeneralTwoBodyDecayer2::partialWidth(inpart,outa,outb);
	}
	}

	double SFFDecayer::threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt) {
	// work out which is the fermion and antifermion
	int ianti(0), iferm(1), iglu(2);
	int itype[2];
	for(unsigned int ix=0;ix<2;++ix) {
	if(decay[ix]->dataPtr()->CC()) itype[ix] = decay[ix]->id()>0 ? 0:1;
	else itype[ix] = 2;
	}
	if(itype[0]==0 && itype[1]!=0) swap(iferm, ianti);
	if(itype[0]==2 && itype[1]==1) swap(iferm, ianti);
	if(itype[0]==0 && itype[1]==0 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
	swap(iferm, ianti);
	if(itype[0]==1 && itype[1]==1 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
	swap(iferm, ianti);

	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);
	SpinorBarWaveFunction::
	constructSpinInfo(wavebar3_ ,decay[iferm],outgoing,true);
	SpinorWaveFunction::
	constructSpinInfo(wave3_ ,decay[ianti],outgoing,true);
	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::Spin1Half,
	PDT::Spin1Half, PDT::Spin1)));
	// create wavefunctions
	SpinorBarWaveFunction::
	calculateWaveFunctions(wavebar3_, decay[iferm],outgoing);
	SpinorWaveFunction::
	calculateWaveFunctions(wave3_ , decay[ianti],outgoing);
	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

	// identify fermion and/or anti-fermion vertex
	AbstractFFVVertexPtr outgoingVertexF;
	AbstractFFVVertexPtr outgoingVertexA;
	identifyVertices(iferm, ianti, inpart, decay, outgoingVertexF, outgoingVertexA,
	inter);

	- const GeneralTwoBodyDecayer::CFlow & colourFlow
	+ const GeneralTwoBodyDecayer2::CFlow & colourFlow
	= colourFlows(inpart, decay);

	Energy2 scale(sqr(inpart.mass()));
	double gs(0.);
	bool couplingSet(false);
	#ifdef GAUGE_CHECK
	double total=0.;
	#endif
	for(unsigned int ifm = 0; ifm < 2; ++ifm) {
	for(unsigned int ia = 0; ia < 2; ++ia) {
	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 =
	incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),
	gluon_[2*ig],swave3_,inpart.mass());

	assert(swave3_.particle()->id()==scalarInter.particle()->id());

	if(!couplingSet) {
	gs = abs(incomingVertex_[inter]->norm());
	couplingSet = true;
	}
	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,wave3_[ia],
	wavebar3_[ifm],scalarInter);
	for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
	(*ME[colourFlow[0][ix].first])(0, ia, ifm, ig) +=
	colourFlow[0][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}
	// radiation from outgoing fermion
	if((decay[iferm]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[iferm]->dataPtr()->charged() && inter==ShowerInteraction::QED)) {
	assert(outgoingVertexF);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[iferm]->dataPtr();
	if(off->CC()) off = off->CC();
	SpinorBarWaveFunction interS =
	outgoingVertexF->evaluate(scale,3,off,wavebar3_[ifm],
	gluon_[2*ig],decay[iferm]->mass());

	assert(wavebar3_[ifm].particle()->id()==interS.particle()->id());

	if(!couplingSet) {
	gs = abs(outgoingVertexF->norm());
	couplingSet = true;
	}
	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,wave3_[ia], interS,swave3_);
	for(unsigned int ix=0;ix<colourFlow[1].size();++ix) {
	(*ME[colourFlow[1][ix].first])(0, ia, ifm, ig) +=
	colourFlow[1][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}

	// radiation from outgoing antifermion
	if((decay[ianti]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[ianti]->dataPtr()->charged() && inter==ShowerInteraction::QED)) {
	assert(outgoingVertexA);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[ianti]->dataPtr();
	if(off->CC()) off = off->CC();
	SpinorWaveFunction interS =
	outgoingVertexA->evaluate(scale,3,off,wave3_[ia],
	gluon_[2*ig],decay[ianti]->mass());

	assert(wave3_[ia].particle()->id()==interS.particle()->id());

	if(!couplingSet) {
	gs = abs(outgoingVertexA->norm());
	couplingSet = true;
	}
	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,interS,wavebar3_[ifm],swave3_);
	for(unsigned int ix=0;ix<colourFlow[2].size();++ix) {
	(*ME[colourFlow[2][ix].first])(0, ia, ifm, ig) +=
	colourFlow[2][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;
	}

	void SFFDecayer::identifyVertices(const int iferm, const int ianti,
	const Particle & inpart, const ParticleVector & decay,
	AbstractFFVVertexPtr & outgoingVertexF,
	AbstractFFVVertexPtr & outgoingVertexA,
	ShowerInteraction inter) {
	// QCD
	if(inter==ShowerInteraction::QCD) {
	// work out which fermion each outgoing vertex corresponds to
	// two outgoing vertices
	if( inpart.dataPtr() ->iColour()==PDT::Colour0 &&
	((decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar) \|\|
	(decay[iferm]->dataPtr()->iColour()==PDT::Colour8 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour8))) {
	if(outgoingVertex1_[inter]==outgoingVertex2_[inter]) {
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else if (outgoingVertex2_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	}
	else if(inpart.dataPtr() ->iColour()==PDT::Colour8 &&
	decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar) {
	if(outgoingVertex1_[inter]==outgoingVertex2_[inter]) {
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else if (outgoingVertex2_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	}

	// one outgoing vertex
	else if(inpart.dataPtr()->iColour()==PDT::Colour3){
	if(decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour0){
	if (outgoingVertex1_[inter]) outgoingVertexF = outgoingVertex1_[inter];
	else if(outgoingVertex2_[inter]) outgoingVertexF = outgoingVertex2_[inter];
	}
	else if (decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour8) {
	if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[ianti]->dataPtr()))) {
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	else {
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	}
	else if(decay[iferm]->dataPtr()->iColour()==PDT::Colour3bar &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar) {
	if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else {
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	}
	}
	else if(inpart.dataPtr()->iColour()==PDT::Colour3bar){
	if(decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar &&
	decay[iferm]->dataPtr()->iColour()==PDT::Colour0){
	if (outgoingVertex1_[inter]) outgoingVertexA = outgoingVertex1_[inter];
	else if(outgoingVertex2_[inter]) outgoingVertexA = outgoingVertex2_[inter];
	}
	else if (decay[iferm]->dataPtr()->iColour()==PDT::Colour8 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar){
	if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else {
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	}
	else if(decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour3) {
	if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else {
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	}
	}
	else if(inpart.dataPtr()->iColour()==PDT::Colour6 \|\|
	inpart.dataPtr()->iColour()==PDT::Colour6bar) {
	if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else {
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	}

	if (! ((incomingVertex_[inter] && (outgoingVertexF \|\| outgoingVertexA)) \|\|
	( outgoingVertexF && outgoingVertexA))) {
	throw Exception()
	<< "Invalid vertices for QCD radiation in SFF decay in SFFDecayer::identifyVertices"
	<< Exception::runerror;
	}
	}
	// QED
	else {
	if(decay[iferm]->dataPtr()->charged()) {
	if (outgoingVertex1_[inter] &&
	outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr())))
	outgoingVertexF = outgoingVertex1_[inter];
	else
	outgoingVertexF = outgoingVertex2_[inter];
	}
	if(decay[ianti]->dataPtr()->charged()) {
	if (outgoingVertex1_[inter] &&
	outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[ianti]->dataPtr())))
	outgoingVertexA = outgoingVertex1_[inter];
	else
	outgoingVertexA = outgoingVertex2_[inter];
	}
	}
	}
	diff --git a/Decay/General/SFFDecayer.h b/Decay/General/SFFDecayer.h
	--- a/Decay/General/SFFDecayer.h
	+++ b/Decay/General/SFFDecayer.h
	@@ -1,234 +1,241 @@
	// -- C++ --
	//
	// SFFDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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.
	//
	#ifndef HERWIG_SFFDecayer_H
	#define HERWIG_SFFDecayer_H
	//
	// This is the declaration of the SFFDecayer class.
	//

	-#include "GeneralTwoBodyDecayer.h"
	+#include "GeneralTwoBodyDecayer2.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Helicity/Vertex/Scalar/FFSVertex.h"
	#include "ThePEG/Helicity/Vertex/Scalar/VSSVertex.h"
	#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"

	namespace Herwig {
	using namespace ThePEG;
	using Helicity::FFSVertexPtr;

	/** \ingroup Decay
	* The SFFDecayer class implements the decay of a scalar to 2
	* fermions in a general model. It holds an FFSVertex pointer that
	* must be typecast from the VertexBase pointer held in
	- * GeneralTwoBodyDecayer. It implents the virtual functions me2() and
	+ * GeneralTwoBodyDecayer2. It implents the virtual functions me2() and
	* partialWidth().
	*
	- * @see GeneralTwoBodyDecayer
	+ * @see GeneralTwoBodyDecayer2
	*/
	-class SFFDecayer: public GeneralTwoBodyDecayer {
	+class SFFDecayer: public GeneralTwoBodyDecayer2 {

	public:

	/**
	* The default constructor.
	*/
	SFFDecayer() {}
	-
	- /** @name Virtual functions required by the Decayer class. */
	- //@{
	- /**
	+
	+ /**
	* Return the matrix element squared for a given mode and phase-space channel.
	- * @param ichan The channel we are calculating the matrix element for.
	+ * @param ichan The channel we are calculating the matrix element for.
	* @param part The decaying Particle.
	- * @param decay The particles produced in the decay.
	+ * @param outgoing The particles produced in the decay
	+ * @param momenta The momenta of the particles produced in the decay
	* @param meopt Option for the calculation of the matrix element
	* @return The matrix element squared for the phase-space configuration.
	*/
	- virtual double me2(const int ichan, const Particle & part,
	- const ParticleVector & decay, MEOption meopt) const;
	+ double me2(const int ichan,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	+ MEOption meopt) const;
	+
	+ /**
	+ * Construct the SpinInfos for the particles produced in the decay
	+ */
	+ virtual void constructSpinInfo(const Particle & part,
	+ ParticleVector outgoing) const;

	/**
	* Function to return partial Width
	* @param inpart The decaying particle.
	* @param outa One of the decay products.
	* @param outb The other decay product.
	*/
	virtual Energy partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const;

	/**
	* Has a POWHEG style correction
	*/
	virtual POWHEGType hasPOWHEGCorrection() {
	POWHEGType output = FSR;
	for(auto vertex : vertex_) {
	if(vertex->orderInAllCouplings()!=1) {
	output = No;
	break;
	}
	}
	return output;
	}

	/**
	* Three-body matrix element including additional QCD radiation
	*/
	virtual double threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter,
	MEOption meopt);

	/**
	* Indentify outgoing vertices for the fermion and antifermion
	*/
	void identifyVertices(const int iferm, const int ianti,
	const Particle & inpart, const ParticleVector & decay,
	AbstractFFVVertexPtr & outgoingVertexF,
	AbstractFFVVertexPtr & outgoingVertexA,
	ShowerInteraction inter);

	/**
	* Set the information on the decay
	*/
	virtual void setDecayInfo(PDPtr incoming, PDPair outgoing,
	vector<VertexBasePtr>,
	map<ShowerInteraction,VertexBasePtr> &,
	const vector<map<ShowerInteraction,VertexBasePtr> > &,
	map<ShowerInteraction,VertexBasePtr>);
	//@}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	SFFDecayer & operator=(const SFFDecayer &);

	private:

	/**
	* Abstract pointer to AbstractFFSVertex
	*/
	vector<AbstractFFSVertexPtr> vertex_;

	/**
	* Pointer to the perturbative vertex
	*/
	vector<FFSVertexPtr> perturbativeVertex_;

	/**
	* Abstract pointer to AbstractVSSVertex for QCD radiation from incoming scalar
	*/
	map<ShowerInteraction,AbstractVSSVertexPtr> incomingVertex_;

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
	*/
	map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex1_;

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
	*/
	map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex2_;

	/**
	* Spin density matrix
	*/
	mutable RhoDMatrix rho_;

	/**
	* Scalar wavefunction
	*/
	mutable ScalarWaveFunction swave_;

	/**
	* Spinor wavefunction
	*/
	mutable vector<SpinorWaveFunction> wave_;

	/**
	* Barred spinor wavefunction
	*/
	mutable vector<SpinorBarWaveFunction> wavebar_;

	/**
	* Spin density matrix for 3 body decay
	*/
	mutable RhoDMatrix rho3_;

	/**
	* Scalar wavefunction for 3 body decay
	*/
	mutable ScalarWaveFunction swave3_;

	/**
	* Spinor wavefunction for 3 body decay
	*/
	mutable vector<SpinorWaveFunction> wave3_;

	/**
	* Barred spinor wavefunction for 3 body decay
	*/
	mutable vector<SpinorBarWaveFunction> wavebar3_;

	/**
	* Vector wavefunction for 3 body decay
	*/
	mutable vector<VectorWaveFunction> gluon_;


	};

	}

	#endif /* HERWIG_SFFDecayer_H */
	diff --git a/Decay/General/SSVDecayer.cc b/Decay/General/SSVDecayer.cc
	--- a/Decay/General/SSVDecayer.cc
	+++ b/Decay/General/SSVDecayer.cc
	@@ -1,368 +1,372 @@
	// -- C++ --
	//
	// SSVDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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));
	perturbativeVertex_.push_back(dynamic_ptr_cast<VSSVertexPtr> (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();
	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
	outgoingVertexV_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[1].at(inter));
	}
	}
	}

	void SSVDecayer::persistentOutput(PersistentOStream & os) const {
	os << vertex_ << perturbativeVertex_
	<< incomingVertex_ << outgoingVertexS_
	<< outgoingVertexV_ << fourPointVertex_;
	}

	void SSVDecayer::persistentInput(PersistentIStream & is, int) {
	is >> vertex_ >> perturbativeVertex_
	>> incomingVertex_ >> outgoingVertexS_
	>> outgoingVertexV_ >> fourPointVertex_;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	-DescribeClass<SSVDecayer,GeneralTwoBodyDecayer>
	+DescribeClass<SSVDecayer,GeneralTwoBodyDecayer2>
	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");

	}

	-double SSVDecayer::me2(const int , const Particle & inpart,
	- const ParticleVector & decay,
	+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(decay[0]->dataPtr()->iSpin() != PDT::Spin0) swap(isc,ivec);
	+ 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>(&inpart),incoming);
	- swave_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),incoming);
	+ calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&part),incoming);
	+ swave_ = ScalarWaveFunction(part.momentum(),part.dataPtr(),incoming);
	// fix rho if no correlations
	fixRho(rho_);
	}
	- 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,false);
	- }
	VectorWaveFunction::
	- calculateWaveFunctions(vector_,decay[ivec],outgoing,false);
	- 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.;
	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.;
	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(inpart.dataPtr(),decay[0]->dataPtr(),
	- decay[1]->dataPtr());
	+ output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]);
	// return the answer
	return output;
	}

	Energy SSVDecayer:: partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const {
	if( inpart.second < outa.second + outb.second ) return ZERO;
	if(perturbativeVertex_.size()==1 &&
	perturbativeVertex_[0]) {
	double mu1sq(sqr(outa.second/inpart.second)),
	mu2sq(sqr(outb.second/inpart.second));
	tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
	if(outa.first->iSpin() == PDT::Spin0) {
	perturbativeVertex_[0]->setCoupling(sqr(inpart.second), outb.first, outa.first,in);
	}
	else {
	swap(mu1sq,mu2sq);
	perturbativeVertex_[0]->setCoupling(sqr(inpart.second), outa.first, outb.first,in);
	}
	double me2(0.);
	if(mu2sq == 0.)
	me2 = -2.*mu1sq - 2.;
	else
	me2 = ( sqr(mu2sq - mu1sq) - 2.*(mu2sq + mu1sq) + 1. )/mu2sq;
	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second, outa.second,
	outb.second);
	Energy output = pcmme2norm(perturbativeVertex_[0]->norm())/8./Constants::pi;
	// colour factor
	output *= colourFactor(inpart.first,outa.first,outb.first);
	// return the answer
	return output;
	}
	else {
	- return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
	+ return GeneralTwoBodyDecayer2::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 &&
	decay[ivect]->dataPtr()->iColour()==PDT::Colour8)
	swap(S,V);
	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
	+ const GeneralTwoBodyDecayer2::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 =
	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;
	}
	#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 =
	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) +=
	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 =
	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) +=
	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) +=
	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;
	}
	diff --git a/Decay/General/SSVDecayer.h b/Decay/General/SSVDecayer.h
	--- a/Decay/General/SSVDecayer.h
	+++ b/Decay/General/SSVDecayer.h
	@@ -1,224 +1,233 @@
	// -- C++ --
	//
	// SSVDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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.
	//
	#ifndef HERWIG_SSVDecayer_H
	#define HERWIG_SSVDecayer_H
	//
	// This is the declaration of the SSVDecayer class.
	//

	-#include "GeneralTwoBodyDecayer.h"
	+#include "GeneralTwoBodyDecayer2.h"
	#include "ThePEG/Helicity/Vertex/Scalar/VSSVertex.h"
	#include "ThePEG/Helicity/Vertex/Vector/VVVVertex.h"
	#include "ThePEG/Helicity/Vertex/Scalar/VVSSVertex.h"
	#include "ThePEG/Repository/EventGenerator.h"

	namespace Herwig {
	using namespace ThePEG;
	using Helicity::VSSVertexPtr;

	/** \ingroup Decay
	* The SSVDecayer class implements the decay of a scalar to a vector
	* and a scalar in a general model. It holds an VSSVertex pointer
	* that must be typecast from the VertexBase pointer held in
	- * GeneralTwoBodyDecayer. It implents the virtual functions me2() and
	+ * GeneralTwoBodyDecayer2. It implents the virtual functions me2() and
	* partialWidth().
	*
	- * @see GeneralTwoBodyDecayer
	+ * @see GeneralTwoBodyDecayer2
	*/
	-class SSVDecayer: public GeneralTwoBodyDecayer {
	+class SSVDecayer: public GeneralTwoBodyDecayer2 {

	public:

	/**
	* The default constructor.
	*/
	SSVDecayer() {}

	/** @name Virtual functions required by the Decayer class. */
	//@{
	/**
	- * Return the matrix element squared for a given mode and phase-space channel
	- * @param ichan The channel we are calculating the matrix element for.
	+ * Return the matrix element squared for a given mode and phase-space channel.
	+ * @param ichan The channel we are calculating the matrix element for.
	* @param part The decaying Particle.
	- * @param decay The particles produced in the decay.
	+ * @param outgoing The particles produced in the decay
	+ * @param momenta The momenta of the particles produced in the decay
	* @param meopt Option for the calculation of the matrix element
	* @return The matrix element squared for the phase-space configuration.
	*/
	- virtual double me2(const int ichan, const Particle & part,
	- const ParticleVector & decay, MEOption meopt) const;
	+ double me2(const int ichan,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	+ MEOption meopt) const;
	+
	+ /**
	+ * Construct the SpinInfos for the particles produced in the decay
	+ */
	+ virtual void constructSpinInfo(const Particle & part,
	+ ParticleVector outgoing) const;

	/**
	* Function to return partial Width
	* @param inpart The decaying particle.
	* @param outa One of the decay products.
	* @param outb The other decay product.
	*/
	virtual Energy partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const;

	/**
	* Has a POWHEG style correction
	*/
	virtual POWHEGType hasPOWHEGCorrection() {
	POWHEGType output = FSR;
	for(auto vertex : vertex_) {
	if(vertex->orderInAllCouplings()!=1) {
	output = No;
	break;
	}
	}
	return output;
	}

	/**
	* Three-body matrix element including additional QCD radiation
	*/
	virtual double threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter,
	MEOption meopt);

	/**
	* Set the information on the decay
	*/
	virtual void setDecayInfo(PDPtr incoming, PDPair outgoing,
	vector<VertexBasePtr>,
	map<ShowerInteraction,VertexBasePtr> &,
	const vector<map<ShowerInteraction,VertexBasePtr> > &,
	map<ShowerInteraction,VertexBasePtr>);
	//@}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	SSVDecayer & operator=(const SSVDecayer &);

	private:

	/**
	* Abstract pointer to AbstractFFVVertex
	*/
	vector<AbstractVSSVertexPtr> vertex_;

	/**
	* Pointer to the perturbative vertex
	*/
	vector<VSSVertexPtr> perturbativeVertex_;

	/**
	* Abstract pointer to AbstractVSSVertex for QCD radiation from incoming scalar
	*/
	map<ShowerInteraction,AbstractVSSVertexPtr> incomingVertex_;

	/**
	* Abstract pointer to AbstractVSSVertex for QCD radiation from outgoing scalar
	*/
	map<ShowerInteraction,AbstractVSSVertexPtr> outgoingVertexS_;

	/**
	* Abstract pointer to AbstractVVVVertex for QCD radiation from outgoing vector
	*/
	map<ShowerInteraction,AbstractVVVVertexPtr> outgoingVertexV_;

	/**
	* Abstract pointer to AbstractVVSSVertex for QCD radiation from 4 point vertex
	*/
	map<ShowerInteraction,AbstractVVSSVertexPtr> fourPointVertex_;

	/**
	* Spinor density matrix
	*/
	mutable RhoDMatrix rho_;

	/**
	* Scalar wavefunction
	*/
	mutable Helicity::ScalarWaveFunction swave_;

	/**
	* Vector wavefunction
	*/
	mutable vector<Helicity::VectorWaveFunction> vector_;

	/**
	* Spin density matrix for 3 body decay
	*/
	mutable RhoDMatrix rho3_;

	/**
	* Scalar wavefunction for 3 body decay
	*/
	mutable Helicity::ScalarWaveFunction swave3_;

	/**
	* Scalar wavefunction for 3 body decay
	*/
	mutable Helicity::ScalarWaveFunction scal_;

	/**
	* Vector wavefunction for 3 body decay
	*/
	mutable vector<Helicity::VectorWaveFunction> vector3_;

	/**
	* Vector wavefunction for 3 body decay
	*/
	mutable vector<Helicity::VectorWaveFunction> gluon_;

	};

	}

	#endif /* HERWIG_SSVDecayer_H */
	diff --git a/Decay/General/TFFDecayer.cc b/Decay/General/TFFDecayer.cc
	--- a/Decay/General/TFFDecayer.cc
	+++ b/Decay/General/TFFDecayer.cc
	@@ -1,353 +1,356 @@
	// -- C++ --
	//
	// TFFDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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 TFFDecayer class.
	//

	#include "TFFDecayer.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 "ThePEG/Helicity/WaveFunction/TensorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
	#include "Herwig/Utilities/Kinematics.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"

	using namespace Herwig;
	using namespace ThePEG::Helicity;
	IBPtr TFFDecayer::clone() const {
	return new_ptr(*this);
	}

	IBPtr TFFDecayer::fullclone() const {
	return new_ptr(*this);
	}

	void TFFDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
	vector<VertexBasePtr> vertex,
	map<ShowerInteraction,VertexBasePtr> &,
	const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
	map<ShowerInteraction,VertexBasePtr> fourV) {
	decayInfo(incoming,outgoing);
	for(auto vert : vertex) {
	vertex_ .push_back(dynamic_ptr_cast<AbstractFFTVertexPtr>(vert));
	perturbativeVertex_.push_back(dynamic_ptr_cast<FFTVertexPtr> (vert));
	}
	vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
	for(auto & inter : itemp) {
	fourPointVertex_[inter] = dynamic_ptr_cast<AbstractFFVTVertexPtr>(fourV.at(inter));
	outgoingVertex1_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr> (outV[0].at(inter));
	outgoingVertex2_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr> (outV[1].at(inter));
	}
	}

	void TFFDecayer::persistentOutput(PersistentOStream & os) const {
	os << vertex_ << perturbativeVertex_
	<< outgoingVertex1_ << outgoingVertex2_
	<< fourPointVertex_;
	}

	void TFFDecayer::persistentInput(PersistentIStream & is, int) {
	is >> vertex_ >> perturbativeVertex_
	>> outgoingVertex1_ >> outgoingVertex2_
	>> fourPointVertex_;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	-DescribeClass<TFFDecayer,GeneralTwoBodyDecayer>
	+DescribeClass<TFFDecayer,GeneralTwoBodyDecayer2>
	describeHerwigTFFDecayer("Herwig::TFFDecayer", "Herwig.so");

	void TFFDecayer::Init() {

	static ClassDocumentation<TFFDecayer> documentation
	("The TFFDecayer class implements the decay of a tensor particle "
	"to 2 fermions ");

	}

	-double TFFDecayer::me2(const int , const Particle & inpart,
	- const ParticleVector & decay,
	+void TFFDecayer::
	+constructSpinInfo(const Particle & part, ParticleVector decay) const {
	+ unsigned int iferm(0),ianti(1);
	+ if(decay[0]->id()>=0) swap(iferm,ianti);
	+ TensorWaveFunction::
	+ constructSpinInfo(tensors_,const_ptr_cast<tPPtr>(&part),
	+ incoming,true,false);
	+ SpinorBarWaveFunction::
	+ constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
	+ SpinorWaveFunction::
	+ constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
	+}
	+
	+double TFFDecayer::me2(const int,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	MEOption meopt) const {
	unsigned int iferm(0),ianti(1);
	if(!ME())
	ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin2,PDT::Spin1Half,PDT::Spin1Half)));
	- if(decay[0]->id()>=0) swap(iferm,ianti);
	+ if(outgoing[0]->id()>=0) swap(iferm,ianti);
	if(meopt==Initialize) {
	TensorWaveFunction::
	- calculateWaveFunctions(tensors_,rho_,const_ptr_cast<tPPtr>(&inpart),
	+ calculateWaveFunctions(tensors_,rho_,const_ptr_cast<tPPtr>(&part),
	incoming,false);
	// fix rho if no correlations
	fixRho(rho_);
	}
	- if(meopt==Terminate) {
	- TensorWaveFunction::
	- constructSpinInfo(tensors_,const_ptr_cast<tPPtr>(&inpart),
	- incoming,true,false);
	- SpinorBarWaveFunction::
	- constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
	- SpinorWaveFunction::
	- constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
	- return 0.;
	- }
	SpinorBarWaveFunction::
	- calculateWaveFunctions(wavebar_,decay[iferm],outgoing);
	+ calculateWaveFunctions(wavebar_,momenta[iferm],outgoing[iferm],Helicity::outgoing);
	SpinorWaveFunction::
	- calculateWaveFunctions(wave_ ,decay[ianti],outgoing);
	- Energy2 scale(sqr(inpart.mass()));
	+ calculateWaveFunctions(wave_ ,momenta[ianti],outgoing[ianti],Helicity::outgoing);
	+ Energy2 scale(sqr(part.mass()));
	unsigned int thel,fhel,ahel;
	for(thel=0;thel<5;++thel) {
	for(fhel=0;fhel<2;++fhel) {
	for(ahel=0;ahel<2;++ahel) {
	if(iferm > ianti) {
	(*ME())(thel,fhel,ahel) = 0.;
	for(auto vert : vertex_)
	(*ME())(thel,fhel,ahel) +=
	vert->evaluate(scale,wave_[ahel],
	wavebar_[fhel],tensors_[thel]);
	}
	else {
	(*ME())(thel,ahel,fhel) = 0.;
	for(auto vert : vertex_)
	(*ME())(thel,ahel,fhel) +=
	vert->evaluate(scale,wave_[ahel],
	wavebar_[fhel],tensors_[thel]);
	}
	}
	}
	}
	double output = (ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
	// colour and identical particle factors
	- output *= colourFactor(inpart.dataPtr(),decay[0]->dataPtr(),
	- decay[1]->dataPtr());
	+ output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]);
	// return the answer
	return output;
	}

	Energy TFFDecayer::partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const {
	if( inpart.second < outa.second + outb.second ) return ZERO;
	if(perturbativeVertex_.size()==1 &&
	perturbativeVertex_[0]) {
	Energy2 scale = sqr(inpart.second);
	tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
	perturbativeVertex_[0]->setCoupling(scale, in, outa.first, outb.first);
	double musq = sqr(outa.second/inpart.second);
	double b = sqrt(1- 4.*musq);
	double me2 = bb(5-2bb)scale/120.UnitRemoval::InvE2;
	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second,outa.second,
	outb.second);
	Energy output = norm(perturbativeVertex_[0]->norm())me2pcm/(8.*Constants::pi);
	// colour factor
	output *= colourFactor(inpart.first,outa.first,outb.first);
	// return the answer
	return output;
	}
	else {
	- return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
	+ return GeneralTwoBodyDecayer2::partialWidth(inpart,outa,outb);
	}
	}


	double TFFDecayer::threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt) {
	// work out which is the fermion and antifermion
	int ianti(0), iferm(1), iglu(2);
	int itype[2];
	for(unsigned int ix=0;ix<2;++ix) {
	if(decay[ix]->dataPtr()->CC()) itype[ix] = decay[ix]->id()>0 ? 0:1;
	else itype[ix] = 2;
	}
	if(itype[0]==0 && itype[1]!=0) swap(iferm, ianti);
	if(itype[0]==2 && itype[1]==1) swap(iferm, ianti);
	if(itype[0]==0 && itype[1]==0 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
	swap(iferm, ianti);
	if(itype[0]==1 && itype[1]==1 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
	swap(iferm, ianti);

	if(meopt==Initialize) {
	// create tensor wavefunction for decaying particle
	TensorWaveFunction::
	calculateWaveFunctions(tensors3_, rho3_, const_ptr_cast<tPPtr>(&inpart), incoming, false);
	}
	// setup spin information when needed
	if(meopt==Terminate) {
	TensorWaveFunction::
	constructSpinInfo(tensors3_, const_ptr_cast<tPPtr>(&inpart),incoming,true, false);
	SpinorBarWaveFunction::
	constructSpinInfo(wavebar3_ ,decay[iferm],outgoing,true);
	SpinorWaveFunction::
	constructSpinInfo(wave3_ ,decay[ianti],outgoing,true);
	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::Spin2, PDT::Spin1Half,
	PDT::Spin1Half, PDT::Spin1)));
	// create wavefunctions
	SpinorBarWaveFunction::
	calculateWaveFunctions(wavebar3_, decay[iferm],outgoing);
	SpinorWaveFunction::
	calculateWaveFunctions(wave3_ , decay[ianti],outgoing);
	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 (! (outgoingVertex1_[inter] && outgoingVertex2_[inter]))
	throw Exception()
	<< "Invalid vertices for QCD radiation in TFF decay in TFFDecayer::threeBodyME"
	<< Exception::runerror;

	// identify fermion and/or anti-fermion vertex
	AbstractFFVVertexPtr outgoingVertexF = outgoingVertex1_[inter];
	AbstractFFVVertexPtr outgoingVertexA = outgoingVertex2_[inter];

	if(outgoingVertex1_[inter]!=outgoingVertex2_[inter] &&
	outgoingVertex1_[inter]->isIncoming(getParticleData(decay[ianti]->id())))
	swap (outgoingVertexF, outgoingVertexA);

	if(! (inpart.dataPtr()->iColour()==PDT::Colour0)){
	throw Exception()
	<< "Invalid vertices for QCD radiation in TFF decay in TFFDecayer::threeBodyME"
	<< Exception::runerror;
	}

	Energy2 scale(sqr(inpart.mass()));

	- const GeneralTwoBodyDecayer::CFlow & colourFlow
	+ const GeneralTwoBodyDecayer2::CFlow & colourFlow
	= colourFlows(inpart, decay);
	double gs(0.);
	bool couplingSet(false);
	#ifdef GAUGE_CHECK
	double total=0.;
	#endif
	for(unsigned int it = 0; it < 5; ++it) {
	for(unsigned int ifm = 0; ifm < 2; ++ifm) {
	for(unsigned int ia = 0; ia < 2; ++ia) {
	for(unsigned int ig = 0; ig < 2; ++ig) {

	// radiation from outgoing fermion
	if((decay[iferm]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[iferm]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertexF);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[iferm]->dataPtr();
	if(off->CC()) off = off->CC();
	SpinorBarWaveFunction interS =
	outgoingVertexF->evaluate(scale,3,off,wavebar3_[ifm],
	gluon_[2*ig],decay[iferm]->mass());

	assert(wavebar3_[ifm].particle()->id()==interS.particle()->id());

	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,wave3_[ia], interS,tensors3_[it]);
	if(!couplingSet) {
	gs = abs(outgoingVertexF->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[1].size();++ix) {
	(*ME[colourFlow[1][ix].first])(it, ifm, ia, ig) +=
	colourFlow[1][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}

	// radiation from outgoing antifermion
	if((decay[ianti]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[ianti]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertexA);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[ianti]->dataPtr();
	if(off->CC()) off = off->CC();
	SpinorWaveFunction interS =
	outgoingVertexA->evaluate(scale,3,off,wave3_[ia],
	gluon_[2*ig],decay[ianti]->mass());

	assert(wave3_[ia].particle()->id()==interS.particle()->id());

	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,interS,wavebar3_[ifm],tensors3_[it]);
	if(!couplingSet) {
	gs = abs(outgoingVertexA->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[2].size();++ix) {
	(*ME[colourFlow[2][ix].first])(it, ifm, ia, ig) +=
	colourFlow[2][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}

	// radiation from 4 point vertex
	if (fourPointVertex_[inter]) {
	Complex diag = fourPointVertex_[inter]->evaluate(scale, wave3_[ia], wavebar3_[ifm],
	gluon_[2*ig], tensors3_[it]);
	for(unsigned int ix=0;ix<colourFlow[3].size();++ix) {
	(*ME[colourFlow[3][ix].first])(it, ifm, ia, 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;
	}

	diff --git a/Decay/General/TFFDecayer.h b/Decay/General/TFFDecayer.h
	--- a/Decay/General/TFFDecayer.h
	+++ b/Decay/General/TFFDecayer.h
	@@ -1,223 +1,232 @@
	// -- C++ --
	//
	// TFFDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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.
	//
	#ifndef HERWIG_TFFDecayer_H
	#define HERWIG_TFFDecayer_H
	//
	// This is the declaration of the TFFDecayer class.
	//

	-#include "GeneralTwoBodyDecayer.h"
	+#include "GeneralTwoBodyDecayer2.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Helicity/Vertex/Tensor/FFTVertex.h"
	#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
	#include "ThePEG/Helicity/Vertex/Tensor/FFVTVertex.h"

	namespace Herwig {
	using namespace ThePEG;
	using Helicity::FFTVertexPtr;

	/** \ingroup Decay
	* The TFFDecayer class implements the decay of a tensor
	* to 2 fermions in a general model. It holds an FFTVertex pointer
	* that must be typecast from the VertexBase pointer held in
	- * GeneralTwoBodyDecayer. It implents the virtual functions me2() and
	+ * GeneralTwoBodyDecayer2. It implents the virtual functions me2() and
	* partialWidth().
	*
	- * @see GeneralTwoBodyDecayer
	+ * @see GeneralTwoBodyDecayer2
	*/
	-class TFFDecayer: public GeneralTwoBodyDecayer {
	+class TFFDecayer: public GeneralTwoBodyDecayer2 {

	public:

	/**
	* The default constructor.
	*/
	TFFDecayer() {}

	/** @name Virtual functions required by the Decayer class. */
	//@{
	/**
	* Return the matrix element squared for a given mode and phase-space channel.
	- * @param ichan The channel we are calculating the matrix element for.
	+ * @param ichan The channel we are calculating the matrix element for.
	* @param part The decaying Particle.
	- * @param decay The particles produced in the decay.
	- * @param meopt Option for the matrix element
	+ * @param outgoing The particles produced in the decay
	+ * @param momenta The momenta of the particles produced in the decay
	+ * @param meopt Option for the calculation of the matrix element
	* @return The matrix element squared for the phase-space configuration.
	*/
	- virtual double me2(const int ichan, const Particle & part,
	- const ParticleVector & decay, MEOption meopt) const;
	+ double me2(const int ichan,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	+ MEOption meopt) const;
	+
	+ /**
	+ * Construct the SpinInfos for the particles produced in the decay
	+ */
	+ virtual void constructSpinInfo(const Particle & part,
	+ ParticleVector outgoing) const;

	/**
	* Function to return partial Width
	* @param inpart The decaying particle.
	* @param outa One of the decay products.
	* @param outb The other decay product.
	*/
	virtual Energy partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const;

	/**
	* Has a POWHEG style correction
	*/
	virtual POWHEGType hasPOWHEGCorrection() {
	POWHEGType output = FSR;
	for(auto vertex : vertex_) {
	if(vertex->orderInAllCouplings()!=1) {
	output = No;
	break;
	}
	}
	return output;
	}

	/**
	* Three-body matrix element including additional QCD radiation
	*/
	virtual double threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt);

	/**
	* Set the information on the decay
	*/
	virtual void setDecayInfo(PDPtr incoming, PDPair outgoing,
	vector<VertexBasePtr>,
	map<ShowerInteraction,VertexBasePtr> &,
	const vector<map<ShowerInteraction,VertexBasePtr> > &,
	map<ShowerInteraction,VertexBasePtr>);
	//@}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	TFFDecayer & operator=(const TFFDecayer &);

	private:

	/**
	* Abstract pointer to AbstractFFTVertex
	*/
	vector<AbstractFFTVertexPtr> vertex_;

	/**
	* Pointer to the perturbative vertex
	*/
	vector<FFTVertexPtr> perturbativeVertex_;

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
	*/
	map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex1_;

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
	*/
	map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex2_;

	/**
	* Abstract pointer to AbstractFFVTVertex for QCD radiation from 4 point vertex
	*/
	map<ShowerInteraction,AbstractFFVTVertexPtr> fourPointVertex_;

	/**
	* Spin density matrix
	*/
	mutable RhoDMatrix rho_;

	/**
	* Polarization tensors for the decaying particle
	*/
	mutable vector<TensorWaveFunction> tensors_;

	/**
	* Spinors for the decay products
	*/
	mutable vector<SpinorWaveFunction> wave_;

	/**
	* Barred spinors for the decay products
	*/
	mutable vector<SpinorBarWaveFunction> wavebar_;

	/**
	* Spin density matrix for 3 body decay
	*/
	mutable RhoDMatrix rho3_;

	/**
	* Tensor wavefunction for 3 body decay
	*/
	mutable vector<TensorWaveFunction> tensors3_;

	/**
	* Spinor wavefunction for 3 body decay
	*/
	mutable vector<SpinorWaveFunction> wave3_;

	/**
	* Barred spinor wavefunction for 3 body decay
	*/
	mutable vector<SpinorBarWaveFunction> wavebar3_;

	/**
	* Vector wavefunction for 3 body decay
	*/
	mutable vector<VectorWaveFunction> gluon_;

	};

	}

	#endif /* HERWIG_TFFDecayer_H */
	diff --git a/Decay/General/VFFDecayer.cc b/Decay/General/VFFDecayer.cc
	--- a/Decay/General/VFFDecayer.cc
	+++ b/Decay/General/VFFDecayer.cc
	@@ -1,472 +1,475 @@
	// -- C++ --
	//
	// VFFDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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 VFFDecayer class.
	//

	#include "VFFDecayer.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 "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
	#include "Herwig/Utilities/Kinematics.h"
	#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"
	using namespace Herwig;
	using namespace ThePEG::Helicity;
	IBPtr VFFDecayer::clone() const {
	return new_ptr(*this);
	}

	IBPtr VFFDecayer::fullclone() const {
	return new_ptr(*this);
	}

	void VFFDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
	vector<VertexBasePtr> vertex,
	map<ShowerInteraction,VertexBasePtr> & inV,
	const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
	map<ShowerInteraction,VertexBasePtr> ) {
	decayInfo(incoming,outgoing);
	for(auto vert : vertex) {
	vertex_ .push_back(dynamic_ptr_cast<AbstractFFVVertexPtr>(vert));
	perturbativeVertex_.push_back(dynamic_ptr_cast<FFVVertexPtr> (vert));
	}
	vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
	for(auto & inter : itemp) {
	incomingVertex_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(inV.at(inter));
	outgoingVertex1_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[0].at(inter));
	outgoingVertex2_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[1].at(inter));
	}
	}

	void VFFDecayer::persistentOutput(PersistentOStream & os) const {
	os << vertex_ << perturbativeVertex_
	<< incomingVertex_ << outgoingVertex1_
	<< outgoingVertex2_;
	}

	void VFFDecayer::persistentInput(PersistentIStream & is, int) {
	is >> vertex_ >> perturbativeVertex_
	>> incomingVertex_ >> outgoingVertex1_
	>> outgoingVertex2_;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	-DescribeClass<VFFDecayer,GeneralTwoBodyDecayer>
	+DescribeClass<VFFDecayer,GeneralTwoBodyDecayer2>
	describeHerwigVFFDecayer("Herwig::VFFDecayer", "Herwig.so");

	void VFFDecayer::Init() {

	static ClassDocumentation<VFFDecayer> documentation
	("The VFFDecayer implements the matrix element for the"
	" decay of a vector to fermion-antifermion pair");

	}

	-double VFFDecayer::me2(const int , const Particle & inpart,
	- const ParticleVector & decay,
	+void VFFDecayer::
	+constructSpinInfo(const Particle & part, ParticleVector decay) const {
	+ int iferm(1),ianti(0);
	+ if(decay[0]->id()>0) swap(iferm,ianti);
	+ VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast<tPPtr>(&part),
	+ incoming,true,false);
	+ SpinorBarWaveFunction::
	+ constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
	+ SpinorWaveFunction::
	+ constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
	+}
	+
	+double VFFDecayer::me2(const int,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	MEOption meopt) const {
	int iferm(1),ianti(0);
	- if(decay[0]->id()>0) swap(iferm,ianti);
	+ if(outgoing[0]->id()>0) swap(iferm,ianti);
	if(!ME())
	ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,PDT::Spin1Half,PDT::Spin1Half)));
	if(meopt==Initialize) {
	VectorWaveFunction::calculateWaveFunctions(vectors_,rho_,
	- const_ptr_cast<tPPtr>(&inpart),
	- incoming,false);
	+ const_ptr_cast<tPPtr>(&part),
	+ incoming,false);
	// fix rho if no correlations
	fixRho(rho_);
	}
	- if(meopt==Terminate) {
	- VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast<tPPtr>(&inpart),
	- incoming,true,false);
	- SpinorBarWaveFunction::
	- constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
	- SpinorWaveFunction::
	- constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
	- return 0.;
	- }
	SpinorBarWaveFunction::
	- calculateWaveFunctions(wavebar_,decay[iferm],outgoing);
	+ calculateWaveFunctions(wavebar_,momenta[iferm],outgoing[iferm],Helicity::outgoing);
	SpinorWaveFunction::
	- calculateWaveFunctions(wave_ ,decay[ianti],outgoing);
	+ calculateWaveFunctions(wave_ ,momenta[ianti],outgoing[ianti],Helicity::outgoing);
	// compute the matrix element
	- Energy2 scale(inpart.mass()*inpart.mass());
	+ Energy2 scale(part.mass()*part.mass());
	for(unsigned int ifm = 0; ifm < 2; ++ifm) { //loop over fermion helicities
	for(unsigned int ia = 0; ia < 2; ++ia) {// loop over antifermion helicities
	for(unsigned int vhel = 0; vhel < 3; ++vhel) {//loop over vector helicities
	- if(iferm > ianti) {
	- (*ME())(vhel, ia, ifm) = 0.;
	- for(auto vert : vertex_)
	- (*ME())(vhel, ia, ifm) +=
	- vert->evaluate(scale,wave_[ia],
	- wavebar_[ifm],vectors_[vhel]);
	- }
	- else {
	- (*ME())(vhel,ifm,ia)= 0.;
	- for(auto vert : vertex_)
	- (*ME())(vhel,ifm,ia) +=
	- vert->evaluate(scale,wave_[ia],
	- wavebar_[ifm],vectors_[vhel]);
	- }
	+ if(iferm > ianti) {
	+ (*ME())(vhel, ia, ifm) = 0.;
	+ for(auto vert : vertex_)
	+ (*ME())(vhel, ia, ifm) +=
	+ vert->evaluate(scale,wave_[ia],
	+ wavebar_[ifm],vectors_[vhel]);
	+ }
	+ else {
	+ (*ME())(vhel,ifm,ia)= 0.;
	+ for(auto vert : vertex_)
	+ (*ME())(vhel,ifm,ia) +=
	+ vert->evaluate(scale,wave_[ia],
	+ wavebar_[ifm],vectors_[vhel]);
	+ }
	}
	}
	}
	double output=(ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
	// colour and identical particle factors
	- output *= colourFactor(inpart.dataPtr(),decay[0]->dataPtr(),
	- decay[1]->dataPtr());
	+ output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]);
	// return the answer
	return output;
	}

	Energy VFFDecayer::partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const {
	if( inpart.second < outa.second + outb.second ) return ZERO;
	if(perturbativeVertex_.size()==1 &&
	perturbativeVertex_[0]) {
	double mu1(outa.second/inpart.second), mu2(outb.second/inpart.second);
	tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
	perturbativeVertex_[0]->setCoupling(sqr(inpart.second), outa.first, outb.first,in);
	Complex cl(perturbativeVertex_[0]->left()), cr(perturbativeVertex_[0]->right());
	double me2 = (norm(cl) + norm(cr))*( sqr(sqr(mu1) - sqr(mu2))
	+ sqr(mu1) + sqr(mu2) - 2.)
	- 6.(clconj(cr) + crconj(cl)).real()mu1*mu2;
	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second,outa.second,
	outb.second);
	Energy output = -norm(perturbativeVertex_[0]->norm())me2pcm /
	(24.*Constants::pi);
	// colour factor
	output *= colourFactor(inpart.first,outa.first,outb.first);
	// return the answer
	return output;
	}
	else {
	- return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
	+ return GeneralTwoBodyDecayer2::partialWidth(inpart,outa,outb);
	}
	}


	double VFFDecayer::threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt) {
	// work out which is the fermion and antifermion
	int ianti(0), iferm(1), iglu(2);
	int itype[2];
	for(unsigned int ix=0;ix<2;++ix) {
	if(decay[ix]->dataPtr()->CC()) itype[ix] = decay[ix]->id()>0 ? 0:1;
	else itype[ix] = 2;
	}
	if(itype[0]==0 && itype[1]!=0) swap(iferm, ianti);
	if(itype[0]==2 && itype[1]==1) swap(iferm, ianti);
	if(itype[0]==0 && itype[1]==0 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
	swap(iferm, ianti);
	if(itype[0]==1 && itype[1]==1 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
	swap(iferm, ianti);

	if(meopt==Initialize) {
	// create vector wavefunction for decaying particle
	VectorWaveFunction::calculateWaveFunctions(vector3_, rho3_, const_ptr_cast<tPPtr>(&inpart),
	incoming, false);
	}
	// setup spin information when needed
	if(meopt==Terminate) {
	VectorWaveFunction::
	constructSpinInfo(vector3_ ,const_ptr_cast<tPPtr>(&inpart),outgoing,true,false);
	SpinorBarWaveFunction::
	constructSpinInfo(wavebar3_,decay[iferm],outgoing,true);
	SpinorWaveFunction::
	constructSpinInfo(wave3_ ,decay[ianti],outgoing,true);
	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::Spin1, PDT::Spin1Half,
	PDT::Spin1Half, PDT::Spin1)));
	// create wavefunctions
	SpinorBarWaveFunction::
	calculateWaveFunctions(wavebar3_, decay[iferm],outgoing);
	SpinorWaveFunction::
	calculateWaveFunctions(wave3_ , decay[ianti],outgoing);
	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

	// identify fermion and/or anti-fermion vertex
	AbstractFFVVertexPtr outgoingVertexF;
	AbstractFFVVertexPtr outgoingVertexA;
	identifyVertices(iferm, ianti, inpart, decay, outgoingVertexF, outgoingVertexA,inter);

	Energy2 scale(sqr(inpart.mass()));

	- const GeneralTwoBodyDecayer::CFlow & colourFlow
	+ const GeneralTwoBodyDecayer2::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 ifm = 0; ifm < 2; ++ifm) {
	for(unsigned int ia = 0; ia < 2; ++ia) {
	for(unsigned int ig = 0; ig < 2; ++ig) {
	// radiation from the incoming vector
	if((inpart.dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(inpart.dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(incomingVertex_[inter]);

	VectorWaveFunction vectorInter =
	incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),vector3_[iv],
	gluon_[2*ig],inpart.mass());

	assert(vector3_[iv].particle()->id()==vectorInter.particle()->id());

	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,wave3_[ia],wavebar3_[ifm],vectorInter);
	if(!couplingSet) {
	gs = abs(incomingVertex_[inter]->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
	(*ME[colourFlow[0][ix].first])(iv, ia, ifm, ig) +=
	colourFlow[0][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}

	// radiation from outgoing fermion
	if((decay[iferm]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[iferm]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertexF);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[iferm]->dataPtr();
	if(off->CC()) off = off->CC();
	SpinorBarWaveFunction interS =
	outgoingVertexF->evaluate(scale,3,off,wavebar3_[ifm],
	gluon_[2*ig],decay[iferm]->mass());

	assert(wavebar3_[ifm].particle()->id()==interS.particle()->id());

	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,wave3_[ia], interS,vector3_[iv]);
	if(!couplingSet) {
	gs = abs(outgoingVertexF->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[1].size();++ix) {
	(*ME[colourFlow[1][ix].first])(iv, ia, ifm, ig) +=
	colourFlow[1][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}

	// radiation from outgoing antifermion
	if((decay[ianti]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[ianti]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertexA);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[ianti]->dataPtr();
	if(off->CC()) off = off->CC();
	SpinorWaveFunction interS =
	outgoingVertexA->evaluate(scale,3,off,wave3_[ia],
	gluon_[2*ig],decay[ianti]->mass());

	assert(wave3_[ia].particle()->id()==interS.particle()->id());

	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,interS,wavebar3_[ifm],vector3_[iv]);
	if(!couplingSet) {
	gs = abs(outgoingVertexA->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[2].size();++ix) {
	(*ME[colourFlow[2][ix].first])(iv, ia, ifm, ig) +=
	colourFlow[2][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 output
	return output;
	}


	void VFFDecayer::identifyVertices(const int iferm, const int ianti,
	const Particle & inpart, const ParticleVector & decay,
	AbstractFFVVertexPtr & outgoingVertexF,
	AbstractFFVVertexPtr & outgoingVertexA,
	ShowerInteraction inter){
	// QCD vertices
	if(inter==ShowerInteraction::QCD) {
	// work out which fermion each outgoing vertex corresponds to
	// two outgoing vertices
	if( inpart.dataPtr() ->iColour()==PDT::Colour0 &&
	((decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar) \|\|
	(decay[iferm]->dataPtr()->iColour()==PDT::Colour8 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour8))){
	if(outgoingVertex1_[inter]==outgoingVertex2_[inter]){
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))){
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else if (outgoingVertex2_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))){
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	}
	else if(inpart.dataPtr() ->iColour()==PDT::Colour8 &&
	decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar){
	if(outgoingVertex1_[inter]==outgoingVertex2_[inter]){
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))){
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else if (outgoingVertex2_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))){
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	}

	// one outgoing vertex
	else if(inpart.dataPtr()->iColour()==PDT::Colour3){
	if(decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour0){
	if (outgoingVertex1_[inter]) outgoingVertexF = outgoingVertex1_[inter];
	else if(outgoingVertex2_[inter]) outgoingVertexF = outgoingVertex2_[inter];
	}
	else if (decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour8){
	if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[ianti]->dataPtr()))){
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	else {
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	}
	}
	else if(inpart.dataPtr()->iColour()==PDT::Colour3bar){
	if(decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar &&
	decay[iferm]->dataPtr()->iColour()==PDT::Colour0){
	if (outgoingVertex1_[inter]) outgoingVertexA = outgoingVertex1_[inter];
	else if(outgoingVertex2_[inter]) outgoingVertexA = outgoingVertex2_[inter];
	}
	else if (decay[iferm]->dataPtr()->iColour()==PDT::Colour8 &&
	decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar){
	if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))){
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else {
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	}
	}
	else if(inpart.dataPtr()->iColour()==PDT::Colour6 \|\|
	inpart.dataPtr()->iColour()==PDT::Colour6bar) {
	if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
	outgoingVertexF = outgoingVertex1_[inter];
	outgoingVertexA = outgoingVertex2_[inter];
	}
	else {
	outgoingVertexF = outgoingVertex2_[inter];
	outgoingVertexA = outgoingVertex1_[inter];
	}
	}

	if (! ((incomingVertex_[inter] && (outgoingVertexF \|\| outgoingVertexA)) \|\|
	( outgoingVertexF && outgoingVertexA)))
	throw Exception()
	<< "Invalid vertices for QCD radiation in VFF decay in VFFDecayer::identifyVertices"
	<< Exception::runerror;
	}
	// QED
	else {
	if(decay[iferm]->dataPtr()->charged()) {
	if (outgoingVertex1_[inter] &&
	outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr())))
	outgoingVertexF = outgoingVertex1_[inter];
	else
	outgoingVertexF = outgoingVertex2_[inter];
	}
	if(decay[ianti]->dataPtr()->charged()) {
	if (outgoingVertex1_[inter] &&
	outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[ianti]->dataPtr())))
	outgoingVertexA = outgoingVertex1_[inter];
	else
	outgoingVertexA = outgoingVertex2_[inter];
	}
	}
	}

	diff --git a/Decay/General/VFFDecayer.h b/Decay/General/VFFDecayer.h
	--- a/Decay/General/VFFDecayer.h
	+++ b/Decay/General/VFFDecayer.h
	@@ -1,232 +1,241 @@
	// -- C++ --
	//
	// VFFDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 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.
	//
	#ifndef HERWIG_VFFDecayer_H
	#define HERWIG_VFFDecayer_H
	//
	// This is the declaration of the VFFDecayer class.
	//

	-#include "GeneralTwoBodyDecayer.h"
	+#include "GeneralTwoBodyDecayer2.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
	#include "ThePEG/Helicity/Vertex/Vector/VVVVertex.h"

	namespace Herwig {
	using namespace ThePEG;
	using Helicity::FFVVertexPtr;

	/** \ingroup Decay
	* The VFFDecayer class implements the decay of a vector
	* to 2 fermions in a general model. It holds an FFVVertex pointer
	* that must be typecast from the VertexBase pointer held in
	- * GeneralTwoBodyDecayer. It implents the virtual functions me2() and
	+ * GeneralTwoBodyDecayer2. It implents the virtual functions me2() and
	* partialWidth().
	*
	- * @see GeneralTwoBodyDecayer
	+ * @see GeneralTwoBodyDecayer2
	*/
	-class VFFDecayer: public GeneralTwoBodyDecayer {
	+class VFFDecayer: public GeneralTwoBodyDecayer2 {

	public:

	/**
	* The default constructor.
	*/
	VFFDecayer() {}

	public:

	/** @name Virtual functions required by the Decayer class. */
	//@{
	/**
	- * Return the matrix element squared for a given mode and phase-space channel.
	- * @param ichan The channel we are calculating the matrix element for.
	+ * Return the matrix element squared for a given mode and phase-space channel.
	+ * @param ichan The channel we are calculating the matrix element for.
	* @param part The decaying Particle.
	- * @param decay The particles produced in the decay.
	- * @param meopt Option for the matrix element
	+ * @param outgoing The particles produced in the decay
	+ * @param momenta The momenta of the particles produced in the decay
	+ * @param meopt Option for the calculation of the matrix element
	* @return The matrix element squared for the phase-space configuration.
	*/
	- virtual double me2(const int ichan, const Particle & part,
	- const ParticleVector & decay, MEOption meopt) const;
	+ double me2(const int ichan,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	+ MEOption meopt) const;
	+
	+ /**
	+ * Construct the SpinInfos for the particles produced in the decay
	+ */
	+ virtual void constructSpinInfo(const Particle & part,
	+ ParticleVector outgoing) const;

	/**
	* Function to return partial Width
	* @param inpart The decaying particle.
	* @param outa One of the decay products.
	* @param outb The other decay product.
	*/
	virtual Energy partialWidth(PMPair inpart, PMPair outa,
	PMPair outb) const;

	/**
	* Has a POWHEG style correction
	*/
	virtual POWHEGType hasPOWHEGCorrection() {
	POWHEGType output = FSR;
	for(auto vertex : vertex_) {
	if(vertex->orderInAllCouplings()!=1) {
	output = No;
	break;
	}
	}
	return output;
	}

	/**
	* Three-body matrix element including additional QCD radiation
	*/
	virtual double threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt);

	/**
	* Indentify outgoing vertices for the fermion and antifermion
	*/
	void identifyVertices(const int iferm, const int ianti,
	const Particle & inpart, const ParticleVector & decay,
	AbstractFFVVertexPtr & abstractOutgoingVertexF,
	AbstractFFVVertexPtr & abstractOutgoingVertexA,
	ShowerInteraction inter);

	/**
	* Set the information on the decay
	*/
	virtual void setDecayInfo(PDPtr incoming, PDPair outgoing,
	vector<VertexBasePtr>,
	map<ShowerInteraction,VertexBasePtr> &,
	const vector<map<ShowerInteraction,VertexBasePtr> > &,
	map<ShowerInteraction,VertexBasePtr>);
	//@}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	VFFDecayer & operator=(const VFFDecayer &);

	private:

	/**
	* Abstract pointer to AbstractFFVVertex
	*/
	vector<AbstractFFVVertexPtr> vertex_;

	/**
	* Pointer to the perturbative vertex
	*/
	vector<FFVVertexPtr> perturbativeVertex_;

	/**
	* Abstract pointer to AbstractVVVVertex for QCD radiation from incoming vector
	*/
	map<ShowerInteraction,AbstractVVVVertexPtr> incomingVertex_;

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
	*/
	map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex1_;

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
	*/
	map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex2_;

	/**
	* Spin density matrix
	*/
	mutable RhoDMatrix rho_;

	/**
	* Polarization vectors for the decaying particle
	*/
	mutable vector<VectorWaveFunction> vectors_;

	/**
	* Spinors for the decay products
	*/
	mutable vector<SpinorWaveFunction> wave_;

	/**
	* Barred spinors for the decay products
	*/
	mutable vector<SpinorBarWaveFunction> wavebar_;
	/**
	* Spin density matrix for 3 body decay
	*/
	mutable RhoDMatrix rho3_;

	/**
	* Scalar wavefunction for 3 body decay
	*/
	mutable vector<VectorWaveFunction> vector3_;

	/**
	* Spinor wavefunction for 3 body decay
	*/
	mutable vector<SpinorWaveFunction> wave3_;

	/**
	* Barred spinor wavefunction for 3 body decay
	*/
	mutable vector<SpinorBarWaveFunction> wavebar3_;

	/**
	* Vector wavefunction for 3 body decay
	*/
	mutable vector<VectorWaveFunction> gluon_;

	};

	}

	#endif /* HERWIG_VFFDecayer_H */

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 8:16 PM (11 h, 49 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4226349
Default Alt Text: (95 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions