No OneTemporary
Actions

Size

117 KB

Subscribers

None

View Options

	diff --git a/Decay/General/GeneralTwoBodyDecayer.h b/Decay/General/GeneralTwoBodyDecayer.h
	--- a/Decay/General/GeneralTwoBodyDecayer.h
	+++ b/Decay/General/GeneralTwoBodyDecayer.h
	@@ -1,306 +1,271 @@
	// -- C++ --
	//
	// GeneralTwoBodyDecayer.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_GeneralTwoBodyDecayer_H
	#define HERWIG_GeneralTwoBodyDecayer_H
	//
	// This is the declaration of the GeneralTwoBodyDecayer class.
	//

	#include "Herwig/Decay/PerturbativeDecayer.h"
	#include "Herwig/Decay/DecayPhaseSpaceMode.h"
	#include "ThePEG/Helicity/Vertex/VertexBase.h"
	#include "GeneralTwoBodyDecayer.fh"
	+#include "GeneralTwoBodyDecayer2.h"

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

	/** \ingroup Decay
	* The GeneralTwoBodyDecayer class is designed to be the base class
	* for 2 body decays for some general model. It inherits from
	* PerturbativeDecayer and implements the modeNumber() virtual function
	* that is the same for all of the decays. A decayer for
	* a specific spin configuration should inherit from this and implement
	* the me2() and partialWidth() member functions. The colourConnections()
	* member should be called from inside me2() in the inheriting decayer
	* to set up the colour lines.
	*
	* @see \ref GeneralTwoBodyDecayerInterfaces "The interfaces"
	* defined for GeneralTwoBodyDecayer.
	* @see PerturbativeDecayer
	*/
	class GeneralTwoBodyDecayer: public PerturbativeDecayer {

	public:

	/** A ParticleData ptr and (possible) mass pair.*/
	typedef pair<tcPDPtr, Energy> PMPair;

	public:

	/**
	* The default constructor.
	*/
	GeneralTwoBodyDecayer() : maxWeight_(1.), colour_(1,DVector(1,1.))
	{}


	/** @name Virtual functions required by the Decayer class. */
	//@{
	/**
	* For a given decay mode and a given particle instance, perform the
	* decay and return the decay products. As this is the base class this
	* is not implemented.
	* @return The vector of particles produced in the decay.
	*/
	virtual ParticleVector decay(const Particle & parent,
	const tPDVector & children) const;

	/**
	* Which of the possible decays is required
	* @param cc Is this mode the charge conjugate
	* @param parent The decaying particle
	* @param children The decay products
	*/
	virtual int modeNumber(bool & cc, tcPDPtr parent,const tPDVector & children) const;

	/**
	* 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 calculation of the matrix element
	* @return The matrix element squared for the phase-space configuration.
	*/
	virtual double me2(const int , const Particle & part,
	const ParticleVector & decay, MEOption meopt) const = 0;

	/**
	* 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;

	/**
	* Specify the \f$1\to2\f$ matrix element to be used in the running width
	* calculation.
	* @param dm The DecayMode
	* @param mecode The code for the matrix element as described
	* in the GenericWidthGenerator class.
	* @param coupling The coupling for the matrix element.
	* @return True if the the order of the particles in the
	* decayer is the same as the DecayMode tag.
	*/
	virtual bool twoBodyMEcode(const DecayMode & dm, int & mecode,
	double & coupling) const;

	/**
	* An overidden member to calculate a branching ratio for a certain
	* particle instance.
	* @param dm The DecayMode of the particle
	* @param p The particle object
	* @param oldbrat The branching fraction given in the DecayMode object
	*/
	virtual double brat(const DecayMode & dm, const Particle & p,
	double oldbrat) const;
	//@}

	/**
	* 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>) =0;

	protected:

	/** @name Functions used by inheriting decayers. */
	//@{
	/**
	* Set integration weight
	* @param wgt Maximum integration weight
	*/
	void setWeight(double wgt) { maxWeight_ = wgt; }

	/**
	* Set colour connections
	* @param parent Parent particle
	* @param out Particle vector containing particles to
	* connect colour lines
	*/
	void colourConnections(const Particle & parent,
	const ParticleVector & out) const;

	/**
	* Compute the spin and colour factor
	*/
	double colourFactor(tcPDPtr in, tcPDPtr out1, tcPDPtr out2) const;

	/**
	* Calculate matrix element ratio R/B
	*/
	double matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
	const ParticleVector & decay3, MEOption meopt,
	ShowerInteraction inter);

	/**
	* Set the information on the decay
	*/
	void decayInfo(PDPtr incoming, PDPair outgoing);
	//@}

	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 Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit();

	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();
	//@}

	protected:

	/**
	* Member for the generation of additional hard radiation
	*/
	//@{
	/**
	* Return the matrix of colour factors
	*/
	typedef vector<pair<int,double > > CFlowPairVec;
	typedef vector<CFlowPairVec> CFlow;

	const vector<DVector> & getColourFactors(const Particle & inpart,
	const ParticleVector & decay,
	unsigned int & nflow);

	const CFlow & colourFlows(const Particle & inpart,
	const ParticleVector & decay);

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

	//@}

	private:

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

	private:

	/**
	* Store the incoming particle
	*/
	PDPtr incoming_;

	/**
	* Outgoing particles
	*/
	vector<PDPtr> outgoing_;

	/**
	* Maximum weight for integration
	*/
	double maxWeight_;

	/**
	* Store colour factors for ME calc.
	*/
	vector<DVector> colour_;

	};

	-
	-/**
	- * Write a map with ShowerInteraction as the key
	- */
	-template<typename T, typename Cmp, typename A>
	-inline PersistentOStream & operator<<(PersistentOStream & os,
	- const map<ShowerInteraction,T,Cmp,A> & m) {
	- os << m.size();
	- if(m.find(ShowerInteraction::QCD)!=m.end()) {
	- os << 0 << m.at(ShowerInteraction::QCD);
	- }
	- if(m.find(ShowerInteraction::QED)!=m.end()) {
	- os << 1 << m.at(ShowerInteraction::QED);
	- }
	- return os;
	}

	-/**
	- * Read a map with ShowerInteraction as the key
	- */
	-template <typename T, typename Cmp, typename A>
	-inline PersistentIStream & operator>>(PersistentIStream & is, map<ShowerInteraction,T,Cmp,A> & m) {
	- m.clear();
	- long size;
	- int k;
	- is >> size;
	- while ( size-- && is ) {
	- is >> k;
	- if(k==0)
	- is >> m[ShowerInteraction::QCD];
	- else if(k==1)
	- is >> m[ShowerInteraction::QED];
	- else
	- assert(false);
	- }
	- return is;
	-}
	-}

	#endif /* HERWIG_GeneralTwoBodyDecayer_H */
	diff --git a/Decay/General/GeneralTwoBodyDecayer2.cc b/Decay/General/GeneralTwoBodyDecayer2.cc
	--- a/Decay/General/GeneralTwoBodyDecayer2.cc
	+++ b/Decay/General/GeneralTwoBodyDecayer2.cc
	@@ -1,819 +1,823 @@

	// -- C++ --
	//
	// GeneralTwoBodyDecayer2.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 GeneralTwoBodyDecayer2 class.
	//

	#include "GeneralTwoBodyDecayer2.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/Utilities/Exception.h"
	#include "Herwig/Shower/RealEmissionProcess.h"

	using namespace Herwig;

	ParticleVector GeneralTwoBodyDecayer2::decay(const Particle & parent,
	const tPDVector & children) const {

	// return empty vector if products heavier than parent
	Energy mout(ZERO);
	for(tPDVector::const_iterator it=children.begin();
	it!=children.end();++it) mout+=(**it).massMin();
	if(mout>parent.mass()) return ParticleVector();
	// generate the decay
	bool cc;
	int imode=modeNumber(cc,parent.dataPtr(),children);
	// generate the kinematics
	ParticleVector decay=generate(generateIntermediates(),cc,imode,parent);
	// make the colour connections
	colourConnections(parent, decay);
	// return the answer
	return decay;
	}

	void GeneralTwoBodyDecayer2::doinit() {
	PerturbativeDecayer2::doinit();
	assert( incoming_ && outgoing_.size()==2);
	//create phase space mode
	addMode(new_ptr(PhaseSpaceMode(incoming_,{outgoing_[0],outgoing_[1]},
	maxWeight_)));
	}

	int GeneralTwoBodyDecayer2::modeNumber(bool & cc, tcPDPtr parent,
	const tPDVector & children) const {
	long parentID = parent->id();
	long id1 = children[0]->id();
	long id2 = children[1]->id();
	cc = false;
	long out1 = outgoing_[0]->id();
	long out2 = outgoing_[1]->id();
	if( parentID == incoming_->id() &&
	((id1 == out1 && id2 == out2) \|\|
	(id1 == out2 && id2 == out1)) ) {
	return 0;
	}
	else if(incoming_->CC() && parentID == incoming_->CC()->id()) {
	cc = true;
	if( outgoing_[0]->CC()) out1 = outgoing_[0]->CC()->id();
	if( outgoing_[1]->CC()) out2 = outgoing_[1]->CC()->id();
	if((id1 == out1 && id2 == out2) \|\|
	(id1 == out2 && id2 == out1)) return 0;
	}
	return -1;
	}

	void GeneralTwoBodyDecayer2::
	colourConnections(const Particle & parent,
	const ParticleVector & out) const {
	PDT::Colour incColour(parent.data().iColour());
	PDT::Colour outaColour(out[0]->data().iColour());
	PDT::Colour outbColour(out[1]->data().iColour());
	//incoming colour singlet
	if(incColour == PDT::Colour0) {
	// colour triplet-colourantitriplet
	if((outaColour == PDT::Colour3 && outbColour == PDT::Colour3bar) \|\|
	(outaColour == PDT::Colour3bar && outbColour == PDT::Colour3)) {
	bool ac(out[0]->id() < 0);
	out[0]->colourNeighbour(out[1],!ac);
	}
	//colour octet
	else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour8) {
	out[0]->colourNeighbour(out[1]);
	out[0]->antiColourNeighbour(out[1]);
	}
	// colour singlets
	else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour0) {
	}
	// unknown
	else
	throw Exception() << "Unknown outgoing colours for decaying "
	<< "colour singlet in "
	<< "GeneralTwoBodyDecayer2::colourConnections "
	<< outaColour << " " << outbColour
	<< Exception::runerror;
	}
	//incoming colour triplet
	else if(incColour == PDT::Colour3) {
	// colour triplet + singlet
	if(outaColour == PDT::Colour3 && outbColour == PDT::Colour0) {
	out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	}
	//opposite order
	else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour3) {
	out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	}
	// octet + triplet
	else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour3) {
	out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	out[1]->antiColourNeighbour(out[0]);
	}
	//opposite order
	else if(outaColour == PDT::Colour3 && outbColour == PDT::Colour8) {
	out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	out[0]->antiColourNeighbour(out[1]);
	}
	else if(outaColour == PDT::Colour3bar && outaColour == PDT::Colour3bar) {
	tColinePtr col[2] = {ColourLine::create(out[0],true),
	ColourLine::create(out[1],true)};
	parent.colourLine()->setSinkNeighbours(col[0],col[1]);
	}
	else
	throw Exception() << "Unknown outgoing colours for decaying "
	<< "colour triplet in "
	<< "GeneralTwoBodyDecayer2::colourConnections() "
	<< outaColour << " " << outbColour
	<< Exception::runerror;
	}
	// incoming colour anti triplet
	else if(incColour == PDT::Colour3bar) {
	// colour antitriplet +singlet
	if(outaColour == PDT::Colour3bar && outbColour == PDT::Colour0) {
	out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	}
	//opposite order
	else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour3bar) {
	out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	}
	//octet + antitriplet
	else if(outaColour == PDT::Colour3bar && outbColour == PDT::Colour8) {
	out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	out[0]->colourNeighbour(out[1]);
	}
	//opposite order
	else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour3bar) {
	out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	out[1]->colourNeighbour(out[0]);
	}
	else if(outaColour == PDT::Colour3 && outbColour == PDT::Colour3) {
	tColinePtr col[2] = {ColourLine::create(out[0]),
	ColourLine::create(out[1])};
	parent.antiColourLine()->setSourceNeighbours(col[0],col[1]);
	}
	else
	throw Exception() << "Unknown outgoing colours for decaying "
	<< "colour antitriplet "
	<< "in GeneralTwoBodyDecayer2::colourConnections() "
	<< outaColour << " " << outbColour
	<< Exception::runerror;
	}
	//incoming colour octet
	else if(incColour == PDT::Colour8) {
	// triplet-antitriplet
	if(outaColour == PDT::Colour3&&outbColour == PDT::Colour3bar) {
	out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	}
	// opposite order
	else if(outbColour == PDT::Colour3&&outaColour == PDT::Colour3bar) {
	out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	}
	// neutral octet
	else if(outaColour == PDT::Colour0&&outbColour == PDT::Colour8) {
	out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	}
	else if(outbColour == PDT::Colour0&&outaColour == PDT::Colour8) {
	out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	}
	else
	throw Exception() << "Unknown outgoing colours for decaying "
	<< "colour octet "
	<< "in GeneralTwoBodyDecayer2::colourConnections() "
	<< outaColour << " " << outbColour
	<< Exception::runerror;
	}
	else if(incColour == PDT::Colour6) {
	if(outaColour == PDT::Colour3 && outbColour == PDT::Colour3) {
	tPPtr tempParent = const_ptr_cast<tPPtr>(&parent);
	Ptr<MultiColour>::pointer parentColour =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>
	(tempParent->colourInfo());

	tColinePtr line1 = const_ptr_cast<tColinePtr>(parentColour->colourLines()[0]);
	line1->addColoured(dynamic_ptr_cast<tPPtr>(out[0]));

	tColinePtr line2 = const_ptr_cast<tColinePtr>(parentColour->colourLines()[1]);
	line2->addColoured(dynamic_ptr_cast<tPPtr>(out[1]));
	}
	else
	throw Exception() << "Unknown outgoing colours for decaying "
	<< "colour sextet "
	<< "in GeneralTwoBodyDecayer2::colourConnections() "
	<< outaColour << " " << outbColour
	<< Exception::runerror;
	}
	else if(incColour == PDT::Colour6bar) {
	if(outaColour == PDT::Colour3bar && outbColour == PDT::Colour3bar) {
	tPPtr tempParent = const_ptr_cast<tPPtr>(&parent);
	Ptr<MultiColour>::pointer parentColour =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>
	(tempParent->colourInfo());

	tColinePtr line1 = const_ptr_cast<tColinePtr>(parentColour->antiColourLines()[0]);
	line1->addAntiColoured(dynamic_ptr_cast<tPPtr>(out[0]));

	tColinePtr line2 = const_ptr_cast<tColinePtr>(parentColour->antiColourLines()[1]);
	line2->addAntiColoured(dynamic_ptr_cast<tPPtr>(out[1]));
	}
	else
	throw Exception() << "Unknown outgoing colours for decaying "
	<< "colour anti-sextet "
	<< "in GeneralTwoBodyDecayer2::colourConnections() "
	<< outaColour << " " << outbColour
	<< Exception::runerror;
	}
	else
	throw Exception() << "Unknown incoming colour in "
	<< "GeneralTwoBodyDecayer2::colourConnections() "
	<< incColour
	<< Exception::runerror;
	}

	bool GeneralTwoBodyDecayer2::twoBodyMEcode(const DecayMode & dm, int & mecode,
	double & coupling) const {
	assert(dm.parent()->id() == incoming_->id());
	ParticleMSet::const_iterator pit = dm.products().begin();
	long id1 = (*pit)->id();
	++pit;
	long id2 = (*pit)->id();
	long id1t(outgoing_[0]->id()), id2t(outgoing_[1]->id());
	mecode = -1;
	coupling = 1.;
	if( id1 == id1t && id2 == id2t ) {
	return true;
	}
	else if( id1 == id2t && id2 == id1t ) {
	return false;
	}
	else
	assert(false);
	return false;
	}


	void GeneralTwoBodyDecayer2::persistentOutput(PersistentOStream & os) const {
	os << incoming_ << outgoing_ << maxWeight_;
	}

	void GeneralTwoBodyDecayer2::persistentInput(PersistentIStream & is, int) {
	is >> incoming_ >> outgoing_ >> maxWeight_;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeAbstractClass<GeneralTwoBodyDecayer2,PerturbativeDecayer2>
	describeHerwigGeneralTwoBodyDecayer2("Herwig::GeneralTwoBodyDecayer2", "Herwig.so");

	void GeneralTwoBodyDecayer2::Init() {

	static ClassDocumentation<GeneralTwoBodyDecayer2> documentation
	("This class is designed to be a base class for all 2 body decays"
	"in a general model");

	}

	double GeneralTwoBodyDecayer2::brat(const DecayMode &, const Particle & p,
	double oldbrat) const {
	ParticleVector children = p.children();
	if( children.size() != 2 \|\| !p.data().widthGenerator() )
	return oldbrat;

	// partial width for this mode
	Energy scale = p.mass();
	Energy pwidth =
	partialWidth( make_pair(p.dataPtr(), scale),
	make_pair(children[0]->dataPtr(), children[0]->mass()),
	make_pair(children[1]->dataPtr(), children[1]->mass()) );
	Energy width = p.data().widthGenerator()->width(p.data(), scale);
	return pwidth/width;
	}

	void GeneralTwoBodyDecayer2::doinitrun() {
	PerturbativeDecayer2::doinitrun();
	for(unsigned int ix=0;ix<numberModes();++ix) {
	double fact = pow(1.5,int(mode(ix)->incoming().first->iSpin())-1);
	mode(ix)->maxWeight(fact*mode(ix)->maxWeight());
	}
	}

	double GeneralTwoBodyDecayer2::colourFactor(tcPDPtr in, tcPDPtr out1,
	tcPDPtr out2) const {
	// identical particle symmetry factor
	double output = out1->id()==out2->id() ? 0.5 : 1.;
	// colour neutral incoming particle
	if(in->iColour()==PDT::Colour0) {
	// both colour neutral
	if(out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour0)
	output *= 1.;
	// colour triplet/ antitriplet
	else if((out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour3bar) \|\|
	(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour3 ) ) {
	output *= 3.;
	}
	// colour octet colour octet
	else if(out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour8 ) {
	output *= 8.;
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for decay colour neutral particle in "
	<< "GeneralTwoBodyDecayer2::colourFactor() for "
	<< in->PDGName() << " -> "
	<< out1->PDGName() << " " << out2->PDGName()
	<< Exception::runerror;
	}
	// triplet
	else if(in->iColour()==PDT::Colour3) {
	// colour triplet + neutral
	if((out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour3) \|\|
	(out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour0) ) {
	output *= 1.;
	}
	// colour triplet + octet
	else if((out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour3) \|\|
	(out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour8) ) {
	output *= 4./3.;
	}
	// colour anti triplet anti triplet
	else if(out1->iColour()==PDT::Colour3bar &&
	out2->iColour()==PDT::Colour3bar) {
	output *= 2.;
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for decay colour triplet particle in "
	<< "GeneralTwoBodyDecayer2::colourFactor() for "
	<< in->PDGName() << " -> "
	<< out1->PDGName() << " " << out2->PDGName()
	<< Exception::runerror;
	}
	// anti triplet
	else if(in->iColour()==PDT::Colour3bar) {
	// colour anti triplet + neutral
	if((out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour3bar ) \|\|
	(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour0 ) ) {
	output *= 1.;
	}
	// colour anti triplet + octet
	else if((out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour3bar ) \|\|
	(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour8 ) ) {
	output *= 4./3.;
	}
	// colour triplet triplet
	else if(out1->iColour()==PDT::Colour3 &&
	out2->iColour()==PDT::Colour3) {
	output *= 2.;
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for decay colour anti triplet particle in "
	<< "GeneralTwoBodyDecayer2::colourFactor() for "
	<< in->PDGName() << " -> "
	<< out1->PDGName() << " " << out2->PDGName()
	<< Exception::runerror;
	}
	else if(in->iColour()==PDT::Colour8) {
	// colour octet + neutral
	if((out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour8 ) \|\|
	(out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour0 ) ) {
	output *= 1.;
	}
	// colour triplet/antitriplet
	else if((out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour3bar) \|\|
	(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour3 ) ) {
	output *= 0.5;
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for decay colour octet particle in "
	<< "GeneralTwoBodyDecayer2::colourFactor() for "
	<< in->PDGName() << " -> "
	<< out1->PDGName() << " " << out2->PDGName()
	<< Exception::runerror;
	}
	else if(in->iColour()==PDT::Colour6) {
	// colour sextet -> triplet triplet
	if( out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour3 ) {
	output *= 1.;
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for decay colour sextet particle in "
	<< "GeneralTwoBodyDecayer2::colourFactor() for "
	<< in->PDGName() << " -> "
	<< out1->PDGName() << " " << out2->PDGName()
	<< Exception::runerror;
	}
	else if(in->iColour()==PDT::Colour6bar) {
	// colour sextet -> triplet triplet
	if( out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour3bar ) {
	output *= 1.;
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for decay colour anti-sextet particle in "
	<< "GeneralTwoBodyDecayer2::colourFactor() for "
	<< in->PDGName() << " -> "
	<< out1->PDGName() << " " << out2->PDGName()
	<< Exception::runerror;
	}
	else
	throw Exception() << "Unknown colour "
	<< in->iColour() << " for the decaying particle in "
	<< "GeneralTwoBodyDecayer2::colourFactor() for "
	<< in->PDGName() << " -> "
	<< out1->PDGName() << " " << out2->PDGName()
	<< Exception::runerror;
	return output;
	}

	Energy GeneralTwoBodyDecayer2::partialWidth(PMPair inpart, PMPair outa,
	- PMPair outb) const {
	+ PMPair outb) const {
	// select the number of the mode
	tPDVector children;
	children.push_back(const_ptr_cast<PDPtr>(outa.first));
	children.push_back(const_ptr_cast<PDPtr>(outb.first));
	bool cc;
	int nmode=modeNumber(cc,inpart.first,children);
	tcPDPtr newchild[2] = {mode(nmode)->outgoing()[0],
	mode(nmode)->outgoing()[1]};
	// make the particles
	Lorentz5Momentum pparent = Lorentz5Momentum(inpart.second);
	PPtr parent = inpart.first->produceParticle(pparent);
	- Lorentz5Momentum pout[2];
	+ vector<Lorentz5Momentum> pout(2);
	double ctheta,phi;
	Kinematics::generateAngles(ctheta,phi);
	Kinematics::twoBodyDecay(pparent, outa.second, outb.second,
	ctheta, phi,pout[0],pout[1]);
	if( ( !cc && outa.first!=newchild[0]) \|\|
	( cc && !(( outa.first->CC() && outa.first->CC() == newchild[0])\|\|
	( !outa.first->CC() && outa.first == newchild[0]) )))
	swap(pout[0],pout[1]);
	- ParticleVector decay;
	- decay.push_back(newchild[0]->produceParticle(pout[0]));
	- decay.push_back(newchild[1]->produceParticle(pout[1]));
	- double me = me2(-1,*parent,decay,Initialize);
	+ tPDVector out = {const_ptr_cast<tPDPtr>(outa.first),
	+ const_ptr_cast<tPDPtr>(outb.first)};
	+ double me = me2(-1,*parent,out,pout,Initialize);
	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second,
	outa.second, outb.second);
	return me/(8.Constants::pi)pcm;
	}

	void GeneralTwoBodyDecayer2::decayInfo(PDPtr incoming, PDPair outgoing) {
	incoming_=incoming;
	outgoing_.clear();
	outgoing_.push_back(outgoing.first );
	outgoing_.push_back(outgoing.second);
	}

	double GeneralTwoBodyDecayer2::matrixElementRatio(const Particle & inpart,
	- const ParticleVector & decay2,
	- const ParticleVector & decay3,
	- MEOption meopt,
	- ShowerInteraction inter) {
	+ const ParticleVector & decay2,
	+ const ParticleVector & decay3,
	+ MEOption meopt,
	+ ShowerInteraction inter) {
	// calculate R/B
	- double B = me2 (0, inpart, decay2, meopt);
	+ tPDVector outgoing = {const_ptr_cast<tPDPtr>(decay2[0]->dataPtr()),
	+ const_ptr_cast<tPDPtr>(decay2[1]->dataPtr())};
	+ const vector<Lorentz5Momentum> mom = {decay2[0]->momentum(),
	+ decay2[1]->momentum()};
	+
	+ double B = me2 (0, inpart, outgoing,mom, meopt);
	double R = threeBodyME(0, inpart, decay3, inter, meopt);
	return R/B;

	}

	const vector<DVector> & GeneralTwoBodyDecayer2::getColourFactors(const Particle & inpart,
	const ParticleVector & decay,
	unsigned int & nflow) {
	// calculate the colour factors for the three-body decay
	vector<int> sing,trip,atrip,oct,sex,asex;
	for(unsigned int it=0;it<decay.size();++it) {
	if (decay[it]->dataPtr()->iColour() == PDT::Colour0 ) sing. push_back(it);
	else if(decay[it]->dataPtr()->iColour() == PDT::Colour3 ) trip. push_back(it);
	else if(decay[it]->dataPtr()->iColour() == PDT::Colour3bar ) atrip.push_back(it);
	else if(decay[it]->dataPtr()->iColour() == PDT::Colour8 ) oct. push_back(it);
	else if(decay[it]->dataPtr()->iColour() == PDT::Colour6 ) sex. push_back(it);
	else if(decay[it]->dataPtr()->iColour() == PDT::Colour6bar ) asex. push_back(it);
	}

	// identical particle symmetry factor
	double symFactor=1.;
	if (( sing.size()==2 && decay[ sing[0]]->id()==decay[ sing[1]]->id()) \|\|
	( trip.size()==2 && decay[ trip[0]]->id()==decay[ trip[1]]->id()) \|\|
	(atrip.size()==2 && decay[atrip[0]]->id()==decay[atrip[1]]->id()) \|\|
	( oct.size()==2 && decay[ oct[0]]->id()==decay[ oct[1]]->id()) \|\|
	( sex.size()==2 && decay[ sex[0]]->id()==decay[ sex[1]]->id()) \|\|
	( asex.size()==2 && decay[ asex[0]]->id()==decay[ asex[1]]->id()))
	symFactor /= 2.;
	else if (oct.size()==3 &&
	decay[oct[0]]->id()==decay[oct[1]]->id() &&
	decay[oct[0]]->id()==decay[oct[2]]->id())
	symFactor /= 6.;

	colour_ = vector<DVector>(1,DVector(1,symFactor*1.));

	// decaying colour singlet
	if(inpart.dataPtr()->iColour() == PDT::Colour0) {
	if(trip.size()==1 && atrip.size()==1 && oct.size()==1) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor*4.));
	}
	else if(trip.size()==1 && atrip.size()==1 && sing.size()==1) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor*3.));
	}
	else if (oct.size()==3){
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor*24.));
	}
	else if(sing.size()==3) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor));
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for decay colour scalar particle in "
	<< "GeneralTwoBodyDecayer2::getColourFactors() for "
	<< inpart. dataPtr()->PDGName() << " -> "
	<< decay[0]->dataPtr()->PDGName() << " "
	<< decay[1]->dataPtr()->PDGName() << " "
	<< decay[2]->dataPtr()->PDGName()
	<< Exception::runerror;
	}
	// decaying colour triplet
	else if(inpart.dataPtr()->iColour() == PDT::Colour3) {
	if(trip.size()==1 && sing.size()==1 && oct.size()==1) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor*4./3.));
	}
	else if(trip.size()==1 && sing.size()==2) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor));
	}
	else if(trip.size()==1 && oct.size()==2) {
	nflow = 2;
	colour_.clear();
	colour_.resize(2,DVector(2,0.));
	colour_[0][0] = symFactor16./9.; colour_[0][1] = -symFactor2./9.;
	colour_[1][0] = -symFactor2./9.; colour_[1][1] = symFactor16./9.;
	}
	else if(atrip.size()==2 && sing.size()==1) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor*2.));
	}
	else if(atrip.size()==2 && oct.size()==1) {
	nflow = 2;
	colour_.clear();
	colour_.resize(2,DVector(2,0.));
	colour_[0][0] = 8./3.symFactor; colour_[0][1] =-4./3.symFactor;
	colour_[1][0] = -4./3.symFactor; colour_[1][1] = 8./3.symFactor;
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for decay colour triplet particle in "
	<< "GeneralTwoBodyDecayer2::getColourFactors() for "
	<< inpart. dataPtr()->PDGName() << " -> "
	<< decay[0]->dataPtr()->PDGName() << " "
	<< decay[1]->dataPtr()->PDGName() << " "
	<< decay[2]->dataPtr()->PDGName()
	<< Exception::runerror;
	}
	// decaying colour anti-triplet
	else if(inpart.dataPtr()->iColour() == PDT::Colour3bar) {
	if(atrip.size()==1 && sing.size()==1 && oct.size()==1) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor*4./3.));
	}
	else if(atrip.size()==1 && sing.size()==2) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor));
	}
	else if(atrip.size()==1 && oct.size()==2){
	nflow = 2;
	colour_.clear();
	colour_ .resize(2,DVector(2,0.));
	colour_[0][0] = symFactor16./9.; colour_[0][1] = -symFactor2./9.;
	colour_[1][0] = -symFactor2./9.; colour_[1][1] = symFactor16./9.;
	}
	else if(trip.size()==2 && sing.size()==1) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor*2.));
	}
	else if(trip.size()==2 && oct.size()==1) {
	nflow = 2;
	colour_.clear();
	colour_.resize(2,DVector(2,0.));
	colour_[0][0] = 8./3.symFactor; colour_[0][1] =-4./3.symFactor;
	colour_[1][0] = -4./3.symFactor; colour_[1][1] = 8./3.symFactor;
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for decay colour anti-triplet particle in "
	<< "GeneralTwoBodyDecayer2::getColourFactors() for "
	<< inpart. dataPtr()->PDGName() << " -> "
	<< decay[0]->dataPtr()->PDGName() << " "
	<< decay[1]->dataPtr()->PDGName() << " "
	<< decay[2]->dataPtr()->PDGName()
	<< Exception::runerror;
	}
	// decaying colour octet
	else if(inpart.dataPtr()->iColour() == PDT::Colour8) {
	if(oct.size()==1 && trip.size()==1 && atrip.size()==1) {
	nflow = 2;
	colour_.clear();
	colour_.resize(2,DVector(2,0.));
	colour_[0][0] = symFactor2./3. ; colour_[0][1] = -symFactor1./12.;
	colour_[1][0] = -symFactor1./12.; colour_[1][1] = symFactor2./3. ;
	}
	else if (sing.size()==1 && trip.size()==1 && atrip.size()==1) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor*0.5));
	}
	else if (oct.size()==2 && sing.size()==1) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor*3.));
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for a decaying colour octet particle in "
	<< "GeneralTwoBodyDecayer2::getColourFactors() for "
	<< inpart. dataPtr()->PDGName() << " -> "
	<< decay[0]->dataPtr()->PDGName() << " "
	<< decay[1]->dataPtr()->PDGName() << " "
	<< decay[2]->dataPtr()->PDGName()
	<< Exception::runerror;
	}
	// Sextet
	else if(inpart.dataPtr()->iColour() == PDT::Colour6) {
	if(trip.size()==2 && sing.size()==1) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor));
	}
	else if(trip.size()==2 && oct.size()==1) {
	nflow = 2;
	colour_.clear();
	colour_.resize(2,DVector(2,0.));
	colour_[0][0] = 4./3.symFactor; colour_[0][1] = 1./3.symFactor;
	colour_[1][0] = 1./3.symFactor; colour_[1][1] = 4./3.symFactor;
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for a decaying colour sextet particle in "
	<< "GeneralTwoBodyDecayer2::getColourFactors() for "
	<< inpart. dataPtr()->PDGName() << " -> "
	<< decay[0]->dataPtr()->PDGName() << " "
	<< decay[1]->dataPtr()->PDGName() << " "
	<< decay[2]->dataPtr()->PDGName()
	<< Exception::runerror;
	}
	// anti Sextet
	else if(inpart.dataPtr()->iColour() == PDT::Colour6bar) {
	if(atrip.size()==2 && sing.size()==1) {
	nflow = 1;
	colour_ = vector<DVector>(1,DVector(1,symFactor));
	}
	else if(atrip.size()==2 && oct.size()==1) {
	nflow = 2;
	colour_.clear();
	colour_.resize(2,DVector(2,0.));
	colour_[0][0] = 4./3.symFactor; colour_[0][1] = 1./3.symFactor;
	colour_[1][0] = 1./3.symFactor; colour_[1][1] = 4./3.symFactor;
	}
	else
	throw Exception() << "Unknown colour for the outgoing particles"
	<< " for a decaying colour anti-sextet particle in "
	<< "GeneralTwoBodyDecayer2::getColourFactors() for "
	<< inpart. dataPtr()->PDGName() << " -> "
	<< decay[0]->dataPtr()->PDGName() << " "
	<< decay[1]->dataPtr()->PDGName() << " "
	<< decay[2]->dataPtr()->PDGName()
	<< Exception::runerror;
	}
	else
	throw Exception() << "Unknown colour for the decaying particle in "
	<< "GeneralTwoBodyDecayer2::getColourFactors() for "
	<< inpart. dataPtr()->PDGName() << " -> "
	<< decay[0]->dataPtr()->PDGName() << " "
	<< decay[1]->dataPtr()->PDGName() << " "
	<< decay[2]->dataPtr()->PDGName()
	<< Exception::runerror;
	return colour_;
	}

	const GeneralTwoBodyDecayer2::CFlow &
	GeneralTwoBodyDecayer2::colourFlows(const Particle & inpart,
	const ParticleVector & decay) {
	// static initialization of commonly used colour structures
	static const CFlow init = CFlow(3, CFlowPairVec(1, make_pair(0, 1.)));
	static CFlow tripflow = init;
	static CFlow atripflow = init;
	static CFlow octflow = init;
	static CFlow sexflow = init;
	static CFlow epsflow = init;
	static const CFlow fpflow = CFlow(4, CFlowPairVec(1, make_pair(0, 1.)));

	static bool initialized = false;

	if (! initialized) {
	tripflow[2].resize(2, make_pair(0,1.));
	tripflow[2][0] = make_pair(0, 1.);
	tripflow[2][1] = make_pair(1,-1.);
	tripflow[1][0] = make_pair(1, 1.);

	atripflow[1].resize(2, make_pair(0,1.));
	atripflow[1][0] = make_pair(0, 1.);
	atripflow[1][1] = make_pair(1,-1.);
	atripflow[2][0] = make_pair(1, 1.);

	octflow[0].resize(2, make_pair(0,1.));
	octflow[0][0] = make_pair(0,-1.);
	octflow[0][1] = make_pair(1, 1.);
	octflow[2][0] = make_pair(1, 1.);

	sexflow[0].resize(2, make_pair(0,1.));
	sexflow[0][1] = make_pair(1,1.);
	sexflow[2][0] = make_pair(1,1.);

	epsflow[0].resize(2, make_pair(0,-1.));
	epsflow[0][0] = make_pair(0,-1.);
	epsflow[0][1] = make_pair(1,-1.);
	epsflow[2][0] = make_pair(1,1.);
	initialized = true;
	}


	// main function body
	int sing=0,trip=0,atrip=0,oct=0,sex=0,asex=0;
	for (size_t it=0; it<decay.size(); ++it) {
	switch ( decay[it]->dataPtr()->iColour() ) {
	case PDT::Colour0: ++sing; break;
	case PDT::Colour3: ++trip; break;
	case PDT::Colour3bar: ++atrip; break;
	case PDT::Colour8: ++oct; break;
	case PDT::Colour6: ++sex; break;
	case PDT::Colour6bar: ++asex; break;
	/// @todo: handle these better
	case PDT::ColourUndefined: break;
	case PDT::Coloured: break;
	}
	}

	const CFlow * retval = 0;
	// decaying colour triplet
	if(inpart.dataPtr()->iColour() == PDT::Colour3 &&
	trip==1 && oct==2) {
	retval = &tripflow;
	}
	// decaying colour anti-triplet
	else if(inpart.dataPtr()->iColour() == PDT::Colour3bar &&
	atrip==1 && oct==2){
	retval = &atripflow;
	}
	// decaying colour octet
	else if(inpart.dataPtr()->iColour() == PDT::Colour8 &&
	oct==1 && trip==1 && atrip==1) {
	retval = &octflow;
	}
	// decaying colour sextet
	else if(inpart.dataPtr()->iColour() ==PDT::Colour6 &&
	oct==1 && trip==2) {
	retval = &sexflow;
	}
	// decaying anti colour sextet
	else if(inpart.dataPtr()->iColour() ==PDT::Colour6bar &&
	oct==1 && atrip==2) {
	retval = &sexflow;
	}
	// decaying colour triplet (eps)
	else if(inpart.dataPtr()->iColour() == PDT::Colour3 &&
	atrip==2 && oct==1) {
	retval = &epsflow;
	}
	// decaying colour anti-triplet (eps)
	else if(inpart.dataPtr()->iColour() == PDT::Colour3bar &&
	trip==2 && oct==1){
	retval = &epsflow;
	}
	else {
	retval = &fpflow;
	}

	return *retval;
	}

	double GeneralTwoBodyDecayer2::threeBodyME(const int , const Particle &,
	const ParticleVector &,
	ShowerInteraction, MEOption) {
	throw Exception() << "Base class GeneralTwoBodyDecayer2::threeBodyME() "
	<< "called, should have an implementation in the inheriting class"
	<< Exception::runerror;
	return 0.;
	}
	diff --git a/Decay/General/GeneralTwoBodyDecayer2.h b/Decay/General/GeneralTwoBodyDecayer2.h
	--- a/Decay/General/GeneralTwoBodyDecayer2.h
	+++ b/Decay/General/GeneralTwoBodyDecayer2.h
	@@ -1,306 +1,295 @@
	// -- C++ --
	//
	// GeneralTwoBodyDecayer.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_GeneralTwoBodyDecayer2_H
	#define HERWIG_GeneralTwoBodyDecayer2_H
	//
	// This is the declaration of the GeneralTwoBodyDecayer2 class.
	//

	#include "Herwig/Decay/PerturbativeDecayer2.h"
	#include "Herwig/Decay/PhaseSpaceMode.h"
	#include "ThePEG/Helicity/Vertex/VertexBase.h"
	#include "GeneralTwoBodyDecayer2.fh"

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

	/** \ingroup Decay
	* The GeneralTwoBodyDecayer2 class is designed to be the base class
	* for 2 body decays for some general model. It inherits from
	* PerturbativeDecayer2 and implements the modeNumber() virtual function
	* that is the same for all of the decays. A decayer for
	* a specific spin configuration should inherit from this and implement
	* the me2() and partialWidth() member functions. The colourConnections()
	* member should be called from inside me2() in the inheriting decayer
	* to set up the colour lines.
	*
	* @see \ref GeneralTwoBodyDecayer2Interfaces "The interfaces"
	* defined for GeneralTwoBodyDecayer2.
	* @see PerturbativeDecayer2
	*/
	class GeneralTwoBodyDecayer2: public PerturbativeDecayer2 {

	public:

	/** A ParticleData ptr and (possible) mass pair.*/
	typedef pair<tcPDPtr, Energy> PMPair;

	public:

	/**
	* The default constructor.
	*/
	GeneralTwoBodyDecayer2() : maxWeight_(1.), colour_(1,DVector(1,1.))
	{}


	/** @name Virtual functions required by the Decayer2 class. */
	//@{
	/**
	* For a given decay mode and a given particle instance, perform the
	* decay and return the decay products. As this is the base class this
	* is not implemented.
	* @return The vector of particles produced in the decay.
	*/
	virtual ParticleVector decay(const Particle & parent,
	const tPDVector & children) const;

	/**
	* Which of the possible decays is required
	* @param cc Is this mode the charge conjugate
	* @param parent The decaying particle
	* @param children The decay products
	*/
	virtual int modeNumber(bool & cc, tcPDPtr parent,const tPDVector & children) const;
	-
	- /**
	- * 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 calculation of the matrix element
	- * @return The matrix element squared for the phase-space configuration.
	- */
	- virtual double me2(const int , const Particle & part,
	- const ParticleVector & decay, MEOption meopt) const = 0;

	/**
	* 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;

	/**
	* Specify the \f$1\to2\f$ matrix element to be used in the running width
	* calculation.
	* @param dm The DecayMode
	* @param mecode The code for the matrix element as described
	* in the GenericWidthGenerator class.
	* @param coupling The coupling for the matrix element.
	* @return True if the the order of the particles in the
	* decayer is the same as the DecayMode tag.
	*/
	virtual bool twoBodyMEcode(const DecayMode & dm, int & mecode,
	double & coupling) const;

	/**
	* An overidden member to calculate a branching ratio for a certain
	* particle instance.
	* @param dm The DecayMode of the particle
	* @param p The particle object
	* @param oldbrat The branching fraction given in the DecayMode object
	*/
	virtual double brat(const DecayMode & dm, const Particle & p,
	double oldbrat) const;
	//@}

	/**
	* 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>) =0;

	protected:

	/** @name Functions used by inheriting decayers. */
	//@{
	/**
	* Set integration weight
	* @param wgt Maximum integration weight
	*/
	void setWeight(double wgt) { maxWeight_ = wgt; }

	/**
	* Set colour connections
	* @param parent Parent particle
	* @param out Particle vector containing particles to
	* connect colour lines
	*/
	void colourConnections(const Particle & parent,
	const ParticleVector & out) const;

	/**
	* Compute the spin and colour factor
	*/
	double colourFactor(tcPDPtr in, tcPDPtr out1, tcPDPtr out2) const;

	/**
	* Calculate matrix element ratio R/B
	*/
	double matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
	const ParticleVector & decay3, MEOption meopt,
	ShowerInteraction inter);

	/**
	* Set the information on the decay
	*/
	void decayInfo(PDPtr incoming, PDPair outgoing);
	//@}

	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 Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit();

	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();
	//@}

	protected:

	/**
	* Member for the generation of additional hard radiation
	*/
	//@{
	/**
	* Return the matrix of colour factors
	*/
	typedef vector<pair<int,double > > CFlowPairVec;
	typedef vector<CFlowPairVec> CFlow;

	const vector<DVector> & getColourFactors(const Particle & inpart,
	const ParticleVector & decay,
	unsigned int & nflow);

	const CFlow & colourFlows(const Particle & inpart,
	const ParticleVector & decay);

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

	//@}

	private:

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

	private:

	/**
	* Store the incoming particle
	*/
	PDPtr incoming_;

	/**
	* Outgoing particles
	*/
	vector<PDPtr> outgoing_;

	/**
	* Maximum weight for integration
	*/
	double maxWeight_;

	/**
	* Store colour factors for ME calc.
	*/
	vector<DVector> colour_;

	};


	/**
	* Write a map with ShowerInteraction as the key
	*/
	template<typename T, typename Cmp, typename A>
	inline PersistentOStream & operator<<(PersistentOStream & os,
	const map<ShowerInteraction,T,Cmp,A> & m) {
	os << m.size();
	if(m.find(ShowerInteraction::QCD)!=m.end()) {
	os << 0 << m.at(ShowerInteraction::QCD);
	}
	if(m.find(ShowerInteraction::QED)!=m.end()) {
	os << 1 << m.at(ShowerInteraction::QED);
	}
	return os;
	}

	/**
	* Read a map with ShowerInteraction as the key
	*/
	template <typename T, typename Cmp, typename A>
	inline PersistentIStream & operator>>(PersistentIStream & is, map<ShowerInteraction,T,Cmp,A> & m) {
	m.clear();
	long size;
	int k;
	is >> size;
	while ( size-- && is ) {
	is >> k;
	if(k==0)
	is >> m[ShowerInteraction::QCD];
	else if(k==1)
	is >> m[ShowerInteraction::QED];
	else
	assert(false);
	}
	return is;
	}
	}

	#endif /* HERWIG_GeneralTwoBodyDecayer2_H */
	diff --git a/Decay/General/SVVDecayer.cc b/Decay/General/SVVDecayer.cc
	--- a/Decay/General/SVVDecayer.cc
	+++ b/Decay/General/SVVDecayer.cc
	@@ -1,434 +1,440 @@
	// -- C++ --
	//
	// SVVDecayer.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 SVVDecayer class.
	//

	#include "SVVDecayer.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/Vertex/Scalar/VVSVertex.h"
	#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
	#include "Herwig/Utilities/Kinematics.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"

	using namespace Herwig;
	using namespace ThePEG::Helicity;

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

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

	void SVVDecayer::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<AbstractVVSVertexPtr>(vert));
	perturbativeVertex_.push_back(dynamic_ptr_cast<VVSVertexPtr> (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<AbstractVVVVertexPtr>(outV[0].at(inter));
	outgoingVertex2_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[1].at(inter));
	}
	}

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

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

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

	void SVVDecayer::Init() {

	static ClassDocumentation<SVVDecayer> documentation
	("This implements the decay of a scalar to 2 vector bosons.");

	}

	-double SVVDecayer::me2(const int , const Particle & inpart,
	- const ParticleVector& decay,
	+
	+void SVVDecayer::
	+constructSpinInfo(const Particle & part, ParticleVector decay) const {
	+ bool photon[2];
	+ for(unsigned int ix=0;ix<2;++ix)
	+ photon[ix] = decay[ix]->mass()==ZERO;
	+ ScalarWaveFunction::
	+ constructSpinInfo(const_ptr_cast<tPPtr>(&part),incoming,true);
	+ for(unsigned int ix=0;ix<2;++ix)
	+ VectorWaveFunction::
	+ constructSpinInfo(vectors_[ix],decay[ix],outgoing,true,photon[ix]);
	+}
	+
	+double SVVDecayer::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::Spin1,PDT::Spin1)));
	bool photon[2];
	for(unsigned int ix=0;ix<2;++ix)
	- photon[ix] = decay[ix]->mass()==ZERO;
	+ photon[ix] = outgoing[ix]->mass()==ZERO;
	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);
	- for(unsigned int ix=0;ix<2;++ix)
	- VectorWaveFunction::
	- constructSpinInfo(vectors_[ix],decay[ix],outgoing,true,photon[ix]);
	- }
	for(unsigned int ix=0;ix<2;++ix)
	VectorWaveFunction::
	- calculateWaveFunctions(vectors_[ix],decay[ix],outgoing,photon[ix]);
	+ calculateWaveFunctions(vectors_[ix],momenta[ix],outgoing[ix],Helicity::outgoing,photon[ix]);


	- Energy2 scale(sqr(inpart.mass()));
	+ Energy2 scale(sqr(part.mass()));
	unsigned int iv1,iv2;
	for(iv2 = 0; iv2 < 3; ++iv2) {
	if( photon[1] && iv2 == 1 ) ++iv2;
	for(iv1=0;iv1<3;++iv1) {
	if( photon[0] && iv1 == 1) ++iv1;
	(*ME())(0, iv1, iv2) = 0.;
	for(auto vert : vertex_)
	(*ME())(0, iv1, iv2) += vert->evaluate(scale,vectors_[0][iv1],
	vectors_[1][iv2],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 SVVDecayer::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, outa.first ,
	outb.first, in);
	double mu1sq = sqr(outa.second/inpart.second);
	double mu2sq = sqr(outb.second/inpart.second);
	double m1pm2 = mu1sq + mu2sq;
	double me2(0.);
	if( mu1sq > 0. && mu2sq > 0.)
	me2 = ( m1pm2(m1pm2 - 2.) + 8.mu1sq*mu2sq + 1.)/4./mu1sq/mu2sq;
	else if( mu1sq == 0. \|\| mu2sq == 0. )
	me2 = 3.;
	else
	me2 = 4.;

	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second,outa.second,
	outb.second);
	Energy output = norm(perturbativeVertex_[0]->norm())*
	me2pcm/(8Constants::pi)/scale*UnitRemoval::E2;
	// 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 SVVDecayer::threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt) {
	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);
	}
	if(meopt==Terminate) {
	ScalarWaveFunction::
	constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),incoming,true);
	VectorWaveFunction::
	constructSpinInfo(vectors3_[0],decay[0],outgoing,true,false);
	VectorWaveFunction::
	constructSpinInfo(vectors3_[1],decay[1],outgoing,true,false);
	VectorWaveFunction::
	constructSpinInfo(gluon_ ,decay[2],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::Spin1,
	PDT::Spin1, PDT::Spin1)));
	bool massless[2];
	for(unsigned int ix=0;ix<2;++ix)
	massless[ix] = decay[ix]->mass()!=ZERO;
	// create wavefunctions
	VectorWaveFunction::calculateWaveFunctions(vectors3_[0],decay[0],outgoing,massless[0]);
	VectorWaveFunction::calculateWaveFunctions(vectors3_[1],decay[1],outgoing,massless[1]);
	VectorWaveFunction::calculateWaveFunctions(gluon_ ,decay[2],outgoing,true);

	// gauge 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[2]->momentum(),
	decay[2]->dataPtr(),10,
	outgoing));
	}
	}
	#endif

	// get the outgoing vertices
	AbstractVVVVertexPtr outgoingVertex1;
	AbstractVVVVertexPtr outgoingVertex2;
	identifyVertices(inpart,decay, outgoingVertex1, outgoingVertex2,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 iv1 = 0; iv1 < 3; ++iv1) {
	if(massless[0] && iv1==1) continue;
	for(unsigned int iv2 = 0; iv2 < 3; ++iv2) {
	if(massless[1] && iv2==1) continue;
	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]);
	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,vectors3_[0][iv1],
	vectors3_[1][iv2],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, iv1, iv2, ig) +=
	colourFlow[0][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}
	// radiation from the 1st outgoing vector
	if((decay[0]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[0]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertex1);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[0]->dataPtr();
	if(off->CC()) off = off->CC();
	VectorWaveFunction vectorInter =
	outgoingVertex1->evaluate(scale,3,off,gluon_[2*ig],vectors3_[0][iv1],decay[0]->mass());

	assert(vectors3_[0][iv1].particle()->id()==vectorInter.particle()->id());

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

	if((decay[1]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[1]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertex2);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[1]->dataPtr();
	if(off->CC()) off = off->CC();
	VectorWaveFunction vectorInter =
	outgoingVertex2->evaluate(scale,3,off, gluon_[2*ig],vectors3_[1][iv2],decay[1]->mass());

	assert(vectors3_[1][iv2].particle()->id()==vectorInter.particle()->id());

	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,vectors3_[0][iv1],vectorInter,swave3_);
	if(!couplingSet) {
	gs = abs(outgoingVertex2->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[2].size();++ix) {
	(*ME[colourFlow[2][ix].first])(0, iv1, iv2, 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 SVVDecayer::identifyVertices(const Particle & inpart, const ParticleVector & decay,
	AbstractVVVVertexPtr & outgoingVertex1,
	AbstractVVVVertexPtr & outgoingVertex2,
	ShowerInteraction inter) {
	if(inter==ShowerInteraction::QCD) {
	// work out which scalar each outgoing vertex corresponds to
	// two outgoing vertices
	if( inpart.dataPtr() ->iColour()==PDT::Colour0 &&
	((decay[0]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour3bar) \|\|
	(decay[0]->dataPtr()->iColour()==PDT::Colour8 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour8))){
	if(outgoingVertex1_[inter]==outgoingVertex2_[inter]){
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	else if (outgoingVertex1_[inter]->isIncoming(getParticleData(decay[0]->id()))){
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	else if (outgoingVertex2_[inter]->isIncoming(getParticleData(decay[0]->id()))){
	outgoingVertex1 = outgoingVertex2_[inter];
	outgoingVertex2 = outgoingVertex1_[inter];
	}
	}
	else if(inpart.dataPtr() ->iColour()==PDT::Colour8 &&
	decay[0]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour3bar){
	if(outgoingVertex1_[inter]==outgoingVertex2_[inter]){
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	else if (outgoingVertex1_[inter]->isIncoming(getParticleData(decay[0]->id()))){
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	else if (outgoingVertex2_[inter]->isIncoming(getParticleData(decay[0]->id()))){
	outgoingVertex1 = outgoingVertex2_[inter];
	outgoingVertex2 = outgoingVertex1_[inter];
	}
	}

	// one outgoing vertex
	else if(inpart.dataPtr()->iColour()==PDT::Colour3){
	if(decay[0]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour0){
	if (outgoingVertex1_[inter]) outgoingVertex1 = outgoingVertex1_[inter];
	else if(outgoingVertex2_[inter]) outgoingVertex1 = outgoingVertex2_[inter];
	}
	else if (decay[0]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour8){
	if (outgoingVertex1_[inter]->isIncoming(getParticleData(decay[1]->dataPtr()->id()))){
	outgoingVertex1 = outgoingVertex2_[inter];
	outgoingVertex2 = outgoingVertex1_[inter];
	}
	else {
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	}
	}
	else if(inpart.dataPtr()->iColour()==PDT::Colour3bar){
	if(decay[1]->dataPtr()->iColour()==PDT::Colour3bar &&
	decay[0]->dataPtr()->iColour()==PDT::Colour0){
	if (outgoingVertex1_[inter]) outgoingVertex2 = outgoingVertex1_[inter];
	else if(outgoingVertex2_[inter]) outgoingVertex2 = outgoingVertex2_[inter];
	}
	else if (decay[0]->dataPtr()->iColour()==PDT::Colour8 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour3bar){
	if (outgoingVertex1_[inter]->isIncoming(getParticleData(decay[0]->dataPtr()->id()))){
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	else {
	outgoingVertex1 = outgoingVertex2_[inter];
	outgoingVertex2 = outgoingVertex1_[inter];
	}
	}
	}

	if (! ((incomingVertex_[inter] && (outgoingVertex1 \|\| outgoingVertex2)) \|\|
	( outgoingVertex1 && outgoingVertex2)))
	throw Exception()
	<< "Invalid vertices for QCD radiation in SVV decay in SVVDecayer::identifyVertices"
	<< Exception::runerror;
	}
	else {
	if(decay[0]->dataPtr()->charged()) {
	if (outgoingVertex1_[inter] &&
	outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[0]->dataPtr())))
	outgoingVertex1 = outgoingVertex1_[inter];
	else
	outgoingVertex1 = outgoingVertex2_[inter];
	}
	if(decay[1]->dataPtr()->charged()) {
	if (outgoingVertex1_[inter] &&
	outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[1]->dataPtr())))
	outgoingVertex2 = outgoingVertex1_[inter];
	else
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	}
	}
	diff --git a/Decay/General/SVVDecayer.h b/Decay/General/SVVDecayer.h
	--- a/Decay/General/SVVDecayer.h
	+++ b/Decay/General/SVVDecayer.h
	@@ -1,230 +1,239 @@
	// -- C++ --
	//
	// SVVDecayer.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_SVVDecayer_H
	#define HERWIG_SVVDecayer_H
	//
	// This is the declaration of the SVVDecayer class.
	//

	-#include "GeneralTwoBodyDecayer.h"
	+#include "GeneralTwoBodyDecayer2.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Helicity/Vertex/Scalar/VVSVertex.h"
	#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"

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

	/** \ingroup Decay
	* This SVVDecayer class implements the decay of a scalar to
	* 2 vector bosons using either the tree level VVSVertex or the loop vertex.
	* It inherits from
	- * GeneralTwoBodyDecayer and implements the virtual member functions me2()
	+ * GeneralTwoBodyDecayer2 and implements the virtual member functions me2()
	* and partialWidth(). It also stores a pointer to the VVSVertex.
	*
	- * @see GeneralTwoBodyDecayer
	+ * @see GeneralTwoBodyDecayer2
	*
	*/
	-class SVVDecayer: public GeneralTwoBodyDecayer {
	+class SVVDecayer: public GeneralTwoBodyDecayer2 {

	public:

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

	/** @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;

	/**
	* 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>);

	/**
	* 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);
	//@}

	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;
	//@}

	protected:

	/**
	* Find the vertices for the decay
	*/
	void identifyVertices(const Particle & inpart, const ParticleVector & decay,
	AbstractVVVVertexPtr & outgoingVertex1,
	AbstractVVVVertexPtr & outgoingVertex2,
	ShowerInteraction inter);

	private:

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

	private:

	/**
	* Abstract pointer to general VVS vertex
	*/
	vector<AbstractVVSVertexPtr> vertex_;

	/**
	* Pointer to the perturbative form
	*/
	vector<VVSVertexPtr> perturbativeVertex_;

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

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from the 1st outgoing vector
	*/
	map<ShowerInteraction,AbstractVVVVertexPtr> outgoingVertex1_;

	/**
	* Abstract pointer to AbstractFFVVertex for QCD radiation from the 2nd outgoing vector
	*/
	map<ShowerInteraction,AbstractVVVVertexPtr> outgoingVertex2_;

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

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

	/**
	* Vector wavefunctions
	*/
	mutable vector<Helicity::VectorWaveFunction> vectors_[2];

	private:

	/**
	* Member for the POWHEG correction
	*/
	//@{
	/**
	* Spin density matrix for 3 body decay
	*/
	mutable RhoDMatrix rho3_;

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

	/**
	* Vector wavefunctions
	*/
	mutable vector<Helicity::VectorWaveFunction> vectors3_[2];

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

	}

	#endif /* HERWIG_SVVDecayer_H */
	diff --git a/Decay/General/TVVDecayer.cc b/Decay/General/TVVDecayer.cc
	--- a/Decay/General/TVVDecayer.cc
	+++ b/Decay/General/TVVDecayer.cc
	@@ -1,340 +1,344 @@
	// -- C++ --
	//
	// TVVDecayer.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 TVVDecayer class.
	//

	#include "TVVDecayer.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/TensorWaveFunction.h"
	#include "Herwig/Utilities/Kinematics.h"
	#include "ThePEG/Helicity/LorentzTensor.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"
	using namespace Herwig;
	using namespace ThePEG::Helicity;

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

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

	void TVVDecayer::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<AbstractVVTVertexPtr>(vert));
	perturbativeVertex_.push_back(dynamic_ptr_cast<VVTVertexPtr> (vert));
	}
	vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
	for(auto & inter : itemp) {
	fourPointVertex_[inter] = dynamic_ptr_cast<AbstractVVVTVertexPtr>(fourV.at(inter));
	outgoingVertex1_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr> (outV[0].at(inter));
	outgoingVertex2_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr> (outV[1].at(inter));
	}
	}

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

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

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

	void TVVDecayer::Init() {

	static ClassDocumentation<TVVDecayer> documentation
	("This class implements the decay of a tensor to 2 vector bosons");

	}

	-double TVVDecayer::me2(const int , const Particle & inpart,
	- const ParticleVector & decay,
	+void TVVDecayer::
	+constructSpinInfo(const Particle & part, ParticleVector decay) const {
	+ bool photon[2];
	+ for(unsigned int ix=0;ix<2;++ix)
	+ photon[ix] = decay[ix]->mass()==ZERO;
	+ TensorWaveFunction::
	+ constructSpinInfo(tensors_,const_ptr_cast<tPPtr>(&part),
	+ incoming,true,false);
	+ for(unsigned int ix=0;ix<2;++ix)
	+ VectorWaveFunction::
	+ constructSpinInfo(vectors_[ix],decay[ix],outgoing,true,photon[ix]);
	+}
	+
	+double TVVDecayer::me2(const int,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	MEOption meopt) const {
	if(!ME())
	ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin2,PDT::Spin1,PDT::Spin1)));
	bool photon[2];
	for(unsigned int ix=0;ix<2;++ix)
	- photon[ix] = decay[ix]->mass()==ZERO;
	+ photon[ix] = outgoing[ix]->mass()==ZERO;
	if(meopt==Initialize) {
	TensorWaveFunction::
	- calculateWaveFunctions(tensors_,rho_,const_ptr_cast<tPPtr>(&inpart),
	- incoming,false);
	+ 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);
	- for(unsigned int ix=0;ix<2;++ix)
	- VectorWaveFunction::
	- constructSpinInfo(vectors_[ix],decay[ix],outgoing,true,photon[ix]);
	- return 0.;
	- }
	for(unsigned int ix=0;ix<2;++ix)
	VectorWaveFunction::
	- calculateWaveFunctions(vectors_[ix],decay[ix],outgoing,photon[ix]);
	- Energy2 scale(sqr(inpart.mass()));
	+ calculateWaveFunctions(vectors_[ix],momenta[ix],outgoing[ix],Helicity::outgoing,photon[ix]);
	+ Energy2 scale(sqr(part.mass()));
	unsigned int thel,v1hel,v2hel;
	for(thel=0;thel<5;++thel) {
	for(v1hel=0;v1hel<3;++v1hel) {
	for(v2hel=0;v2hel<3;++v2hel) {
	- (*ME())(thel,v1hel,v2hel) = 0.;
	- for(auto vert : vertex_)
	- (*ME())(thel,v1hel,v2hel) += vert->evaluate(scale,
	- vectors_[0][v1hel],
	- vectors_[1][v2hel],
	- tensors_[thel]);
	- if(photon[1]) ++v2hel;
	+ (*ME())(thel,v1hel,v2hel) = 0.;
	+ for(auto vert : vertex_)
	+ (*ME())(thel,v1hel,v2hel) += vert->evaluate(scale,
	+ vectors_[0][v1hel],
	+ vectors_[1][v2hel],
	+ tensors_[thel]);
	+ if(photon[1]) ++v2hel;
	}
	if(photon[0]) ++v1hel;
	}
	}
	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 TVVDecayer::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, outa.first, outb.first, in);
	double mu2 = sqr(outa.second/inpart.second);
	double b = sqrt(1 - 4.*mu2);
	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second,outa.second,
	outb.second);
	Energy2 me2;
	if(outa.second > ZERO && outb.second > ZERO)
	me2 = scale(30 - 20.bb + 3.pow(b,4))/120.;
	else
	me2 = scale/10.;

	Energy output = norm(perturbativeVertex_[0]->norm())me2pcm
	/(8.Constants::pi)UnitRemoval::InvE2;
	// 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 TVVDecayer::threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt) {
	bool massless[2];
	for(unsigned int ix=0;ix<2;++ix)
	massless[ix] = decay[ix]->mass()==ZERO;
	int iglu(2);
	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);
	for(unsigned int ix=0;ix<2;++ix)
	VectorWaveFunction::
	constructSpinInfo(vectors3_[ix],decay[ix ],outgoing,true, massless[ix]);
	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::Spin1,
	PDT::Spin1, PDT::Spin1)));
	// create wavefunctions
	for(unsigned int ix=0;ix<2;++ix)
	VectorWaveFunction::
	calculateWaveFunctions(vectors3_[ix],decay[ix ],outgoing,massless[ix]);
	VectorWaveFunction::
	calculateWaveFunctions(gluon_ ,decay[iglu ],outgoing,true);

	// gauge 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

	// work out which vector each outgoing vertex corresponds to
	if(outgoingVertex1_[inter]!=outgoingVertex2_[inter] &&
	outgoingVertex1_[inter]->isIncoming(getParticleData(decay[1]->id())))
	swap(outgoingVertex1_[inter], outgoingVertex2_[inter]);

	if (! (outgoingVertex1_[inter] && outgoingVertex2_[inter]))
	throw Exception()
	<< "Invalid vertices for radiation in TVV decay in TVVDecayer::threeBodyME"
	<< Exception::runerror;

	if( !(!inpart.dataPtr()->coloured() && inter ==ShowerInteraction::QCD) &&
	!(!inpart.dataPtr()->charged() && inter ==ShowerInteraction::QED))
	throw Exception()
	<< "Invalid vertices for radiation in TVV decay in TVVDecayer::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 iv0 = 0; iv0 < 3; ++iv0) {
	for(unsigned int iv1 = 0; iv1 < 3; ++iv1) {
	for(unsigned int ig = 0; ig < 2; ++ig) {

	// radiation from first outgoing vector
	if((decay[0]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[0]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertex1_[inter]);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[0]->dataPtr();
	if(off->CC()) off = off->CC();
	VectorWaveFunction vectInter =
	outgoingVertex1_[inter]->evaluate(scale,3,off,gluon_[2*ig],
	vectors3_[0][iv0],decay[0]->mass());

	assert(vectors3_[0][iv0].particle()->PDGName()==vectInter.particle()->PDGName());

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

	// radiation from second outgoing vector
	if((decay[1]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[1]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertex2_[inter]);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[1]->dataPtr();
	if(off->CC()) off = off->CC();
	VectorWaveFunction vectInter =
	outgoingVertex2_[inter]->evaluate(scale,3,off,vectors3_[1][iv1],
	gluon_[2*ig],decay[1]->mass());

	assert(vectors3_[1][iv1].particle()->PDGName()==vectInter.particle()->PDGName());

	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,vectInter,vectors3_[0][iv0],
	tensors3_[it]);
	if(!couplingSet) {
	gs = abs(outgoingVertex2_[inter]->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[2].size();++ix) {
	(*ME[colourFlow[2][ix].first])(it, iv0, iv1, 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, vectors3_[0][iv0],
	vectors3_[1][iv1],gluon_[2*ig],
	tensors3_[it]);
	for(unsigned int ix=0;ix<colourFlow[3].size();++ix) {
	(*ME[colourFlow[3][ix].first])(it, iv0, iv1, ig) +=
	colourFlow[3][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}
	}
	if(massless[1]) ++iv1;
	}
	if(massless[0]) ++iv0;
	}
	}

	// 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/TVVDecayer.h b/Decay/General/TVVDecayer.h
	--- a/Decay/General/TVVDecayer.h
	+++ b/Decay/General/TVVDecayer.h
	@@ -1,213 +1,222 @@
	// -- C++ --
	//
	// TVVDecayer.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_TVVDecayer_H
	#define HERWIG_TVVDecayer_H
	//
	// This is the declaration of the TVVDecayer class.
	//

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

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

	/** \ingroup Decay
	* The TVVDecayer class implements the decay of a tensor
	* to 2 vector bosons in a general model. It holds a VVTVertex 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 TVVDecayer: public GeneralTwoBodyDecayer {
	+class TVVDecayer: public GeneralTwoBodyDecayer2 {

	public:

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

	/** @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.
	*/
	TVVDecayer & operator=(const TVVDecayer &);

	private:

	/**
	* Abstract pointer to AbstractVVTVertex
	*/
	vector<AbstractVVTVertexPtr> vertex_;

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

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

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

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

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

	/**
	* Polarization tensors of decaying particle
	*/
	mutable vector<Helicity::TensorWaveFunction> tensors_;

	/**
	* Polarization vectors of outgoing vector bosons
	*/
	mutable vector<Helicity::VectorWaveFunction> vectors_[2];

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

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

	/**
	* Polarization vectors of outgoing vector bosons
	*/
	mutable vector<Helicity::VectorWaveFunction> vectors3_[2];

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

	};

	}

	#endif /* HERWIG_TVVDecayer_H */
	diff --git a/Decay/General/VVVDecayer.cc b/Decay/General/VVVDecayer.cc
	--- a/Decay/General/VVVDecayer.cc
	+++ b/Decay/General/VVVDecayer.cc
	@@ -1,443 +1,448 @@
	// -- C++ --
	//
	// VVVDecayer.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 VVVDecayer class.
	//

	#include "VVVDecayer.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 "Herwig/Utilities/Kinematics.h"
	#include "Herwig/Decay/GeneralDecayMatrixElement.h"

	using namespace Herwig;
	using namespace ThePEG::Helicity;

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

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

	void VVVDecayer::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<AbstractVVVVertexPtr>(vert));
	perturbativeVertex_ .push_back(dynamic_ptr_cast<VVVVertexPtr> (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<AbstractVVVVertexPtr>(outV[0].at(inter));
	outgoingVertex2_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[1].at(inter));
	fourPointVertex_[inter] = dynamic_ptr_cast<AbstractVVVVVertexPtr>(fourV.at(inter));
	}
	}

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

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

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

	void VVVDecayer::Init() {

	static ClassDocumentation<VVVDecayer> documentation
	("The VVVDecayer class implements the decay of a vector boson "
	"into 2 vector bosons");

	}

	-double VVVDecayer::me2(const int , const Particle & inpart,
	- const ParticleVector & decay,
	+void VVVDecayer::
	+constructSpinInfo(const Particle & part, ParticleVector decay) const {
	+ bool massless[2];
	+ for(unsigned int ix=0;ix<2;++ix)
	+ massless[ix] = (decay[ix]->id()==ParticleID::gamma \|\|
	+ decay[ix]->id()==ParticleID::g);
	+ VectorWaveFunction::constructSpinInfo(vectors_[0],const_ptr_cast<tPPtr>(&part),
	+ incoming,true,false);
	+ for(unsigned int ix=0;ix<2;++ix)
	+ VectorWaveFunction::
	+ constructSpinInfo(vectors_[ix+1],decay[ix],outgoing,true,massless[ix]);
	+}
	+
	+double VVVDecayer::me2(const int,const Particle & part,
	+ const tPDVector & outgoing,
	+ const vector<Lorentz5Momentum> & momenta,
	MEOption meopt) const {
	if(!ME())
	ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,PDT::Spin1,PDT::Spin1)));
	bool massless[2];
	for(unsigned int ix=0;ix<2;++ix)
	- massless[ix] = (decay[ix]->id()==ParticleID::gamma \|\|
	- decay[ix]->id()==ParticleID::g);
	+ massless[ix] = (outgoing[ix]->id()==ParticleID::gamma \|\|
	+ outgoing[ix]->id()==ParticleID::g);
	if(meopt==Initialize) {
	VectorWaveFunction::calculateWaveFunctions(vectors_[0],rho_,
	- const_ptr_cast<tPPtr>(&inpart),
	+ const_ptr_cast<tPPtr>(&part),
	incoming,false);
	// fix rho if no correlations
	fixRho(rho_);
	}
	- if(meopt==Terminate) {
	- VectorWaveFunction::constructSpinInfo(vectors_[0],const_ptr_cast<tPPtr>(&inpart),
	- incoming,true,false);
	- for(unsigned int ix=0;ix<2;++ix)
	- VectorWaveFunction::
	- constructSpinInfo(vectors_[ix+1],decay[ix],outgoing,true,massless[ix]);
	- return 0.;
	- }
	for(unsigned int ix=0;ix<2;++ix)
	VectorWaveFunction::
	- calculateWaveFunctions(vectors_[ix+1],decay[ix],outgoing,massless[ix]);
	- Energy2 scale(sqr(inpart.mass()));
	+ calculateWaveFunctions(vectors_[ix],momenta[ix],outgoing[ix],Helicity::outgoing,massless[ix]);
	+ Energy2 scale(sqr(part.mass()));
	for(unsigned int iv3=0;iv3<3;++iv3) {
	for(unsigned int iv2=0;iv2<3;++iv2) {
	for(unsigned int iv1=0;iv1<3;++iv1) {
	(*ME())(iv1,iv2,iv3) = 0.;
	for(auto vert : vertex_) {
	(*ME())(iv1,iv2,iv3) += vert->
	evaluate(scale,vectors_[1][iv2],vectors_[2][iv3],vectors_[0][iv1]);
	}
	}
	}
	}
	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 VVVDecayer::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), in,
	outa.first, outb.first);
	double mu1(outa.second/inpart.second), mu1sq(sqr(mu1)),
	mu2(outb.second/inpart.second), mu2sq(sqr(mu2));
	double vn = norm(perturbativeVertex_[0]->norm());
	if(vn == ZERO \|\| mu1sq == ZERO \|\| mu2sq == ZERO) return ZERO;
	double me2 =
	(mu1 - mu2 - 1.)(mu1 - mu2 + 1.)(mu1 + mu2 - 1.)*(mu1 + mu2 + 1.)
	* (sqr(mu1sq) + sqr(mu2sq) + 10.(mu1sqmu2sq + mu1sq + mu2sq) + 1.)
	/4./mu1sq/mu2sq;
	Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second,outa.second,
	outb.second);
	Energy pWidth = vnme2pcm/24./Constants::pi;
	// colour factor
	pWidth *= colourFactor(inpart.first,outa.first,outb.first);
	// return the answer
	return pWidth;
	}
	else {
	- return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
	+ return GeneralTwoBodyDecayer2::partialWidth(inpart,outa,outb);
	}
	}

	double VVVDecayer::threeBodyME(const int , const Particle & inpart,
	const ParticleVector & decay,
	ShowerInteraction inter, MEOption meopt) {
	if(meopt==Initialize) {
	// create vector wavefunction for decaying particle
	VectorWaveFunction::calculateWaveFunctions(vector3_, rho3_, const_ptr_cast<tPPtr>(&inpart),
	incoming, false);
	}
	if(meopt==Terminate) {
	VectorWaveFunction::
	constructSpinInfo(vector3_ ,const_ptr_cast<tPPtr>(&inpart),outgoing,true,false);
	VectorWaveFunction::
	constructSpinInfo(vectors3_[0],decay[0],outgoing,true,false);
	VectorWaveFunction::
	constructSpinInfo(vectors3_[1],decay[1],outgoing,true,false);
	VectorWaveFunction::
	constructSpinInfo(gluon_ ,decay[2],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::Spin1,
	PDT::Spin1, PDT::Spin1)));
	bool massless[2];
	for(unsigned int ix=0;ix<2;++ix)
	massless[ix] = decay[ix]->mass()!=ZERO;
	// create wavefunctions
	VectorWaveFunction::calculateWaveFunctions(vectors3_[0],decay[0],outgoing,massless[0]);
	VectorWaveFunction::calculateWaveFunctions(vectors3_[1],decay[1],outgoing,massless[1]);
	VectorWaveFunction::calculateWaveFunctions(gluon_ ,decay[2],outgoing,true);

	// gauge 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[2]->momentum(),
	decay[2]->dataPtr(),10,
	outgoing));
	}
	}
	#endif

	// get the outgoing vertices
	AbstractVVVVertexPtr outgoingVertex1;
	AbstractVVVVertexPtr outgoingVertex2;
	identifyVertices(inpart,decay, outgoingVertex1, outgoingVertex2,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 iv0 = 0; iv0 < 3; ++iv0) {
	for(unsigned int iv1 = 0; iv1 < 3; ++iv1) {
	if(massless[0] && iv1==1) continue;
	for(unsigned int iv2 = 0; iv2 < 3; ++iv2) {
	if(massless[1] && iv2==1) continue;
	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_[iv0],
	gluon_[2*ig],inpart.mass());

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

	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,vectorInter,vectors3_[0][iv1],
	vectors3_[1][iv2]);
	if(!couplingSet) {
	gs = abs(incomingVertex_[inter]->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
	(*ME[colourFlow[0][ix].first])(iv0, iv1, iv2, ig) +=
	colourFlow[0][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}
	// radiation from the 1st outgoing vector
	if((decay[0]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[0]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertex1);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[0]->dataPtr();
	if(off->CC()) off = off->CC();
	VectorWaveFunction vectorInter =
	outgoingVertex1->evaluate(scale,3,off,gluon_[2*ig],vectors3_[0][iv1],decay[0]->mass());

	assert(vectors3_[0][iv1].particle()->id()==vectorInter.particle()->id());

	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,vector3_[iv0],vectorInter,vectors3_[1][iv2]);
	if(!couplingSet) {
	gs = abs(outgoingVertex1->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[1].size();++ix) {
	(*ME[colourFlow[1][ix].first])(iv0, iv1, iv2, ig) +=
	colourFlow[1][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}
	// radiation from the second outgoing vector
	if((decay[1]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) \|\|
	(decay[1]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
	assert(outgoingVertex2);
	// ensure you get correct outgoing particle from first vertex
	tcPDPtr off = decay[1]->dataPtr();
	if(off->CC()) off = off->CC();
	VectorWaveFunction vectorInter =
	outgoingVertex2->evaluate(scale,3,off, gluon_[2*ig],vectors3_[1][iv2],decay[1]->mass());

	assert(vectors3_[1][iv2].particle()->id()==vectorInter.particle()->id());

	Complex diag = 0.;
	for(auto vertex : vertex_)
	diag += vertex->evaluate(scale,vector3_[iv0],vectors3_[0][iv1],vectorInter);
	if(!couplingSet) {
	gs = abs(outgoingVertex2->norm());
	couplingSet = true;
	}
	for(unsigned int ix=0;ix<colourFlow[2].size();++ix) {
	(*ME[colourFlow[2][ix].first])(iv0, iv1, iv2, ig) +=
	colourFlow[2][ix].second*diag;
	}
	#ifdef GAUGE_CHECK
	total+=norm(diag);
	#endif
	}
	// 4-point vertex
	if (fourPointVertex_[inter]) {
	Complex diag = fourPointVertex_[inter]->evaluate(scale,0, vector3_[iv0],
	vectors3_[0][iv1],
	vectors3_[1][iv2],gluon_[2*ig]);
	for(unsigned int ix=0;ix<colourFlow[2].size();++ix) {
	(*ME[colourFlow[3][ix].first])(iv0, iv1, iv2, 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;
	}

	void VVVDecayer::identifyVertices(const Particle & inpart, const ParticleVector & decay,
	AbstractVVVVertexPtr & outgoingVertex1,
	AbstractVVVVertexPtr & outgoingVertex2,
	ShowerInteraction inter) {
	if(inter==ShowerInteraction::QCD) {
	// work out which scalar each outgoing vertex corresponds to
	// two outgoing vertices
	if( inpart.dataPtr() ->iColour()==PDT::Colour0 &&
	((decay[0]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour3bar) \|\|
	(decay[0]->dataPtr()->iColour()==PDT::Colour8 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour8))){
	if(outgoingVertex1_[inter]==outgoingVertex2_[inter]){
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	else if (outgoingVertex1_[inter]->isIncoming(getParticleData(decay[0]->id()))){
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	else if (outgoingVertex2_[inter]->isIncoming(getParticleData(decay[0]->id()))){
	outgoingVertex1 = outgoingVertex2_[inter];
	outgoingVertex2 = outgoingVertex1_[inter];
	}
	}
	else if(inpart.dataPtr() ->iColour()==PDT::Colour8 &&
	decay[0]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour3bar){
	if(outgoingVertex1_[inter]==outgoingVertex2_[inter]){
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	else if (outgoingVertex1_[inter]->isIncoming(getParticleData(decay[0]->id()))){
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	else if (outgoingVertex2_[inter]->isIncoming(getParticleData(decay[0]->id()))){
	outgoingVertex1 = outgoingVertex2_[inter];
	outgoingVertex2 = outgoingVertex1_[inter];
	}
	}

	// one outgoing vertex
	else if(inpart.dataPtr()->iColour()==PDT::Colour3){
	if(decay[0]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour0){
	if (outgoingVertex1_[inter]) outgoingVertex1 = outgoingVertex1_[inter];
	else if(outgoingVertex2_[inter]) outgoingVertex1 = outgoingVertex2_[inter];
	}
	else if (decay[0]->dataPtr()->iColour()==PDT::Colour3 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour8){
	if (outgoingVertex1_[inter]->isIncoming(getParticleData(decay[1]->dataPtr()->id()))){
	outgoingVertex1 = outgoingVertex2_[inter];
	outgoingVertex2 = outgoingVertex1_[inter];
	}
	else {
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	}
	}
	else if(inpart.dataPtr()->iColour()==PDT::Colour3bar){
	if(decay[1]->dataPtr()->iColour()==PDT::Colour3bar &&
	decay[0]->dataPtr()->iColour()==PDT::Colour0){
	if (outgoingVertex1_[inter]) outgoingVertex2 = outgoingVertex1_[inter];
	else if(outgoingVertex2_[inter]) outgoingVertex2 = outgoingVertex2_[inter];
	}
	else if (decay[0]->dataPtr()->iColour()==PDT::Colour8 &&
	decay[1]->dataPtr()->iColour()==PDT::Colour3bar){
	if (outgoingVertex1_[inter]->isIncoming(getParticleData(decay[0]->dataPtr()->id()))){
	outgoingVertex1 = outgoingVertex1_[inter];
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	else {
	outgoingVertex1 = outgoingVertex2_[inter];
	outgoingVertex2 = outgoingVertex1_[inter];
	}
	}
	}

	if (! ((incomingVertex_[inter] && (outgoingVertex1 \|\| outgoingVertex2)) \|\|
	( outgoingVertex1 && outgoingVertex2)))
	throw Exception()
	<< "Invalid vertices for QCD radiation in VVV decay in VVVDecayer::identifyVertices"
	<< Exception::runerror;
	}
	else {
	if(decay[0]->dataPtr()->charged()) {
	if (outgoingVertex1_[inter] &&
	outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[0]->dataPtr())))
	outgoingVertex1 = outgoingVertex1_[inter];
	else
	outgoingVertex1 = outgoingVertex2_[inter];
	}
	if(decay[1]->dataPtr()->charged()) {
	if (outgoingVertex1_[inter] &&
	outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[1]->dataPtr())))
	outgoingVertex2 = outgoingVertex1_[inter];
	else
	outgoingVertex2 = outgoingVertex2_[inter];
	}
	}
	}
	diff --git a/Decay/General/VVVDecayer.h b/Decay/General/VVVDecayer.h
	--- a/Decay/General/VVVDecayer.h
	+++ b/Decay/General/VVVDecayer.h
	@@ -1,228 +1,237 @@
	// -- C++ --
	//
	// VVVDecayer.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_VVVDecayer_H
	#define HERWIG_VVVDecayer_H
	//
	// This is the declaration of the VVVDecayer class.
	//

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

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

	/** \ingroup Decay
	* The VVVDecayer class implements the decay of a vector
	* to 2 vectors in a general model. It holds an VVVVertex 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 VVVDecayer: public GeneralTwoBodyDecayer {
	+class VVVDecayer: public GeneralTwoBodyDecayer2 {

	public:

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

	/** @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;

	/**
	* 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>);

	/**
	* 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);
	//@}

	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;
	//@}

	protected:

	/**
	* Find the vertices for the decay
	*/
	void identifyVertices(const Particle & inpart, const ParticleVector & decay,
	AbstractVVVVertexPtr & outgoingVertex1,
	AbstractVVVVertexPtr & outgoingVertex2,
	ShowerInteraction inter);

	private:

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

	private:

	/**
	* Abstract pointer to AbstractVVVVertex
	*/
	vector<AbstractVVVVertexPtr> vertex_;

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

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

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

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

	/**
	* Abstract pointer to AbstractVVVVertex for QCD radiation from the 4-point vertex
	*/
	map<ShowerInteraction,AbstractVVVVVertexPtr> fourPointVertex_;

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

	/**
	* Vector wavefunctions
	*/
	mutable vector<Helicity::VectorWaveFunction> vectors_[3];

	private:

	/**
	* Members for the POWHEG correction
	*/
	//@{
	/**
	* Spin density matrix for 3 body decay
	*/
	mutable RhoDMatrix rho3_;

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

	/**
	* Vector wavefunctions
	*/
	mutable vector<Helicity::VectorWaveFunction> vectors3_[2];

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

	}

	#endif /* HERWIG_VVVDecayer_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 3:15 PM (1 d, 16 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3804927
Default Alt Text: (117 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions