No OneTemporary
Actions

Size

61 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,1010 +1,1025 @@
	// -- 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 << _intOpt << _relerr;
	}

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

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

	static Switch<GeneralThreeBodyDecayer,unsigned int> interfacePartialWidthIntegration
	("PartialWidthIntegration",
	"Switch to control the partial width integration",
	&GeneralThreeBodyDecayer::_intOpt, 0, false, false);
	static SwitchOption interfacePartialWidthIntegrationAllPoles
	(interfacePartialWidthIntegration,
	"AllPoles",
	"Include all potential poles",
	0);
	static SwitchOption interfacePartialWidthIntegrationShallowestPole
	(interfacePartialWidthIntegration,
	"ShallowestPole",
	"Only include the shallowest pole",
	1);

	static Parameter<GeneralThreeBodyDecayer,double> interfaceRelativeError
	("RelativeError",
	"The relative error for the GQ integration of the partial width",
	&GeneralThreeBodyDecayer::_relerr, 1e-2, 1e-10, 1.,
	false, false, Interface::limited);

	}

	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()) != abs(_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,
	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(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));
	int jac = 0;
	double power = 0.0;
	if ( _diagrams[ix].intermediate->mass() == ZERO ) {
	jac = 1;
	power = -2.0;
	}
	if(_diagrams[ix].channelType==TBDiagram::channel23) {
	newchannel->addIntermediate(extpart[0],0,0.0,-1,1);
	newchannel->addIntermediate(_diagrams[ix].intermediate,jac,power, 2,3);
	}
	else if(_diagrams[ix].channelType==TBDiagram::channel13) {
	newchannel->addIntermediate(extpart[0],0,0.0,-1,2);
	newchannel->addIntermediate(_diagrams[ix].intermediate,jac,power, 1,3);
	}
	else if(_diagrams[ix].channelType==TBDiagram::channel12) {
	newchannel->addIntermediate(extpart[0],0,0.0,-1,3);
	newchannel->addIntermediate(_diagrams[ix].intermediate,jac,power, 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 if (singlet.size()==1&&antitriplet.size()==2) {
	tColinePtr col[2] = {ColourLine::create(outgoing[antitriplet[0]],true),
	ColourLine::create(outgoing[antitriplet[1]],true)};
	parent.colourLine()->setSinkNeighbours(col[0],col[1]);
	}
	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 if (singlet.size()==1&&triplet.size()==2) {
	tColinePtr col[2] = {ColourLine::create(outgoing[triplet[0]]),
	ColourLine::create(outgoing[triplet[1]])};
	parent.antiColourLine()->setSourceNeighbours(col[0],col[1]);
	}
	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 if(triplet.size()==3) {
	tColinePtr col[2];
	if(colourFlow()==0) {
	outgoing[0]->incomingColour (const_ptr_cast<tPPtr>(&parent));
	col[0] = ColourLine::create(outgoing[1]);
	col[1] = ColourLine::create(outgoing[2]);
	}
	else if(colourFlow()==1) {
	outgoing[1]->incomingColour (const_ptr_cast<tPPtr>(&parent));
	col[0] = ColourLine::create(outgoing[0]);
	col[1] = ColourLine::create(outgoing[2]);
	}
	else if(colourFlow()==2) {
	outgoing[2]->incomingColour (const_ptr_cast<tPPtr>(&parent));
	col[0] = ColourLine::create(outgoing[0]);
	col[1] = ColourLine::create(outgoing[1]);
	}
	else
	assert(false);
	parent.antiColourLine()->setSourceNeighbours(col[0],col[1]);
	}
	else if(antitriplet.size()==3) {
	tColinePtr col[2];
	if(colourFlow()==0) {
	outgoing[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	col[0] = ColourLine::create(outgoing[1],true);
	col[1] = ColourLine::create(outgoing[2],true);
	}
	else if(colourFlow()==1) {
	outgoing[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	col[0] = ColourLine::create(outgoing[0],true);
	col[1] = ColourLine::create(outgoing[2],true);
	}
	else if(colourFlow()==2) {
	outgoing[2]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	col[0] = ColourLine::create(outgoing[0],true);
	col[1] = ColourLine::create(outgoing[1],true);
	}
	else
	assert(false);
	parent.colourLine()->setSinkNeighbours(col[0],col[1]);
	}
	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)};
	Energy absmin = -1e20*GeV;
	int minType = -1;
	for(unsigned int iy=0;iy<_diagmap.size();++iy) {
	unsigned int ix=_diagmap[iy];
	if(getProcessInfo()[ix].channelType==TBDiagram::fourPoint) continue;
	Energy dm1(min-getProcessInfo()[ix].intermediate->mass());
	Energy dm2(getProcessInfo()[ix].intermediate->mass());
	int itype(0);
	if (getProcessInfo()[ix].channelType==TBDiagram::channel23) {
	dm1 -= outgoing()[0]->mass();
	dm2 -= outgoing()[1]->mass()+outgoing()[2]->mass();
	itype = 3;
	}
	else if(getProcessInfo()[ix].channelType==TBDiagram::channel13) {
	dm1 -= outgoing()[1]->mass();
	dm2 -= outgoing()[0]->mass()+outgoing()[2]->mass();
	itype = 2;
	}
	else if(getProcessInfo()[ix].channelType==TBDiagram::channel12) {
	dm1 -= outgoing()[2]->mass();
	dm2 -= outgoing()[0]->mass()+outgoing()[1]->mass();
	itype = 1;
	}
	if((dm1<ZERO\|\|dm2<ZERO)&&getProcessInfo()[ix].intermediate->id()!=ParticleID::gamma) {
	if (imin[itype].first < 0 \|\|
	(dm1<ZERO && imin[itype].second < dm1) ) {
	imin[itype] = make_pair(ix,dm1);
	if(dm1<ZERO&&absmin<dm1) {
	absmin = dm1;
	minType = itype;
	}
	}
	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(widthOption() ==3 ? ZERO : 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;
	}
	}
	// physical poles, use them and return
	if(nchannel>0) {
	inweights = vector<double>(nchannel,1./double(nchannel));
	return;
	}
	// use shallowest pole
	else if(_intOpt==1&&minType>0) {
	intype.push_back(minType);
	if(getProcessInfo()[imin[minType].first].intermediate->id()!=ParticleID::gamma) {
	inpow.push_back(0.);
	inmass.push_back(getProcessInfo()[imin[minType].first].intermediate->mass());
	inwidth.push_back(widthOption() ==3 ? ZERO : getProcessInfo()[imin[minType].first].intermediate->width());
	}
	else {
	inpow.push_back(-2.);
	inmass.push_back(-1.*GeV);
	inwidth.push_back(-1.*GeV);
	}
	inweights = vector<double>(1,1.);
	return;
	}
	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(widthOption() ==3 ? ZERO : 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(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 GeneralThreeBodyDecayer::"
	<< "setColourFactors(), 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 "
	<< "GeneralThreeBodyDecayer::setColourFactors()"
	<< " 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 if(sng.size()==1&&atrip.size()==2) {
	_nflow = 1;
	_colour = vector<DVector>(1,DVector(1,2.));
	_colourLargeNC = vector<DVector>(1,DVector(1,2.));
	}
	else {
	generator()->log() << "Unknown colour structure for "
	<< "triplet decay in "
	<< "GeneralThreeBodyDecayer::setColourFactors()"
	<< " 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 "
	<< "GeneralThreeBodyDecayer::setColourFactors()"
	<< " 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 if(sng.size()==1&&trip.size()==2) {
	_nflow = 1;
	_colour = vector<DVector>(1,DVector(1,2.));
	_colourLargeNC = vector<DVector>(1,DVector(1,2.));
	}
	else {
	generator()->log() << "Unknown colour antitriplet decay in "
	<< "GeneralThreeBodyDecayer::setColourFactors()"
	<< " for " << name << " omitting decay\n";
	return false;
	}
	}
	// colour octet particle
	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));
	}
	// three (anti)triplets
	else if(trip.size()==3\|\|atrip.size()==3) {
	_nflow = 3;
	_colour = vector<DVector>(3,DVector(3,0.));
	_colourLargeNC = vector<DVector>(3,DVector(3,0.));
	_colour[0][0] = 1.; _colour[1][1] = 1.; _colour[2][2] = 1.;
	_colour[0][1] = -0.5; _colour[1][0] = -0.5;
	_colour[0][2] = -0.5; _colour[2][0] = -0.5;
	_colour[1][2] = -0.5; _colour[2][1] = -0.5;
	_colourLargeNC = vector<DVector>(3,DVector(3,0.));
	_colourLargeNC[0][0] = 1.; _colourLargeNC[1][1] = 1.; _colourLargeNC[2][2] = 1.;
	// sett the factors for the diagrams
	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(_diagrams[ix].channelType+1,sign));
	_diagrams[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(_diagrams[ix].channelType+1,sign));
	}
	}
	// unknown
	else {
	generator()->log() << "Unknown colour octet decay in "
	<< "GeneralThreeBodyDecayer::setColourFactors()"
	<< " for " << name << " omitting decay\n";
	return false;
	}
	}
	+ else if (_incoming->iColour() == PDT::Colour6 ) {
	+ generator()->log() << "Unknown colour sextet decay in "
	+ << "GeneralThreeBodyDecayer::setColourFactors()"
	+ << " for " << name << " omitting decay\n";
	+ return false;
	+ }
	+ else if (_incoming->iColour() == PDT::Colour6bar ) {
	+ generator()->log() << "Unknown colour anti-sextet decay in "
	+ << "GeneralThreeBodyDecayer::setColourFactors()"
	+ << " for " << name << " omitting decay\n";
	+ return false;
	+ }
	+
	+ assert(_nflow != 999);
	+
	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";
	}
	for(unsigned int ix=0;ix<_diagrams.size();++ix) {
	generator()->log() << "colour flow for diagram : " << ix;
	for(unsigned int iy=0;iy<_diagrams[ix].colourFlow.size();++iy)
	generator()->log() << "(" << _diagrams[ix].colourFlow[iy].first << ","
	<< _diagrams[ix].colourFlow[iy].second << "); ";
	generator()->log() << "\n";
	}
	}
	return true;
	}
	diff --git a/Models/General/NBodyDecayConstructorBase.cc b/Models/General/NBodyDecayConstructorBase.cc
	--- a/Models/General/NBodyDecayConstructorBase.cc
	+++ b/Models/General/NBodyDecayConstructorBase.cc
	@@ -1,627 +1,627 @@
	// -- 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/Repository/Repository.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_
	<< createModes_ << minReleaseFraction_ << maxBoson_ << maxList_
	<< removeOnShell_ << excludeEffective_ << includeTopOnShell_
	<< excludedVerticesVector_ << excludedVerticesSet_
	<< excludedParticlesVector_ << excludedParticlesSet_;
	}

	void NBodyDecayConstructorBase::persistentInput(PersistentIStream & is , int) {
	is >> init_ >> iteration_ >> points_ >> info_ >> decayConstructor_
	>> createModes_ >> minReleaseFraction_ >> maxBoson_ >> maxList_
	>> 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;
	if ( Debug::level > 0 )
	Repository::cout() << "Constructing N-body decays for "
	<< parent->PDGName() << '\n';
	// first create prototype 1->2 decays
	std::stack<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 stack
	PrototypeVertexPtr proto = prototypes.top();
	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(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<newDiagrams.size();++iy) {
	generator()->log() << "Diagram: " << iy << "\n";
	generator()->log() << newDiagrams[iy] << "\n";
	}
	}
	+ createDecayMode(newDiagrams,possibleOnShell,symfac);
	}
	}
	}

	void NBodyDecayConstructorBase::createDecayMode(vector<NBDiagram> &,
	bool, double) {
	throw Exception() << "In NBodyDecayConstructorBase::createDecayMode() which"
	<< " should have be overridden in an inheriting class"
	<< Exception::abortnow;
	}

File Metadata

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

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions