No OneTemporary
Actions

Size

147 KB

Subscribers

None

View Options

	diff --git a/Decay/General/GeneralThreeBodyDecayer.cc b/Decay/General/GeneralThreeBodyDecayer.cc
	--- a/Decay/General/GeneralThreeBodyDecayer.cc
	+++ b/Decay/General/GeneralThreeBodyDecayer.cc
	@@ -1,830 +1,837 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the GeneralThreeBodyDecayer class.
	//

	#include "GeneralThreeBodyDecayer.h"
	#include "Herwig++/Decay/DecayPhaseSpaceMode.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig++/PDT/ThreeBodyAllOnCalculator.h"

	using namespace Herwig;

	void GeneralThreeBodyDecayer::persistentOutput(PersistentOStream & os) const {
	os << _incoming << _outgoing << _diagrams << _diagmap << _colour << _colourLargeNC
	<< _nflow << _widthopt << _reftag << _reftagcc;
	}

	void GeneralThreeBodyDecayer::persistentInput(PersistentIStream & is, int) {
	is >> _incoming >> _outgoing >> _diagrams >> _diagmap >> _colour >> _colourLargeNC
	>> _nflow >> _widthopt >> _reftag >> _reftagcc;
	}

	AbstractClassDescription<GeneralThreeBodyDecayer>
	GeneralThreeBodyDecayer::initGeneralThreeBodyDecayer;
	// Definition of the static class description member.

	void GeneralThreeBodyDecayer::Init() {

	static ClassDocumentation<GeneralThreeBodyDecayer> documentation
	("The GeneralThreeBodyDecayer class is the base class for the implementation of"
	" all three body decays based on spin structures in Herwig++.");

	static Switch<GeneralThreeBodyDecayer,unsigned int> interfaceWidthOption
	("WidthOption",
	"Option for the treatment of the widths of the intermediates",
	&GeneralThreeBodyDecayer::_widthopt, 1, false, false);
	static SwitchOption interfaceWidthOptionFixed
	(interfaceWidthOption,
	"Fixed",
	"Use fixed widths",
	1);
	static SwitchOption interfaceWidthOptionRunning
	(interfaceWidthOption,
	"Running",
	"Use running widths",
	2);
	static SwitchOption interfaceWidthOptionZero
	(interfaceWidthOption,
	"Zero",
	"Set the widths to zero",
	3);

	}

	ParticleVector GeneralThreeBodyDecayer::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;
	}

	int GeneralThreeBodyDecayer::
	modeNumber(bool & cc, tcPDPtr in, const tPDVector & outin) const {
	assert( !_reftag.empty() && !_reftagcc.empty() );
	// check number of outgoing particles
	if( outin.size() != 3 \|\| abs(in->id()) != _incoming->id() ) return -1;
	OrderedParticles testmode(outin.begin(), outin.end());
	OrderedParticles::const_iterator dit = testmode.begin();
	string testtag(in->name() + "->");
	for( unsigned int i = 1; dit != testmode.end(); ++dit, ++i) {
	testtag += (**dit).name();
	if( i != 3 ) testtag += string(",");
	}
	if( testtag == _reftag ) {
	cc = false;
	return 0;
	}
	else if ( testtag == _reftagcc ) {
	cc = true;
	return 0;
	}
	else return -1;
	}

	bool GeneralThreeBodyDecayer::setDecayInfo(PDPtr incoming,
	vector<PDPtr> outgoing,
	- const vector<TBDiagram> & process) {
	+ const vector<TBDiagram> & process,
	+ double symfac) {
	// set the member variables from the info supplied
	_incoming = incoming;
	_outgoing = outgoing;
	_diagrams = process;
	assert( _outgoing.size() == 3 );
	// Construct reference tags for testing in modeNumber function
	OrderedParticles refmode(_outgoing.begin(), _outgoing.end());
	OrderedParticles::const_iterator dit = refmode.begin();
	_reftag = _incoming->name() + "->";
	for( unsigned int i = 1; dit != refmode.end(); ++dit, ++i) {
	_reftag += (**dit).name();
	if( i != 3 ) _reftag += string(",");
	}
	//CC-mode
	refmode.clear();
	_reftagcc = _incoming->CC() ? _incoming->CC()->name() :
	_incoming->name();
	_reftagcc += "->";
	for( unsigned int i = 0; i < 3; ++i ) {
	if( _outgoing[i]->CC() ) refmode.insert( _outgoing[i]->CC() );
	else refmode.insert( _outgoing[i] );
	}
	dit = refmode.begin();
	for( unsigned int i = 1; dit != refmode.end(); ++dit , ++i) {
	_reftagcc += (**dit).name();
	if( i != 3 ) _reftagcc += string(",");
	}
	// set the colour factors and return the answer
	- return setColourFactors();
	+ return setColourFactors(symfac);
	}

	void GeneralThreeBodyDecayer::doinit() {
	DecayIntegrator::doinit();
	// create the phase space integrator
	tPDVector extpart(1,_incoming);
	extpart.insert(extpart.end(),_outgoing.begin(),_outgoing.end());
	// create the integration channels for the decay
	DecayPhaseSpaceModePtr mode(new_ptr(DecayPhaseSpaceMode(extpart,this,true)));
	DecayPhaseSpaceChannelPtr newchannel;
	// create the phase-space channels for the integration
	unsigned int nmode(0);
	for(unsigned int ix=0;ix<_diagrams.size();++ix) {
	if(_diagrams[ix].channelType==TBDiagram::fourPoint\|\|
	_diagrams[ix].channelType==TBDiagram::UNDEFINED) continue;
	// create the new channel
	newchannel=new_ptr(DecayPhaseSpaceChannel(mode));
	if(_diagrams[ix].channelType==TBDiagram::channel23) {
	newchannel->addIntermediate(extpart[0],0,0.0,-1,1);
	newchannel->addIntermediate(_diagrams[ix].intermediate,0,0.0, 2,3);
	}
	else if(_diagrams[ix].channelType==TBDiagram::channel13) {
	newchannel->addIntermediate(extpart[0],0,0.0,-1,2);
	newchannel->addIntermediate(_diagrams[ix].intermediate,0,0.0, 1,3);
	}
	else if(_diagrams[ix].channelType==TBDiagram::channel12) {
	newchannel->addIntermediate(extpart[0],0,0.0,-1,3);
	newchannel->addIntermediate(_diagrams[ix].intermediate,0,0.0, 1,2);
	}
	_diagmap.push_back(ix);
	mode->addChannel(newchannel);
	++nmode;
	}
	if(nmode==0) {
	string mode = extpart[0]->PDGName() + "->";
	for(unsigned int ix=1;ix<extpart.size();++ix) mode += extpart[ix]->PDGName() + " ";
	throw Exception() << "No decay channels in GeneralThreeBodyDecayer::"
	<< "doinit() for " << mode << "\n" << Exception::runerror;
	}
	// add the mode
	vector<double> wgt(nmode,1./double(nmode));
	addMode(mode,1.,wgt);
	}

	double GeneralThreeBodyDecayer::
	threeBodyMatrixElement(const int imode, const Energy2 q2,
	const Energy2 s3, const Energy2 s2,
	const Energy2 s1, const Energy m1,
	const Energy m2, const Energy m3) const {
	// calculate the momenta of the outgoing particles
	Energy m0=sqrt(q2);
	// energies
	Energy eout[3] = {0.5*(q2+sqr(m1)-s1)/m0,
	0.5*(q2+sqr(m2)-s2)/m0,
	0.5*(q2+sqr(m3)-s3)/m0};
	// magnitudes of the momenta
	Energy pout[3] = {sqrt(sqr(eout[0])-sqr(m1)),
	sqrt(sqr(eout[1])-sqr(m2)),
	sqrt(sqr(eout[2])-sqr(m3))};
	double cos2 = 0.5*(sqr(pout[0])+sqr(pout[1])-sqr(pout[2]))/pout[0]/pout[1];
	double cos3 = 0.5*(sqr(pout[0])-sqr(pout[1])+sqr(pout[2]))/pout[0]/pout[2];
	double sin2 = sqrt(1.-sqr(cos2)), sin3 = sqrt(1.-sqr(cos3));
	Lorentz5Momentum out[3]=
	{Lorentz5Momentum( ZERO , ZERO , pout[0] , eout[0] , m1),
	Lorentz5Momentum( pout[1]sin2 , ZERO , -pout[1]cos2 , eout[1] , m2),
	Lorentz5Momentum( -pout[2]sin3 , ZERO , -pout[2]cos3 , eout[2] , m3)};
	// create the incoming
	PPtr inpart=mode(imode)->externalParticles(0)->
	produceParticle(Lorentz5Momentum(sqrt(q2)));
	// and outgoing particles
	ParticleVector decay;
	for(unsigned int ix=1;ix<4;++ix)
	decay.push_back(mode(imode)->externalParticles(ix)->produceParticle(out[ix-1]));
	// return the matrix element
	return me2(-1,*inpart,decay,Initialize);
	}

	double GeneralThreeBodyDecayer::brat(const DecayMode &, const Particle & p,
	double oldbrat) const {
	ParticleVector children = p.children();
	if( children.size() != 3 \|\| !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()),
	make_pair(children[2]->dataPtr(), children[2]->mass()) );
	Energy width = p.data().widthGenerator()->width(p.data(), scale);
	return pwidth/width;
	}

	Energy GeneralThreeBodyDecayer::partialWidth(PMPair inpart, PMPair outa,
	PMPair outb, PMPair outc) const {
	if(inpart.second<outa.second+outb.second+outc.second) return ZERO;
	// create the object to calculate the width if it doesn't all ready exist
	if(!_widthcalc) {
	string tag = _incoming->name() + "->";
	tag += _outgoing[0]->name() + "," + _outgoing[1]->name() + ","
	+ _outgoing[2]->name() + ";";
	DMPtr dm = generator()->findDecayMode(tag);
	_widthcalc = threeBodyMEIntegrator(*dm);
	}
	return _widthcalc->partialWidth(sqr(inpart.second));
	}

	void GeneralThreeBodyDecayer::
	colourConnections(const Particle & parent,
	const ParticleVector & out) const {
	// first extract the outgoing particles and intermediate
	PPtr inter;
	ParticleVector outgoing;
	if(!generateIntermediates()) {
	outgoing=out;
	}
	else {
	// find the diagram
	unsigned int idiag = diagramMap()[mode(imode())->selectedChannel()];
	PPtr child;
	for(unsigned int ix=0;ix<out.size();++ix) {
	if(out[ix]->children().empty()) child = out[ix];
	else inter = out[ix];
	}
	outgoing.resize(3);
	switch(_diagrams[idiag].channelType) {
	case TBDiagram::channel23:
	outgoing[0] = child;
	outgoing[1] = inter->children()[0];
	outgoing[2] = inter->children()[1];
	break;
	case TBDiagram::channel13:
	outgoing[0] = inter->children()[0];
	outgoing[1] = child;
	outgoing[2] = inter->children()[1];
	break;
	case TBDiagram::channel12:
	outgoing[0] = inter->children()[0];
	outgoing[1] = inter->children()[1];
	outgoing[2] = child;
	break;
	default:
	throw Exception() << "unknown diagram type in GeneralThreeBodyDecayer::"
	<< "colourConnections()" << Exception::runerror;
	}
	}
	// extract colour of the incoming and outgoing particles
	PDT::Colour inColour(parent.data().iColour());
	vector<PDT::Colour> outColour;
	vector<int> singlet,octet,triplet,antitriplet;
	for(unsigned int ix=0;ix<outgoing.size();++ix) {
	outColour.push_back(outgoing[ix]->data().iColour());
	switch(outColour.back()) {
	case PDT::Colour0 :
	singlet.push_back(ix);
	break;
	case PDT::Colour3 :
	triplet.push_back(ix);
	break;
	case PDT::Colour3bar:
	antitriplet.push_back(ix);
	break;
	case PDT::Colour8 :
	octet.push_back(ix);
	break;
	default:
	throw Exception() << "Unknown colour for particle in GeneralThreeBodyDecayer::"
	<< "colourConnections()" << Exception::runerror;
	}
	}
	// colour neutral decaying particle
	if ( inColour == PDT::Colour0) {
	// options are all neutral or triplet/antitriplet+ neutral
	if(singlet.size()==3) return;
	else if(singlet.size()==1&&triplet.size()==1&&antitriplet.size()==1) {
	outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
	// add intermediate if needed
	if(inter&&inter->coloured()) {
	if(inter->dataPtr()->iColour()==PDT::Colour3)
	outgoing[triplet[0]]->colourLine()->addColoured(inter);
	else if(inter->dataPtr()->iColour()==PDT::Colour3bar)
	outgoing[triplet[0]]->colourLine()->addAntiColoured(inter);
	}
	}
	else if(octet.size()==1&&triplet.size()==1&&antitriplet.size()==1) {
	outgoing[ triplet[0]]->antiColourNeighbour(outgoing[octet[0]]);
	outgoing[antitriplet[0]]-> colourNeighbour(outgoing[octet[0]]);
	if(inter&&inter->coloured()) {
	if(inter->dataPtr()->iColour()==PDT::Colour3)
	outgoing[antitriplet[0]]->antiColourLine()->addColoured(inter);
	else if(inter->dataPtr()->iColour()==PDT::Colour3bar)
	outgoing[ triplet[0]]-> colourLine()->addAntiColoured(inter);
	else if(inter->dataPtr()->iColour()==PDT::Colour8) {
	outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter);
	outgoing[ triplet[0]]-> colourLine()->addColoured(inter);
	}
	}
	}
	else if(triplet.size()==3) {
	tColinePtr col[3] = {ColourLine::create(outgoing[0]),
	ColourLine::create(outgoing[1]),
	ColourLine::create(outgoing[2])};
	col[0]->setSourceNeighbours(col[1],col[2]);
	}
	else if(antitriplet.size()==3) {
	tColinePtr col[3] = {ColourLine::create(outgoing[0],true),
	ColourLine::create(outgoing[1],true),
	ColourLine::create(outgoing[2],true)};
	col[0]->setSinkNeighbours(col[1],col[2]);
	}
	else {
	string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
	+ out[1]->PDGName() + " " + out[2]->PDGName();
	throw Exception()
	<< "Unknown colour structure in GeneralThreeBodyDecayer::"
	<< "colourConnections() for singlet decaying particle "
	<< mode << Exception::runerror;
	}
	}
	// colour triplet decaying particle
	else if( inColour == PDT::Colour3) {
	if(singlet.size()==2&&triplet.size()==1) {
	outgoing[triplet[0]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	if(inter&&inter->coloured())
	outgoing[triplet[0]]->colourLine()->addColoured(inter);
	}
	else if(antitriplet.size()==1&&triplet.size()==2) {
	if(colourFlow()==0) {
	outgoing[triplet[0]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[antitriplet[0]]->colourNeighbour(outgoing[triplet[1]]);
	if(inter&&inter->coloured()) {
	switch (inter->dataPtr()->iColour()) {
	case PDT::Colour8:
	inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[triplet[1]]->colourLine()->addAntiColoured(inter);
	break;
	default:
	string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
	+ out[1]->PDGName() + " " + out[2]->PDGName();
	throw Exception() << "Unknown colour for intermediate in "
	<< "GeneralThreeBodyDecayer::"
	<< "colourConnections() for "
	<< "decaying colour triplet "
	<< mode << Exception::runerror;
	}
	}
	}
	else {
	outgoing[triplet[1]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[antitriplet[0]]->colourNeighbour(outgoing[triplet[0]]);
	if(inter&&inter->coloured()) {
	switch (inter->dataPtr()->iColour()) {
	case PDT::Colour8:
	inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[triplet[0]]->colourLine()->addAntiColoured(inter);
	break;
	default:
	string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
	+ out[1]->PDGName() + " " + out[2]->PDGName();
	throw Exception() << "Unknown colour for intermediate in "
	<< "GeneralThreeBodyDecayer::"
	<< "colourConnections() for "
	<< "decaying colour triplet "
	<< mode << Exception::runerror;
	}
	}
	}
	}
	else if (singlet.size()==1&&triplet.size()==1&&octet.size()==1) {
	if(inter) {
	if(inter->children()[0]->dataPtr()->iColour()==PDT::Colour8 \|\|
	inter->children()[1]->dataPtr()->iColour()==PDT::Colour8) {
	inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[octet[0]]->incomingColour(inter);
	outgoing[octet[0]]->colourNeighbour(outgoing[triplet[0]]);
	}
	else {
	outgoing[octet[0]]->incomingColour(inter);
	outgoing[octet[0]]->colourNeighbour(inter);
	outgoing[triplet[0]]->incomingColour(inter);
	}
	}
	else {
	outgoing[octet[0]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[octet[0]]->colourNeighbour(outgoing[triplet[0]]);
	}
	}
	else {
	string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
	+ out[1]->PDGName() + " " + out[2]->PDGName();
	throw Exception()
	<< "Unknown colour structure in GeneralThreeBodyDecayer::"
	<< "colourConnections() for triplet decaying particle "
	<< mode << Exception::runerror;
	}
	}
	else if( inColour == PDT::Colour3bar) {
	if(singlet.size()==2&&antitriplet.size()==1) {
	outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	}
	else if(antitriplet.size()==2&&triplet.size()==1) {
	if(colourFlow()==0) {
	outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[1]]);
	if(inter&&inter->coloured()) {
	switch (inter->dataPtr()->iColour()) {
	case PDT::Colour8:
	inter->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[antitriplet[1]]->antiColourLine()->addAntiColoured(inter);
	break;
	default:
	string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
	+ out[1]->PDGName() + " " + out[2]->PDGName();
	throw Exception() << "Unknown colour for intermediate in"
	<< " GeneralThreeBodyDecayer::"
	<< "colourConnections() for "
	<< "decaying colour antitriplet "
	<< mode << Exception::runerror;
	}
	}
	}
	else {
	outgoing[antitriplet[1]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
	if(inter&&inter->coloured()) {
	switch (inter->dataPtr()->iColour()) {
	case PDT::Colour8:
	inter->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter);
	break;
	default:
	string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
	+ out[1]->PDGName() + " " + out[2]->PDGName();
	throw Exception() << "Unknown colour for intermediate in "
	<< "GeneralThreeBodyDecayer::"
	<< "colourConnections() for "
	<< "decaying colour antitriplet "
	<< mode << Exception::runerror;
	}
	}
	}
	}
	else if (singlet.size()==1&&antitriplet.size()==1&&octet.size()==1) {
	if(inter) {
	if(inter->children()[0]->dataPtr()->iColour()==PDT::Colour8 \|\|
	inter->children()[1]->dataPtr()->iColour()==PDT::Colour8) {
	inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[octet[0]]->incomingAntiColour(inter);
	outgoing[octet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
	}
	else {
	outgoing[octet[0]]->incomingAntiColour(inter);
	outgoing[octet[0]]->antiColourNeighbour(inter);
	outgoing[antitriplet[0]]->incomingAntiColour(inter);
	}
	}
	else {
	outgoing[octet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	outgoing[octet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
	}
	}
	else {
	string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
	+ out[1]->PDGName() + " " + out[2]->PDGName();
	throw Exception()
	<< "Unknown colour structure in GeneralThreeBodyDecayer::"
	<< "colourConnections() for anti-triplet decaying particle"
	<< mode << Exception::runerror;
	}
	}
	else if( inColour == PDT::Colour8) {
	if(triplet.size()==1&&antitriplet.size()==1&&singlet.size()==1) {
	outgoing[ triplet[0]]->incomingColour (const_ptr_cast<tPPtr>(&parent));
	outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	if(inter&&inter->coloured()) {
	switch (inter->dataPtr()->iColour()) {
	case PDT::Colour3:
	outgoing[triplet[0]]->colourLine()->addColoured(inter);
	break;
	case PDT::Colour3bar:
	outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter);
	break;
	default:
	string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
	+ out[1]->PDGName() + " " + out[2]->PDGName();
	throw Exception() << "Unknown colour for intermediate"
	<< " in GeneralThreeBodyDecayer::"
	<< "colourConnections() for "
	<< "decaying colour octet "
	<< mode << Exception::runerror;
	}
	}
	}
	else {
	string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
	+ out[1]->PDGName() + " " + out[2]->PDGName();
	throw Exception()
	<< "Unknown colour structure in GeneralThreeBodyDecayer::"
	<< "colourConnections() for octet decaying particle"
	<< mode << Exception::runerror;
	}
	}
	}

	void GeneralThreeBodyDecayer::
	constructIntegratorChannels(vector<int> & intype, vector<Energy> & inmass,
	vector<Energy> & inwidth, vector<double> & inpow,
	vector<double> & inweights) const {
	Energy min = incoming()->mass();
	int nchannel(0);
	pair<int,Energy> imin[4]={make_pair(-1,-1.GeV),make_pair(-1,-1.GeV),
	make_pair(-1,-1.GeV),make_pair(-1,-1.GeV)};
	for(unsigned int iy=0;iy<_diagmap.size();++iy) {
	unsigned int ix=_diagmap[iy];
	Energy deltam(min);
	if(getProcessInfo()[ix].channelType==TBDiagram::fourPoint) continue;
	int itype(0);
	if (getProcessInfo()[ix].channelType==TBDiagram::channel23) {
	deltam -= outgoing()[0]->mass();
	itype = 3;
	}
	else if(getProcessInfo()[ix].channelType==TBDiagram::channel13) {
	deltam -= outgoing()[1]->mass();
	itype = 2;
	}
	else if(getProcessInfo()[ix].channelType==TBDiagram::channel12) {
	deltam -= outgoing()[2]->mass();
	itype = 1;
	}
	deltam -= getProcessInfo()[ix].intermediate->mass();
	// if(deltam<ZERO&&getProcessInfo()[ix].intermediate->width()>ZERO) {
	if(deltam<ZERO) {
	if (imin[itype].first < 0 ) imin[itype] = make_pair(ix,deltam);
	else if (imin[itype].second<deltam) imin[itype] = make_pair(ix,deltam);
	}
	if(deltam<ZERO) continue;
	if(getProcessInfo()[ix].intermediate->id()!=ParticleID::gamma) {
	intype.push_back(itype);
	inpow.push_back(0.);
	inmass.push_back(getProcessInfo()[ix].intermediate->mass());
	inwidth.push_back(getProcessInfo()[ix].intermediate->width());
	++nchannel;
	}
	else if(getProcessInfo()[ix].intermediate->id()==ParticleID::gamma) {
	intype.push_back(itype);
	inpow.push_back(-2.);
	inmass.push_back(-1.*GeV);
	inwidth.push_back(-1.*GeV);
	++nchannel;
	}
	}
	for(unsigned int ix=1;ix<4;++ix) {
	if(imin[ix].first>=0) {
	intype.push_back(ix);
	if(getProcessInfo()[imin[ix].first].intermediate->id()!=ParticleID::gamma) {
	inpow.push_back(0.);
	inmass.push_back(getProcessInfo()[imin[ix].first].intermediate->mass());
	inwidth.push_back(getProcessInfo()[imin[ix].first].intermediate->width());
	}
	else {
	inpow.push_back(-2.);
	inmass.push_back(-1.*GeV);
	inwidth.push_back(-1.*GeV);
	}
	++nchannel;
	}
	}
	inweights = vector<double>(nchannel,1./double(nchannel));
	}

	-bool GeneralThreeBodyDecayer::setColourFactors() {
	+bool GeneralThreeBodyDecayer::setColourFactors(double symfac) {
	string name = _incoming->PDGName() + "->";
	vector<int> sng,trip,atrip,oct;
	unsigned int iloc(0);
	for(vector<PDPtr>::const_iterator it = _outgoing.begin();
	it != _outgoing.end();++it) {
	name += (**it).PDGName() + " ";
	if ((**it).iColour() == PDT::Colour0 ) sng.push_back(iloc) ;
	else if((**it).iColour() == PDT::Colour3 ) trip.push_back(iloc) ;
	else if((**it).iColour() == PDT::Colour3bar ) atrip.push_back(iloc);
	else if((**it).iColour() == PDT::Colour8 ) oct.push_back(iloc) ;
	++iloc;
	}
	// colour neutral decaying particle
	if ( _incoming->iColour() == PDT::Colour0) {
	// options are all neutral or triplet/antitriplet+ neutral
	if(sng.size()==3) {
	_nflow = 1;
	_colour = vector<DVector>(1,DVector(1,1.));
	_colourLargeNC = vector<DVector>(1,DVector(1,1.));
	}
	else if(sng.size()==1&&trip.size()==1&&atrip.size()==1) {
	_nflow = 1;
	_colour = vector<DVector>(1,DVector(1,3.));
	_colourLargeNC = vector<DVector>(1,DVector(1,3.));
	}
	else if(trip.size()==1&&atrip.size()==1&&oct.size()==1) {
	_nflow = 1;
	_colour = vector<DVector>(1,DVector(1,4.));
	_colourLargeNC = vector<DVector>(1,DVector(1,4.));
	}
	else if( trip.size() == 3 \|\| atrip.size() == 3 ) {
	_nflow = 1;
	_colour = vector<DVector>(1,DVector(1,6.));
	_colourLargeNC = vector<DVector>(1,DVector(1,6.));
	for(unsigned int ix=0;ix<_diagrams.size();++ix) {
	tPDPtr inter = _diagrams[ix].intermediate;
	if(inter->CC()) inter = inter->CC();
	unsigned int io[2]={1,2};
	double sign = _diagrams[ix].channelType == TBDiagram::channel13 ? -1. : 1.;
	for(unsigned int iy=0;iy<3;++iy) {
	if (iy==1) io[0]=0;
	else if(iy==2) io[1]=1;
	tPDVector decaylist = _diagrams[ix].vertices.second->search(iy, inter);
	if(decaylist.empty()) continue;
	bool found=false;
	for(unsigned int iz=0;iz<decaylist.size();iz+=3) {
	if(decaylist[iz+io[0]]->id()==_diagrams[ix].outgoingPair.first &&
	decaylist[iz+io[1]]->id()==_diagrams[ix].outgoingPair.second) {
	sign *= 1.;
	found = true;
	}
	else if(decaylist[iz+io[0]]->id()==_diagrams[ix].outgoingPair.second &&
	decaylist[iz+io[1]]->id()==_diagrams[ix].outgoingPair.first ) {
	sign *= -1.;
	found = true;
	}
	}
	if(found) {
	if(iy==1) sign *=-1.;
	break;
	}
	}
	_diagrams[ix]. colourFlow = vector<CFPair>(1,make_pair(1,sign));
	_diagrams[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(1,sign));
	}
	}
	else {
	generator()->log() << "Unknown colour flow structure for "
	<< "colour neutral decay "
	<< name
	<< " in getColourFactors(), omitting decay\n";
	return false;
	}
	}
	// colour triplet decaying particle
	else if( _incoming->iColour() == PDT::Colour3) {
	if(sng.size()==2&&trip.size()==1) {
	_nflow = 1;
	_colour = vector<DVector>(1,DVector(1,1.));
	_colourLargeNC = vector<DVector>(1,DVector(1,1.));
	}
	else if(trip.size()==2&&atrip.size()==1) {
	_nflow = 2;
	_colour.resize(2,DVector(2,0.));
	_colour[0][0] = 3.; _colour[0][1] = 1.;
	_colour[1][0] = 1.; _colour[1][1] = 3.;
	_colourLargeNC.resize(2,DVector(2,0.));
	_colourLargeNC[0][0] = 3.; _colourLargeNC[1][1] = 3.;
	// sort out the contribution of the different diagrams to the colour
	// flows
	for(unsigned int ix=0;ix<_diagrams.size();++ix) {
	// colour singlet intermediate
	if(_diagrams[ix].intermediate->iColour()==PDT::Colour0) {
	if(_diagrams[ix].channelType==trip[0]) {
	_diagrams[ix]. colourFlow = vector<CFPair>(1,make_pair(1,1.));
	_diagrams[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
	}
	else {
	_diagrams[ix].colourFlow = vector<CFPair>(1,make_pair(2,1.));
	_diagrams[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
	}
	}
	// colour octet intermediate
	else if(_diagrams[ix].intermediate->iColour()==PDT::Colour8) {
	if(_diagrams[ix].channelType==trip[0]) {
	vector<CFPair> flow(1,make_pair(2, 0.5 ));
	_diagrams[ix].largeNcColourFlow = flow;
	flow.push_back( make_pair(1,-1./6.));
	_diagrams[ix].colourFlow=flow;
	}
	else {
	vector<CFPair> flow(1,make_pair(1, 0.5 ));
	_diagrams[ix].largeNcColourFlow = flow;
	flow.push_back( make_pair(2,-1./6.));
	_diagrams[ix].colourFlow=flow;
	}
	}
	else {
	generator()->log() << "Unknown colour for the intermediate in "
	<< "triplet -> triplet triplet antitriplet in "
	<< "ThreeBodyDecayConstructor::getColourFactors()"
	<< " for " << name << " omitting decay\n";
	return false;
	}
	}
	}
	else if(trip.size()==1&&oct.size()==1&&sng.size()==1) {
	_nflow = 1;
	_colour = vector<DVector>(1,DVector(1,4./3.));
	_colourLargeNC = vector<DVector>(1,DVector(1,4./3.));
	}
	else {
	generator()->log() << "Unknown colour structure for "
	<< "triplet decay in "
	<< "ThreeBodyDecayConstructor::getColourFactors()"
	<< " for " << name << " omitting decay\n";
	return false;
	}
	}
	// colour antitriplet decaying particle
	else if( _incoming->iColour() == PDT::Colour3bar) {
	if(sng.size()==2&&atrip.size()==1) {
	_nflow = 1;
	_colour = vector<DVector>(1,DVector(1,1.));
	_colourLargeNC = vector<DVector>(1,DVector(1,1.));
	}
	else if(atrip.size()==2&&trip.size()==1) {
	_nflow = 2;
	_colour.resize(2,DVector(2,0.));
	_colour[0][0] = 3.; _colour[0][1] = 1.;
	_colour[1][0] = 1.; _colour[1][1] = 3.;
	_colourLargeNC.resize(2,DVector(2,0.));
	_colourLargeNC[0][0] = 3.; _colourLargeNC[1][1] = 3.;
	// sort out the contribution of the different diagrams to the colour
	// flows
	for(unsigned int ix=0;ix<_diagrams.size();++ix) {
	// colour singlet intermediate
	if(_diagrams[ix].intermediate->iColour()==PDT::Colour0) {
	if(_diagrams[ix].channelType==atrip[0]) {
	_diagrams[ix]. colourFlow = vector<CFPair>(1,make_pair(1,1.));
	_diagrams[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
	}
	else {
	_diagrams[ix].colourFlow = vector<CFPair>(1,make_pair(2,1.));
	_diagrams[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
	}
	}
	// colour octet intermediate
	else if(_diagrams[ix].intermediate->iColour()==PDT::Colour8) {
	if(_diagrams[ix].channelType==atrip[0]) {
	vector<CFPair> flow(1,make_pair(2, 0.5 ));
	_diagrams[ix].largeNcColourFlow = flow;
	flow.push_back( make_pair(1,-1./6.));
	_diagrams[ix].colourFlow=flow;
	}
	else {
	vector<CFPair> flow(1,make_pair(1, 0.5 ));
	_diagrams[ix].largeNcColourFlow = flow;
	flow.push_back( make_pair(2,-1./6.));
	_diagrams[ix].colourFlow=flow;
	}
	}
	else {
	generator()->log() << "Unknown colour for the intermediate in "
	<< "antitriplet -> antitriplet antitriplet triplet in "
	<< "ThreeBodyDecayConstructor::getColourFactors()"
	<< " for " << name << " omitting decay\n";
	return false;
	}
	}
	}
	else if(atrip.size()==1&&oct.size()==1&&sng.size()==1) {
	_nflow = 1;
	_colour = vector<DVector>(1,DVector(1,4./3.));
	_colourLargeNC = vector<DVector>(1,DVector(1,4./3.));
	}
	else {
	generator()->log() << "Unknown colour antitriplet decay in "
	<< "ThreeBodyDecayConstructor::getColourFactors()"
	<< " for " << name << " omitting decay\n";
	return false;
	}
	}
	else if( _incoming->iColour() == PDT::Colour8) {
	// triplet antitriplet
	if(trip.size() == 1 && atrip.size() == 1 && sng.size() == 1) {
	_nflow = 1;
	_colour = vector<DVector>(1,DVector(1,0.5));
	_colourLargeNC = vector<DVector>(1,DVector(1,0.5));
	}
	else {
	generator()->log() << "Unknown colour octet decay in "
	<< "ThreeBodyDecayConstructor::getColourFactors()"
	<< " for " << name << " omitting decay\n";
	return false;
	}
	}
	+ for(unsigned int ix=0;ix<_nflow;++ix) {
	+ for(unsigned int iy=0;iy<_nflow;++iy) {
	+ _colour [ix][iy] /= symfac;
	+ _colourLargeNC[ix][iy] /= symfac;
	+ }
	+ }
	if( Debug::level > 1 ) {
	generator()->log() << "Mode: " << name << " has colour factors\n";
	for(unsigned int ix=0;ix<_nflow;++ix) {
	for(unsigned int iy=0;iy<_nflow;++iy) {
	generator()->log() << _colour[ix][iy] << " ";
	}
	generator()->log() << "\n";
	}
	}
	return true;
	}
	diff --git a/Decay/General/GeneralThreeBodyDecayer.h b/Decay/General/GeneralThreeBodyDecayer.h
	--- a/Decay/General/GeneralThreeBodyDecayer.h
	+++ b/Decay/General/GeneralThreeBodyDecayer.h
	@@ -1,338 +1,339 @@
	// -- C++ --
	#ifndef HERWIG_GeneralThreeBodyDecayer_H
	#define HERWIG_GeneralThreeBodyDecayer_H
	//
	// This is the declaration of the GeneralThreeBodyDecayer class.
	//

	#include "Herwig++/Decay/DecayIntegrator.h"
	#include "Herwig++/Models/General/TBDiagram.h"
	#include "GeneralThreeBodyDecayer.fh"

	namespace Herwig {
	using namespace ThePEG;

	/**
	* Here is the documentation of the GeneralThreeBodyDecayer class.
	*
	* @see \ref GeneralThreeBodyDecayerInterfaces "The interfaces"
	* defined for GeneralThreeBodyDecayer.
	*/
	class GeneralThreeBodyDecayer: public DecayIntegrator {

	public:

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


	public:

	/**
	* The default constructor.
	*/
	GeneralThreeBodyDecayer() : _nflow(999), _widthopt(1),
	_reftag(), _reftagcc(), _iflow(999)
	{}

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

	/**
	* The matrix element to be integrated for the three-body decays as a function
	* of the invariant masses of pairs of the outgoing particles.
	* @param imode The mode for which the matrix element is needed.
	* @param q2 The scale, \e i.e. the mass squared of the decaying particle.
	* @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
	* @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
	* @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
	* @param m1 The mass of the first outgoing particle.
	* @param m2 The mass of the second outgoing particle.
	* @param m3 The mass of the third outgoing particle.
	* @return The matrix element
	*/
	virtual double threeBodyMatrixElement(const int imode, const Energy2 q2,
	const Energy2 s3, const Energy2 s2,
	const Energy2 s1, const Energy m1,
	const Energy m2, const Energy m3) const;

	/**
	* Function to return partial Width
	* @param inpart The decaying particle.
	* @param outa First decay product.
	* @param outb Second decay product.
	* @param outc Third decay product.
	*/
	virtual Energy partialWidth(PMPair inpart, PMPair outa,
	PMPair outb, PMPair outc) 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 diagrams
	*/
	bool setDecayInfo(PDPtr incoming,vector<PDPtr> outgoing,
	- const vector<TBDiagram> & process);
	+ const vector<TBDiagram> & process,
	+ double symfac);

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

	protected:

	/**
	* Access the TBDiagrams that store the required information
	* to create the diagrams
	*/
	const vector<TBDiagram> & getProcessInfo() const {
	return _diagrams;
	}

	/**
	* Incoming particle
	*/
	PDPtr incoming() const { return _incoming; }

	/**
	* Outgoing particles
	*/
	const vector<PDPtr> & outgoing() const { return _outgoing; }

	/**
	* Number of colour flows
	*/
	unsigned int numberOfFlows() const { return _nflow; }

	/**
	* Set up the colour factors
	*/
	- bool setColourFactors();
	+ bool setColourFactors(double symfac);

	/**
	* Return the matrix of colour factors
	*/
	const vector<DVector> & getColourFactors() const { return _colour; }

	/**
	* Return the matrix of colour factors
	*/
	const vector<DVector> & getLargeNcColourFactors() const {
	return _colourLargeNC;
	}

	/**
	* Get the mapping between the phase-space channel and the diagram
	*/
	const vector<unsigned int> & diagramMap() const {
	return _diagmap;
	}

	/**
	* Option for the handling of the widths of the intermediate particles
	*/
	unsigned int widthOption() const { return _widthopt; }

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

	/**
	* Method to construct the channels for the integrator to give the partial width
	* @param intype Types of the channels
	* @param inmass Mass for the channels
	* @param inwidth Width for the channels
	* @param inpow Power for the channels
	* @param inweights Weights for the channels
	*/
	void constructIntegratorChannels(vector<int> & intype, vector<Energy> & inmass,
	vector<Energy> & inwidth, vector<double> & inpow,
	vector<double> & inweights) const;

	/**
	* Set the colour flow
	* @param flow The value for the colour flow
	*/
	void colourFlow(unsigned int flow) const { _iflow = flow; }

	/**
	* Set the colour flow
	*/
	unsigned int const & colourFlow() const { return _iflow; }

	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is an abstract class with persistent data.
	*/
	static AbstractClassDescription<GeneralThreeBodyDecayer> initGeneralThreeBodyDecayer;

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

	private:

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

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

	/**
	* Store the diagrams for the decay
	*/
	vector<TBDiagram> _diagrams;

	/**
	* Map between the diagrams and the phase-space channels
	*/
	vector<unsigned int> _diagmap;

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

	/**
	* Store cololur factors for ME calc at large N_c
	*/
	vector<DVector> _colourLargeNC;

	/**
	* The number of colourflows.
	*/
	unsigned int _nflow;

	/**
	* Reference to object to calculate the partial width
	*/
	mutable WidthCalculatorBasePtr _widthcalc;

	/**
	* Option for the treatment of the widths
	*/
	unsigned int _widthopt;

	/**
	* Store a decay tag for this mode that can be tested when
	* trying to determine whether it can be generated by
	* this Decayer
	*/
	string _reftag;

	/**
	* Store a decay tag for the cc-mode that can be tested when
	* trying to determine whether it can be generated by
	* this Decayer
	*/
	string _reftagcc;

	/**
	* The colour flow
	*/
	mutable unsigned int _iflow;
	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of GeneralThreeBodyDecayer. */
	template <>
	struct BaseClassTrait<Herwig::GeneralThreeBodyDecayer,1> {
	/** Typedef of the first base class of GeneralThreeBodyDecayer. */
	typedef Herwig::DecayIntegrator NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the GeneralThreeBodyDecayer class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::GeneralThreeBodyDecayer>
	: public ClassTraitsBase<Herwig::GeneralThreeBodyDecayer> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::GeneralThreeBodyDecayer"; }
	};

	/** @endcond */

	}


	#endif /* HERWIG_GeneralThreeBodyDecayer_H */
	diff --git a/Models/General/FourBodyDecayConstructor.cc b/Models/General/FourBodyDecayConstructor.cc
	--- a/Models/General/FourBodyDecayConstructor.cc
	+++ b/Models/General/FourBodyDecayConstructor.cc
	@@ -1,418 +1,243 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the FourBodyDecayConstructor class.
	//

	#include "FourBodyDecayConstructor.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig++/Models/StandardModel/StandardModel.h"
	#include "Herwig++/Decay/General/GeneralFourBodyDecayer.h"
	#include "Herwig++/Decay/DecayPhaseSpaceMode.h"
	#include "DecayConstructor.h"
	#include <queue>

	using namespace Herwig;

	FourBodyDecayConstructor::~FourBodyDecayConstructor() {}

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

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

	void FourBodyDecayConstructor::persistentOutput(PersistentOStream & os) const {
	os << interOpt_ << widthOpt_;
	}

	void FourBodyDecayConstructor::persistentInput(PersistentIStream & is, int) {
	is >> interOpt_ >> widthOpt_;
	}

	DescribeClass<FourBodyDecayConstructor,NBodyDecayConstructorBase>
	describeFourBodyDecayConstructor("Herwig::FourBodyDecayConstructor","Herwig.so");

	void FourBodyDecayConstructor::Init() {

	static ClassDocumentation<FourBodyDecayConstructor> documentation
	("The FourBodyDecayConstructor class implements a small number"
	" of 4-body decays in general models");

	static Switch<FourBodyDecayConstructor,unsigned int> interfaceWidthOption
	("WidthOption",
	"Option for the treatment of the widths of the intermediates",
	&FourBodyDecayConstructor::widthOpt_, 1, false, false);
	static SwitchOption interfaceWidthOptionFixed
	(interfaceWidthOption,
	"Fixed",
	"Use fixed widths",
	1);
	static SwitchOption interfaceWidthOptionRunning
	(interfaceWidthOption,
	"Running",
	"Use running widths",
	2);
	static SwitchOption interfaceWidthOptionZero
	(interfaceWidthOption,
	"Zero",
	"Set the widths to zero",
	3);

	static Switch<FourBodyDecayConstructor,unsigned int> interfaceIntermediateOption
	("IntermediateOption",
	"Option for the inclusion of intermediates in the event",
	&FourBodyDecayConstructor::interOpt_, 0, false, false);
	static SwitchOption interfaceIntermediateOptionAlways
	(interfaceIntermediateOption,
	"Always",
	"Always include the intermediates",
	1);
	static SwitchOption interfaceIntermediateOptionNever
	(interfaceIntermediateOption,
	"Never",
	"Never include the intermediates",
	2);
	static SwitchOption interfaceIntermediateOptionOnlyIfOnShell
	(interfaceIntermediateOption,
	"OnlyIfOnShell",
	"Only if there are on-shell diagrams",
	0);

	}

	void FourBodyDecayConstructor::DecayList(const set<PDPtr> & particles) {
	if( particles.empty() ) return;
	NBodyDecayConstructorBase::DecayList(particles);
	}

	void FourBodyDecayConstructor::
	-createDecayMode(vector<PrototypeVertexPtr> & diagrams) {
	+createDecayMode(vector<NBDiagram> & diagrams,
	+ bool possibleOnShell, double symfac) {
	// some basic checks for the modes we are interested in
	// only looking at scalars
	- if(diagrams[0]->incoming->iSpin()!=PDT::Spin0) return;
	+ if(diagrams[0].incoming->iSpin()!=PDT::Spin0) return;
	// which decay to 4 fermions
	unsigned int nferm=0;
	- for(OrderedParticles::const_iterator it=diagrams[0]->outPart.begin();
	- it!=diagrams[0]->outPart.end();++it) {
	+ for(OrderedParticles::const_iterator it=diagrams[0].outgoing.begin();
	+ it!=diagrams[0].outgoing.end();++it) {
	if((**it).iSpin()==PDT::Spin1Half) ++nferm;
	}
	if(nferm!=4) return;
	// check for on-shell intermediates
	- bool possibleOnShell=false;
	- for(unsigned int iy=0;iy<diagrams.size();++iy)
	- possibleOnShell \|= diagrams[iy]->possibleOnShell;
	bool inter = interOpt_ == 1 \|\| (interOpt_ == 0 && possibleOnShell);
	// incoming particle
	- tPDPtr inpart = diagrams[0]->incoming;
	+ tPDPtr inpart = diagrams[0].incoming;
	// outgoing particles
	- OrderedParticles outgoing=diagrams[0]->outPart;
	+ OrderedParticles outgoing=diagrams[0].outgoing;
	// incoming particle is now unstable
	inpart->stable(false);
	// construct the tag for the decay mode
	string tag = inpart->name() + "->";
	for(OrderedParticles::const_iterator it = outgoing.begin();
	it != outgoing.end(); ++it) {
	if(it!=outgoing.begin()) tag += ",";
	tag += (**it).name();
	}
	tag += ";";
	tDMPtr dm = generator()->findDecayMode(tag);
	// if mode disabled zero BR and return
	if( decayConstructor()->disableDecayMode(tag) ) {
	// If mode alread exists, ie has been read from file,
	// disable it
	if( dm ) {
	generator()->preinitInterface(dm, "BranchingRatio", "set", "0.0");
	generator()->preinitInterface(dm, "OnOff", "set", "Off");
	}
	return;
	}
	// create mode if needed
	if( createDecayModes() && (!dm \|\| inpart->id() == ParticleID::h0) ) {
	// create the decayer
	- GeneralFourBodyDecayerPtr decayer = createDecayer(diagrams,inter);
	+ GeneralFourBodyDecayerPtr decayer = createDecayer(diagrams,inter,symfac);
	if(!decayer) {
	if(Debug::level > 1 ) generator()->log() << "Can't create the decayer for "
	<< tag << " so mode not created\n";
	return;
	}
	// create the decay mode
	tDMPtr ndm = generator()->preinitCreateDecayMode(tag);
	if(ndm) {
	string test = generator()->preinitInterface(ndm, "Decayer", "set",
	decayer->fullName());
	generator()->preinitInterface(ndm, "OnOff", "set", "On");
	Energy width =
	decayer->partialWidth(inpart,outgoing);
	setBranchingRatio(ndm, width);
	}
	else
	throw NBodyDecayConstructorError()
	<< "FourBodyDecayConstructor::createDecayMode - Needed to create "
	<< "new decaymode but one could not be created for the tag "
	<< tag << Exception::warning;
	}
	// otherwise
	else if (dm && (dm->decayer()->fullName()).find("Mambo") != string::npos) {
	// create the decayer
	- GeneralFourBodyDecayerPtr decayer = createDecayer(diagrams,inter);
	+ GeneralFourBodyDecayerPtr decayer = createDecayer(diagrams,inter,symfac);
	if(!decayer) {
	if(Debug::level > 1 ) generator()->log() << "Can't create the decayer for "
	<< tag << " so mode not created\n";
	return;
	}
	generator()->preinitInterface(dm, "Decayer", "set",
	decayer->fullName());
	}
	//update CC mode if it exists
	if( inpart->CC() )
	inpart->CC()->synchronize();
	}

	-
	-namespace {
	-
	-double factorial(const int i) {
	- if(i>1) return i*factorial(i-1);
	- else return 1.;
	-}
	-
	-void constructIdenticalSwaps(unsigned int depth,
	- vector<vector<unsigned int> > identical,
	- vector<unsigned int> channelType,
	- list<vector<unsigned int> > & swaps) {
	- if(depth==0) {
	- unsigned int size = identical.size();
	- for(unsigned ix=0;ix<size;++ix) {
	- for(unsigned int iy=2;iy<identical[ix].size();++iy)
	- identical.push_back(identical[ix]);
	- }
	- }
	- if(depth+1!=identical.size()) {
	- constructIdenticalSwaps(depth+1,identical,channelType,swaps);
	- }
	- else {
	- swaps.push_back(channelType);
	- }
	- for(unsigned int ix=0;ix<identical[depth].size();++ix) {
	- for(unsigned int iy=ix+1;iy<identical[depth].size();++iy) {
	- vector<unsigned int> newType=channelType;
	- swap(newType[identical[depth][ix]],newType[identical[depth][iy]]);
	- // at bottom of chain
	- if(depth+1==identical.size()) {
	- swaps.push_back(newType);
	- }
	- else {
	- constructIdenticalSwaps(depth+1,identical,newType,swaps);
	- }
	- }
	- }
	-}
	-
	-void identicalFromSameDecay(unsigned int & loc, const NBVertex & vertex,
	- vector<vector<unsigned int> > & sameDecay) {
	- list<pair<PDPtr,NBVertex> >::const_iterator it = vertex.vertices.begin();
	- while(it!=vertex.vertices.end()) {
	- if(it->second.incoming) {
	- identicalFromSameDecay(loc,it->second,sameDecay);
	- ++it;
	- continue;
	- }
	- ++loc;
	- long id = it->first->id();
	- ++it;
	- if(it->second.incoming) continue;
	- if(it->first->id()!=id) continue;
	- sameDecay.push_back(vector<unsigned int>());
	- sameDecay.back().push_back(loc-1);
	- while(!it->second.incoming&&it->first->id()==id) {
	- ++loc;
	- ++it;
	- sameDecay.back().push_back(loc-1);
	- }
	- };
	-}
	-
	-}
	-
	GeneralFourBodyDecayerPtr
	-FourBodyDecayConstructor::createDecayer(vector<PrototypeVertexPtr> & diagrams,
	- bool inter) const {
	+FourBodyDecayConstructor::createDecayer(vector<NBDiagram> & diagrams,
	+ bool inter, double symfac) const {
	if(diagrams.empty()) return GeneralFourBodyDecayerPtr();
	// extract the external particles for the process
	- PDPtr incoming = diagrams[0]->incoming;
	+ PDPtr incoming = diagrams[0].incoming;
	// outgoing particles
	- vector<PDPtr> outgoing(diagrams[0]->outPart.begin(),
	- diagrams[0]->outPart.end());
	- vector<NBDiagram> newDiagrams;
	- double symfac(1.);
	- // convert the diagrams
	- for(unsigned int ix=0;ix<diagrams.size();++ix) {
	- symfac = 1.;
	- NBDiagram templateDiagram = NBDiagram(diagrams[ix]);
	- // extract the ordering
	- vector<int> order(templateDiagram.channelType.size());
	- for(unsigned int iz=0;iz<order.size();++iz) {
	- order[templateDiagram.channelType[iz]-1]=iz;
	- }
	- // find any identical particles
	- vector<vector<unsigned int> > identical;
	- OrderedParticles::const_iterator it=templateDiagram.outgoing.begin();
	- unsigned int iloc=0;
	- while(it!=templateDiagram.outgoing.end()) {
	- OrderedParticles::const_iterator lt = templateDiagram.outgoing.lower_bound(*it);
	- OrderedParticles::const_iterator ut = templateDiagram.outgoing.upper_bound(*it);
	- unsigned int nx=0;
	- for(OrderedParticles::const_iterator jt=lt;jt!=ut;++jt) {++nx;}
	- if(nx==1) {
	- ++it;
	- ++iloc;
	- }
	- else {
	- symfac *= factorial(nx);
	- identical.push_back(vector<unsigned int>());
	- for(OrderedParticles::const_iterator jt=lt;jt!=ut;++jt) {
	- identical.back().push_back(order[iloc]);
	- ++iloc;
	- }
	- it = ut;
	- }
	- }
	- // that's it if there outgoing are unqiue
	- if(identical.empty()) {
	- newDiagrams.push_back(templateDiagram);
	- }
	- else {
	- // find any identical particles which shouldn't be swapped
	- unsigned int loc=0;
	- vector<vector<unsigned int> > sameDecay;
	- identicalFromSameDecay(loc,templateDiagram,sameDecay);
	- // compute the swaps
	- list<vector<unsigned int> > swaps;
	- constructIdenticalSwaps(0,identical,templateDiagram.channelType,swaps);
	- // special if identical from same decay
	- if(!sameDecay.empty()) {
	- for(vector<vector<unsigned int> >::const_iterator st=sameDecay.begin();
	- st!=sameDecay.end();++st) {
	- for(list<vector<unsigned int> >::iterator it=swaps.begin();
	- it!=swaps.end();++it) {
	- for(unsigned int iy=0;iy<(*st).size();++iy) {
	- for(unsigned int iz=iy+1;iz<(*st).size();++iz) {
	- if((it)[(st)[iy]]>(it)[(st)[iz]])
	- swap((it)[(st)[iy]],(it)[(st)[iz]]);
	- }
	- }
	- }
	- }
	- }
	- // remove any dupliciates
	- for(list<vector<unsigned int> >::iterator it=swaps.begin();
	- it!=swaps.end();++it) {
	- for(list<vector<unsigned int> >::iterator jt=it;
	- jt!=swaps.end();++jt) {
	- if(it==jt) continue;
	- bool different=false;
	- for(unsigned int iz=0;iz<(*it).size();++iz) {
	- if((it)[iz]!=(jt)[iz]) {
	- different = true;
	- break;
	- }
	- }
	- if(!different) {
	- jt = swaps.erase(jt);
	- --jt;
	- }
	- }
	- }
	- // special for identical decay products
	- if(templateDiagram.vertices.begin()->second.incoming) {
	- if( templateDiagram.vertices.begin() ->first ==
	- (++templateDiagram.vertices.begin())->first) {
	- if(*( diagrams[ix]->outgoing.begin() ->second) ==
	- *((++diagrams[ix]->outgoing.begin())->second)) {
	- for(list<vector<unsigned int> >::iterator it=swaps.begin();
	- it!=swaps.end();++it) {
	- for(list<vector<unsigned int> >::iterator jt=it;
	- jt!=swaps.end();++jt) {
	- if(it==jt) continue;
	- if((it)[0]==(jt)[2]&&(it)[1]==(jt)[3]) {
	- jt = swaps.erase(jt);
	- --jt;
	- }
	- }
	- }
	- }
	- }
	- }
	- for(list<vector<unsigned int> >::iterator it=swaps.begin();
	- it!=swaps.end();++it) {
	- newDiagrams.push_back(templateDiagram);
	- newDiagrams.back().channelType = *it;
	- }
	- }
	- }
	+ vector<PDPtr> outgoing(diagrams[0].outgoing.begin(),
	+ diagrams[0].outgoing.end());
	// get the name for the object
	string objectname ("/Herwig/Decays/");
	- string classname = DecayerClassName(incoming, diagrams[0]->outPart, objectname);
	+ string classname = DecayerClassName(incoming, diagrams[0].outgoing, objectname);
	if(classname=="") return GeneralFourBodyDecayerPtr();
	// create the object
	GeneralFourBodyDecayerPtr decayer =
	dynamic_ptr_cast<GeneralFourBodyDecayerPtr>
	(generator()->preinitCreate(classname, objectname));
	// set up the decayer and return if doesn't work
	- if(!decayer->setDecayInfo(incoming,outgoing,newDiagrams,symfac))
	+ if(!decayer->setDecayInfo(incoming,outgoing,diagrams,symfac))
	return GeneralFourBodyDecayerPtr();
	// set decayer options from base class
	setDecayerInterfaces(objectname);
	// set the width option
	ostringstream value;
	value << widthOpt_;
	generator()->preinitInterface(objectname, "WidthOption", "set", value.str());
	// set the intermediates option
	ostringstream value2;
	value2 << inter;
	generator()->preinitInterface(objectname, "GenerateIntermediates", "set",
	value2.str());
	// initialize the decayer
	decayer->init();
	// return the decayer
	return decayer;
	}

	string FourBodyDecayConstructor::DecayerClassName(tcPDPtr incoming,
	const OrderedParticles & outgoing,
	string & objname) const {
	string classname("Herwig::");
	// spins of the outgoing particles
	unsigned int ns(0),nf(0),nv(0);
	objname += incoming->PDGName() + "2";
	for(OrderedParticles::const_iterator it=outgoing.begin();
	it!=outgoing.end();++it) {
	if ((**it).iSpin()==PDT::Spin0 ) ++ns;
	else if((**it).iSpin()==PDT::Spin1Half) ++nf;
	else if((**it).iSpin()==PDT::Spin1 ) ++nv;
	objname += (**it).PDGName();
	}
	objname += "Decayer";
	if(incoming->iSpin()==PDT::Spin0) {
	if(nf==4) classname += "StoFFFFDecayer";
	else classname = "";
	}
	else {
	classname="";
	}
	return classname;
	}
	diff --git a/Models/General/FourBodyDecayConstructor.h b/Models/General/FourBodyDecayConstructor.h
	--- a/Models/General/FourBodyDecayConstructor.h
	+++ b/Models/General/FourBodyDecayConstructor.h
	@@ -1,141 +1,141 @@
	// -- C++ --
	#ifndef THEPEG_FourBodyDecayConstructor_H
	#define THEPEG_FourBodyDecayConstructor_H
	//
	// This is the declaration of the FourBodyDecayConstructor class.
	//

	#include "NBodyDecayConstructorBase.h"
	#include "ThePEG/Helicity/Vertex/VertexBase.h"
	#include "Herwig++/Decay/General/GeneralFourBodyDecayer.fh"
	#include "PrototypeVertex.h"

	namespace Herwig {

	using namespace ThePEG;
	using Helicity::VertexBasePtr;

	/**
	* Here is the documentation of the FourBodyDecayConstructor class.
	*
	* @see \ref FourBodyDecayConstructorInterfaces "The interfaces"
	* defined for FourBodyDecayConstructor.
	*/
	class FourBodyDecayConstructor: public NBodyDecayConstructorBase {

	public:

	/**
	* The default constructor.
	*/
	FourBodyDecayConstructor() :
	interOpt_(0), widthOpt_(1) {}

	/**
	* Destructor
	*/
	~FourBodyDecayConstructor();

	/**
	* Function used to determine allowed decaymodes, to be implemented
	* in derived class.
	* @param particles vector of ParticleData pointers containing
	* particles in model
	*/
	virtual void DecayList(const set<PDPtr> & particles);

	/**
	* Number of outgoing lines. Required for correct ordering.
	*/
	virtual unsigned int numBodies() const {return 4;}

	/**
	* Create a decay mode
	*/
	- void createDecayMode(vector<PrototypeVertexPtr> &);
	+ void createDecayMode(vector<NBDiagram> &,bool,double);

	/**
	* Create the decayer
	* @param diagrams The diagrams for the decay
	* @param inter Option for intermediates
	*/
	- GeneralFourBodyDecayerPtr createDecayer(vector<PrototypeVertexPtr> & diagrams,
	- bool inter) const;
	+ GeneralFourBodyDecayerPtr createDecayer(vector<NBDiagram> & diagrams,
	+ bool inter, double symfac) const;

	/**
	* Contruct the classname and object name for the Decayer
	* @param incoming The incoming particle
	* @param outgoing The decay products
	* @param objname a string containing the default path of the Decayer object
	*/
	string DecayerClassName(tcPDPtr incoming, const OrderedParticles & outgoing,
	string & objname) const;

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

	private:

	/**
	* Option for the inclusion of intermediates
	*/
	unsigned int interOpt_;

	/**
	* How to treat the widths of the intermediate particles
	*/
	unsigned int widthOpt_;

	};

	}

	#endif /* THEPEG_FourBodyDecayConstructor_H */
	diff --git a/Models/General/NBodyDecayConstructorBase.cc b/Models/General/NBodyDecayConstructorBase.cc
	--- a/Models/General/NBodyDecayConstructorBase.cc
	+++ b/Models/General/NBodyDecayConstructorBase.cc
	@@ -1,447 +1,643 @@
	// -- C++ --
	//
	// NBodyDecayConstructorBase.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 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 NBodyDecayConstructorBase class.
	//

	#include "NBodyDecayConstructorBase.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "DecayConstructor.h"
	#include "Herwig++/Models/StandardModel/StandardModel.h"

	using namespace Herwig;
	using namespace ThePEG;

	void NBodyDecayConstructorBase::persistentOutput(PersistentOStream & os ) const {
	os << init_ << iteration_ << points_ << info_ << decayConstructor_
	<< removeOnShell_ << excludeEffective_ << includeTopOnShell_
	<< excludedVerticesVector_ << excludedVerticesSet_
	<< excludedParticlesVector_ << excludedParticlesSet_;
	}

	void NBodyDecayConstructorBase::persistentInput(PersistentIStream & is , int) {
	is >> init_ >> iteration_ >> points_ >> info_ >> decayConstructor_
	>> removeOnShell_ >> excludeEffective_ >> includeTopOnShell_
	>> excludedVerticesVector_ >> excludedVerticesSet_
	>> excludedParticlesVector_ >> excludedParticlesSet_;
	}

	AbstractClassDescription<NBodyDecayConstructorBase>
	NBodyDecayConstructorBase::initNBodyDecayConstructorBase;
	// Definition of the static class description member.

	void NBodyDecayConstructorBase::Init() {

	static ClassDocumentation<NBodyDecayConstructorBase> documentation
	("The NBodyDecayConstructorBase class is the base class for the automatic"
	"construction of the decay modes");

	static Switch<NBodyDecayConstructorBase,bool> interfaceInitializeDecayers
	("InitializeDecayers",
	"Initialize new decayers",
	&NBodyDecayConstructorBase::init_, true, false, false);
	static SwitchOption interfaceInitializeDecayersInitializeDecayersOn
	(interfaceInitializeDecayers,
	"Yes",
	"Initialize new decayers to find max weights",
	true);
	static SwitchOption interfaceInitializeDecayersoff
	(interfaceInitializeDecayers,
	"No",
	"Use supplied weights for integration",
	false);

	static Parameter<NBodyDecayConstructorBase,int> interfaceInitIteration
	("InitIteration",
	"Number of iterations to optimise integration weights",
	&NBodyDecayConstructorBase::iteration_, 1, 0, 10,
	false, false, true);

	static Parameter<NBodyDecayConstructorBase,int> interfaceInitPoints
	("InitPoints",
	"Number of points to generate when optimising integration",
	&NBodyDecayConstructorBase::points_, 1000, 100, 100000000,
	false, false, true);

	static Switch<NBodyDecayConstructorBase,bool> interfaceOutputInfo
	("OutputInfo",
	"Whether to output information about the decayers",
	&NBodyDecayConstructorBase::info_, false, false, false);
	static SwitchOption interfaceOutputInfoOff
	(interfaceOutputInfo,
	"No",
	"Do not output information regarding the created decayers",
	false);
	static SwitchOption interfaceOutputInfoOn
	(interfaceOutputInfo,
	"Yes",
	"Output information regarding the decayers",
	true);

	static Switch<NBodyDecayConstructorBase,bool> interfaceCreateDecayModes
	("CreateDecayModes",
	"Whether to create the ThePEG::DecayMode objects as well as the decayers",
	&NBodyDecayConstructorBase::createModes_, true, false, false);
	static SwitchOption interfaceCreateDecayModesOn
	(interfaceCreateDecayModes,
	"Yes",
	"Create the ThePEG::DecayMode objects",
	true);
	static SwitchOption interfaceCreateDecayModesOff
	(interfaceCreateDecayModes,
	"No",
	"Only create the Decayer objects",
	false);

	static Switch<NBodyDecayConstructorBase,unsigned int> interfaceRemoveOnShell
	("RemoveOnShell",
	"Remove on-shell diagrams as should be treated as a sequence of 1->2 decays",
	&NBodyDecayConstructorBase::removeOnShell_, 1, false, false);
	static SwitchOption interfaceRemoveOnShellYes
	(interfaceRemoveOnShell,
	"Yes",
	"Remove the diagrams if neither the production of decay or the intermediate"
	" can happen",
	1);
	static SwitchOption interfaceRemoveOnShellNo
	(interfaceRemoveOnShell,
	"No",
	"Never remove the intermediate",
	0);
	static SwitchOption interfaceRemoveOnShellProduction
	(interfaceRemoveOnShell,
	"Production",
	"Remove the diagram if the on-shell production of the intermediate is allowed",
	2);

	static RefVector<NBodyDecayConstructorBase,VertexBase> interfaceExcludedVertices
	("ExcludedVertices",
	"Vertices which are not included in the three-body decayers",
	&NBodyDecayConstructorBase::excludedVerticesVector_,
	-1, false, false, true, true, false);

	static RefVector<NBodyDecayConstructorBase,ParticleData> interfaceExcludedIntermediates
	("ExcludedIntermediates",
	"Excluded intermediate particles",
	&NBodyDecayConstructorBase::excludedParticlesVector_,
	-1, false, false, true, true, false);

	static Switch<NBodyDecayConstructorBase,bool> interfaceExcludeEffectiveVertices
	("ExcludeEffectiveVertices",
	"Exclude effectice vertices",
	&NBodyDecayConstructorBase::excludeEffective_, true, false, false);
	static SwitchOption interfaceExcludeEffectiveVerticesYes
	(interfaceExcludeEffectiveVertices,
	"Yes",
	"Exclude the effective vertices",
	true);
	static SwitchOption interfaceExcludeEffectiveVerticesNo
	(interfaceExcludeEffectiveVertices,
	"No",
	"Don't exclude the effective vertices",
	false);

	static Parameter<NBodyDecayConstructorBase,double> interfaceMinReleaseFraction
	("MinReleaseFraction",
	"The minimum energy release for a three-body decay, as a "
	"fraction of the parent mass.",
	&NBodyDecayConstructorBase::minReleaseFraction_, 1e-3, 0.0, 1.0,
	false, false, Interface::limited);

	static Switch<NBodyDecayConstructorBase,unsigned int> interfaceMaximumGaugeBosons
	("MaximumGaugeBosons",
	"Maximum number of electroweak gauge bosons"
	" to be produced as decay products",
	&NBodyDecayConstructorBase::maxBoson_, 1, false, false);
	static SwitchOption interfaceMaximumGaugeBosonsNone
	(interfaceMaximumGaugeBosons,
	"None",
	"Produce no W/Zs",
	0);
	static SwitchOption interfaceMaximumGaugeBosonsSingle
	(interfaceMaximumGaugeBosons,
	"Single",
	"Produce at most one W/Zs",
	1);
	static SwitchOption interfaceMaximumGaugeBosonsDouble
	(interfaceMaximumGaugeBosons,
	"Double",
	"Produce at most two W/Zs",
	2);
	static SwitchOption interfaceMaximumGaugeBosonsTriple
	(interfaceMaximumGaugeBosons,
	"Triple",
	"Produce at most three W/Zs",
	3);

	static Switch<NBodyDecayConstructorBase,unsigned int> interfaceMaximumNewParticles
	("MaximumNewParticles",
	"Maximum number of particles from the list of "
	"decaying particles to be allowed as decay products",
	&NBodyDecayConstructorBase::maxList_, 0, false, false);
	static SwitchOption interfaceMaximumNewParticlesNone
	(interfaceMaximumNewParticles,
	"None",
	"No particles from the list",
	0);
	static SwitchOption interfaceMaximumNewParticlesOne
	(interfaceMaximumNewParticles,
	"One",
	"A single particle from the list",
	1);
	static SwitchOption interfaceMaximumNewParticlesTwo
	(interfaceMaximumNewParticles,
	"Two",
	"Two particles from the list",
	2);
	static SwitchOption interfaceMaximumNewParticlesThree
	(interfaceMaximumNewParticles,
	"Three",
	"Three particles from the list",
	3);
	static SwitchOption interfaceMaximumNewParticlesFour
	(interfaceMaximumNewParticles,
	"Four",
	"Four particles from the list",
	4);

	static Switch<NBodyDecayConstructorBase,bool> interfaceIncludeOnShellTop
	("IncludeOnShellTop",
	"Include the on-shell diagrams involving t -> bW",
	&NBodyDecayConstructorBase::includeTopOnShell_, false, false, false);
	static SwitchOption interfaceIncludeOnShellTopYes
	(interfaceIncludeOnShellTop,
	"Yes",
	"Inlude them",
	true);
	static SwitchOption interfaceIncludeOnShellTopNo
	(interfaceIncludeOnShellTop,
	"No",
	"Don't include them",
	true);
	}

	void NBodyDecayConstructorBase::setBranchingRatio(tDMPtr dm, Energy pwidth) {
	//Need width and branching ratios for all currently created decay modes
	PDPtr parent = const_ptr_cast<PDPtr>(dm->parent());
	DecaySet modes = parent->decayModes();
	if( modes.empty() ) return;
	double dmbrat(0.);
	if( modes.size() == 1 ) {
	parent->width(pwidth);
	if( pwidth > ZERO ) parent->cTau(hbarc/pwidth);
	dmbrat = 1.;
	}
	else {
	Energy currentwidth(parent->width());
	Energy newWidth(currentwidth + pwidth);
	parent->width(newWidth);
	if( newWidth > ZERO ) parent->cTau(hbarc/newWidth);
	//need to reweight current branching fractions if there are any
	double factor = newWidth > ZERO ? double(currentwidth/newWidth) : 0.;
	for(DecaySet::const_iterator dit = modes.begin();
	dit != modes.end(); ++dit) {
	if( *dit == dm \|\| !(**dit).on() ) continue;
	double newbrat = (*dit).brat()factor;
	ostringstream brf;
	brf << setprecision(13)<< newbrat;
	generator()->preinitInterface(*dit, "BranchingRatio",
	"set", brf.str());
	}
	//set brat for current mode
	dmbrat = newWidth > ZERO ? double(pwidth/newWidth) : 0.;
	}
	ostringstream br;
	br << setprecision(13) << dmbrat;
	generator()->preinitInterface(dm, "BranchingRatio",
	"set", br.str());
	}

	void NBodyDecayConstructorBase::setDecayerInterfaces(string fullname) const {
	if( initialize() ) {
	ostringstream value;
	value << initialize();
	generator()->preinitInterface(fullname, "Initialize", "set",
	value.str());
	value.str("");
	value << iteration();
	generator()->preinitInterface(fullname, "Iteration", "set",
	value.str());
	value.str("");
	value << points();
	generator()->preinitInterface(fullname, "Points", "set",
	value.str());
	}
	// QED stuff if needed
	if(decayConstructor()->QEDGenerator())
	generator()->preinitInterface(fullname, "PhotonGenerator", "set",
	decayConstructor()->QEDGenerator()->fullName());
	string outputmodes;
	if( info() ) outputmodes = string("Output");
	else outputmodes = string("NoOutput");
	generator()->preinitInterface(fullname, "OutputModes", "set",
	outputmodes);
	}

	void NBodyDecayConstructorBase::doinit() {
	Interfaced::doinit();
	excludedVerticesSet_ = set<VertexBasePtr>(excludedVerticesVector_.begin(),
	excludedVerticesVector_.end());
	excludedParticlesSet_ = set<PDPtr>(excludedParticlesVector_.begin(),
	excludedParticlesVector_.end());
	if(removeOnShell_==0&&numBodies()>2)
	generator()->log() << "Warning: Including diagrams with on-shell "
	<< "intermediates in " << numBodies() << "-body BSM decays, this"
	<< " can lead to double counting and is not"
	<< " recommended unless you really know what you are doing\n"
	<< "This can be switched off using\n set "
	<< fullName() << ":RemoveOnShell Yes\n";
	}

	+namespace {
	+
	+double factorial(const int i) {
	+ if(i>1) return i*factorial(i-1);
	+ else return 1.;
	+}
	+
	+void constructIdenticalSwaps(unsigned int depth,
	+ vector<vector<unsigned int> > identical,
	+ vector<unsigned int> channelType,
	+ list<vector<unsigned int> > & swaps) {
	+ if(depth==0) {
	+ unsigned int size = identical.size();
	+ for(unsigned ix=0;ix<size;++ix) {
	+ for(unsigned int iy=2;iy<identical[ix].size();++iy)
	+ identical.push_back(identical[ix]);
	+ }
	+ }
	+ if(depth+1!=identical.size()) {
	+ constructIdenticalSwaps(depth+1,identical,channelType,swaps);
	+ }
	+ else {
	+ swaps.push_back(channelType);
	+ }
	+ for(unsigned int ix=0;ix<identical[depth].size();++ix) {
	+ for(unsigned int iy=ix+1;iy<identical[depth].size();++iy) {
	+ vector<unsigned int> newType=channelType;
	+ swap(newType[identical[depth][ix]],newType[identical[depth][iy]]);
	+ // at bottom of chain
	+ if(depth+1==identical.size()) {
	+ swaps.push_back(newType);
	+ }
	+ else {
	+ constructIdenticalSwaps(depth+1,identical,newType,swaps);
	+ }
	+ }
	+ }
	+}
	+
	+void identicalFromSameDecay(unsigned int & loc, const NBVertex & vertex,
	+ vector<vector<unsigned int> > & sameDecay) {
	+ list<pair<PDPtr,NBVertex> >::const_iterator it = vertex.vertices.begin();
	+ while(it!=vertex.vertices.end()) {
	+ if(it->second.incoming) {
	+ identicalFromSameDecay(loc,it->second,sameDecay);
	+ ++it;
	+ continue;
	+ }
	+ ++loc;
	+ long id = it->first->id();
	+ ++it;
	+ if(it->second.incoming) continue;
	+ if(it->first->id()!=id) continue;
	+ sameDecay.push_back(vector<unsigned int>());
	+ sameDecay.back().push_back(loc-1);
	+ while(!it->second.incoming&&it->first->id()==id) {
	+ ++loc;
	+ ++it;
	+ sameDecay.back().push_back(loc-1);
	+ }
	+ };
	+}
	+
	+}
	+
	void NBodyDecayConstructorBase::DecayList(const set<PDPtr> & particles) {
	if( particles.empty() ) return;
	// cast the StandardModel to the Hw++ one to get the vertices
	tHwSMPtr model = dynamic_ptr_cast<tHwSMPtr>(generator()->standardModel());
	unsigned int nv(model->numberOfVertices());
	// loop over the particles and create the modes
	for(set<PDPtr>::const_iterator ip=particles.begin();
	ip!=particles.end();++ip) {
	// get the decaying particle
	tPDPtr parent = *ip;
	// first create prototype 1->2 decays
	std::queue<PrototypeVertexPtr> prototypes;
	for(unsigned int iv = 0; iv < nv; ++iv) {
	VertexBasePtr vertex = model->vertex(iv);
	if(excluded(vertex)) continue;
	PrototypeVertex::createPrototypes(parent, vertex, prototypes);
	}
	// now expand them into potential decay modes
	list<vector<PrototypeVertexPtr> > modes;
	while(!prototypes.empty()) {
	// get the first prototype from the queue
	PrototypeVertexPtr proto = prototypes.front();
	prototypes.pop();
	// multiplcity too low
	if(proto->npart<numBodies()) {
	// loop over all vertices and expand
	for(unsigned int iv = 0; iv < nv; ++iv) {
	VertexBasePtr vertex = model->vertex(iv);
	if(excluded(vertex) \|\|
	proto->npart+vertex->getNpoint()>numBodies()+2) continue;
	PrototypeVertex::expandPrototypes(proto,vertex,prototypes,
	excludedParticlesSet_);
	}
	}
	// multiplcity too high disgard
	else if(proto->npart>numBodies()) {
	continue;
	}
	// right multiplicity
	else {
	// check it's kinematical allowed for physical masses
	if( proto->incomingMass() < proto->outgoingMass() ) continue;
	// and for constituent masses
	Energy outgoingMass = proto->outgoingConstituentMass();
	if( proto->incomingMass() < proto->outgoingConstituentMass() ) continue;
	// remove processes which aren't kinematically allowed within
	// the release fraction
	if( proto->incomingMass() - outgoingMass <=
	minReleaseFraction_ * proto->incomingMass() ) continue;
	// check the external particles
	if(!proto->checkExternal()) continue;
	// check if first piece on-shell
	bool onShell = proto->canBeOnShell(removeOnShell(),proto->incomingMass(),true);
	// special treatment for three-body top decays
	if(onShell) {
	if(includeTopOnShell_ &&
	abs(proto->incoming->id())==ParticleID::t && proto->npart==3) {
	unsigned int nprimary=0;
	bool foundW=false, foundQ=false;
	for(OrderedVertices::const_iterator it = proto->outgoing.begin();
	it!=proto->outgoing.end();++it) {
	if(abs(it->first->id())==ParticleID::Wplus)
	foundW = true;
	if(abs(it->first->id())==ParticleID::b \|\|
	abs(it->first->id())==ParticleID::s \|\|
	abs(it->first->id())==ParticleID::d)
	foundQ = true;
	++nprimary;
	}
	if(!(foundW&&foundQ&&nprimary==2)) continue;
	}
	else continue;
	}
	// check if should be added to an existing decaymode
	bool added = false;
	for(list<vector<PrototypeVertexPtr> >::iterator it=modes.begin();
	it!=modes.end();++it) {
	// is the same ?
	if(!(it)[0]->sameDecay(proto)) continue;
	// it is the same
	added = true;
	// check if it is a duplicate
	bool already = false;
	for(unsigned int iz = 0; iz < (*it).size(); ++iz) {
	if( (it)[iz] == *proto) {
	already = true;
	break;
	}
	}
	if(!already) (*it).push_back(proto);
	break;
	}
	if(!added) modes.push_back(vector<PrototypeVertexPtr>(1,proto));
	}
	}
	// now look at the decay modes
	for(list<vector<PrototypeVertexPtr> >::iterator mit = modes.begin();
	mit!=modes.end();++mit) {
	// count the number of gauge bosons and particles from the list
	unsigned int nlist(0),nbos(0);
	for(OrderedParticles::const_iterator it=(*mit)[0]->outPart.begin();
	it!=(*mit)[0]->outPart.end();++it) {
	if(abs((**it).id()) == ParticleID::Wplus \|\|
	abs((**it).id()) == ParticleID::Z0) ++nbos;
	if(particles.find(*it)!=particles.end()) ++nlist;
	if((it).CC() && particles.find((it).CC())!=particles.end()) ++nlist;
	}
	// if too many ignore the mode
	if(nbos > maxBoson_ \|\| nlist > maxList_) continue;
	+ // translate the prototypes into diagrams
	+ vector<NBDiagram> newDiagrams;
	+ double symfac(1.);
	+ bool possibleOnShell=false;
	+ for(unsigned int ix=0;ix<(*mit).size();++ix) {
	+ symfac = 1.;
	+ possibleOnShell \|= (*mit)[ix]->possibleOnShell;
	+ NBDiagram templateDiagram = NBDiagram((*mit)[ix]);
	+ // extract the ordering
	+ vector<int> order(templateDiagram.channelType.size());
	+ for(unsigned int iz=0;iz<order.size();++iz) {
	+ order[templateDiagram.channelType[iz]-1]=iz;
	+ }
	+ // find any identical particles
	+ vector<vector<unsigned int> > identical;
	+ OrderedParticles::const_iterator it=templateDiagram.outgoing.begin();
	+ unsigned int iloc=0;
	+ while(it!=templateDiagram.outgoing.end()) {
	+ OrderedParticles::const_iterator lt = templateDiagram.outgoing.lower_bound(*it);
	+ OrderedParticles::const_iterator ut = templateDiagram.outgoing.upper_bound(*it);
	+ unsigned int nx=0;
	+ for(OrderedParticles::const_iterator jt=lt;jt!=ut;++jt) {++nx;}
	+ if(nx==1) {
	+ ++it;
	+ ++iloc;
	+ }
	+ else {
	+ symfac *= factorial(nx);
	+ identical.push_back(vector<unsigned int>());
	+ for(OrderedParticles::const_iterator jt=lt;jt!=ut;++jt) {
	+ identical.back().push_back(order[iloc]);
	+ ++iloc;
	+ }
	+ it = ut;
	+ }
	+ }
	+ // that's it if there outgoing are unqiue
	+ if(identical.empty()) {
	+ newDiagrams.push_back(templateDiagram);
	+ continue;
	+ }
	+ // find any identical particles which shouldn't be swapped
	+ unsigned int loc=0;
	+ vector<vector<unsigned int> > sameDecay;
	+ identicalFromSameDecay(loc,templateDiagram,sameDecay);
	+ // compute the swaps
	+ list<vector<unsigned int> > swaps;
	+ constructIdenticalSwaps(0,identical,templateDiagram.channelType,swaps);
	+ // special if identical from same decay
	+ if(!sameDecay.empty()) {
	+ for(vector<vector<unsigned int> >::const_iterator st=sameDecay.begin();
	+ st!=sameDecay.end();++st) {
	+ for(list<vector<unsigned int> >::iterator it=swaps.begin();
	+ it!=swaps.end();++it) {
	+ for(unsigned int iy=0;iy<(*st).size();++iy) {
	+ for(unsigned int iz=iy+1;iz<(*st).size();++iz) {
	+ if((it)[(st)[iy]]>(it)[(st)[iz]])
	+ swap((it)[(st)[iy]],(it)[(st)[iz]]);
	+ }
	+ }
	+ }
	+ }
	+ }
	+ // remove any dupliciates
	+ for(list<vector<unsigned int> >::iterator it=swaps.begin();
	+ it!=swaps.end();++it) {
	+ for(list<vector<unsigned int> >::iterator jt=it;
	+ jt!=swaps.end();++jt) {
	+ if(it==jt) continue;
	+ bool different=false;
	+ for(unsigned int iz=0;iz<(*it).size();++iz) {
	+ if((it)[iz]!=(jt)[iz]) {
	+ different = true;
	+ break;
	+ }
	+ }
	+ if(!different) {
	+ jt = swaps.erase(jt);
	+ --jt;
	+ }
	+ }
	+ }
	+ // special for identical decay products
	+ if(templateDiagram.vertices.begin()->second.incoming) {
	+ if( templateDiagram.vertices.begin() ->first ==
	+ (++templateDiagram.vertices.begin())->first) {
	+ if(( (mit)[ix]->outgoing.begin() ->second) ==
	+ ((++(mit)[ix]->outgoing.begin())->second)) {
	+ for(list<vector<unsigned int> >::iterator it=swaps.begin();
	+ it!=swaps.end();++it) {
	+ for(list<vector<unsigned int> >::iterator jt=it;
	+ jt!=swaps.end();++jt) {
	+ if(it==jt) continue;
	+ if((it)[0]==(jt)[2]&&(it)[1]==(jt)[3]) {
	+ jt = swaps.erase(jt);
	+ --jt;
	+ }
	+ }
	+ }
	+ }
	+ }
	+ }
	+ for(list<vector<unsigned int> >::iterator it=swaps.begin();
	+ it!=swaps.end();++it) {
	+ newDiagrams.push_back(templateDiagram);
	+ newDiagrams.back().channelType = *it;
	+ }
	+ }
	// create the decay
	- createDecayMode(*mit);
	+ createDecayMode(newDiagrams,possibleOnShell,symfac);
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	+
	if( Debug::level > 1 ) {
	generator()->log() << "Mode: ";
	generator()->log() << (*mit)[0]->incoming->PDGName() << " -> ";
	for(OrderedParticles::const_iterator it=(*mit)[0]->outPart.begin();
	it!=(*mit)[0]->outPart.end();++it)
	generator()->log() << (**it).PDGName() << " ";
	generator()->log() << "\n";
	generator()->log() << "There are " << (*mit).size() << " diagrams\n";
	for(unsigned int iy=0;iy<(*mit).size();++iy) {
	generator()->log() << "Diagram: " << iy << "\n";
	generator()->log() << (mit)[iy] << "\n";
	}
	}
	}
	}
	}

	-void NBodyDecayConstructorBase::createDecayMode(vector<PrototypeVertexPtr> &) {
	+void NBodyDecayConstructorBase::createDecayMode(vector<NBDiagram> &,
	+ bool, double) {
	throw Exception() << "In NBodyDecayConstructorBase::createDecayMode() which"
	<< " should have be overridden in an inheriting class"
	<< Exception::abortnow;
	}
	diff --git a/Models/General/NBodyDecayConstructorBase.h b/Models/General/NBodyDecayConstructorBase.h
	--- a/Models/General/NBodyDecayConstructorBase.h
	+++ b/Models/General/NBodyDecayConstructorBase.h
	@@ -1,334 +1,336 @@
	// -- C++ --
	//
	// NBodyDecayConstructorBase.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_NBodyDecayConstructorBase_H
	#define HERWIG_NBodyDecayConstructorBase_H
	//
	// This is the declaration of the NBodyDecayConstructorBase class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "ThePEG/Utilities/Exception.h"
	#include "ThePEG/PDT/ParticleData.h"
	#include "NBodyDecayConstructorBase.fh"
	#include "PrototypeVertex.h"
	#include "DecayConstructor.fh"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* This is the base class for NBodyDecayConstructors. An N-body
	* decay constructor should inherit from this and implement the
	* DecayList virtual funtcion to create the decays and decayers.
	*
	* @see \ref NBodyDecayConstructorBaseInterfaces "The interfaces"
	* defined for NBodyDecayConstructor.
	*/
	class NBodyDecayConstructorBase: public Interfaced {

	public:

	/**
	* The default constructor.
	*/
	NBodyDecayConstructorBase() :
	init_(true),iteration_(1), points_(1000), info_(false),
	createModes_(true), removeOnShell_(1), excludeEffective_(true),
	minReleaseFraction_(1e-3), maxBoson_(1), maxList_(1),
	includeTopOnShell_(false) {}

	/**
	* Function used to determine allowed decaymodes, to be implemented
	* in derived class.
	* @param particles vector of ParticleData pointers containing
	* particles in model
	*/
	virtual void DecayList(const set<PDPtr> & particles);

	/**
	* Number of outgoing lines. Required for correct ordering.
	*/
	virtual unsigned int numBodies() const = 0;

	/**
	* Set the pointer to the DecayConstrcutor
	*/
	void decayConstructor(tDecayConstructorPtr d) {
	decayConstructor_ = d;
	}

	protected:

	/**
	* Method to set up the decay mode, should be overidden in inheriting class
	*/
	- virtual void createDecayMode(vector<PrototypeVertexPtr> & mode);
	+ virtual void createDecayMode(vector<NBDiagram> & mode,
	+ bool possibleOnShell,
	+ double symfac);


	/**
	* Set the branching ratio of this mode. This requires
	* calculating a new width for the decaying particle and reweighting
	* the current branching fractions.
	* @param dm The decaymode for which to set the branching ratio
	* @param pwidth The calculated width of the mode
	*/
	void setBranchingRatio(tDMPtr dm, Energy pwidth);

	/**
	* Set the interfaces of the decayers depending on the flags stored.
	* @param name Fullname of the decayer in the EventGenerator
	* including the path
	*/
	void setDecayerInterfaces(string name) const;

	/**
	* Whether to initialize decayers or not
	*/
	bool initialize() const { return init_; }

	/**
	* Number of iterations if initializing (default 1)
	*/
	int iteration() const { return iteration_; }

	/**
	* Number of points to do in initialization
	*/
	int points() const { return points_; }

	/**
	* Whether to output information on the decayers
	*/
	bool info() const { return info_; }

	/**
	* Whether to create the DecayModes as well as the Decayer objects
	*/
	bool createDecayModes() const { return createModes_; }

	/**
	* Maximum number of electroweak gauge bosons
	*/
	unsigned int maximumGaugeBosons() const { return maxBoson_;}

	/**
	* Maximum number of particles from the list whose decays we are calculating
	*/
	unsigned int maximumList() const { return maxList_;}

	/**
	* Minimum energy release fraction
	*/
	double minimumReleaseFraction() const {return minReleaseFraction_;}

	/**
	* Get the pointer to the DecayConstructor object
	*/
	tDecayConstructorPtr decayConstructor() const {
	return decayConstructor_;
	}

	/**
	* Option for on-shell particles
	*/
	unsigned int removeOnShell() const { return removeOnShell_; }

	/**
	* Check if a vertex is excluded
	*/
	bool excluded(VertexBasePtr vertex) const {
	// skip an effective vertex
	if( excludeEffective_ &&
	int(vertex->orderInGs() + vertex->orderInGem()) != int(vertex->getNpoint())-2)
	return true;
	// check if explicitly forbidden
	return excludedVerticesSet_.find(vertex)!=excludedVerticesSet_.end();
	}

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

	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is an abstract class with persistent data.
	*/
	static AbstractClassDescription<NBodyDecayConstructorBase>
	initNBodyDecayConstructorBase;

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

	private:

	/**
	* Whether to initialize decayers or not
	*/
	bool init_;

	/**
	* Number of iterations if initializing (default 1)
	*/
	int iteration_;

	/**
	* Number of points to do in initialization
	*/
	int points_;

	/**
	* Whether to output information on the decayers
	*/
	bool info_;

	/**
	* Whether to create the DecayModes as well as the Decayer objects
	*/
	bool createModes_;

	/**
	* Whether or not to remove on-shell diagrams
	*/
	unsigned int removeOnShell_;

	/**
	* Excluded Vertices
	*/
	vector<VertexBasePtr> excludedVerticesVector_;

	/**
	* Excluded Vertices
	*/
	set<VertexBasePtr> excludedVerticesSet_;

	/**
	* Excluded Particles
	*/
	vector<PDPtr> excludedParticlesVector_;

	/**
	* Excluded Particles
	*/
	set<PDPtr> excludedParticlesSet_;

	/**
	* Whether or not to exclude effective vertices
	*/
	bool excludeEffective_;

	/**
	* A pointer to the DecayConstructor object
	*/
	tDecayConstructorPtr decayConstructor_;

	/**
	* The minimum energy release for a three-body decay as a
	* fraction of the parent mass
	*/
	double minReleaseFraction_;

	/**
	* Maximum number of EW gauge bosons
	*/
	unsigned int maxBoson_;

	/**
	* Maximum number of particles from the decaying particle list
	*/
	unsigned int maxList_;

	/**
	* Include on-shell for \f$t\to b W\f$
	*/
	bool includeTopOnShell_;
	};

	/** An Exception class that can be used by all inheriting classes to
	* indicate a setup problem. */
	class NBodyDecayConstructorError : public Exception {

	public:

	NBodyDecayConstructorError() : Exception() {}

	NBodyDecayConstructorError(const string & str,
	Severity sev) : Exception(str,sev)
	{}
	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of NBodyDecayConstructorBase. */
	template <>
	struct BaseClassTrait<Herwig::NBodyDecayConstructorBase,1> {
	/** Typedef of the first base class of NBodyDecayConstructorBase. */
	typedef Interfaced NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the NBodyDecayConstructorBase class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::NBodyDecayConstructorBase>
	: public ClassTraitsBase<Herwig::NBodyDecayConstructorBase> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::NBodyDecayConstructorBase"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_NBodyDecayConstructorBase_H */
	diff --git a/Models/General/TBDiagram.h b/Models/General/TBDiagram.h
	--- a/Models/General/TBDiagram.h
	+++ b/Models/General/TBDiagram.h
	@@ -1,216 +1,212 @@
	// -- C++ --
	//
	// TBDiagram.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_TBDiagram_H
	#define HERWIG_TBDiagram_H
	//
	// This is the declaration of the TBDiagram struct.
	//

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Helicity/Vertex/VertexBase.h"
	#include "ThePEG/Utilities/EnumIO.h"
	#include "PrototypeVertex.h"

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

	/** Pair of particle ids. */
	typedef pair<long, long> IDPair;

	/** Pair of bool's*/
	typedef pair<bool, bool> BPair;

	/** Convenient typedef of VertexBasePtr */
	typedef VertexBasePtr VBPtr;

	/** Pair of VertexBasePtrs */
	typedef pair<VBPtr, VBPtr> VBPair;

	/** Pair of int,double */
	typedef pair<unsigned int, double> CFPair;

	/**
	* The TBDiagram struct contains information about a \f$1\to3\f$ decay that
	* has been automatically generated by ThreeBodyDecayConstructor.
	*/
	struct TBDiagram {

	/** Enumeration for channel type */
	enum Channel {UNDEFINED = -1, channel23=0, channel13=1, channel12=2, fourPoint=3};

	/** Standard Constructor */
	TBDiagram()
	: incoming(0), outgoing(0), outgoingPair(make_pair(0, 0)),
	channelType(UNDEFINED), colourFlow(0), ids(4, 0)
	{}

	/** Constructor taking ids as arguments.*/
	TBDiagram(long a, long b, IDPair c)
	: incoming(a), outgoing(b), outgoingPair(c),
	channelType(UNDEFINED), colourFlow(0), ids(4, 0) {
	ids[0] = a;
	ids[1] = b;
	ids[2] = c.first;
	ids[3] = c.second;
	}

	/** Incoming particle id's */
	long incoming;

	/** Outgoing particle from first vertex */
	long outgoing;

	/** Outgoing particle id's fropm resonance*/
	IDPair outgoingPair;

	/** ParticleData pointer to intermediate, null for 4-point vertices */
	PDPtr intermediate;

	/** The two vertices for the diagram */
	VBPair vertices;

	/** Enum of channel type */
	Channel channelType;

	/** Store colour flow at \f$N_c=3\f$ information */
	mutable vector<CFPair> colourFlow;

	/** Store colour flow at \f$N_c=\infty\f$ information */
	mutable vector<CFPair> largeNcColourFlow;

	/** Store the ids in a vector for easy use of comparison operator. */
	vector<long> ids;

	/**
	* Test whether this and x are the same decay
	* @param x The other process to check
	*/
	bool sameDecay(const TBDiagram & x) const {
	if(ids[0] != x.ids[0]) return false;
	bool used[4]={false,false,false,false};
	for(unsigned int ix=1;ix<4;++ix) {
	bool found=false;
	for(unsigned int iy=1;iy<4;++iy) {
	if(used[iy]) continue;
	if(ids[ix]==x.ids[iy]) {
	used[iy]=true;
	found=true;
	break;
	}
	}
	if(!found) return false;
	}
	return true;
	}

	- /** Constructor from PrototypeVertex */
	- TBDiagram(PrototypeVertexPtr vertex)
	+ /** Constructor from NBDiagram */
	+ TBDiagram(const NBDiagram & diagram)
	: channelType(UNDEFINED),
	- colourFlow (vector<CFPair>(1,make_pair(1,1.))),
	- largeNcColourFlow(vector<CFPair>(1,make_pair(1,1.))), ids(4,0) {
	- incoming = vertex->incoming->id();
	- PrototypeVertexPtr child;
	- if(vertex->outgoing.size()==2) {
	- for(OrderedVertices::const_iterator it = vertex->outgoing.begin();
	- it!=vertex->outgoing.end();++it) {
	- if(it->second) child=it->second;
	- else outgoing = it->first->id();
	- }
	+ colourFlow (diagram.colourFlow),
	+ largeNcColourFlow(diagram.largeNcColourFlow), ids(4,0) {
	+ incoming = diagram.incoming->id();
	+ if(diagram.vertices.size()==2) {
	+ const NBVertex & child = (++diagram.vertices.begin())->second;
	+ outgoing = diagram.vertices.begin()->first->id();
	unsigned int iloc=0;
	- for(OrderedVertices::const_iterator it = child->outgoing.begin();
	- it!=child->outgoing.end();++it) {
	- if (iloc==0) outgoingPair.first = it->first->id();
	- else if (iloc==1) outgoingPair.second = it->first->id();
	+ for(OrderedParticles::const_iterator it = child.outgoing.begin();
	+ it!=child.outgoing.end();++it) {
	+ if (iloc==0) outgoingPair.first = (**it).id();
	+ else if (iloc==1) outgoingPair.second = (**it).id();
	++iloc;
	}
	- intermediate = child->incoming;
	- vertices.second = child->vertex;
	+ intermediate = child.incoming;
	+ vertices.second = child.vertex;
	}
	else {
	unsigned int iloc=0;
	- for(OrderedVertices::const_iterator it = vertex->outgoing.begin();
	- it!=vertex->outgoing.end();++it) {
	- if(iloc==0) outgoing = it->first->id();
	- else if (iloc==1) outgoingPair.first = it->first->id();
	- else if (iloc==2) outgoingPair.second = it->first->id();
	+ for(OrderedParticles::const_iterator it = diagram.outgoing.begin();
	+ it!=diagram.outgoing.end();++it) {
	+ if(iloc==0) outgoing = (**it).id();
	+ else if (iloc==1) outgoingPair.first = (**it).id();
	+ else if (iloc==2) outgoingPair.second = (**it).id();
	++iloc;
	}
	}
	- vertices.first = vertex->vertex;
	+ vertices.first = diagram.vertex;
	ids[0] = incoming;
	ids[1] = outgoing;
	ids[2] = outgoingPair.first;
	ids[3] = outgoingPair.second;
	}

	};

	/**
	* Test whether two diagrams are identical.
	*/
	inline bool operator==(const TBDiagram & x, const TBDiagram & y) {
	if( x.incoming != y.incoming) return false;
	if( x.outgoing != y.outgoing) return false;
	if( x.intermediate != y.intermediate) return false;
	if( x.vertices != y.vertices ) return false;
	if( (x.outgoingPair.first == y.outgoingPair.first &&
	x.outgoingPair.second == y.outgoingPair.second ) \|\|
	(x.outgoingPair.first == y.outgoingPair.second &&
	x.outgoingPair.second == y.outgoingPair.first ) ) return true;
	else return false;
	}

	/**
	* Output to a stream
	*/
	inline ostream & operator<<(ostream & os, const TBDiagram & diag) {
	os << diag.incoming << " -> ";
	os << diag.outgoing << " + ( ";
	if(diag.intermediate) os << diag.intermediate->id();
	os << " ) -> ";
	os << diag.outgoingPair.first << " " << diag.outgoingPair.second
	<< " channel: " << diag.channelType << " ";
	for(size_t cf = 0; cf < diag.colourFlow.size(); ++cf)
	os << "(" << diag.colourFlow[cf].first << ","
	<<diag.colourFlow[cf].second << ")";
	os << '\n';
	return os;
	}

	/**
	* Output operator to allow the structure to be persistently written
	* @param os The output stream
	* @param x The TBDiagram
	*/
	inline PersistentOStream & operator<<(PersistentOStream & os,
	const TBDiagram & x) {
	os << x.incoming << x.outgoing << x.outgoingPair << x.intermediate
	<< x.vertices << oenum(x.channelType) << x.colourFlow
	<< x.largeNcColourFlow << x.ids;
	return os;
	}

	/**
	* Input operator to allow persistently written data to be read in
	* @param is The input stream
	* @param x The TBDiagram
	*/
	inline PersistentIStream & operator>>(PersistentIStream & is,
	TBDiagram & x) {
	is >> x.incoming >> x.outgoing >> x.outgoingPair >> x.intermediate
	>> x.vertices >> ienum(x.channelType) >> x.colourFlow
	>> x.largeNcColourFlow >> x.ids;
	return is;
	}

	}

	#endif /* HERWIG_TBDiagram_H */
	diff --git a/Models/General/ThreeBodyDecayConstructor.cc b/Models/General/ThreeBodyDecayConstructor.cc
	--- a/Models/General/ThreeBodyDecayConstructor.cc
	+++ b/Models/General/ThreeBodyDecayConstructor.cc
	@@ -1,288 +1,289 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ThreeBodyDecayConstructor class.
	//

	#include "ThreeBodyDecayConstructor.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig++/Models/StandardModel/StandardModel.h"
	#include "Herwig++/Decay/General/GeneralThreeBodyDecayer.h"
	#include "Herwig++/Decay/DecayPhaseSpaceMode.h"
	#include "Herwig++/PDT/ThreeBodyAllOnCalculator.h"
	#include "ThePEG/PDT/StandardMatchers.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "DecayConstructor.h"
	#include "WeakCurrentDecayConstructor.h"

	using namespace Herwig;

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

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

	void ThreeBodyDecayConstructor::persistentOutput(PersistentOStream & os) const {
	os << interOpt_ << widthOpt_;
	}

	void ThreeBodyDecayConstructor::persistentInput(PersistentIStream & is, int) {
	is >> interOpt_ >> widthOpt_;
	}

	ClassDescription<ThreeBodyDecayConstructor>
	ThreeBodyDecayConstructor::initThreeBodyDecayConstructor;
	// Definition of the static class description member.

	void ThreeBodyDecayConstructor::Init() {

	static ClassDocumentation<ThreeBodyDecayConstructor> documentation
	("The ThreeBodyDecayConstructor class constructs the three body decay modes");

	static Switch<ThreeBodyDecayConstructor,unsigned int> interfaceWidthOption
	("WidthOption",
	"Option for the treatment of the widths of the intermediates",
	&ThreeBodyDecayConstructor::widthOpt_, 1, false, false);
	static SwitchOption interfaceWidthOptionFixed
	(interfaceWidthOption,
	"Fixed",
	"Use fixed widths",
	1);
	static SwitchOption interfaceWidthOptionRunning
	(interfaceWidthOption,
	"Running",
	"Use running widths",
	2);
	static SwitchOption interfaceWidthOptionZero
	(interfaceWidthOption,
	"Zero",
	"Set the widths to zero",
	3);

	static Switch<ThreeBodyDecayConstructor,unsigned int> interfaceIntermediateOption
	("IntermediateOption",
	"Option for the inclusion of intermediates in the event",
	&ThreeBodyDecayConstructor::interOpt_, 0, false, false);
	static SwitchOption interfaceIntermediateOptionAlways
	(interfaceIntermediateOption,
	"Always",
	"Always include the intermediates",
	1);
	static SwitchOption interfaceIntermediateOptionNever
	(interfaceIntermediateOption,
	"Never",
	"Never include the intermediates",
	2);
	static SwitchOption interfaceIntermediateOptionOnlyIfOnShell
	(interfaceIntermediateOption,
	"OnlyIfOnShell",
	"Only if there are on-shell diagrams",
	0);

	}

	void ThreeBodyDecayConstructor::DecayList(const set<PDPtr> & particles) {
	if( particles.empty() ) return;
	// special for weak decays
	for(unsigned int ix=0;ix<decayConstructor()->decayConstructors().size();++ix) {
	Ptr<Herwig::WeakCurrentDecayConstructor>::pointer
	weak = dynamic_ptr_cast<Ptr<Herwig::WeakCurrentDecayConstructor>::pointer >
	(decayConstructor()->decayConstructors()[ix]);
	if(!weak) continue;
	weakMassCut_ = max(weakMassCut_,weak->massCut());
	}
	NBodyDecayConstructorBase::DecayList(particles);
	}

	GeneralThreeBodyDecayerPtr ThreeBodyDecayConstructor::
	-createDecayer(vector<TBDiagram> & diagrams, bool inter) const {
	+createDecayer(vector<TBDiagram> & diagrams, bool inter,
	+ double symfac) const {
	if(diagrams.empty()) return GeneralThreeBodyDecayerPtr();
	// extract the external particles for the process
	PDPtr incoming = getParticleData(diagrams[0].incoming);
	// outgoing particles
	OrderedParticles outgoing;
	outgoing.insert(getParticleData(diagrams[0].outgoing ));
	outgoing.insert(getParticleData(diagrams[0].outgoingPair.first ));
	outgoing.insert(getParticleData(diagrams[0].outgoingPair.second));
	// sort out ordering and labeling of diagrams
	vector<PDPtr> outVector(outgoing.begin(),outgoing.end());
	for(unsigned int ix=0;ix<diagrams.size();++ix) {
	unsigned int iy=0;
	for(;iy<3;++iy)
	if(diagrams[ix].outgoing == outVector[iy]->id()) break;
	if(diagrams[ix].channelType == TBDiagram::UNDEFINED) {
	diagrams[ix].channelType = TBDiagram::Channel(iy);
	if( ( iy == 0 && outVector[1]->id() != diagrams[ix].outgoingPair.first)\|\|
	( iy == 1 && outVector[0]->id() != diagrams[ix].outgoingPair.first)\|\|
	( iy == 2 && outVector[0]->id() != diagrams[ix].outgoingPair.first) )
	swap(diagrams[ix].outgoingPair.first, diagrams[ix].outgoingPair.second);
	}
	}
	// get the object name
	string objectname ("/Herwig/Decays/");
	string classname = DecayerClassName(incoming, outgoing, objectname);
	if(classname=="") return GeneralThreeBodyDecayerPtr();
	// create the object
	GeneralThreeBodyDecayerPtr decayer =
	dynamic_ptr_cast<GeneralThreeBodyDecayerPtr>
	(generator()->preinitCreate(classname, objectname));
	// set up the decayer and return if doesn't work
	- if(!decayer->setDecayInfo(incoming,outVector,diagrams))
	+ if(!decayer->setDecayInfo(incoming,outVector,diagrams,symfac))
	return GeneralThreeBodyDecayerPtr();
	// set decayer options from base class
	setDecayerInterfaces(objectname);
	// set the width option
	ostringstream value;
	value << widthOpt_;
	generator()->preinitInterface(objectname, "WidthOption", "set", value.str());
	// set the intermediates option
	ostringstream value2;
	value2 << inter;
	generator()->preinitInterface(objectname, "GenerateIntermediates", "set",
	value2.str());
	// initialize the decayer
	decayer->init();
	// return the decayer
	return decayer;
	}

	string ThreeBodyDecayConstructor::
	DecayerClassName(tcPDPtr incoming, const OrderedParticles & outgoing,
	string & objname) const {
	string classname("Herwig::");
	// spins of the outgoing particles
	unsigned int ns(0),nf(0),nv(0);
	objname += incoming->PDGName() + "2";
	for(OrderedParticles::const_iterator it=outgoing.begin();
	it!=outgoing.end();++it) {
	if ((**it).iSpin()==PDT::Spin0 ) ++ns;
	else if((**it).iSpin()==PDT::Spin1Half) ++nf;
	else if((**it).iSpin()==PDT::Spin1 ) ++nv;
	objname += (**it).PDGName();
	}
	objname += "Decayer";
	if(incoming->iSpin()==PDT::Spin0) {
	if(ns==1&&nf==2) classname += "StoSFFDecayer";
	else if(nf==2&&nv==1) classname += "StoFFVDecayer";
	else classname = "";
	}
	else if(incoming->iSpin()==PDT::Spin1Half) {
	if(nf==3) classname += "FtoFFFDecayer";
	else if(nf==1&&nv==2) classname += "FtoFVVDecayer";
	else classname = "";
	}
	else if(incoming->iSpin()==PDT::Spin1) {
	if(nf==2&&nv==1) classname += "VtoFFVDecayer";
	else classname = "";
	}
	else {
	classname="";
	}
	return classname;
	}

	void ThreeBodyDecayConstructor::
	-createDecayMode(vector<PrototypeVertexPtr> & mode) {
	+createDecayMode(vector<NBDiagram> & mode,
	+ bool possibleOnShell,
	+ double symfac) {
	// convert the diagrams from the general to the three body structure
	vector<TBDiagram> diagrams;
	- bool possibleOnShell=false;
	for(unsigned int iy=0;iy<mode.size();++iy) {
	diagrams.push_back(TBDiagram(mode[iy]));
	// remove weak processes simulated using the weak current
	if(weakMassCut_>ZERO && diagrams.back().intermediate &&
	abs(diagrams.back().intermediate->id())==ParticleID::Wplus) {
	Energy deltaM =
	getParticleData(diagrams.back().incoming)->mass() -
	getParticleData(diagrams.back().outgoing)->mass();
	if(deltaM<weakMassCut_) diagrams.pop_back();
	}
	- possibleOnShell \|= mode[iy]->possibleOnShell;
	}
	if(diagrams.empty()) return;
	// check if possible on-shell internal particles
	bool inter = interOpt_ == 1 \|\| (interOpt_ == 0 && possibleOnShell);
	// incoming particle
	tPDPtr inpart = getParticleData(diagrams[0].incoming);
	// outgoing particles
	OrderedParticles outgoing;
	outgoing.insert(getParticleData(diagrams[0].outgoing));
	outgoing.insert(getParticleData(diagrams[0].outgoingPair.first ));
	outgoing.insert(getParticleData(diagrams[0].outgoingPair.second));
	// incoming particle is now unstable
	inpart->stable(false);
	// construct the tag for the decay mode
	string tag = inpart->name() + "->";
	unsigned int iprod=0;
	for(OrderedParticles::const_iterator it = outgoing.begin();
	it != outgoing.end(); ++it) {
	++iprod;
	tag += (**it).name();
	if(iprod != 3) tag += ",";
	}
	tag += ";";
	tDMPtr dm = generator()->findDecayMode(tag);
	if( decayConstructor()->disableDecayMode(tag) ) {
	// If mode alread exists, ie has been read from file,
	// disable it
	if( dm ) {
	generator()->preinitInterface(dm, "BranchingRatio", "set", "0.0");
	generator()->preinitInterface(dm, "OnOff", "set", "Off");
	}
	return;
	}
	// create mode if needed
	if( createDecayModes() && (!dm \|\| inpart->id() == ParticleID::h0) ) {
	// create the decayer
	- GeneralThreeBodyDecayerPtr decayer = createDecayer(diagrams,inter);
	+ GeneralThreeBodyDecayerPtr decayer = createDecayer(diagrams,inter,symfac);
	if(!decayer) {
	if(Debug::level > 1 ) generator()->log() << "Can't create the decayer for "
	<< tag << " so mode not created\n";
	return;
	}
	tDMPtr ndm = generator()->preinitCreateDecayMode(tag);
	if(ndm) {
	generator()->preinitInterface(ndm, "Decayer", "set",
	decayer->fullName());
	generator()->preinitInterface(ndm, "OnOff", "set", "On");
	OrderedParticles::const_iterator pit=outgoing.begin();
	tPDPtr pa = *pit; ++pit;
	tPDPtr pb = *pit; ++pit;
	tPDPtr pc = *pit;
	Energy width =
	decayer->partialWidth(make_pair(inpart,inpart->mass()),
	make_pair(pa,pa->mass()) ,
	make_pair(pb,pb->mass()) ,
	make_pair(pc,pc->mass()));
	setBranchingRatio(ndm, width);
	}
	else
	throw NBodyDecayConstructorError()
	<< "ThreeBodyDecayConstructor::createDecayMode - Needed to create "
	<< "new decaymode but one could not be created for the tag "
	<< tag << Exception::warning;
	}
	else if( dm ) {
	if((dm->decayer()->fullName()).find("Mambo") != string::npos) {
	// create the decayer
	- GeneralThreeBodyDecayerPtr decayer = createDecayer(diagrams,inter);
	+ GeneralThreeBodyDecayerPtr decayer = createDecayer(diagrams,inter,symfac);
	if(!decayer) {
	if(Debug::level > 1 ) generator()->log() << "Can't create the decayer for "
	<< tag << " so mode not created\n";
	return;
	}
	generator()->preinitInterface(dm, "Decayer", "set",
	decayer->fullName());
	}
	}
	//update CC mode if it exists
	if( inpart->CC() ) inpart->CC()->synchronize();
	}
	diff --git a/Models/General/ThreeBodyDecayConstructor.h b/Models/General/ThreeBodyDecayConstructor.h
	--- a/Models/General/ThreeBodyDecayConstructor.h
	+++ b/Models/General/ThreeBodyDecayConstructor.h
	@@ -1,181 +1,183 @@
	// -- C++ --
	#ifndef HERWIG_ThreeBodyDecayConstructor_H
	#define HERWIG_ThreeBodyDecayConstructor_H
	//
	// This is the declaration of the ThreeBodyDecayConstructor class.
	//

	#include "NBodyDecayConstructorBase.h"
	#include "ThePEG/Helicity/Vertex/VertexBase.h"
	#include "TBDiagram.h"
	#include "PrototypeVertex.h"
	#include "Herwig++/Decay/General/GeneralThreeBodyDecayer.fh"

	namespace Herwig {
	using namespace ThePEG;

	using Helicity::VertexBasePtr;

	/**
	* The ThreeBodyDecayConstructor class inherits from the dummy base class
	* NBodyDecayConstructorBase and implements the necessary functions in
	* order to create the 3 body decaymodes for a given set of vertices
	* stored in a Model class.
	*
	* @see \ref ThreeBodyDecayConstructorInterfaces "The interfaces"
	* defined for ThreeBodyDecayConstructor.
	* @see NBodyDecayConstructor
	*/
	class ThreeBodyDecayConstructor: public NBodyDecayConstructorBase {

	public:

	/**
	* The default constructor.
	*/
	ThreeBodyDecayConstructor() :
	interOpt_(0), widthOpt_(1), weakMassCut_(-GeV) {}

	/**
	* Function used to determine allowed decaymodes, to be implemented
	* in derived class.
	*@param part vector of ParticleData pointers containing particles in model
	*/
	virtual void DecayList(const set<PDPtr> & part);

	/**
	* Number of outgoing lines. Required for correct ordering.
	*/
	virtual unsigned int numBodies() const { return 3; }

	protected:

	/**
	* Create the decayer
	* @param diagrams The diagrams for the decay
	* @param inter Option for intermediates
	*/
	GeneralThreeBodyDecayerPtr createDecayer(vector<TBDiagram> & diagrams,
	- bool inter) const;
	+ bool inter,double symfac) const;

	/**
	* Contruct the classname and object name for the Decayer
	* @param incoming The incoming particle
	* @param outgoing The decay products
	* @param objname a string containing the default path of the Decayer object
	*/
	string DecayerClassName(tcPDPtr incoming, const OrderedParticles & outgoing,
	string & objname) const;

	/**
	* Create the DecayMode from the diagrams
	* @param diagrams The diagrams
	* @param inter Option for intermediates
	*/
	- virtual void createDecayMode(vector<PrototypeVertexPtr> & mode);
	+ virtual void createDecayMode(vector<NBDiagram> & mode,
	+ bool possibleOnShell,
	+ double symfac);

	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 static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<ThreeBodyDecayConstructor> initThreeBodyDecayConstructor;

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

	private:

	/**
	* Option for the inclusion of intermediates
	*/
	unsigned int interOpt_;

	/**
	* How to treat the widths of the intermediate particles
	*/
	unsigned int widthOpt_;

	/**
	* Cut off or decays via the weak current
	*/
	Energy weakMassCut_;
	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of ThreeBodyDecayConstructor. */
	template <>
	struct BaseClassTrait<Herwig::ThreeBodyDecayConstructor,1> {
	/** Typedef of the first base class of ThreeBodyDecayConstructor. */
	typedef Herwig::NBodyDecayConstructorBase NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the ThreeBodyDecayConstructor class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::ThreeBodyDecayConstructor>
	: public ClassTraitsBase<Herwig::ThreeBodyDecayConstructor> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::ThreeBodyDecayConstructor"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_ThreeBodyDecayConstructor_H */
	diff --git a/Models/Susy/SusyBase.cc b/Models/Susy/SusyBase.cc
	--- a/Models/Susy/SusyBase.cc
	+++ b/Models/Susy/SusyBase.cc
	@@ -1,912 +1,922 @@
	// -- C++ --
	//
	// SusyBase.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 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 SusyBase class.
	//

	#include "SusyBase.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Repository/Repository.h"
	#include "ThePEG/Utilities/StringUtils.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/PDT/MassGenerator.h"
	#include "ThePEG/PDT/WidthGenerator.h"
	#include "ThePEG/PDT/DecayMode.h"

	using namespace Herwig;

	SusyBase::SusyBase() : readFile_(false), topModesFromFile_(false),
	tolerance_(1e-6), MPlanck_(2.4e18*GeV),
	gravitino_(false),
	tanBeta_(0), mu_(ZERO),
	M1_(ZERO), M2_(ZERO), M3_(ZERO),
	mH12_(ZERO),mH22_(ZERO),
	meL_(ZERO),mmuL_(ZERO),mtauL_(ZERO),
	meR_(ZERO),mmuR_(ZERO),mtauR_(ZERO),
	mq1L_(ZERO),mq2L_(ZERO),mq3L_(ZERO),
	mdR_(ZERO),muR_(ZERO),msR_(ZERO),
	mcR_(ZERO),mbR_(ZERO),mtR_(ZERO),
	gluinoPhase_(1.)
	{}

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

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

	void SusyBase::doinit() {
	addVertex(WSFSFVertex_);
	addVertex(NFSFVertex_);
	addVertex(GFSFVertex_);
	addVertex(HSFSFVertex_);
	addVertex(CFSFVertex_);
	addVertex(GSFSFVertex_);
	addVertex(GGSQSQVertex_);
	addVertex(GSGSGVertex_);
	addVertex(NNZVertex_);
	if(NNPVertex_) addVertex(NNPVertex_);
	if(GNGVertex_) addVertex(GNGVertex_);
	addVertex(CCZVertex_);
	addVertex(CNWVertex_);
	addVertex(GOGOHVertex_);
	addVertex(WHHVertex_);
	addVertex(NCTVertex_);
	if(gravitino_) {
	if(GVNHVertex_) addVertex(GVNHVertex_);
	if(GVNVVertex_) addVertex(GVNVVertex_);
	if(GVFSVertex_) addVertex(GVFSVertex_);
	}
	StandardModel::doinit();
	}

	void SusyBase::persistentOutput(PersistentOStream & os) const {
	os << readFile_ << topModesFromFile_ << gravitino_
	<< NMix_ << UMix_ << VMix_ << WSFSFVertex_
	<< NFSFVertex_ << GFSFVertex_ << HSFSFVertex_ << CFSFVertex_
	<< GSFSFVertex_ << GGSQSQVertex_ << GSGSGVertex_
	<< NNZVertex_ << NNPVertex_ << CCZVertex_ << CNWVertex_
	<< GOGOHVertex_ << WHHVertex_ << GNGVertex_ << NCTVertex_
	<< GVNHVertex_ << GVNVVertex_ << GVFSVertex_
	<< tanBeta_ << ounit(mu_,GeV)
	<< ounit(M1_,GeV) << ounit(M2_,GeV) << ounit(M3_,GeV)
	<< ounit(mH12_,GeV2) << ounit(mH22_,GeV2)
	<< ounit(meL_,GeV) << ounit(mmuL_,GeV) << ounit(mtauL_,GeV)
	<< ounit(meR_,GeV) << ounit(mmuR_,GeV) << ounit(mtauR_,GeV)
	<< ounit(mq1L_,GeV) << ounit(mq2L_,GeV) << ounit(mq3L_,GeV)
	<< ounit(mdR_,GeV) << ounit(muR_,GeV) << ounit(msR_,GeV)
	<< ounit(mcR_,GeV) << ounit(mbR_,GeV) << ounit(mtR_,GeV)
	<< gluinoPhase_ << tolerance_ << ounit(MPlanck_,GeV);
	}

	void SusyBase::persistentInput(PersistentIStream & is, int) {
	is >> readFile_ >> topModesFromFile_ >> gravitino_
	>> NMix_ >> UMix_ >> VMix_ >> WSFSFVertex_
	>> NFSFVertex_ >> GFSFVertex_ >> HSFSFVertex_ >> CFSFVertex_
	>> GSFSFVertex_ >> GGSQSQVertex_ >> GSGSGVertex_
	>> NNZVertex_ >> NNPVertex_ >> CCZVertex_ >> CNWVertex_
	>> GOGOHVertex_ >> WHHVertex_ >> GNGVertex_ >> NCTVertex_
	>> GVNHVertex_ >> GVNVVertex_ >> GVFSVertex_
	>> tanBeta_ >> iunit(mu_,GeV)
	>> iunit(M1_,GeV) >> iunit(M2_,GeV) >> iunit(M3_,GeV)
	>> iunit(mH12_,GeV2) >> iunit(mH22_,GeV2)
	>> iunit(meL_,GeV) >> iunit(mmuL_,GeV) >> iunit(mtauL_,GeV)
	>> iunit(meR_,GeV) >> iunit(mmuR_,GeV) >> iunit(mtauR_,GeV)
	>> iunit(mq1L_,GeV) >> iunit(mq2L_,GeV) >> iunit(mq3L_,GeV)
	>> iunit(mdR_,GeV) >> iunit(muR_,GeV) >> iunit(msR_,GeV)
	>> iunit(mcR_,GeV) >> iunit(mbR_,GeV) >> iunit(mtR_,GeV)
	>> gluinoPhase_ >> tolerance_ >> iunit(MPlanck_,GeV);
	}

	ClassDescription<SusyBase> SusyBase::initSusyBase;
	// Definition of the static class description member.

	void SusyBase::Init() {

	static ClassDocumentation<SusyBase> documentation
	("This is the base class for any SUSY model.",
	"SUSY spectrum files follow the Les Houches accord"
	" \\cite{Skands:2003cj,Allanach:2008qq}.",
	" %\\cite{Skands:2003cj}\n"
	"\\bibitem{Skands:2003cj}\n"
	" P.~Skands {\\it et al.},\n"
	" ``SUSY Les Houches accord: Interfacing SUSY spectrum calculators, decay\n"
	" %packages, and event generators,''\n"
	" JHEP {\\bf 0407}, 036 (2004)\n"
	" [arXiv:hep-ph/0311123].\n"
	" %%CITATION = JHEPA,0407,036;%%\n"
	"%\\cite{Allanach:2008qq}\n"
	"\\bibitem{Allanach:2008qq}\n"
	" B.~Allanach {\\it et al.},\n"
	" %``SUSY Les Houches Accord 2,''\n"
	" Comput.\\ Phys.\\ Commun.\\ {\\bf 180}, 8 (2009)\n"
	" [arXiv:0801.0045 [hep-ph]].\n"
	" %%CITATION = CPHCB,180,8;%%\n"
	);

	static Switch<SusyBase,bool> interfaceTopModes
	("TopModes",
	"Whether ro use the Herwig++ SM top decays or those from the SLHA file",
	&SusyBase::topModesFromFile_, false, false, false);
	static SwitchOption interfaceTopModesFile
	(interfaceTopModes,
	"File",
	"Take the modes from the files",
	true);
	static SwitchOption interfaceTopModesHerwig
	(interfaceTopModes,
	"Herwig",
	"Use the SM ones",
	false);

	static Reference<SusyBase,Helicity::AbstractVSSVertex> interfaceVertexWSS
	("Vertex/WSFSF",
	"Reference to Susy W SF SF vertex",
	&SusyBase::WSFSFVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractFFSVertex> interfaceVertexNFSF
	("Vertex/NFSF",
	"Reference to the neutralino-fermion-sfermion vertex",
	&SusyBase::NFSFVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractFFSVertex> interfaceVertexGFSF
	("Vertex/GFSF",
	"Reference to the gluino-fermion-sfermion vertex",
	&SusyBase::GFSFVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractSSSVertex> interfaceVertexHSFSF
	("Vertex/HSFSF",
	"Reference to the Higgs-fermion-sfermion vertex",
	&SusyBase::HSFSFVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractFFSVertex> interfaceVertexCFSF
	("Vertex/CFSF",
	"Reference to the chargino-fermion-sfermion vertex",
	&SusyBase::CFSFVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractVSSVertex> interfaceVertexGSFSF
	("Vertex/GSFSF",
	"Reference to the gluon-sfermion-sfermion vertex",
	&SusyBase::GSFSFVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractVVSSVertex> interfaceVertexGGSS
	("Vertex/GGSQSQ",
	"Reference to the gluon-gluon-squark-squark vertex",
	&SusyBase::GGSQSQVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractFFVVertex> interfaceVertexGSGSG
	("Vertex/GSGSG",
	"Reference to the gluon-gluino-gluino vertex",
	&SusyBase::GSGSGVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractFFVVertex> interfaceVertexNNZ
	("Vertex/NNZ",
	"Reference to Z-~chi_i0-~chi_i0 vertex",
	&SusyBase::NNZVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractFFVVertex> interfaceVertexNNP
	("Vertex/NNP",
	"Reference to photon-~chi_i0-~chi_i0 vertex",
	&SusyBase::NNPVertex_, false, false, true, true, false);

	static Reference<SusyBase,Helicity::AbstractFFVVertex> interfaceVertexGNG
	("Vertex/GNG",
	"Reference to gluon-~chi_i0-gluino vertex",
	&SusyBase::GNGVertex_, false, false, true, true, false);

	static Reference<SusyBase,Helicity::AbstractFFVVertex> interfaceVertexCCZ
	("Vertex/CCZ",
	"Reference to ~chi_i+-~chi_i-Z vertex",
	&SusyBase::CCZVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractFFVVertex> interfaceVertexCNW
	("Vertex/CNW",
	"Reference to ~chi_i+-chi_i0-W vertex",
	&SusyBase::CNWVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractFFSVertex> interfaceVertexGOGOH
	("Vertex/GOGOH",
	"Reference to the gaugino-gaugino-higgs vertex",
	&SusyBase::GOGOHVertex_, false, false, true, false);

	static Reference<SusyBase,Helicity::AbstractVSSVertex> interfaceVertexWHH
	("Vertex/SSWHH",
	"Reference to Susy WHHVertex",
	&SusyBase::WHHVertex_, false, false, true, false);

	static Reference<SusyBase,AbstractFFSVertex> interfaceVertexNCT
	("Vertex/NCT",
	"Vertex for the flavour violating coupling of the top squark "
	"to the neutralino and charm quark.",
	&SusyBase::NCTVertex_, false, false, true, false, false);

	static Reference<SusyBase,AbstractRFSVertex> interfaceVertexGVNH
	("Vertex/GVNH",
	"Vertex for the interfaction of the gravitino-neutralino"
	" and Higgs bosons",
	&SusyBase::GVNHVertex_, false, false, true, true, false);

	static Reference<SusyBase,AbstractRFVVertex> interfaceVertexGVNV
	("Vertex/GVNV",
	"Vertex for the interfaction of the gravitino-neutralino"
	" and vector bosons",
	&SusyBase::GVNVVertex_, false, false, true, true, false);

	static Reference<SusyBase,AbstractRFSVertex> interfaceVertexGVFS
	("Vertex/GVFS",
	"Vertex for the interfaction of the gravitino-fermion"
	" and sfermion",
	&SusyBase::GVFSVertex_, false, false, true, true, false);

	static Parameter<SusyBase,double> interfaceBRTolerance
	("BRTolerance",
	"Tolerance for the sum of branching ratios to be difference from one.",
	&SusyBase::tolerance_, 1e-6, 1e-8, 0.01,
	false, false, Interface::limited);

	static Parameter<SusyBase,Energy> interfaceMPlanck
	("MPlanck",
	"The Planck mass for GMSB models",
	&SusyBase::MPlanck_, GeV, 2.4e18GeV, 1.e16GeV, 1.e20*GeV,
	false, false, Interface::limited);

	}

	void SusyBase::readSetup(istream & is) {
	string filename = dynamic_ptr_cast<istringstream*>(&is)->str();
	if(readFile_)
	throw SetupException()
	<< "A second SLHA file " << filename << " has been opened."
	<< "This is probably unintended and as it can cause crashes"
	<< " and other unpredictable behaviour it is not allowed."
	<< Exception::runerror;
	CFileLineReader cfile;
	cfile.open(filename);
	if( !cfile ) throw SetupException()
	<< "SusyBase::readSetup - An error occurred in opening the "
	<< "spectrum file \"" << filename << "\". A SUSY model cannot be "
	<< "run without this."
	<< Exception::runerror;
	useMe();
	//Before reading the spectrum/decay files the SM higgs
	//decay modes, mass and width generators need to be turned off.
	PDPtr h0 = getParticleData(ParticleID::h0);
	h0->widthGenerator(WidthGeneratorPtr());
	h0->massGenerator(MassGenPtr());
	h0->width(ZERO);
	DecaySet::const_iterator dit = h0->decayModes().begin();
	DecaySet::const_iterator dend = h0->decayModes().end();
	for( ; dit != dend; ++dit ) {
	const InterfaceBase * ifb =
	BaseRepository::FindInterface(*dit, "BranchingRatio");
	ifb->exec(**dit, "set", "0.0");
	ifb = BaseRepository::FindInterface(*dit, "OnOff");
	ifb->exec(**dit, "set", "Off");
	}
	// if taking the top modes from the file
	// delete the SM stuff
	if(topModesFromFile_) {
	PDPtr top = getParticleData(ParticleID::t);
	top->widthGenerator(WidthGeneratorPtr());
	top->massGenerator(MassGenPtr());
	DecaySet::const_iterator dit = top->decayModes().begin();
	DecaySet::const_iterator dend = top->decayModes().end();
	for( ; dit != dend; ++dit ) {
	const InterfaceBase * ifb =
	BaseRepository::FindInterface(*dit, "BranchingRatio");
	ifb->exec(**dit, "set", "0.0");
	ifb = BaseRepository::FindInterface(*dit, "OnOff");
	ifb->exec(**dit, "set", "Off");
	}
	}
	// read first line and check if this is a Les Houches event file
	cfile.readline();
	bool lesHouches = cfile.find("<LesHouchesEvents");
	bool reading = !lesHouches;
	if(lesHouches) cfile.readline();
	//function pointer for putting all characters to lower case.
	int (*pf)(int) = tolower;
	do {
	string line = cfile.getline();
	// check for start of slha block in SLHA files
	if(lesHouches && !reading) {
	if(line.find("<slha")==0) reading = true;
	if(!cfile.readline()) break;
	continue;
	}
	// ignore comment lines
	if(line[0] == '#') {
	if(!cfile.readline()) break;
	continue;
	}
	// make everything lower case
	transform(line.begin(), line.end(), line.begin(), pf);
	// start of a block
	if(line.find("block") == 0) {
	string name = StringUtils::car(StringUtils::cdr(line), " #");
	name = StringUtils::stripws(name);
	// mixing matrix
	if((name.find("mix") != string::npos &&
	name.find("hmix") != 0)) {
	unsigned int row(0),col(0);
	MixingVector vals = readMatrix(cfile,row,col);
	mixings_[name] = make_pair(make_pair(row,col),vals);
	}
	else if(name.find("au") == 0 \|\| name.find("ad") == 0 \|\|
	name.find("ae") == 0 ) {
	string test = StringUtils::car(line, "#");
	while (test.find("=")!= string::npos) {
	test = StringUtils::cdr(test, "=");
	}
	istringstream is(test);
	double scale;
	is >> scale;
	unsigned int row(0),col(0);
	MixingVector vals = readMatrix(cfile,row,col);
	if(scale>1e10) continue;
	mixings_[name] = make_pair(make_pair(row,col),vals);
	}
	else if( name.find("info") == string::npos) {
	readBlock(cfile,name,line);
	}
	else {
	if(!cfile.readline()) break;
	}
	continue;
	}
	// decays
	else if( line.find("decay") == 0 ) {
	readDecay(cfile, line);
	continue;
	}
	else if( lesHouches && line.find("</slha") == 0 ) {
	reading = false;
	break;
	}
	if(!cfile.readline()) break;
	}
	while(true);
	// extract the relevant parameters
	extractParameters();
	// create the mixing matrices we need
	createMixingMatrices();
	// set the masses, this has to be done after the
	// mixing matrices have been created
	resetRepositoryMasses();
	// have now read the file
	readFile_=true;
	}

	void SusyBase::readBlock(CFileLineReader & cfile,string name,string linein) {
	if(!cfile)
	throw SetupException()
	<< "SusyBase::readBlock() - The input stream is in a bad state"
	<< Exception::runerror;
	// storage or the parameters
	string test = StringUtils::car(linein, "#");
	ParamMap store;
	bool set = true;
	// special for the alpha block
	if(name.find("alpha") == 0 ) {
	double alpha;
	cfile.readline();
	string line = cfile.getline();
	istringstream iss(line);
	iss >> alpha;
	store.insert(make_pair(1,alpha));
	}
	else {
	// extract the scale from the block if present
	if(test.find("=")!= string::npos) {
	while(test.find("=")!=string::npos)
	test= StringUtils::cdr(test,"=");
	istringstream is(test);
	double scale;
	is >> scale;
	// only store the lowest scale block
	if(parameters_.find(name)!=parameters_.end()) {
	set = scale < parameters_[name][-1];
	}
	else {
	store.insert(make_pair(-1,scale));
	}
	}
	while(cfile.readline()) {
	string line = cfile.getline();
	// skip comments
	if(line[0] == '#') continue;
	// reached the end
	if( line[0] == 'B' \|\| line[0] == 'b' \|\|
	line[0] == 'D' \|\| line[0] == 'd' \|\|
	line[0] == '<' ) {
	cfile.resetline();
	break;
	}
	istringstream is(line);
	long index;
	double value;
	if(name.find("rvlam")!= string::npos\|\|
	name=="rvt" \|\| name=="rvtp" \|\| name=="rvtpp") {
	int i,j,k;
	is >> i >> j >> k >> value;
	index = i100+j10+k;
	}
	else {
	is >> index >> value;
	}
	store.insert(make_pair(index, value));
	}
	}
	if(set) parameters_[name]=store;
	}

	void SusyBase::readDecay(CFileLineReader & cfile,
	string decay) const{
	map<string,ParamMap>::const_iterator fit=parameters_.find("mass");
	if(fit==parameters_.end())
	throw Exception() << "SusyBase::readDecay() "
	<< "BLOCK MASS not found in input file"
	<< "it must be before the decays so the kinematics can be checked"
	<< Exception::runerror;
	ParamMap theMasses = fit->second;
	if(!cfile)
	throw SetupException()
	<<"SusyBase::readDecay - The input stream is in a bad state";
	// extract parent PDG code and width
	long parent(0);
	Energy width(ZERO);
	istringstream iss(decay);
	string dummy;
	iss >> dummy >> parent >> iunit(width, GeV);
	PDPtr inpart = getParticleData(parent);
	if(!topModesFromFile_&&abs(parent)==ParticleID::t) {
	cfile.readline();
	return;
	}
	if(!inpart) throw SetupException()
	<< "SusyBase::readDecay() - "
	<< "A ParticleData object with the PDG code "
	<< parent << " does not exist. "
	<< Exception::runerror;
	inpart->width(width);
	if( width > ZERO ) inpart->cTau(hbarc/width);
	inpart->widthCut(5.*width);
	ParamMap::iterator it = theMasses.find(abs(inpart->id()));
	Energy inMass = it!=theMasses.end() ? abs(it->second*GeV) : inpart->mass();
	string prefix("decaymode " + inpart->name() + "->");
	double brsum(0.);
	unsigned int nmode = 0;
	while(cfile.readline()) {
	string line = cfile.getline();
	// skip comments
	if(line[0] == '#') continue;
	// reached the end
	if( line[0] == 'B' \|\| line[0] == 'b' \|\|
	line[0] == 'D' \|\| line[0] == 'd' \|\|
	line[0] == '<' ) {
	cfile.resetline();
	break;
	}
	// read the mode
	// get the branching ratio and no of decay products
	istringstream is(line);
	double brat(0.);
	unsigned int nda(0),npr(0);
	is >> brat >> nda;
	vector<tcPDPtr> products,bosons;
	Energy mout(ZERO),moutnoWZ(ZERO);
	string tag = prefix;
	while( true ) {
	long t;
	is >> t;
	if( is.fail() ) break;
	if( t == abs(parent) ) {
	throw SetupException()
	<< "An error occurred while read a decay of the "
	<< inpart->PDGName() << ". One of its products has the same PDG code "
	<< "as the parent particle. Please check the SLHA file.\n"
	<< Exception::runerror;
	}
	tcPDPtr p = getParticleData(t);
	if( !p ) {
	throw SetupException()
	<< "SusyBase::readDecay() - An unknown PDG code has been encounterd "
	<< "while reading a decay mode. ID: " << t
	<< Exception::runerror;
	}
	++npr;
	tag += p->name() + ",";
	ParamMap::iterator it = theMasses.find(abs(p->id()));
	Energy mass = it!=theMasses.end() ? abs(it->second*GeV) : p->mass();
	mout += mass;
	if(abs(p->id())==ParticleID::Wplus\|\|p->id()==ParticleID::Z0) {
	bosons.push_back(p);
	}
	else {
	products.push_back(p);
	moutnoWZ += mass;
	}
	}
	if( npr != nda ) {
	throw SetupException()
	<< "SusyBase::readDecay - While reading a decay of the "
	<< inpart->PDGName() << " from an SLHA file, an inconsistency "
	<< "between the number of decay products and the value in "
	<< "the 'NDA' column was found. Please check if the spectrum "
	<< "file is correct.\n"
	<< Exception::warning;
	}
	if( npr > 1 ) {
	tag.replace(tag.size() - 1, 1, ";");
	// normal option
	if(mout<=inMass) {
	inpart->stable(false);
	brsum += brat;
	createDecayMode(tag, brat);
	}
	// no possible off-shell gauge bosons throw it away
	else if(bosons.empty() \|\| bosons.size()>2 \|\|
	moutnoWZ>=inMass) {
	cerr << "SusyBase::readDecay() "
	<< "The decay " << tag << " cannot proceed for on-shell "
	<< "particles, skipping it.\n";
	}
	else {
	Energy maxMass = inMass - moutnoWZ;
	string newTag = prefix;
	for(unsigned int ix=0;ix<products.size();++ix)
	newTag += products[ix]->name() + ",";
	if(bosons.size()==1) {
	cerr << "SusyBase::readDecay() "
	<< "The decay " << tag << " cannot proceed for on-shell\n"
	<< "particles, replacing gauge boson with its decay products\n";
	vector<pair<double,string> > modes =
	createWZDecayModes(newTag,brat,bosons[0],maxMass);
	for(unsigned int ix=0;ix<modes.size();++ix) {
	modes[ix].second.replace(modes[ix].second.size() - 1, 1, ";");
	createDecayMode(modes[ix].second,modes[ix].first);
	brsum += modes[ix].first;
	}
	}
	else if(bosons.size()==2) {
	bool identical = bosons[0]->id()==bosons[1]->id();
	if(maxMass>bosons[0]->mass()&&maxMass>bosons[1]->mass()) {
	cerr << "SusyBase::readDecay() "
	<< "The decay " << tag << " cannot proceed for on-shell\n"
	<< "particles, replacing one of the gauge bosons"
	<< " with its decay products\n";
	unsigned int imax = identical ? 1 : 2;
	if(imax==2) brat *= 0.5;
	for(unsigned int ix=0;ix<imax;++ix) {
	string newTag2 = newTag+bosons[ix]->name()+',';
	vector<pair<double,string> > modes =
	createWZDecayModes(newTag2,brat,bosons[ix],maxMass);
	for(unsigned int ix=0;ix<modes.size();++ix) {
	modes[ix].second.replace(modes[ix].second.size() - 1, 1, ";");
	createDecayMode(modes[ix].second,modes[ix].first);
	brsum += modes[ix].first;
	}
	}
	}
	else {
	cerr << "SusyBase::readDecay() "
	<< "The decay " << tag << " cannot proceed for on-shell\n"
	<< "particles, and has too many off-shell gauge bosons,"
	<< " skipping it.\n";
	}
	}
	else {
	cerr << "SusyBase::readDecay() "
	<< "The decay " << tag << " cannot proceed for on-shell\n"
	<< "particles, and has too many outgoing gauge bosons skipping it.\n";
	}
	}
	}
	}
	if( abs(brsum - 1.) > tolerance_ && nmode!=0 ) {
	cerr << "Warning: The total branching ratio for " << inpart->PDGName()
	<< " from the spectrum file does not sum to 1. The branching fractions"
	<< " will be rescaled.\n";
	cerr << setprecision(13) << abs(brsum - 1.) << "\n";
	}
	}

	const MixingVector
	SusyBase::readMatrix(CFileLineReader & cfile,
	unsigned int & row, unsigned int & col) {
	if(!cfile)
	throw SetupException()
	<< "SusyBase::readMatrix() - The input stream is in a bad state."
	<< Exception::runerror;
	unsigned int rowmax(0), colmax(0);
	MixingVector values;
	while(cfile.readline()) {
	string line = cfile.getline();
	// skip comments
	if(line[0] == '#') continue;
	// reached the end
	if( line[0] == 'B' \|\| line[0] == 'b' \|\|
	line[0] == 'D' \|\| line[0] == 'd' \|\|
	line[0] == '<' ) {
	cfile.resetline();
	break;
	}
	istringstream is(line);
	unsigned int index1, index2;
	double real(0.), imag(0.);
	is >> index1 >> index2 >> real >> imag;
	values.push_back(MixingElement(index1,index2,Complex(real, imag)));
	if(index1 > rowmax) rowmax = index1;
	if(index2 > colmax) colmax = index2;

	}
	col=colmax;
	row=rowmax;
	return values;
	}

	void SusyBase::createDecayMode(string tag, double brat) const {
	ostringstream cmd;
	cmd << tag << string(" ")
	<< setprecision(13) << brat << string(" 1 /Herwig/Decays/Mambo");
	Repository::exec(cmd.str(), cerr);
	}

	void SusyBase::createMixingMatrix(MixingMatrixPtr & matrix,
	string name, const MixingVector & values,
	MatrixSize size) {
	matrix = new_ptr(MixingMatrix(size.first,size.second));
	for(unsigned int ix=0; ix < values.size(); ++ix)
	(*matrix)(values[ix].row-1,values[ix].col-1) = values[ix].value;

	if(name == "nmix") {
	vector<long> ids(4);
	ids[0] = 1000022; ids[1] = 1000023;
	ids[2] = 1000025; ids[3] = 1000035;
	matrix->setIds(ids);
	}
	else if(name == "nmnmix") {
	vector<long> ids(5);
	ids[0] = 1000022; ids[1] = 1000023;
	ids[2] = 1000025; ids[3] = 1000035;
	ids[4] = 1000045;
	matrix->setIds(ids);
	}
	else if(name == "umix") {
	vector<long> ids(2);
	ids[0] = 1000024; ids[1] = 1000037;
	matrix->setIds(ids);
	}
	else if(name == "vmix") {
	vector<long> ids(2);
	ids[0] = 1000024; ids[1] = 1000037;
	matrix->setIds(ids);
	}
	else if(name == "stopmix") {
	vector<long> ids(2);
	ids[0] = 1000006; ids[1] = 2000006;
	matrix->setIds(ids);
	}
	else if(name == "sbotmix") {
	vector<long> ids(2);
	ids[0] = 1000005; ids[1] = 2000005;
	matrix->setIds(ids);
	}
	else if(name == "staumix") {
	vector<long> ids(2);
	ids[0] = 1000015; ids[1] = 2000015;
	matrix->setIds(ids);
	}
	else if(name == "nmhmix") {
	vector<long> ids(3);
	ids[0] = 25; ids[1] = 35; ids[2] = 45;
	matrix->setIds(ids);
	}
	else if(name == "nmamix") {
	vector<long> ids(2);
	ids[0] = 36; ids[1] = 46;
	matrix->setIds(ids);
	}
	else
	throw SetupException() << "Cannot find correct title for mixing matrix "
	<< name << Exception::runerror;
	}

	void SusyBase::resetRepositoryMasses() {
	map<string,ParamMap>::const_iterator fit=parameters_.find("mass");
	if(fit==parameters_.end())
	throw Exception() << "BLOCK MASS not found in input file"
	<< " can't set masses of SUSY particles"
	<< Exception::runerror;
	ParamMap theMasses = fit->second;
	for(ParamMap::iterator it = theMasses.begin(); it != theMasses.end();
	++it) {
	long id = it->first;
	double mass = it->second;
	//a negative mass requires an adjustment to the
	//associated mixing matrix by a factor of i
	if(mass < 0.0) adjustMixingMatrix(id);
	PDPtr part = getParticleData(id);
	if(!part) throw SetupException()
	<< "SusyBase::resetRepositoryMasses() - Particle with PDG code " << id
	<< " not found." << Exception::warning;
	//Find interface nominal mass interface
	const InterfaceBase * ifb = BaseRepository::FindInterface(part, "NominalMass");
	ostringstream os;
	os << abs(it->second);
	ifb->exec(*part, "set", os.str());
	// switch on gravitino interactions?
	gravitino_ \|= id== ParticleID::SUSY_Gravitino;
	}
	theMasses.clear();
	}

	void SusyBase::adjustMixingMatrix(long id) {
	//get correct mixing matrix
	switch(id) {
	case 1000021 :
	gluinoPhase_ = Complex(0.,1.);
	break;
	case 1000022 :
	case 1000023 :
	case 1000025 :
	case 1000035 :
	case 1000045 :
	if(NMix_)
	NMix_->adjustPhase(id);
	else
	throw SetupException() << "SusyBase::adjustMixingMatrix - "
	<< "The neutralino mixing matrix pointer "
	<< "is null!" << Exception::runerror;
	break;
	case 1000024 :
	case 1000037 :
	if(UMix_)
	UMix_->adjustPhase(id);
	else
	throw SetupException() << "SusyBase::adjustMixingMatrix - "
	<< "The U-Type chargino mixing matrix pointer "
	<< "is null!" << Exception::runerror;
	if(VMix_)
	VMix_->adjustPhase(id);
	else
	throw SetupException() << "SusyBase::adjustMixingMatrix - "
	<< "The V-Type chargino mixing matrix pointer "
	<< "is null!" << Exception::runerror;
	break;
	default :
	throw SetupException()
	<< "SusyBase::adjustMixingMatrix - Trying to adjust mixing matrix "
	<< "phase for a particle that does not have a mixing matrix "
	<< "associated with it. " << id << " must have a negative mass in "
	<< "the spectrum file, this should only occur for particles that mix."
	<< Exception::runerror;
	}
	}

	void SusyBase::createMixingMatrices() {
	map<string,pair<MatrixSize, MixingVector > >::const_iterator it;
	for(it=mixings_.begin();it!=mixings_.end();++it) {
	string name=it->first;
	// create the gaugino mixing matrices
	if(name == "nmix" \|\| name == "nmnmix" ){
	createMixingMatrix(NMix_,name,it->second.second,it->second.first);
	}
	else if (name == "umix" ) {
	createMixingMatrix(UMix_,name,it->second.second,it->second.first);
	}
	else if (name == "vmix") {
	createMixingMatrix(VMix_,name,it->second.second,it->second.first);
	}
	}
	}

	void SusyBase::extractParameters(bool checkmodel) {
	map<string,ParamMap>::const_iterator pit;
	ParamMap::const_iterator it;
	// try and get tan beta from extpar first
	pit=parameters_.find("extpar");
	// extract tan beta
	tanBeta_ = -1.;
	- if(pit!=parameters_.end()) {
	- it = pit->second.find(25);
	- if(it!=pit->second.end()) tanBeta_ = it->second;
	+ if(tanBeta_<0.) {
	+ pit=parameters_.find("hmix");
	+ if(pit!=parameters_.end()) {
	+ it = pit->second.find(2);
	+ if(it!=pit->second.end()) tanBeta_ = it->second;
	+ }
	+ }
	+ if(tanBeta_<0.) {
	+ pit=parameters_.find("extpar");
	+ if(pit!=parameters_.end()) {
	+ it = pit->second.find(25);
	+ if(it!=pit->second.end()) tanBeta_ = it->second;
	+ }
	}
	// otherwise from minpar
	if(tanBeta_<0.) {
	pit=parameters_.find("minpar");
	if(pit!=parameters_.end()) {
	it = pit->second.find(3);
	if(it!=pit->second.end()) tanBeta_ = it->second;
	}
	}
	if(tanBeta_<0.)
	throw Exception() << "SusyBase::extractParameters() "
	<< "Can't find tan beta in BLOCK MINPAR"
	<< " or BLOCK EXTPAR " << Exception::runerror;
	// extract parameters from hmix
	pit=parameters_.find("hmix");
	if(pit==parameters_.end()) {
	if(generator())
	generator()->logWarning(Exception("SusyBase::extractParameters() BLOCK HMIX not found setting mu to zero\n",
	Exception::warning));
	else
	cerr << "SusyBase::extractParameters() BLOCK HMIX not found setting mu to zero\n";
	mu_=ZERO;
	}
	else {
	mu_=findValue(pit,1,"HMIX","mu")*GeV;
	}
	pit = parameters_.find("msoft");
	if( pit == parameters_.end() )
	throw Exception() << "SusyBase::extractParameters() "
	<< "BLOCK MSOFT not found in "
	<< "SusyBase::extractParameters()"
	<< Exception::runerror;
	M1_ = findValue(pit,1 ,"MSOFT","M_1" )*GeV;
	M2_ = findValue(pit,2 ,"MSOFT","M_2" )*GeV;
	M3_ = findValue(pit,3 ,"MSOFT","M_3" )*GeV;
	mH12_ = findValue(pit,21,"MSOFT","m_H1^2")*GeV2;
	mH22_ = findValue(pit,22,"MSOFT","m_H2^2")*GeV2;
	meL_ = findValue(pit,31,"MSOFT","M_eL" )*GeV;
	mmuL_ = findValue(pit,32,"MSOFT","M_muL" )*GeV;
	mtauL_ = findValue(pit,33,"MSOFT","M_tauL")*GeV;
	meR_ = findValue(pit,34,"MSOFT","M_eR" )*GeV;
	mmuR_ = findValue(pit,35,"MSOFT","M_muR" )*GeV;
	mtauR_ = findValue(pit,36,"MSOFT","M_tauR")*GeV;
	mq1L_ = findValue(pit,41,"MSOFT","M_q1L" )*GeV;
	mq2L_ = findValue(pit,42,"MSOFT","M_q2L" )*GeV;
	mq3L_ = findValue(pit,43,"MSOFT","M_q3L" )*GeV;
	muR_ = findValue(pit,44,"MSOFT","M_uR" )*GeV;
	mcR_ = findValue(pit,45,"MSOFT","M_cR" )*GeV;
	mtR_ = findValue(pit,46,"MSOFT","M_tR" )*GeV;
	mdR_ = findValue(pit,47,"MSOFT","M_dR" )*GeV;
	msR_ = findValue(pit,48,"MSOFT","M_sR" )*GeV;
	mbR_ = findValue(pit,49,"MSOFT","M_bR" )*GeV;
	if(checkmodel) {
	throw Exception() << "The SusyBase class should not be used as a "
	<< "Model class, use one of the models which inherit"
	<< " from it" << Exception::runerror;
	}
	}

	vector<pair<double,string> >
	SusyBase::createWZDecayModes(string tag, double brat,
	tcPDPtr boson, Energy maxMass) const {
	vector<pair<double,string> > modes;
	double sum(0.);
	for(DecaySet::const_iterator dit=boson->decayModes().begin();
	dit!=boson->decayModes().end();++dit) {
	tcDMPtr mode = *dit;
	if(!mode->on()) continue;
	string extra;
	Energy outMass(ZERO);
	for(ParticleMSet::const_iterator pit=mode->products().begin();
	pit!=mode->products().end();++pit) {
	extra += (**pit).name() + ",";
	outMass += (**pit).mass();
	}
	if(outMass<maxMass) {
	sum += mode->brat();
	modes.push_back(make_pair(mode->brat(),tag+extra));
	}
	}
	for(unsigned int ix=0;ix<modes.size();++ix)
	modes[ix].first *= brat/sum;
	return modes;
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:49 PM (1 d, 7 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805923
Default Alt Text: (147 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions