No OneTemporary
Actions

Size

45 KB

Subscribers

None

View Options

	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::queue<PrototypeVertexPtr> prototypes;
	+ 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 queue
	- PrototypeVertexPtr proto = prototypes.front();
	+ // 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";
	}
	}
	}
	}
	}

	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/PrototypeVertex.cc b/Models/General/PrototypeVertex.cc
	--- a/Models/General/PrototypeVertex.cc
	+++ b/Models/General/PrototypeVertex.cc
	@@ -1,277 +1,277 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the FourBodyDecayConstructor class.
	//

	#include "ThePEG/Helicity/Vertex/VertexBase.h"
	#include "ThePEG/PDT/StandardMatchers.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Utilities/Exception.h"
	#include "PrototypeVertex.h"

	using namespace Herwig;

	void PrototypeVertex::createPrototypes(tPDPtr inpart, VertexBasePtr vertex,
	- std::queue<PrototypeVertexPtr> & prototypes) {
	+ std::stack<PrototypeVertexPtr> & prototypes) {
	int id = inpart->id();
	if(!vertex->isIncoming(inpart)) return;
	for(unsigned int list=0;list<vertex->getNpoint();++list) {
	tPDVector decaylist = vertex->search(list, inpart);
	tPDVector::size_type nd = decaylist.size();
	for( tPDVector::size_type i = 0; i < nd; i += vertex->getNpoint() ) {
	tPDVector pout(decaylist.begin()+i,
	decaylist.begin()+i+vertex->getNpoint());
	OrderedVertices out;
	for(unsigned int ix=1;ix<pout.size();++ix) {
	if(pout[ix]->id() == id ) swap(pout[0], pout[ix]);
	if(pout[ix]->CC()) pout[ix] = pout[ix]->CC();
	out.insert(make_pair(pout[ix],PrototypeVertexPtr()));
	}
	if(vertex->getNpoint()==3) {
	// remove radiation
	if((pout[0]==pout[1] && (pout[2]->id()==ParticleID::gamma\|\|
	pout[2]->id()==ParticleID::g\|\|
	pout[2]->id()==ParticleID::Z0)) \|\|
	(pout[0]==pout[2] && (pout[1]->id()==ParticleID::gamma\|\|
	pout[1]->id()==ParticleID::g\|\|
	pout[1]->id()==ParticleID::Z0)))
	continue;
	}
	prototypes.push(new_ptr(PrototypeVertex(inpart,out,vertex,
	int(vertex->getNpoint())-1)));
	}
	}
	}

	PrototypeVertexPtr PrototypeVertex::replicateTree(PrototypeVertexPtr parent,
	PrototypeVertexPtr oldChild,
	PrototypeVertexPtr & newChild) {
	PrototypeVertexPtr newParent =
	new_ptr(PrototypeVertex(parent->incoming,OrderedVertices(),
	parent->vertex,parent->npart));
	for(OrderedVertices::const_iterator it = parent->outgoing.begin();
	it!=parent->outgoing.end();++it) {
	PrototypeVertexPtr child = it->second ?
	replicateTree(it->second,oldChild,newChild) :
	PrototypeVertexPtr();
	newParent->outgoing.insert(make_pair(it->first,child));
	if(child) child->parent = newParent;
	if(it->second==oldChild) newChild=child;
	}
	if(oldChild==parent) newChild=newParent;
	return newParent;
	}

	void PrototypeVertex::expandPrototypes(PrototypeVertexPtr proto, VertexBasePtr vertex,
	- std::queue<PrototypeVertexPtr> & prototypes,
	+ std::stack<PrototypeVertexPtr> & prototypes,
	const set<PDPtr> & excluded) {
	for(OrderedVertices::const_iterator it = proto->outgoing.begin();
	it!=proto->outgoing.end();++it) {
	if(it->second) {
	expandPrototypes(it->second,vertex,prototypes,excluded);
	}
	else {
	if(!vertex->isIncoming(it->first)) continue;
	if(excluded.find(it->first)!=excluded.end()) continue;
	if(it->first->CC() &&
	excluded.find(it->first->CC())!=excluded.end()) continue;
	int id = it->first->id();
	PrototypeVertexPtr parent=proto;
	while(parent->parent) parent=parent->parent;
	for(unsigned int il = 0; il < vertex->getNpoint(); ++il) {
	tPDVector decaylist = vertex->search(il,it->first );
	tPDVector::size_type nd = decaylist.size();
	for( tPDVector::size_type i = 0; i < nd; i += vertex->getNpoint() ) {
	tPDVector pout(decaylist.begin()+i,
	decaylist.begin()+i+vertex->getNpoint());
	OrderedVertices outgoing;
	for(unsigned int iy=1;iy<pout.size();++iy) {
	if(pout[iy]->id() == id ) swap(pout[0], pout[iy]);
	if(pout[iy]->CC()) pout[iy] = pout[iy]->CC();
	outgoing.insert(make_pair(pout[iy],PrototypeVertexPtr()));
	}
	if(vertex->getNpoint()==3) {
	if((pout[0]==pout[1] && (pout[2]->id()==ParticleID::gamma\|\|
	pout[2]->id()==ParticleID::g\|\|
	pout[2]->id()==ParticleID::Z0)) \|\|
	(pout[0]==pout[2] && (pout[1]->id()==ParticleID::gamma\|\|
	pout[1]->id()==ParticleID::g\|\|
	pout[1]->id()==ParticleID::Z0)))
	continue;
	// remove weak decays of quarks other than top
	if(StandardQCDPartonMatcher::Check(pout[0]->id()) &&
	((StandardQCDPartonMatcher::Check(pout[1]->id()) &&
	abs(pout[2]->id())==ParticleID::Wplus)\|\|
	(StandardQCDPartonMatcher::Check(pout[2]->id())&&
	abs(pout[1]->id())==ParticleID::Wplus))) continue;
	// remove weak decays of leptons
	if((abs(pout[0]->id())>=11&&abs(pout[0]->id())<=16) &&
	(((abs(pout[1]->id())>=11&&abs(pout[1]->id())<=16) &&
	abs(pout[2]->id())==ParticleID::Wplus)\|\|
	((abs(pout[2]->id())>=11&&abs(pout[2]->id())<=16)&&
	abs(pout[1]->id())==ParticleID::Wplus))) continue;
	}
	PrototypeVertexPtr newBranch =
	new_ptr(PrototypeVertex(it->first,
	outgoing,vertex,int(vertex->getNpoint())-1));
	PrototypeVertexPtr newChild;
	PrototypeVertexPtr newVertex = replicateTree(parent,proto,newChild);
	newBranch->parent = newChild;
	OrderedVertices::iterator kt = newChild->outgoing.begin();
	for(OrderedVertices::const_iterator jt = proto->outgoing.begin();
	jt!=it;++jt,++kt) {;}
	pair< tPDPtr, PrototypeVertexPtr > newPair = make_pair(kt->first,newBranch);
	newChild->outgoing.erase(kt);
	newChild->outgoing.insert(newPair);
	newChild->incrementN(newBranch->npart-1);
	prototypes.push(newVertex);
	}
	}
	}
	}
	}

	bool PrototypeVertex::canBeOnShell(unsigned int opt,Energy maxMass,bool first) {
	Energy outMass = outgoingMass();
	if(!first) {
	bool in = maxMass>incomingMass();
	bool out = incomingMass()>outMass;
	if(opt!=0) {
	if( in && ( out \|\| opt==2 ) ) return true;
	}
	else if (incoming->width() == ZERO) {
	tPrototypeVertexPtr original = this;
	while(original->parent) {
	original=original->parent;
	};
	ostringstream name;
	name << original->incoming->PDGName() << " -> ";
	for(OrderedParticles::const_iterator it = original->outPart.begin();
	it!= original->outPart.end();++it)
	name << (**it).PDGName() << " ";
	Throw<InitException>()
	<< "Trying to include on-shell diagram in decay"
	<< name.str() << "including on-shell "
	<< incoming->PDGName() << " with zero width.\n"
	<< "You should make sure that the width for the intermediate is either"
	<< " read from an SLHA file or the intermediate is included in the "
	<< "DecayParticles list of the ModelGenerator.\n"
	<< Exception::runerror;
	}
	}
	else maxMass = incomingMass();
	// check the decay products
	for(OrderedVertices::const_iterator it = outgoing.begin();
	it!=outgoing.end();++it) {
	if(!it->second) continue;
	Energy newMax = maxMass-outMass+it->second->outgoingMass();
	if(it->second->canBeOnShell(opt,newMax,false)) {
	if(first) possibleOnShell = true;
	return true;
	}
	}
	return false;
	}

	bool PrototypeVertex::sameDecay(const PrototypeVertex & x) const {
	if(incoming != x.incoming) return false;
	if(outPart.empty()) setOutgoing();
	if(x.outPart.empty()) x.setOutgoing();
	OrderedParticles::const_iterator cit = outPart.begin();
	OrderedParticles::const_iterator cjt = x.outPart.begin();
	if(x.npart!=npart) return false;
	for(;cit!=outPart.end();++cit,++cjt) {
	if(cit!=cjt) return false;
	}
	return true;
	}

	namespace {
	void constructTag(vector<unsigned int> & order,NBVertex & vertex,
	const OrderedParticles & outgoing,
	vector<bool> & matched) {
	for(list<pair<PDPtr,NBVertex> >::iterator it=vertex.vertices.begin();
	it!=vertex.vertices.end();++it) {
	// search down the tree
	if(it->second.vertex) {
	constructTag(order,it->second,outgoing,matched);
	}
	// identify this particle
	else {
	unsigned int ix=0;
	for(OrderedParticles::const_iterator pit=outgoing.begin();
	pit!=outgoing.end();++pit) {
	if(*pit==it->first&&!matched[ix]) {
	matched[ix] = true;
	order.push_back(ix+1);
	break;
	}
	++ix;
	}
	}
	}
	}
	}

	NBDiagram::NBDiagram(PrototypeVertexPtr proto)
	: NBVertex(proto),
	colourFlow (vector<CFPair>(1,make_pair(1,1.))),
	largeNcColourFlow(vector<CFPair>(1,make_pair(1,1.))) {
	if(!proto) return;
	// finally work out the channel and the ordering
	vector<bool> matched(outgoing.size(),false);
	for(list<pair<PDPtr,NBVertex> >::iterator it=vertices.begin();
	it!=vertices.end();++it) {
	// search down the tree
	if(it->second.vertex) {
	constructTag(channelType,it->second,outgoing,matched);
	}
	// identify this particle
	else {
	unsigned int ix=0;
	for(OrderedParticles::const_iterator pit=outgoing.begin();
	pit!=outgoing.end();++pit) {
	if(*pit==it->first&&!matched[ix]) {
	matched[ix] = true;
	channelType.push_back(ix+1);
	break;
	}
	++ix;
	}
	}
	}
	}

	NBVertex::NBVertex(PrototypeVertexPtr proto) {
	if(!proto) return;
	incoming = proto->incoming;
	outgoing = proto->outPart;
	vertex = proto->vertex;
	// create the vertices
	for(OrderedVertices::const_iterator it=proto->outgoing.begin();
	it!=proto->outgoing.end();++it) {
	vertices.push_back(make_pair(it->first,NBVertex(it->second)));
	}
	// now let's re-order so that branchings are at the end
	for(list<pair<PDPtr,NBVertex> >::iterator it=vertices.begin();
	it!=vertices.end();++it) {
	if(!it->second.incoming) continue;
	list<pair<PDPtr,NBVertex> >::iterator jt=it;
	for( ; jt!=vertices.end();++jt) {
	if(jt==it) continue;
	if(!jt->second.incoming) {
	break;
	}
	}
	if(jt!=vertices.end()) {
	list<pair<PDPtr,NBVertex> >::iterator kt = it;
	while(kt!=jt) {
	list<pair<PDPtr,NBVertex> >::iterator lt = kt;
	++lt;
	swap(kt,lt);
	kt=lt;
	}
	}
	}
	}
	diff --git a/Models/General/PrototypeVertex.h b/Models/General/PrototypeVertex.h
	--- a/Models/General/PrototypeVertex.h
	+++ b/Models/General/PrototypeVertex.h
	@@ -1,474 +1,474 @@
	// -- C++ --
	//
	// PrototypeVertex.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2007 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_PrototypeVertex_H
	#define HERWIG_PrototypeVertex_H
	-#include <queue>
	+#include <stack>
	#include "ThePEG/Helicity/Vertex/VertexBase.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Utilities/EnumIO.h"
	//
	// This is the declaration of the PrototypeVertex class.
	//

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

	class PrototypeVertex;
	ThePEG_DECLARE_POINTERS(Herwig::PrototypeVertex,PrototypeVertexPtr);

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


	/**
	* A struct to order the particles in the same way as in the DecayMode's
	*/
	struct ParticleOrdering {
	/**
	* Operator for the ordering
	* @param p1 The first ParticleData object
	* @param p2 The second ParticleData object
	*/
	bool operator() (PDPtr p1, PDPtr p2) {
	return abs(p1->id()) > abs(p2->id()) \|\|
	( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) \|\|
	( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
	}
	};

	/**
	* A struct to order the particles in the same way as in the DecayMode's
	*/
	struct VertexOrdering {
	/**
	* Operator for the ordering
	* @param p1 The first ParticleData object
	* @param p2 The second ParticleData object
	*/
	bool operator() (const pair< tPDPtr, PrototypeVertexPtr > & p1,
	const pair< tPDPtr, PrototypeVertexPtr > & p2) const {
	return abs(p1.first->id()) > abs(p2.first->id()) \|\|
	( abs(p1.first->id()) == abs(p2.first->id()) && p1.first->id() > p2.first->id() ) \|\|
	( p1.first->id() == p2.first->id() && p1.first->fullName() > p2.first->fullName() );
	}
	};

	typedef multiset<pair< tPDPtr, PrototypeVertexPtr >,VertexOrdering > OrderedVertices;

	/**
	* A set of ParticleData objects ordered as for the DecayMode's
	*/
	typedef multiset<PDPtr,ParticleOrdering> OrderedParticles;

	/**
	* Storage of a potenital n-body decay
	*/
	class PrototypeVertex : public Base {

	public:

	/**
	* Default Constructor
	*/
	PrototypeVertex() : npart(0), possibleOnShell(false) {}

	/**
	* Constructor
	*/
	PrototypeVertex(tPDPtr in, OrderedVertices out,
	VertexBasePtr v, int n) :
	incoming(in), outgoing(out), vertex(v), npart(n),
	possibleOnShell(false) {}

	/**
	* Incoming particle
	*/
	tPDPtr incoming;

	/**
	* Outgoing particles
	*/
	OrderedVertices outgoing;

	/**
	* The vertex for the interaction
	*/
	VertexBasePtr vertex;

	/**
	* The parent of the vertex
	*/
	tPrototypeVertexPtr parent;

	/**
	* Number of particles
	*/
	unsigned int npart;

	/**
	* Outgoing particles
	*/
	mutable OrderedParticles outPart;

	/**
	* Can have on-shell intermediates
	*/
	bool possibleOnShell;

	/**
	* Increment the number of particles
	*/
	void incrementN(int in) {
	npart += in;
	if(parent) parent->incrementN(in);
	}

	/**
	* Mass of the incoming particle
	*/
	Energy incomingMass() {
	return incoming->mass();
	}

	/**
	* Total mass of all the outgoing particles
	*/
	Energy outgoingMass() {
	Energy mass(ZERO);
	for(OrderedVertices::const_iterator it = outgoing.begin();
	it!=outgoing.end();++it) {
	mass += it->second ?
	it->second->outgoingMass() : it->first->mass();
	}
	return mass;
	}

	/**
	* Total constituent mass of all the outgoing particles
	*/
	Energy outgoingConstituentMass() {
	Energy mass(ZERO);
	for(OrderedVertices::const_iterator it = outgoing.begin();
	it!=outgoing.end();++it) {
	mass += it->second ?
	it->second->outgoingConstituentMass() : it->first->constituentMass();
	}
	return mass;
	}

	/**
	* Check the external particles
	*/
	bool checkExternal(bool first=true) {
	if(outPart.empty()) setOutgoing();
	if(first&&outPart.find(incoming)!=outPart.end()) return false;
	bool output = true;
	for(OrderedVertices::const_iterator it = outgoing.begin();
	it!=outgoing.end();++it) {
	if(it->second&& !it->second->checkExternal(false)) output = false;
	}
	return output;
	}

	/**
	* Set the outgoing particles
	*/
	void setOutgoing() const {
	assert(outPart.empty());
	for(OrderedVertices::const_iterator it = outgoing.begin();
	it!=outgoing.end();++it) {
	if(it->second) {
	it->second->setOutgoing();
	outPart.insert(it->second->outPart.begin(),
	it->second->outPart.end());
	}
	else
	outPart.insert(it->first);
	}
	}

	/**
	* Are there potential on-shell intermediates?
	*/
	bool canBeOnShell(unsigned int opt,Energy maxMass,bool first);

	/**
	* Check if same external particles
	*/
	bool sameDecay(const PrototypeVertex & x) const;

	/**
	* Create a \f$1\to2\f$ prototype
	*/
	static void createPrototypes(tPDPtr inpart, VertexBasePtr vertex,
	- std::queue<PrototypeVertexPtr> & prototypes);
	+ std::stack<PrototypeVertexPtr> & prototypes);

	/**
	* Expand the prototypes by adding more legs
	*/
	static void expandPrototypes(PrototypeVertexPtr proto, VertexBasePtr vertex,
	- std::queue<PrototypeVertexPtr> & prototypes,
	+ std::stack<PrototypeVertexPtr> & prototypes,
	const set<PDPtr> & excluded);

	/**
	* Copy the whole structure with a new branching
	*/
	static PrototypeVertexPtr replicateTree(PrototypeVertexPtr parent,
	PrototypeVertexPtr oldChild,
	PrototypeVertexPtr & newChild);

	};

	/**
	* Output to a stream
	*/
	inline ostream & operator<<(ostream & os, const PrototypeVertex & diag) {
	os << diag.incoming->PDGName() << " -> ";
	bool seq=false;
	for(OrderedVertices::const_iterator it = diag.outgoing.begin();
	it!=diag.outgoing.end();++it) {
	os << it->first->PDGName() << " ";
	if(it->second) seq = true;
	}
	os << " decays via "
	<< diag.vertex->fullName() << " in a "
	<< diag.npart << "-body decay\n";
	if(!seq) return os;
	os << "Followed by\n";
	for(OrderedVertices::const_iterator it = diag.outgoing.begin();
	it!=diag.outgoing.end();++it) {
	if(it->second) os << *it->second;
	}
	return os;
	}

	/**
	* Test whether two diagrams are identical.
	*/
	inline bool operator==(const PrototypeVertex & x, const PrototypeVertex & y) {
	if(x.incoming != y.incoming) return false;
	if(x.vertex != y.vertex) return false;
	if(x.npart != y.npart) return false;
	if(x.outgoing.empty()&&y.outgoing.empty()) return true;
	if(x.outgoing.size() != y.outgoing.size()) return false;
	OrderedVertices::const_iterator xt = x.outgoing.begin();
	OrderedVertices::const_iterator yt = y.outgoing.begin();
	while(xt!=x.outgoing.end()) {
	if(xt->first != yt->first) return false;
	// special for identical particles
	OrderedVertices::const_iterator lxt = x.outgoing.lower_bound(*xt);
	OrderedVertices::const_iterator uxt = x.outgoing.upper_bound(*xt);
	--uxt;
	// just one particle
	if(lxt==uxt) {
	if(xt->second && yt->second) {
	if((xt->second)==(yt->second)) {
	++xt;
	++yt;
	continue;
	}
	else return false;
	}
	else if(xt->second \|\| yt->second)
	return false;
	++xt;
	++yt;
	}
	// identical particles
	else {
	++uxt;
	OrderedVertices::const_iterator lyt = y.outgoing.lower_bound(*xt);
	OrderedVertices::const_iterator uyt = y.outgoing.upper_bound(*xt);
	unsigned int nx=0;
	for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) {++nx;}
	unsigned int ny=0;
	for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) {++ny;}
	if(nx!=ny) return false;
	vector<bool> matched(ny,false);
	for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) {
	bool found = false;
	unsigned int iy=0;
	for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) {
	if(matched[iy]) {
	++iy;
	continue;
	}
	if( (!ixt->second &&!iyt->second) \|\|
	( ixt->second&&iyt->second &&
	(ixt->second)==(iyt->second)) ) {
	matched[iy] = true;
	found = true;
	break;
	}
	++iy;
	}
	if(!found) return false;
	}
	xt=uxt;
	yt=uyt;
	}
	}
	return true;
	}

	/**
	* A simple vertex for the N-body diagram
	*/
	struct NBVertex {

	/**
	* Constructor taking a prototype vertex as the arguments
	*/
	NBVertex(PrototypeVertexPtr proto = PrototypeVertexPtr() );

	/**
	* Incoming particle
	*/
	tPDPtr incoming;

	/**
	* Outgoing particles
	*/
	mutable OrderedParticles outgoing;

	/**
	* The vertices
	*/
	list<pair<PDPtr,NBVertex> > vertices;

	/**
	* The vertex
	*/
	VertexBasePtr vertex;
	};

	/**
	* The NBDiagram struct contains information about a \f$1\ton\f$ decay
	* that has been automatically generated.
	*/
	struct NBDiagram : public NBVertex {

	/**
	* Constructor taking a prototype vertex as the arguments*/
	NBDiagram(PrototypeVertexPtr proto=PrototypeVertexPtr());

	/**
	* The type of channel
	*/
	vector<unsigned int> 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;
	};

	/**
	* Output operator to allow the structure to be persistently written
	* @param os The output stream
	* @param x The NBVertex
	*/
	inline PersistentOStream & operator<<(PersistentOStream & os,
	const NBVertex & x) {
	os << x.incoming << x.outgoing << x.vertices << x.vertex;
	return os;
	}

	/**
	* Input operator to allow persistently written data to be read in
	* @param is The input stream
	* @param x The NBVertex
	*/
	inline PersistentIStream & operator>>(PersistentIStream & is,
	NBVertex & x) {
	is >> x.incoming >> x.outgoing >> x.vertices >> x.vertex;
	return is;
	}

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

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

	/**
	* Output a NBVertex to a stream
	*/
	inline ostream & operator<<(ostream & os, const NBVertex & vertex) {
	os << vertex.incoming->PDGName() << " -> ";
	bool seq=false;
	for(list<pair<PDPtr,NBVertex> >::const_iterator it=vertex.vertices.begin();
	it!=vertex.vertices.end();++it) {
	os << it->first->PDGName() << " ";
	if(it->second.incoming) seq = true;
	}
	os << "via vertex " << vertex.vertex->fullName() << "\n";
	if(!seq) return os;
	os << "Followed by\n";
	for(list<pair<PDPtr,NBVertex> >::const_iterator it=vertex.vertices.begin();
	it!=vertex.vertices.end();++it) {
	if(it->second.incoming) os << it->second;
	}
	return os;
	}

	/**
	* Output a NBDiagram to a stream
	*/
	inline ostream & operator<<(ostream & os, const NBDiagram & diag) {
	os << diag.incoming->PDGName() << " -> ";
	for(OrderedParticles::const_iterator it=diag.outgoing.begin();
	it!=diag.outgoing.end();++it) {
	os << (**it).PDGName() << " ";
	}
	os << " has order ";
	for(unsigned int ix=0;ix<diag.channelType.size();++ix)
	os << diag.channelType[ix] << " ";
	os << "\n";
	os << "First decay " << diag.incoming->PDGName() << " -> ";
	bool seq=false;
	for(list<pair<PDPtr,NBVertex> >::const_iterator it=diag.vertices.begin();
	it!=diag.vertices.end();++it) {
	os << it->first->PDGName() << " ";
	if(it->second.incoming) seq = true;
	}
	os << "via vertex " << diag.vertex->fullName() << "\n";
	if(!seq) return os;
	os << "Followed by\n";
	for(list<pair<PDPtr,NBVertex> >::const_iterator it=diag.vertices.begin();
	it!=diag.vertices.end();++it) {
	if(it->second.incoming) os << it->second;
	}
	return os;
	}

	}

	#endif /* HERWIG_PrototypeVertex_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 8:36 PM (1 d, 4 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806082
Default Alt Text: (45 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions