No OneTemporary
Actions

Size

70 KB

Subscribers

None

View Options

	diff --git a/Models/General/TwoToTwoProcessConstructor.cc b/Models/General/TwoToTwoProcessConstructor.cc
	--- a/Models/General/TwoToTwoProcessConstructor.cc
	+++ b/Models/General/TwoToTwoProcessConstructor.cc
	@@ -1,643 +1,643 @@
	// -- C++ --
	//
	// TwoToTwoProcessConstructor.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 TwoToTwoProcessConstructor class.
	//

	#include "TwoToTwoProcessConstructor.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"

	using namespace Herwig;

	namespace {
	// Helper functor for find_if in duplicate function.
	class SameIncomingAs {
	public:
	SameIncomingAs(tPDPair in) : a(in.first->id()), b(in.second->id()) {}
	bool operator()(tPDPair ppair) const {
	long id1(ppair.first->id()), id2(ppair.second->id());
	return ( id1 == a && id2 == b ) \|\| ( id1 == b && id2 == a );
	}
	private:
	long a, b;
	};

	bool duplicateIncoming(tPDPair ppair, const vector<tPDPair> & incPairs ) {
	vector<tPDPair>::const_iterator it =
	find_if( incPairs.begin(), incPairs.end(), SameIncomingAs(ppair) );
	return it != incPairs.end();
	}
	}

	TwoToTwoProcessConstructor::TwoToTwoProcessConstructor() :
	Nout_(0), nv_(0), allDiagrams_(true),
	processOption_(0), scaleChoice_(0), scaleFactor_(1.)
	{}

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

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

	void TwoToTwoProcessConstructor::doinit() {
	HardProcessConstructor::doinit();
	if(processOption_==2&&outgoing_.size()!=2)
	throw InitException()
	<< "Exclusive processes require exactly"
	<< " two outgoing particles but " << outgoing_.size()
	<< "have been inserted in TwoToTwoProcessConstructor::doinit()."
	<< Exception::runerror;
	Nout_ = outgoing_.size();
	PDVector::size_type ninc = incoming_.size();
	// exit if nothing to do
	if(Nout_==0\|\|ninc==0) return;
	//create vector of initial-state pairs
	for(PDVector::size_type i = 0; i < ninc; ++i) {
	for(PDVector::size_type j = 0; j < ninc; ++j) {
	tPDPair inc = make_pair(incoming_[i], incoming_[j]);

	if( (inc.first->iSpin() > inc.second->iSpin()) \|\|
	(inc.first->iSpin() == inc.second->iSpin() &&
	inc.first->id() < inc.second->id()) )
	swap(inc.first, inc.second);

	if( !duplicateIncoming(inc,incPairs_) ) {
	incPairs_.push_back(inc);
	}
	}
	}
	// excluded vertices
	excludedVertexSet_ =
	set<VertexBasePtr>(excludedVertexVector_.begin(),
	excludedVertexVector_.end());
	}


	void TwoToTwoProcessConstructor::persistentOutput(PersistentOStream & os) const {
	os << vertices_ << incoming_ << outgoing_
	<< allDiagrams_ << processOption_
	<< scaleChoice_ << scaleFactor_ << excluded_ << excludedExternal_
	<< excludedVertexVector_ << excludedVertexSet_;
	}

	void TwoToTwoProcessConstructor::persistentInput(PersistentIStream & is, int) {
	is >> vertices_ >> incoming_ >> outgoing_
	>> allDiagrams_ >> processOption_
	>> scaleChoice_ >> scaleFactor_ >> excluded_ >> excludedExternal_
	>> excludedVertexVector_ >> excludedVertexSet_;
	}

	ClassDescription<TwoToTwoProcessConstructor>
	TwoToTwoProcessConstructor::initTwoToTwoProcessConstructor;
	// Definition of the static class description member.

	void TwoToTwoProcessConstructor::Init() {

	static ClassDocumentation<TwoToTwoProcessConstructor> documentation
	("TwoToTwoProcessConstructor constructs the possible diagrams for "
	"a process given the external particles");

	static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceIn
	("Incoming",
	"Pointers to incoming particles",
	&TwoToTwoProcessConstructor::incoming_, -1, false, false, true, false);

	static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceOut
	("Outgoing",
	"Pointers to incoming particles",
	&TwoToTwoProcessConstructor::outgoing_, -1, false, false, true, false);

	static Switch<TwoToTwoProcessConstructor,bool> interfaceIncludeAllDiagrams
	("IncludeEW",
	"Switch to decide which diagrams to include in ME calc.",
	&TwoToTwoProcessConstructor::allDiagrams_, true, false, false);
	static SwitchOption interfaceIncludeAllDiagramsOff
	(interfaceIncludeAllDiagrams,
	"No",
	"Only include QCD diagrams",
	false);
	static SwitchOption interfaceIncludeAllDiagramsOn
	(interfaceIncludeAllDiagrams,
	"Yes",
	"Include EW+QCD.",
	true);

	static Switch<TwoToTwoProcessConstructor,unsigned int> interfaceProcesses
	("Processes",
	"Whether to generate inclusive or exclusive processes",
	&TwoToTwoProcessConstructor::processOption_, 0, false, false);
	static SwitchOption interfaceProcessesSingleParticleInclusive
	(interfaceProcesses,
	"SingleParticleInclusive",
	"Require at least one particle from the list of outgoing particles"
	" in the hard process",
	0);
	static SwitchOption interfaceProcessesTwoParticleInclusive
	(interfaceProcesses,
	"TwoParticleInclusive",
	"Require that both the particles in the hard processes are in the"
	" list of outgoing particles",
	1);
	static SwitchOption interfaceProcessesExclusive
	(interfaceProcesses,
	"Exclusive",
	"Require that both the particles in the hard processes are in the"
	" list of outgoing particles in every hard process",
	2);

	static Switch<TwoToTwoProcessConstructor,unsigned int> interfaceScaleChoice
	("ScaleChoice",
	"&TwoToTwoProcessConstructor::scaleChoice_",
	&TwoToTwoProcessConstructor::scaleChoice_, 0, false, false);
	static SwitchOption interfaceScaleChoiceDefault
	(interfaceScaleChoice,
	"Default",
	"Use if sHat if intermediates all colour neutral, otherwise the transverse mass",
	0);
	static SwitchOption interfaceScaleChoicesHat
	(interfaceScaleChoice,
	"sHat",
	"Always use sHat",
	1);
	static SwitchOption interfaceScaleChoiceTransverseMass
	(interfaceScaleChoice,
	"TransverseMass",
	"Always use the transverse mass",
	2);

	static Parameter<TwoToTwoProcessConstructor,double> interfaceScaleFactor
	("ScaleFactor",
	"The prefactor used in the scale calculation. The scale used is"
	" that defined by scaleChoice multiplied by this prefactor",
	&TwoToTwoProcessConstructor::scaleFactor_, 1.0, 0.0, 10.0,
	false, false, Interface::limited);

	static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceExcluded
	("Excluded",
	"Particles which are not allowed as intermediates",
	&TwoToTwoProcessConstructor::excluded_, -1, false, false, true, false, false);

	static RefVector<TwoToTwoProcessConstructor,ParticleData> interfaceExcludedExternal
	("ExcludedExternal",
	"Particles which are not allowed as outgoing particles",
	&TwoToTwoProcessConstructor::excludedExternal_, -1,
	false, false, true, false, false);

	static RefVector<TwoToTwoProcessConstructor,VertexBase> interfaceExcludedVertices
	("ExcludedVertices",
	"Vertices which are not included in the 2 -> 2 scatterings",
	&TwoToTwoProcessConstructor::excludedVertexVector_, -1, false, false, true, true, false);

	}

	namespace {
	// Helper functor for find_if below.
	class SameProcessAs {
	public:
	SameProcessAs(const HPDiagram & diag) : a(diag) {}
	bool operator()(const HPDiagram & b) const {
	return a.sameProcess(b);
	}
	private:
	HPDiagram a;
	};
	}

	void TwoToTwoProcessConstructor::constructDiagrams() {
	if(incPairs_.empty() \|\| outgoing_.empty()) return;
	nv_ = model()->numberOfVertices();
	//make sure vertices are initialised
	for(unsigned int ix = 0; ix < nv_; ++ix ) {
	VertexBasePtr vertex = model()->vertex(ix);
	if(excludedVertexSet_.find(vertex) !=
	excludedVertexSet_.end()) continue;
	vertices_.push_back(vertex);
	}
	nv_ = vertices_.size();
	//Create necessary diagrams
	vector<tcPDPair>::size_type is;
	PDVector::size_type os;
	for(is = 0; is < incPairs_.size(); ++is) {
	tPDPair ppi = incPairs_[is];
	for(os = 0; os < Nout_; ++os) {
	long fs = outgoing_[os]->id();
	for(size_t iv = 0; iv < nv_; ++iv) {
	tVertexBasePtr vertexA = vertices_[iv];

	//This skips an effective vertex and the EW ones if
	// we only want the strong diagrams
	if( !allDiagrams_ && vertexA->orderInGs() == 0 )
	continue;

	if(vertexA->getNpoint() == 3) {
	//scattering diagrams
	createTChannels(ppi, fs, vertexA);

	//resonance diagrams
	if( vertexA->isIncoming(ppi.first) &&
	vertexA->isIncoming(ppi.second) )
	createSChannels(ppi, fs, vertexA);
	}
	else
	makeFourPointDiagrams(ppi.first->id(), ppi.second->id(),
	fs, vertexA);
	}
	}
	}
	//need to find all of the diagrams that relate to the same process
	//first insert them into a map which uses the '<' operator
	//to sort the diagrams
	multiset<HPDiagram> grouped;
	HPDVector::iterator dit = processes_.begin();
	HPDVector::iterator dend = processes_.end();
	bool abort=false;
	for( ; dit != dend; ++dit) {
	// check for on-shell s-channel
	tPDPtr out1 = getParticleData(dit->outgoing.first );
	tPDPtr out2 = getParticleData(dit->outgoing.second);
	if(dit->channelType == HPDiagram::sChannel &&
	dit->intermediate->width()==ZERO &&
	dit->intermediate->mass() > out1->mass()+ out2->mass()) {
	tPDPtr in1 = getParticleData(dit->incoming.first );
	tPDPtr in2 = getParticleData(dit->incoming.second);
	generator()->log() << dit->intermediate->PDGName()
	<< " can be on-shell in the process "
	<< in1 ->PDGName() << " " << in2->PDGName() << " -> "
	<< out1->PDGName() << " " << out2->PDGName()
	<< " but has zero width.\nEither set the width, enable "
	<< "calculation of its decays, and hence the width,\n"
	<< "or disable it as a potential intermediate using\n"
	<< "insert " << fullName() << ":Excluded 0 "
	<< dit->intermediate->fullName() << "\n---\n";
	abort = true;
	}
	grouped.insert(*dit);
	}
	if(abort) throw Exception() << "One or more processes with zero width"
	<< " resonant intermediates"
	<< Exception::runerror;
	assert( processes_.size() == grouped.size() );
	processes_.clear();
	typedef multiset<HPDiagram>::const_iterator set_iter;
	set_iter it = grouped.begin(), iend = grouped.end();
	while( it != iend ) {
	pair<set_iter,set_iter> range = grouped.equal_range(*it);
	set_iter itb = range.first;
	HPDVector process;
	for( ; itb != range.second; ++itb ) {
	process.push_back(*itb);
	}
	// if inclusive enforce the exclusivity
	if(processOption_==2) {
	if(!((process[0].outgoing. first==outgoing_[0]->id()&&
	process[0].outgoing.second==outgoing_[1]->id())\|\|
	(process[0].outgoing. first==outgoing_[1]->id()&&
	process[0].outgoing.second==outgoing_[0]->id()))) {
	process.clear();
	it = range.second;
	continue;
	}
	}
	if(find(excludedExternal_.begin(),excludedExternal_.end(),
	getParticleData(process[0].outgoing. first))!=excludedExternal_.end()) {
	process.clear();
	it = range.second;
	continue;
	}
	if(find(excludedExternal_.begin(),excludedExternal_.end(),
	getParticleData(process[0].outgoing.second))!=excludedExternal_.end()) {
	process.clear();
	it = range.second;
	continue;
	}
	createMatrixElement(process);
	process.clear();
	it = range.second;
	}
	}

	void TwoToTwoProcessConstructor::
	createSChannels(tcPDPair inpp, long fs, tVertexBasePtr vertex) {
	//Have 2 incoming particle and a vertex, find the possible offshell
	//particles
	pair<long,long> inc = make_pair(inpp.first->id(), inpp.second->id());
	tPDSet offshells = search(vertex, inpp.first->id(), incoming,
	inpp.second->id(), incoming, outgoing);
	tPDSet::const_iterator it;
	for(it = offshells.begin(); it != offshells.end(); ++it) {
	if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
	for(size_t iv = 0; iv < nv_; ++iv) {
	tVertexBasePtr vertexB = vertices_[iv];
	if( vertexB->getNpoint() != 3) continue;
	if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;

	tPDSet final;
	if( vertexB->isOutgoing(getParticleData(fs)) &&
	vertexB->isIncoming(*it) )
	final = search(vertexB, (*it)->id(), incoming, fs,
	outgoing, outgoing);
	//Now make diagrams
	if(!final.empty())
	makeDiagrams(inc, fs, final, *it, HPDiagram::sChannel,
	make_pair(vertex, vertexB), make_pair(true,true));
	}
	}

	}

	void TwoToTwoProcessConstructor::
	createTChannels(tPDPair inpp, long fs, tVertexBasePtr vertex) {
	pair<long,long> inc = make_pair(inpp.first->id(), inpp.second->id());
	//first try a with c
	tPDSet offshells = search(vertex, inpp.first->id(), incoming, fs,
	outgoing, outgoing);
	tPDSet::const_iterator it;
	for(it = offshells.begin(); it != offshells.end(); ++it) {
	if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
	for(size_t iv = 0; iv < nv_; ++iv) {
	tVertexBasePtr vertexB = vertices_[iv];
	if( vertexB->getNpoint() != 3 ) continue;
	if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;
	tPDSet final;
	if( vertexB->isIncoming(inpp.second) )
	final = search(vertexB, inc.second, incoming, (*it)->id(),
	incoming, outgoing);
	if( !final.empty() )
	makeDiagrams(inc, fs, final, *it, HPDiagram::tChannel,
	make_pair(vertex, vertexB), make_pair(true,true));
	}
	}
	//now try b with c
	offshells = search(vertex, inpp.second->id(), incoming, fs,
	outgoing, incoming);
	for(it = offshells.begin(); it != offshells.end(); ++it) {
	if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
	for(size_t iv = 0; iv < nv_; ++iv) {
	tVertexBasePtr vertexB = vertices_[iv];
	if( vertexB->getNpoint() != 3 ) continue;
	if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;

	tPDSet final;
	if( vertexB->isIncoming(inpp.first) )
	final = search(vertexB, inc.first, incoming, (*it)->id(),
	outgoing, outgoing);
	if( !final.empty() )
	makeDiagrams(inc, fs, final, *it, HPDiagram::tChannel,
	make_pair(vertexB, vertex), make_pair(true, false));
	}
	}

	}

	void TwoToTwoProcessConstructor::makeFourPointDiagrams(long parta, long partb,
	long partc, VBPtr vert) {
	if(processOption_>=1) {
	PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
	getParticleData(partc));
	if(loc==outgoing_.end()) return;
	}
	tPDSet ext = search(vert, parta, incoming, partb,incoming, partc, outgoing);
	if( ext.empty() ) return;
	IDPair in(parta, partb);
	for(tPDSet::const_iterator iter=ext.begin(); iter!=ext.end();
	++iter) {
	if(processOption_>=1) {
	PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
	*iter);
	if(loc==outgoing_.end()) continue;
	}
	HPDiagram nhp(in,make_pair(partc, (*iter)->id()));
	nhp.vertices = make_pair(vert, vert);
	nhp.channelType = HPDiagram::fourPoint;
	fixFSOrder(nhp);
	if( !duplicate(nhp, processes_) ) {
	assignToCF(nhp);
	processes_.push_back(nhp);
	}
	}
	}

	void
	TwoToTwoProcessConstructor::makeDiagrams(IDPair in, long out1, const tPDSet & out2,
	PDPtr inter, HPDiagram::Channel chan,
	VBPair vertexpair, BPair cross) {
	if(processOption_>=1) {
	PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
	getParticleData(out1));
	if(loc==outgoing_.end()) return;
	}
	for(tPDSet::const_iterator it = out2.begin(); it != out2.end(); ++it) {
	if(processOption_>=1) {
	PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
	*it);
	if(loc==outgoing_.end()) continue;
	}
	HPDiagram nhp( in,make_pair(out1, (*it)->id()) );
	nhp.intermediate = inter;
	nhp.vertices = vertexpair;
	nhp.channelType = chan;
	nhp.ordered = cross;
	fixFSOrder(nhp);
	if( !duplicate(nhp, processes_) ) {
	assignToCF(nhp);
	processes_.push_back(nhp);
	}
	}
	}

	set<tPDPtr>
	TwoToTwoProcessConstructor::search(VBPtr vertex, long part1, direction d1,
	long part2, direction d2, direction d3) {
	if(vertex->getNpoint() != 3) return tPDSet();
	if(d1 == incoming && getParticleData(part1)->CC()) part1 = -part1;
	if(d2 == incoming && getParticleData(part2)->CC()) part2 = -part2;
	vector<long> ext;
	tPDSet third;
	for(unsigned int ix = 0;ix < 3; ++ix) {
	vector<long> pdlist = vertex->search(ix, part1);
	ext.insert(ext.end(), pdlist.begin(), pdlist.end());
	}
	for(unsigned int ix = 0; ix < ext.size(); ix += 3) {
	long id0 = ext.at(ix);
	long id1 = ext.at(ix+1);
	long id2 = ext.at(ix+2);
	int pos;
	if((id0 == part1 && id1 == part2) \|\|
	(id0 == part2 && id1 == part1))
	pos = ix + 2;
	else if((id0 == part1 && id2 == part2) \|\|
	(id0 == part2 && id2 == part1))
	pos = ix + 1;
	else if((id1 == part1 && id2 == part2) \|\|
	(id1 == part2 && id2 == part1))
	pos = ix;
	else
	pos = -1;
	if(pos >= 0) {
	tPDPtr p = getParticleData(ext[pos]);
	if(d3 == incoming && p->CC()) p = p->CC();
	third.insert(p);
	}
	}

	return third;
	}

	set<tPDPtr>
	TwoToTwoProcessConstructor::search(VBPtr vertex, long part1, direction d1,
	long part2, direction d2, long part3, direction d3,
	direction d4) {
	if(vertex->getNpoint() != 4) return tPDSet();
	if(d1 == incoming && getParticleData(part1)->CC()) part1 = -part1;
	if(d2 == incoming && getParticleData(part2)->CC()) part2 = -part2;
	if(d3 == incoming && getParticleData(part3)->CC()) part3 = -part3;
	vector<long> ext;
	tPDSet fourth;
	for(unsigned int ix = 0;ix < 4; ++ix) {
	vector<long> pdlist = vertex->search(ix, part1);
	ext.insert(ext.end(), pdlist.begin(), pdlist.end());
	}
	for(unsigned int ix = 0;ix < ext.size(); ix += 4) {
	long id0 = ext.at(ix); long id1 = ext.at(ix + 1);
	long id2 = ext.at(ix + 2); long id3 = ext.at(ix + 3);
	int pos;
	if((id0 == part1 && id1 == part2 && id2 == part3) \|\|
	(id0 == part1 && id1 == part3 && id2 == part2) \|\|
	(id0 == part2 && id1 == part1 && id2 == part3) \|\|
	(id0 == part2 && id1 == part3 && id2 == part1) \|\|
	(id0 == part3 && id1 == part1 && id2 == part2) \|\|
	(id0 == part3 && id1 == part2 && id2 == part1))
	pos = ix + 3;
	else if((id0 == part1 && id1 == part2 && id3 == part3) \|\|
	(id0 == part1 && id1 == part3 && id3 == part2) \|\|
	(id0 == part2 && id1 == part1 && id3 == part3) \|\|
	(id0 == part2 && id1 == part3 && id3 == part1) \|\|
	(id0 == part3 && id1 == part1 && id3 == part2) \|\|
	(id0 == part3 && id1 == part2 && id3 == part1))
	pos = ix + 2;
	else if((id0 == part1 && id2 == part2 && id3 == part3) \|\|
	(id0 == part1 && id2 == part3 && id3 == part2) \|\|
	(id0 == part2 && id2 == part1 && id3 == part3) \|\|
	(id0 == part2 && id2 == part3 && id3 == part1) \|\|
	(id0 == part3 && id2 == part1 && id3 == part2) \|\|
	(id0 == part3 && id2 == part2 && id3 == part1))
	pos = ix + 1;
	else if((id1 == part1 && id2 == part2 && id3 == part3) \|\|
	(id1 == part1 && id2 == part3 && id3 == part2) \|\|
	(id1 == part2 && id2 == part1 && id3 == part3) \|\|
	(id1 == part2 && id2 == part3 && id3 == part1) \|\|
	(id1 == part3 && id2 == part1 && id3 == part2) \|\|
	(id1 == part3 && id2 == part2 && id3 == part1))
	pos = ix;
	else
	pos = -1;

	if(pos >= 0) {
	tPDPtr p = getParticleData(ext[pos]);
	if(d4 == incoming && p->CC())
	p = p->CC();
	fourth.insert(p);
	}
	}
	return fourth;
	}

	void
	TwoToTwoProcessConstructor::createMatrixElement(const HPDVector & process) const {
	if ( process.empty() ) return;
	// external particles
	tcPDVector extpart(4);
	extpart[0] = getParticleData(process[0].incoming.first);
	extpart[1] = getParticleData(process[0].incoming.second);
	extpart[2] = getParticleData(process[0].outgoing.first);
	extpart[3] = getParticleData(process[0].outgoing.second);
	// create the object
	string objectname ("/Herwig/MatrixElements/");
	string classname = MEClassname(extpart, objectname);
	GeneralHardMEPtr matrixElement = dynamic_ptr_cast<GeneralHardMEPtr>
	(generator()->preinitCreate(classname, objectname));
	if( !matrixElement ) {
	- stringstream message;
	+ std::stringstream message;
	message << "TwoToTwoProcessConstructor::createMatrixElement "
	<< "- No matrix element object could be created for "
	<< "the process "
	<< extpart[0]->PDGName() << " " << extpart[0]->iSpin() << ","
	<< extpart[1]->PDGName() << " " << extpart[1]->iSpin() << "->"
	<< extpart[2]->PDGName() << " " << extpart[2]->iSpin() << ","
	<< extpart[3]->PDGName() << " " << extpart[3]->iSpin()
	<< ". Constructed class name: \"" << classname << "\"";
	generator()->logWarning(TwoToTwoProcessConstructorError(message.str(),Exception::warning));
	return;
	}
	// choice for the scale
	unsigned int scale;
	if(scaleChoice_==0) {
	// check coloured initial and final state
	bool inColour = ( extpart[0]->coloured() \|\|
	extpart[1]->coloured());
	bool outColour = ( extpart[2]->coloured() \|\|
	extpart[3]->coloured());
	if(inColour&&outColour) {
	bool coloured = false;
	for(unsigned int ix=0;ix<process.size();++ix) {
	if(process[ix].intermediate&&
	process[ix].intermediate->coloured()) {
	coloured = true;
	break;
	}
	}
	scale = coloured ? 1 : 0;
	}
	else {
	scale = 0;
	}
	}
	else {
	scale = scaleChoice_-1;
	}
	// set the information
	matrixElement->setProcessInfo(process, colourFlow(extpart),
	debug(), scale, scaleFactor_);
	// insert it
	generator()->preinitInterface(subProcess(), "MatrixElements",
	subProcess()->MEs().size(),
	"insert", matrixElement->fullName());
	}

	string TwoToTwoProcessConstructor::MEClassname(const vector<tcPDPtr> & extpart,
	string & objname) const {
	string classname("Herwig::ME");
	for(vector<tcPDPtr>::size_type ix = 0; ix < extpart.size(); ++ix) {
	if(ix == 2) classname += "2";
	if(extpart[ix]->iSpin() == PDT::Spin0) classname += "s";
	else if(extpart[ix]->iSpin() == PDT::Spin1) classname += "v";
	else if(extpart[ix]->iSpin() == PDT::Spin1Half) classname += "f";
	else if(extpart[ix]->iSpin() == PDT::Spin2) classname += "t";
	else {
	- stringstream message;
	+ std::stringstream message;
	message << "MEClassname() : Encountered an unknown spin for "
	<< extpart[ix]->PDGName() << " while constructing MatrixElement "
	<< "classname " << extpart[ix]->iSpin();
	generator()->logWarning(TwoToTwoProcessConstructorError(message.str(),Exception::warning));
	}
	}
	objname += "ME" + extpart[0]->PDGName() + extpart[1]->PDGName() + "2"
	+ extpart[2]->PDGName() + extpart[3]->PDGName();
	return classname;
	}
	diff --git a/Models/Leptoquarks/LeptoquarkModelSLQFFVertex.cc b/Models/Leptoquarks/LeptoquarkModelSLQFFVertex.cc
	--- a/Models/Leptoquarks/LeptoquarkModelSLQFFVertex.cc
	+++ b/Models/Leptoquarks/LeptoquarkModelSLQFFVertex.cc
	@@ -1,337 +1,337 @@
	// -- C++ --
	//
	// LeptoquarkModelSLQFFVertex.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 LeptoquarkModelSLQFFVertex class.
	//

	#include "LeptoquarkModelSLQFFVertex.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Config/Constants.h"

	using namespace Herwig;

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

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

	LeptoquarkModelSLQFFVertex::LeptoquarkModelSLQFFVertex()
	{}

	void LeptoquarkModelSLQFFVertex::doinit() {
	//S0
	- addToList(15,6,-9911561);
	- addToList(-15,-6,9911561);
	+ addToList( 15, 6,-9911561);
	+ addToList(-15,-6, 9911561);

	- addToList(-16,-5,9911561);
	- addToList(16,5,-9911561);
	+ addToList(-16,-5, 9911561);
	+ addToList( 16, 5,-9911561);

	//~S0
	- addToList(-15,-5,9921551);
	- addToList(15,5,-9921551);
	+ addToList(-15,-5, 9921551);
	+ addToList( 15, 5,-9921551);

	//S1 triplet
	//S1p
	- addToList(-15,-5,9931551);
	- addToList(15,5,-9931551);
	+ addToList(-15,-5, 9931551);
	+ addToList( 15, 5,-9931551);
	//S1z
	- addToList(-15,-6,9931561);
	- addToList(15,6,-9931561);
	- addToList(-16,-5,9931561);
	- addToList(16,5,-9931561);
	+ addToList(-15,-6, 9931561);
	+ addToList( 15, 6,-9931561);
	+ addToList(-16,-5, 9931561);
	+ addToList( 16, 5,-9931561);
	//S1m
	- addToList(-16,-6,9931661);
	- addToList(16,6,-9931661);
	+ addToList(-16,-6, 9931661);
	+ addToList( 16, 6,-9931661);

	//S1/2 doublet
	- addToList(15,-6,9941561);
	- addToList(-15,6,-9941561);
	+ addToList( 15,-6, 9941561);
	+ addToList(-15, 6,-9941561);

	- addToList(-15,5,-9941551);
	- addToList(-16,-6,-9941551);
	- addToList(15,-5,9941551);
	- addToList(16,-6,9941551);
	+ addToList(-15, 5,-9941551);
	+ addToList(-16, 6,-9941551);
	+ addToList( 15,-5, 9941551);
	+ addToList( 16,-6, 9941551);


	//S1/2 tilde doublet
	- addToList(5,-16,-9951651);
	- addToList(-5, 16,9951651);
	+ addToList( 5,-16,-9951651);
	+ addToList(-5, 16, 9951651);

	- addToList(-5, 15,9951551);
	- addToList(5,-15,-9951551);
	+ addToList(-5, 15, 9951551);
	+ addToList( 5,-15,-9951551);


	//dS0
	- addToList(15,-5,9961551);
	- addToList(-15,5,-9961551);
	+ addToList( 15,-5, 9961551);
	+ addToList(-15, 5,-9961551);

	- addToList(16,-6,9961551);
	- addToList(-16,6,-9961551);
	+ addToList( 16,-6, 9961551);
	+ addToList(-16, 6,-9961551);

	//~dS0
	- addToList(15,-6,9971561);
	- addToList(-15,6,-9971561);
	+ addToList( 15,-6, 9971561);
	+ addToList(-15, 6,-9971561);


	//dS1 triplet

	//dS1p
	- addToList(15,-6,9981561);
	- addToList(-15,6,-9981561);
	+ addToList( 15,-6, 9981561);
	+ addToList(-15, 6,-9981561);

	//dS1z
	- addToList(16,-6,9981551);
	- addToList(-16,6,-9981551);
	+ addToList( 16,-6, 9981551);
	+ addToList(-16, 6,-9981551);

	- addToList(15,-5,9981551);
	- addToList(-15,5,-9981551);
	+ addToList( 15,-5, 9981551);
	+ addToList(-15, 5,-9981551);

	//dS1m
	- addToList(16,-5,9981651);
	- addToList(-16,5,-9981651);
	+ addToList( 16,-5, 9981651);
	+ addToList(-16, 5,-9981651);

	//dS1/2 doublet
	- addToList(-15,-5,9991551);
	- addToList(15,5,-9991551);
	+ addToList(-15,-5, 9991551);
	+ addToList( 15, 5,-9991551);

	- addToList(-15,-6,9991561);
	- addToList(15,6,-9991561);
	+ addToList(-15,-6, 9991561);
	+ addToList( 15, 6,-9991561);

	- addToList(-16,-5,9991561);
	- addToList(16,5,-9991561);
	+ addToList(-16,-5, 9991561);
	+ addToList( 16, 5,-9991561);

	//dS1/2 tilde doublet
	- addToList(-15,-6,9901561);
	- addToList(15,6,-9901561);
	+ addToList(-15,-6, 9901561);
	+ addToList( 15, 6,-9901561);

	- addToList(-16,-6,9901661);
	- addToList(16,6,-9901661);
	+ addToList(-16,-6, 9901661);
	+ addToList( 16, 6,-9901661);


	orderInGem(0);
	orderInGs(0);
	_theModel = generator()->standardModel();
	tcHwLeptoquarkPtr hwLeptoquark=dynamic_ptr_cast<tcHwLeptoquarkPtr>(_theModel);
	if(hwLeptoquark){
	_CFF=hwLeptoquark->cfermion();
	_cL0 =hwLeptoquark->cleft();
	_cR0 =hwLeptoquark->cright();
	_cR0t = hwLeptoquark->crighttilde();
	_cL1 =hwLeptoquark->cleft1();
	_cL12 =hwLeptoquark->cleft12();
	_cR12 =hwLeptoquark->cright12();
	_cL12t =hwLeptoquark->cleft12tilde();


	_derivscale = hwLeptoquark->fscale();
	_dcL0 =hwLeptoquark->dcleft();
	_dcR0 =hwLeptoquark->dcright();
	_dcR0t = hwLeptoquark->dcrighttilde();
	_dcL1 =hwLeptoquark->dcleft1();
	_dcL12 =hwLeptoquark->dcleft12();
	_dcR12 =hwLeptoquark->dcright12();
	_dcL12t =hwLeptoquark->dcleft12tilde();

	}
	FFSVertex::doinit();
	}

	void LeptoquarkModelSLQFFVertex::persistentOutput(PersistentOStream & os) const {
	os << _CFF << _cL0 << _cR0 << _cR0t << _cL1 << _cL12 << _cR12 << _cL12t << _dcL0 << _dcR0 << _dcR0t << _dcL1 << _dcL12 << _dcR12 << _dcL12t << _derivscale;
	}

	void LeptoquarkModelSLQFFVertex::persistentInput(PersistentIStream & is, int) {
	is >> _CFF >> _cL0 >> _cR0 >> _cR0t >> _cL1 >> _cL12 >> _cR12 >> _cL12t >>_dcL0 >> _dcR0 >> _dcR0t >> _dcL1 >> _dcL12 >> _dcR12 >> _dcL12t >> _derivscale;
	}

	ClassDescription<LeptoquarkModelSLQFFVertex>
	LeptoquarkModelSLQFFVertex::initLeptoquarkModelSLQFFVertex;
	// Definition of the static class description member.


	void LeptoquarkModelSLQFFVertex::Init() {

	static ClassDocumentation<LeptoquarkModelSLQFFVertex> documentation
	("The LeptoquarkModelSLQFFVertex class is the implementation"
	" of the helicity amplitude calculation of the Leptoquark"
	" quark-lepton vertex.");
	}


	void LeptoquarkModelSLQFFVertex::setCoupling(Energy2,tcPDPtr aa ,tcPDPtr bb, tcPDPtr cc) {
	long isc(cc->id()), ism(aa->id()),
	ichg(bb->id());
	long lqid = isc;
	long smid_1 = ism;
	long smid_2 = ichg;
	if(fabs(lqid) < 9900000) {
	lqid = ism;
	smid_1 = ichg;
	smid_2 = isc;
	}
	if(fabs(lqid) < 9900000) {
	lqid = ichg;
	smid_1 = ism;
	smid_2 = isc;
	}
	if( fabs(smid_1) > fabs(smid_2) ) { swap(smid_1, smid_2); }

	double mtop = 174.2;
	double mbot = 4.2;
	double mtau = 1.77699;

	//set the couplings to left and right
	//S0
	if( fabs(isc) == 9911561 \|\| fabs(ism) == 9911561 \|\| fabs(ichg) == 9911561 ) {
	if(fabs(isc) == 5 \|\| fabs(ism) == 5 \|\| fabs(ichg) == 5) {
	_cL = -_cL0; _cR = Complex(0.);
	}
	if(fabs(isc) == 6 \|\| fabs(ism) == 6 \|\| fabs(ichg) == 6) {
	_cL = _cL0;
	_cR = _cR0;
	}
	}
	//~S0
	if( fabs(isc) == 9921551 \|\| fabs(ism) == 9921551 \|\| fabs(ichg) == 9921551 ) {
	_cL = Complex(0.); _cR = _cR0t;
	}

	//S1 triplet
	//Q = + 4/3
	if( fabs(isc) == 9931551 \|\| fabs(ism) == 9931551 \|\| fabs(ichg) == 9931551 ) {
	_cL = sqrt(2.)* _cL1; _cR = Complex(0.);
	}
	//Q = + 1/3
	if( fabs(isc) == 9931561 \|\| fabs(ism) == 9931561 \|\| fabs(ichg) == 9931561 ) {
	_cL = - _cL1; _cR = Complex(0.);
	}
	//Q = - 2/3
	if( fabs(isc) == 9931661 \|\| fabs(ism) == 9931661 \|\| fabs(ichg) == 9931661 ) {
	_cL = sqrt(2.) * _cL1; _cR = Complex(0.);
	}

	//S1/2 doublet

	//Q = + 5/3
	if( fabs(isc) == 9941561 \|\| fabs(ism) == 9941561 \|\| fabs(ichg) == 9941561 ) {
	_cR = _cL12; _cL = _cR12;
	}


	//Q = + 2/3
	if( fabs(isc) == 9941551 \|\| fabs(ism) == 9941551 \|\| fabs(ichg) == 9941551 ) {
	if(fabs(isc) == 5 \|\| fabs(ism) == 5 \|\| fabs(ichg) == 5) {
	_cR = Complex(0.); _cL = - _cR12;
	}
	if(fabs(isc) == 6 \|\| fabs(ism) == 6 \|\| fabs(ichg) == 6) {
	_cL = Complex(0.); _cR = _cL12;
	}
	}

	//S1/2 tilde doublet

	//Q = + 2/3
	if( fabs(isc) == 9951551 \|\| fabs(ism) == 9951551 \|\| fabs(ichg) == 9951551 ) {
	_cR = _cL12t; _cL = Complex(0.);
	}


	//Q = - 1/3
	if( fabs(isc) == 9951651 \|\| fabs(ism) == 9951651 \|\| fabs(ichg) == 9951651 ) {
	_cR = _cL12t; _cL = Complex(0.);
	}

	//dS0
	if( fabs(isc) == 9961551 \|\| fabs(ism) == 9961551 \|\| fabs(ichg) == 9961551) {
	if(fabs(isc) == 5 \|\| fabs(ism) == 5 \|\| fabs(ichg) == 5) {
	_cR = _dcL0 * mbot +_dcR0 * mtau;
	_cL = _dcR0 * mbot + _dcL0 * mtau;
	}
	if(fabs(isc) == 6 \|\| fabs(ism) == 6 \|\| fabs(ichg) == 6) {
	_cR = _dcL0 * mtop;
	_cL = Complex(0.);
	}
	_cL /= sqrt(2.) * _derivscale;
	_cR /= sqrt(2.) * _derivscale;
	}

	//d~S0
	if( fabs(isc) == 9971561 \|\| fabs(ism) == 9971561 \|\| fabs(ichg) == 9971561) {
	_cR = _dcR0t * mtau / (sqrt(2.) * _derivscale);
	_cL = _dcR0t * mtop / (sqrt(2.) * _derivscale);
	}

	//dS1 triplet
	if( fabs(isc) == 9981561 \|\| fabs(ism) == 9981561 \|\| fabs(ichg) == 9981561) {
	_cR = sqrt(2.) * _dcL1 * mtop / (sqrt(2.) * _derivscale);
	_cL = sqrt(2.) * _dcL1 * mtau / (sqrt(2.) * _derivscale);
	}
	if( fabs(isc) == 9981551 \|\| fabs(ism) == 9981551 \|\| fabs(ichg) == 9981551) {
	if(fabs(isc) == 5 \|\| fabs(ism) == 5 \|\| fabs(ichg) == 5) {
	_cR = -_dcL1 * mbot;
	_cL = -_dcL1 * mtau;
	}
	if(fabs(isc) == 6 \|\| fabs(ism) == 6 \|\| fabs(ichg) == 6) {
	_cR = _dcL1 * mtop;
	_cL = Complex(0.);
	}
	_cL /= sqrt(2.) * _derivscale;
	_cR /= sqrt(2.) * _derivscale;

	}

	if( fabs(isc) == 9981651 \|\| fabs(ism) == 9981651 \|\| fabs(ichg) == 9981651) {
	_cL = sqrt(2.) * _dcL1 * mbot / (sqrt(2.) * _derivscale);
	_cR = Complex(0.);
	}


	//dS1/2 doublet
	if( fabs(isc) == 9991551 \|\| fabs(ism) == 9991551 \|\| fabs(ichg) == 9991551 ) {
	_cL = _dcL12 * mbot + _dcR12 * mtau;
	_cR = _dcR12 * mbot + _dcL12 * mtau;
	_cL /= sqrt(2.) * _derivscale;
	_cR /= sqrt(2.) * _derivscale;
	}

	if( fabs(isc) == 9991561 \|\| fabs(ism) == 9991561 \|\| fabs(ichg) == 9991561 ) {
	if(fabs(isc) == 6 \|\| fabs(ism) == 6 \|\| fabs(ichg) == 6) {
	_cL = _dcR12 * mtau;
	_cR = _dcR12 * mtop;
	}
	if(fabs(isc) == 5 \|\| fabs(ism) == 5 \|\| fabs(ichg) == 5) {
	_cL = _dcL12 * mbot;
	}

	_cL /= sqrt(2.) * _derivscale;
	_cR /= sqrt(2.) * _derivscale;
	}

	//dS1/2 tilde doublet
	if( fabs(isc) == 9901561 \|\| fabs(ism) == 9901561 \|\| fabs(ichg) == 9901561 ) {
	_cL = _dcL12t * mtop / (sqrt(2.) * _derivscale); _cR = _dcL12t * mtau / (sqrt(2.) * _derivscale);
	}

	if( fabs(isc) == 9901661 \|\| fabs(ism) == 9901661 \|\| fabs(ichg) == 9901661 ) {
	_cL = _dcL12t * mtop / (sqrt(2.) * _derivscale); _cR = Complex(0.);
	}
	if(smid_1 > 0) {
	left(conj(_cR)); right(conj(_cL));
	} else { left(_cL); right(_cR); }



	norm(_CFF);
	}
	diff --git a/Models/Susy/NMSSM/NMSSMWWHHVertex.cc b/Models/Susy/NMSSM/NMSSMWWHHVertex.cc
	--- a/Models/Susy/NMSSM/NMSSMWWHHVertex.cc
	+++ b/Models/Susy/NMSSM/NMSSMWWHHVertex.cc
	@@ -1,173 +1,173 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the NMSSMWWHHVertex class.
	//

	#include "NMSSMWWHHVertex.h"
	#include "NMSSM.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace Herwig;

	NMSSMWWHHVertex::NMSSMWWHHVertex() : couplast_(0.),q2last_(ZERO),
	sw_(0.), cw_(0.), sb_(0.), cb_(0.)
	{}

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

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

	void NMSSMWWHHVertex::persistentOutput(PersistentOStream & os) const {
	os << sw_ << cw_ << sb_ << cb_ << mixS_ << mixP_;
	}

	void NMSSMWWHHVertex::persistentInput(PersistentIStream & is, int) {
	is >> sw_ >> cw_ >> sb_ >> cb_ >> mixS_ >> mixP_;
	}

	ClassDescription<NMSSMWWHHVertex> NMSSMWWHHVertex::initNMSSMWWHHVertex;
	// Definition of the static class description member.

	void NMSSMWWHHVertex::Init() {

	static ClassDocumentation<NMSSMWWHHVertex> documentation
	("The NMSSMWWHHVertex class implements the coupling of"
	" two electroweak and two Higgs bosons in the NMSSM.");

	}

	void NMSSMWWHHVertex::doinit() {
	int scalar[3]={25,35,45};
	int pseudo[2]={36,46};
	// scalar higgs bosons
	for(unsigned int i=0;i<3;++i) {
	// pair of scalars
	for(unsigned int j=0;j<3;++j) {
	addToList( 24,-24,scalar[i],scalar[j]);
	addToList( 23, 23,scalar[i],scalar[j]);
	}
	// scalar charged
	- addToList( 22, 24,scalar[i], 37);
	- addToList( 22,-24,scalar[i],-37);
	- addToList( 23, 24,scalar[i], 37);
	- addToList( 23,-24,scalar[i],-37);
	+ addToList( 22, 24,scalar[i],-37);
	+ addToList( 22,-24,scalar[i], 37);
	+ addToList( 23, 24,scalar[i],-37);
	+ addToList( 23,-24,scalar[i], 37);
	}
	// pair of pseudoscalars
	for(unsigned int i=0;i<2;++i) {
	for(unsigned int j=0;j<2;++j) {
	addToList( 24,-24,pseudo[i],pseudo[j]);
	addToList( 23, 23,pseudo[i],pseudo[j]);
	}
	// pseudo charged
	- addToList( 22, 24,pseudo[i], 37);
	- addToList( 22,-24,pseudo[i],-37);
	- addToList( 23, 24,pseudo[i], 37);
	- addToList( 23,-24,pseudo[i],-37);
	+ addToList( 22, 24,pseudo[i],-37);
	+ addToList( 22,-24,pseudo[i], 37);
	+ addToList( 23, 24,pseudo[i],-37);
	+ addToList( 23,-24,pseudo[i], 37);
	}
	addToList( 24,-24, 37,-37);
	addToList( 23, 23, 37,-37);
	addToList( 22, 23, 37,-37);
	addToList( 22, 22, 37,-37);
	// cast to NMSSM model
	tcNMSSMPtr model=dynamic_ptr_cast<tcNMSSMPtr>(generator()->standardModel());
	if(!model)
	throw InitException() << "Must have the NMSSM Model in NMSSMWWHHVertex::doinit()"
	<< Exception::runerror;
	// sin and cos theta_W
	sw_ = sqrt( sin2ThetaW());
	cw_ = sqrt(1.-sin2ThetaW());
	// get the mixing matrices
	mixS_ = model->CPevenHiggsMix();
	if(!mixS_) throw InitException() << "Mixing matrix for CP-even neutral Higgs"
	<< " bosons is not set in NMSSMWWHHVertex::doinit()"
	<< Exception::runerror;
	mixP_ = model->CPoddHiggsMix();
	if(!mixP_) throw InitException() << "Mixing matrix for CP-odd neutral Higgs"
	<< " bosons is not set in NMSSMWWHHVertex::doinit()"
	<< Exception::runerror;
	// sin and cos beta
	double beta = atan(model->tanBeta());
	sb_ = sin(beta);
	cb_ = cos(beta);
	// order in couplings
	orderInGem(2);
	orderInGs (0);
	// base class
	VVSSVertex::doinit();
	}

	void NMSSMWWHHVertex::setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,
	tcPDPtr part3, tcPDPtr part4) {
	if(q2!=q2last_\|\|couplast_==0.) {
	couplast_ = sqr(electroMagneticCoupling(q2));
	q2last_=q2;
	}
	int ibos1 = part1->id(), ibos2 = part2->id();
	int isca1 = part3->id(), isca2 = part4->id();
	Complex fact(1.);
	if(abs(ibos1)==abs(ibos2)) {
	fact *= 0.5/sqr(sw_);
	if(ibos1==ParticleID::Z0) fact /= sqr(cw_);
	// charged
	if(abs(isca1)==37) {
	if(ibos1==ParticleID::Z0) fact *= sqr(sqr(cw_)-sqr(sw_));
	else if(ibos1==ParticleID::gamma) fact = 2.;
	}
	// pair of scalars
	else if((isca1-5)%10==0) {
	unsigned int i = (isca1-25)/10;
	unsigned int j = (isca2-25)/10;
	fact = (mixS_)(i,0)(mixS_)(j,0)+(mixS_)(i,1)(*mixS_)(j,1);
	}
	// pair of pseudoscalars
	else if((isca1-6)%10==0) {
	unsigned int i = (isca1-36)/10;
	unsigned int j = (isca2-36)/10;
	fact = (mixP_)(i,0)(mixP_)(j,0)+(mixP_)(i,1)(*mixP_)(j,1);
	}
	else
	assert(false);
	}
	else if(abs(ibos1)==ParticleID::Wplus \|\|
	abs(ibos2)==ParticleID::Wplus) {
	if(abs(ibos1)==ParticleID::gamma \|\|
	abs(ibos2)==ParticleID::gamma) {
	fact *= -0.5/sw_;
	}
	else {
	fact *= 0.5/cw_;
	}
	if((isca1-5)%10==0) {
	unsigned int i = (isca1-25)/10;
	fact = sb_(mixS_)(i,0) - cb_(*mixS_)(i,1);
	}
	else if((isca2-5)%10==0) {
	unsigned int i = (isca2-25)/10;
	fact = sb_(mixS_)(i,0) - cb_(*mixS_)(i,1);
	}
	else if((isca1-6)%10==0) {
	unsigned int i = (isca1-36)/10;
	fact = sb_(mixP_)(i,0) + cb_(*mixP_)(i,1);
	fact *= isca2==ParticleID::Hplus ? -Complex(0.,1.) : Complex(0.,1.);
	}
	else if((isca2-6)%10==0) {
	unsigned int i = (isca2-36)/10;
	fact = sb_(mixP_)(i,0) + cb_(*mixP_)(i,1);
	fact *= isca1==ParticleID::Hplus ? -Complex(0.,1.) : Complex(0.,1.);
	}
	else
	assert(false);
	}
	else {
	fact = (sqr(cw_)-sqr(sw_))/cw_/sw_;
	}
	// set the coupling
	norm(couplast_*fact);
	}
	diff --git a/Models/Susy/SSWWHHVertex.cc b/Models/Susy/SSWWHHVertex.cc
	--- a/Models/Susy/SSWWHHVertex.cc
	+++ b/Models/Susy/SSWWHHVertex.cc
	@@ -1,177 +1,177 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the SSWWHHVertex class.
	//

	#include "SSWWHHVertex.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "MSSM.h"

	using namespace Herwig;

	SSWWHHVertex::SSWWHHVertex() : couplast_(0.), q2last_(ZERO)
	{}

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

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

	void SSWWHHVertex::persistentOutput(PersistentOStream & os) const {
	os << coup_;
	}

	void SSWWHHVertex::persistentInput(PersistentIStream & is, int) {
	is >> coup_;
	}

	ClassDescription<SSWWHHVertex> SSWWHHVertex::initSSWWHHVertex;
	// Definition of the static class description member.

	void SSWWHHVertex::Init() {

	static ClassDocumentation<SSWWHHVertex> documentation
	("The SSWWHHVertex class implements the coupling of two Higgs bosons and"
	"two electroweak vector bosons in the MSSM.");

	}

	void SSWWHHVertex::doinit() {
	int id[3]={25,35,36};
	for(unsigned int ix=0;ix<3;++ix) {
	addToList( 24,-24,id[ix],id[ix]);
	addToList( 23, 23,id[ix],id[ix]);
	- addToList( 22, 24,id[ix], 37);
	- addToList( 22,-24,id[ix],-37);
	- addToList( 23, 24,id[ix], 37);
	- addToList( 23,-24,id[ix],-37);
	+ addToList( 22, 24,id[ix],-37);
	+ addToList( 22,-24,id[ix], 37);
	+ addToList( 23, 24,id[ix],-37);
	+ addToList( 23,-24,id[ix], 37);
	}
	addToList( 24,-24, 37,-37);
	addToList( 23, 23, 37,-37);
	addToList( 22, 23, 37,-37);
	addToList( 22, 22, 37,-37);
	VVSSVertex::doinit();
	tMSSMPtr model = dynamic_ptr_cast<tMSSMPtr>(generator()->standardModel());
	if( !model )
	throw Exception()
	<< "SSWWHHVertex::doinit() - The pointer to the MSSM object is null!"
	<< Exception::abortnow;
	coup_.resize(11);
	double sw2 = sin2ThetaW();
	double sw = sqrt(sw2);
	double cw2 = 1.-sw2;
	double cw = sqrt(cw2);
	double c2w = cw2-sw2;
	double sinalp = sin(model->higgsMixingAngle());
	double cosalp = sqrt(1. - sqr(sinalp));
	double tanbeta = model->tanBeta();
	double sinbeta = tanbeta/sqrt(1. + sqr(tanbeta));
	double cosbeta = sqrt( 1. - sqr(sinbeta) );
	double sinbma = sinbetacosalp - cosbetasinalp;
	double cosbma = cosbetacosalp + sinbetasinalp;
	// WWHH
	coup_[0] = 0.5/sw2;
	// ZZH0H0
	coup_[1] = 0.5/sw2/cw2;
	// ZZH+H-
	coup_[2] = 0.5*sqr(c2w)/cw2/sw2;
	// Z W h0 H+
	coup_[3] =-0.5/cw*cosbma;
	// Z W H0 H+
	coup_[4] = 0.5/cw*sinbma;
	// Z W A0 H+
	coup_[5] =-Complex(0.,0.5)/cw;
	// A A H+H-
	coup_[6] = 2.;
	// A Z H+H-
	coup_[7] = c2w/sw/cw;
	// A W h0 H+
	coup_[8] = 0.5*cosbma/sw;
	// A W H0 H+
	coup_[9] =-0.5*sinbma/sw;
	// A W A0 H+
	coup_[10] = Complex(0.,0.5)/sw;
	}

	void SSWWHHVertex::setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,
	tcPDPtr part3, tcPDPtr part4) {
	if(q2!=q2last_\|\|couplast_==0.) {
	couplast_ = sqr(electroMagneticCoupling(q2));
	q2last_=q2;
	}
	int ibos1 = part1->id(), ibos2 = part2->id();
	int isca1 = part3->id(), isca2 = part4->id();
	if(abs(ibos1)==abs(ibos2)) {
	if(abs(ibos1)==ParticleID::Wplus) {
	norm(couplast_*coup_[0]);
	}
	else if(ibos1==ParticleID::Z0) {
	if(abs(isca1)==ParticleID::Hplus)
	norm(couplast_*coup_[2]);
	else
	norm(couplast_*coup_[1]);
	}
	else if(ibos1==ParticleID::gamma) {
	norm(couplast_*coup_[6]);
	}
	else
	assert(false);
	}
	else if(abs(ibos1)==ParticleID::Wplus \|\|
	abs(ibos2)==ParticleID::Wplus) {
	if(abs(ibos1)==ParticleID::gamma \|\|
	abs(ibos2)==ParticleID::gamma) {
	if(abs(isca1)==ParticleID::h0 \|\|
	abs(isca2)==ParticleID::h0) {
	norm(couplast_*coup_[8]);
	}
	else if(abs(isca1)==ParticleID::H0 \|\|
	abs(isca2)==ParticleID::H0) {
	norm(couplast_*coup_[9]);
	}
	else if(abs(isca1)==ParticleID::A0 \|\|
	abs(isca2)==ParticleID::A0) {
	if(isca1==ParticleID::Hplus \|\|
	isca2==ParticleID::Hplus) {
	norm(couplast_* coup_[10] );
	}
	else {
	norm(couplast_*conj(coup_[10]));
	}
	}
	else
	assert(false);
	}
	else {
	if(abs(isca1)==ParticleID::h0 \|\|
	abs(isca2)==ParticleID::h0) {
	norm(couplast_*coup_[3]);
	}
	else if(abs(isca1)==ParticleID::H0 \|\|
	abs(isca2)==ParticleID::H0) {
	norm(couplast_*coup_[4]);
	}
	else if(abs(isca1)==ParticleID::A0 \|\|
	abs(isca2)==ParticleID::A0) {
	if(isca1==ParticleID::Hplus \|\|
	isca2==ParticleID::Hplus) {
	norm(couplast_* coup_[5] );
	}
	else {
	norm(couplast_*conj(coup_[5]));
	}
	}
	else
	assert(false);
	}
	}
	else {
	norm(couplast_*coup_[7]);
	}
	}
	diff --git a/Models/UED/UEDF1F1W0Vertex.cc b/Models/UED/UEDF1F1W0Vertex.cc
	--- a/Models/UED/UEDF1F1W0Vertex.cc
	+++ b/Models/UED/UEDF1F1W0Vertex.cc
	@@ -1,155 +1,155 @@
	// -- C++ --
	//
	// UEDF1F1W0Vertex.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 UEDF1F1W0Vertex class.
	//

	#include "UEDF1F1W0Vertex.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDT/EnumParticles.h"

	using namespace ThePEG::Helicity;
	using namespace Herwig;

	UEDF1F1W0Vertex::UEDF1F1W0Vertex(): theRadius(ZERO), theQ2Last(ZERO),
	theCoupLast(0.),
	thefermALast(0), thefermBLast(0)
	{}

	void UEDF1F1W0Vertex::doinit() {
	//outgoing W+
	for( long i = 2; i < 17; i += 2 ) {
	if( i == 7 ) i += 5;
	addToList(-5100000 - i, 5100000 + i - 1, 24);
	if( i < 7 ) {
	addToList(-6100000 - i, 6100000 + i - 1, 24);
	}
	}
	addToList(-6100006, 5100005, 24);
	- addToList(-5100005, 6100006, 24);
	+ addToList(-5100006, 6100005, 24);
	//outgoing W-
	for( long i = 1; i < 16; i += 2 ) {
	if( i == 6 ) i += 5;
	addToList(-5100000 - i, 5100001 + i, -24);
	if( i < 6 ) {
	addToList(-6100000 - i, 6100001 + i, -24);
	}
	}
	addToList(-6100005, 5100006, -24);
	addToList(-5100005, 6100006, -24);
	FFVVertex::doinit();
	tUEDBasePtr model = dynamic_ptr_cast<tUEDBasePtr>(generator()->standardModel());
	if(!model)
	throw InitException() << "UEDF1F1W0Vertex::doinit() - The pointer to "
	<< "the UEDBase object is null!"
	<< Exception::runerror;
	theRadius = model->compactRadius();
	orderInGs(0);
	orderInGem(0);
	}

	void UEDF1F1W0Vertex::persistentOutput(PersistentOStream & os) const {
	os << ounit(theRadius,1/GeV);
	}

	void UEDF1F1W0Vertex::persistentInput(PersistentIStream & is, int) {
	is >> iunit(theRadius,1/GeV);
	}

	ClassDescription<UEDF1F1W0Vertex> UEDF1F1W0Vertex::initUEDF1F1W0Vertex;
	// Definition of the static class description member.

	void UEDF1F1W0Vertex::Init() {

	static ClassDocumentation<UEDF1F1W0Vertex> documentation
	("This class implements the coupling of a pair of level-1 KK fermions"
	"to an SM W boson");

	}

	void UEDF1F1W0Vertex::setCoupling(Energy2 q2, tcPDPtr part1, tcPDPtr part2,
	#ifndef NDEBUG
	tcPDPtr part3) {
	#else
	tcPDPtr) {
	#endif
	long ianti(abs(part1->id())), iferm(abs(part2->id()));
	assert( abs(part3->id()) == 24 );
	bool ferma = (iferm >= 5100001 && iferm <= 5100006) \|\|
	(iferm >= 6100001 && iferm <= 6100006) \|\|
	(iferm >= 5100011 && iferm <= 5100016) \|\|
	(iferm >= 6100011 && iferm <= 6100016);
	bool fermb = (ianti >= 5100001 && ianti <= 5100006) \|\|
	(ianti >= 6100001 && ianti <= 6100006) \|\|
	(ianti >= 5100011 && ianti <= 5100016) \|\|
	(ianti >= 6100011 && ianti <= 6100016);
	if( !ferma \|\| !fermb )
	throw HelicityLogicalError() << "UEDF1F1W0Vertex::setCoupling - "
	<< "There is an unknown particle(s) in the "
	<< "UED F^(1) F^(1) W^(0) vertex. ID: "
	<< ianti << " " << iferm
	<< Exception::runerror;
	if(q2 != theQ2Last \|\| theCoupLast == 0. ) {
	theQ2Last = q2;
	theCoupLast = sqrt(0.5)*weakCoupling(q2);
	}
	if(iferm != thefermALast \|\| ianti != thefermBLast) {
	thefermALast = iferm;
	thefermBLast = ianti;
	int stateA(ianti/1000000), stateB(iferm/1000000);
	long sma = (stateA == 6) ? ianti - 6100000 : ianti - 5100000;
	long smb = (stateB == 6) ? iferm - 6100000 : iferm - 5100000;
	double afu(0.), afd(0.);
	if( sma % 2 == 0 ) {
	afu = atan(getParticleData(sma)->mass()*theRadius)/2.;
	afd = atan(getParticleData(smb)->mass()*theRadius)/2.;
	}
	else {
	afd = atan(getParticleData(sma)->mass()*theRadius)/2.;
	afu = atan(getParticleData(smb)->mass()*theRadius)/2.;

	}
	if( stateA == stateB ) {
	if( stateA == 5 ) {
	left(cos(afu)*cos(afd));
	right(cos(afu)*cos(afd));
	}
	else {
	left(sin(afu)*sin(afd));
	right(sin(afu)*sin(afd));
	}
	}
	else {
	if( sma % 2 == 0 ) {
	if( stateA == 5 ) {
	left(cos(afu)*sin(afd));
	right(-cos(afu)*sin(afd));
	}
	else {
	left(sin(afu)*cos(afd));
	right(-sin(afu)*cos(afd));
	}
	}
	else {
	if( stateA == 5 ) {
	left(sin(afu)*cos(afd));
	right(-sin(afu)*cos(afd));
	}
	else {
	left(cos(afu)*sin(afd));
	right(-cos(afu)*sin(afd));
	}
	}
	}
	}
	norm(theCoupLast);
	}
	diff --git a/Tests/Makefile.am b/Tests/Makefile.am
	--- a/Tests/Makefile.am
	+++ b/Tests/Makefile.am
	@@ -1,230 +1,230 @@
	AUTOMAKE_OPTIONS = -Wno-portability
	AM_LDFLAGS += -module -avoid-version -rpath /dummy/path/not/used

	EXTRA_DIST = Inputs python Rivet

	dist-hook:
	rm -rf $(distdir)/Inputs/.svn
	rm -rf $(distdir)/python/.svn
	rm -rf $(distdir)/Rivet/.svn

	EXTRA_LTLIBRARIES = LeptonTest.la GammaTest.la HadronTest.la DISTest.la

	if WANT_LIBFASTJET
	EXTRA_LTLIBRARIES += HadronJetTest.la LeptonJetTest.la
	HadronJetTest_la_SOURCES = \
	Hadron/VHTest.h Hadron/VHTest.cc\
	Hadron/VTest.h Hadron/VTest.cc\
	Hadron/HTest.h Hadron/HTest.cc
	HadronJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
	-I$(FASTJETPATH)
	HadronJetTest_la_LIBADD = $(FASTJETLIBS)
	LeptonJetTest_la_SOURCES = \
	Lepton/TopDecay.h Lepton/TopDecay.cc
	LeptonJetTest_la_CPPFLAGS = $(AM_CPPFLAGS) $(FASTJETINCLUDE) \
	-I$(FASTJETPATH)
	LeptonJetTest_la_LIBADD = $(FASTJETLIBS)
	endif

	LeptonTest_la_SOURCES = \
	Lepton/VVTest.h Lepton/VVTest.cc \
	Lepton/VBFTest.h Lepton/VBFTest.cc \
	Lepton/VHTest.h Lepton/VHTest.cc \
	Lepton/FermionTest.h Lepton/FermionTest.cc

	GammaTest_la_SOURCES = \
	Gamma/GammaMETest.h Gamma/GammaMETest.cc \
	Gamma/GammaPMETest.h Gamma/GammaPMETest.cc

	DISTest_la_SOURCES = \
	DIS/DISTest.h DIS/DISTest.cc

	HadronTest_la_SOURCES = \
	Hadron/HadronVVTest.h Hadron/HadronVVTest.cc\
	Hadron/HadronVBFTest.h Hadron/HadronVBFTest.cc\
	Hadron/WHTest.h Hadron/WHTest.cc\
	Hadron/ZHTest.h Hadron/ZHTest.cc\
	Hadron/VGammaTest.h Hadron/VGammaTest.cc\
	Hadron/ZJetTest.h Hadron/ZJetTest.cc\
	Hadron/WJetTest.h Hadron/WJetTest.cc\
	Hadron/QQHTest.h Hadron/QQHTest.cc


	REPO = $(top_builddir)/src/HerwigDefaults.rpo
	HERWIG = $(top_builddir)/src/Herwig++
	HWREAD = $(HERWIG) read -r $(REPO) -L $(builddir)/.libs
	HWRUN = $(HERWIG) run

	tests : tests-LEP tests-DIS tests-LHC tests-Gamma

	if WANT_LIBFASTJET
	tests-LEP : test-LEP-VV test-LEP-VH test-LEP-VBF test-LEP-BB test-LEP-Quarks test-LEP-Leptons \
	test-LEP-default test-LEP-Powheg test-LEP-TopDecay
	else
	tests-LEP : test-LEP-VV test-LEP-VH test-LEP-VBF test-LEP-BB test-LEP-Quarks test-LEP-Leptons
	endif

	tests-DIS : test-DIS-Charged test-DIS-Neutral

	if WANT_LIBFASTJET
	tests-LHC : test-LHC-WW test-LHC-WZ test-LHC-ZZ test-LHC-ZGamma test-LHC-WGamma \
	test-LHC-ZH test-LHC-WH test-LHC-ZJet test-LHC-WJet test-LHC-Z test-LHC-W test-LHC-ZZVBF test-LHC-VBF \
	test-LHC-WWVBF test-LHC-bbH test-LHC-ttH test-LHC-GammaGamma test-LHC-GammaJet test-LHC-Higgs \
	test-LHC-HiggsJet test-LHC-QCDFast test-LHC-QCD test-LHC-Top test-LHC-Bottom \
	test-LHC-WHJet test-LHC-ZHJet test-LHC-HJet test-LHC-ZShower test-LHC-WShower\
	test-LHC-WHJet-Powheg test-LHC-ZHJet-Powheg test-LHC-HJet-Powheg \
	test-LHC-ZShower-Powheg test-LHC-WShower-Powheg
	else
	tests-LHC : test-LHC-WW test-LHC-WZ test-LHC-ZZ test-LHC-ZGamma test-LHC-WGamma \
	test-LHC-ZH test-LHC-WH test-LHC-ZJet test-LHC-WJet test-LHC-Z test-LHC-W test-LHC-ZZVBF test-LHC-VBF \
	test-LHC-WWVBF test-LHC-bbH test-LHC-ttH test-LHC-GammaGamma test-LHC-GammaJet test-LHC-Higgs \
	test-LHC-HiggsJet test-LHC-QCDFast test-LHC-QCD test-LHC-Top
	endif

	tests-Gamma : test-Gamma-FF test-Gamma-WW test-Gamma-P

	if WANT_LIBFASTJET
	test-LEP-% : Inputs/LEP-%.in LeptonTest.la LeptonJetTest.la
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
	else
	test-LEP-% : Inputs/LEP-%.in LeptonTest.la
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
	endif

	Rivet-LEP-% : Rivet/LEP-%.in
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}

	Rivet-TVT-% : Rivet/TVT-%.in
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
	Rivet-TVT-% : Rivet/TVT-%.in
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}

	Rivet-LHC-% : Rivet/LHC-%.in
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}

	Rivet-Star-% : Rivet/Star-%.in
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}

	Rivet-LEP: Rivet-LEP-22 Rivet-LEP-35 Rivet-LEP-44 Rivet-LEP-91 \
	Rivet-LEP-133 Rivet-LEP-161 Rivet-LEP-172 Rivet-LEP-177 \
	Rivet-LEP-183 Rivet-LEP-189 Rivet-LEP-196 Rivet-LEP-197 \
	Rivet-LEP-200 Rivet-LEP-206 Rivet-LEP-14 \
	- Rivet-LEP-Powheg-14 Rivet-LEP-Powheg-22 \
	+ Rivet-LEP-Powheg-14 Rivet-LEP-Powheg-22 \
	Rivet-LEP-Powheg-35 Rivet-LEP-Powheg-44 \
	Rivet-LEP-Powheg-91 Rivet-LEP-Powheg-133 \
	Rivet-LEP-Powheg-161 Rivet-LEP-Powheg-172 \
	Rivet-LEP-Powheg-177 Rivet-LEP-Powheg-183 \
	Rivet-LEP-Powheg-189 Rivet-LEP-Powheg-196 \
	Rivet-LEP-Powheg-197 Rivet-LEP-Powheg-200 \
	Rivet-LEP-Powheg-206
	for i in LEP-*.aida; do rivet-rmgaps $$i; done;
	rivet-rmgaps LEP-91.aida;
	rivet-rmgaps LEP-Powheg-91.aida
	rm -rf Rivet-LEP
	python/merge-LEP LEP
	python/merge-LEP LEP-Powheg
	rivet-mkhtml -o Rivet-LEP LEP.aida:Hw++ LEP-Powheg.aida:Hw++-Powheg

	Rivet-TVT-WZ: Rivet-TVT-Run-I-Z Rivet-TVT-Powheg-Run-I-Z \
	Rivet-TVT-Run-I-W Rivet-TVT-Powheg-Run-I-W \
	Rivet-TVT-Run-I-WZ Rivet-TVT-Powheg-Run-I-WZ\
	Rivet-TVT-Run-II-Z-e Rivet-TVT-Powheg-Run-II-Z-e \
	Rivet-TVT-Run-II-Z-mu Rivet-TVT-Powheg-Run-II-Z-mu \
	Rivet-TVT-Run-II-W Rivet-TVT-Powheg-Run-II-W
	rivet-rmgaps TVT-Run-II-Z-e.aida;
	rivet-rmgaps TVT-Powheg-Run-II-Z-e.aida;
	rm -rf Rivet-TVT-WZ
	python/merge-TVT-EW TVT-Run-II-W.aida TVT-Run-II-Z-{e,mu}.aida\
	TVT-Run-I-{W,Z,WZ}.aida -o TVT-WZ.aida
	python/merge-TVT-EW TVT-Powheg-Run-II-W.aida TVT-Powheg-Run-II-Z-{e,mu}.aida\
	TVT-Powheg-Run-I-{W,Z,WZ}.aida -o TVT-Powheg-WZ.aida
	rivet-mkhtml -o Rivet-TVT-WZ TVT-WZ.aida:Hw++ TVT-Powheg-WZ.aida:Hw++-Powheg

	Rivet-TVT-Photon: Rivet-TVT-Run-II-DiPhoton Rivet-TVT-Run-II-PromptPhoton
	# Rivet-TVT-Run-I-PromptPhoton
	rm -rf Rivet-TVT-Photon
	python/merge-aida TVT-Run-II-DiPhoton.aida TVT-Run-II-PromptPhoton.aida\
	-o TVT-Photon.aida
	rivet-mkhtml -o Rivet-TVT-Photon TVT-Photon.aida:Hw++


	Rivet-TVT-Jets: Rivet-TVT-Run-II-Jets-1 Rivet-TVT-Run-II-Jets-2 \
	Rivet-TVT-Run-II-Jets-3 Rivet-TVT-Run-II-Jets-4 \
	Rivet-TVT-Run-II-Jets-5 Rivet-TVT-Run-II-Jets-6 \
	Rivet-TVT-Run-II-Jets-7 Rivet-TVT-Run-II-Jets-8 \
	Rivet-TVT-Run-II-Jets-9 Rivet-TVT-Run-II-Jets-10\
	Rivet-TVT-Run-II-Jets-11 Rivet-TVT-Run-II-UE \
	Rivet-TVT-Run-I-Jets-1 Rivet-TVT-Run-I-Jets-2 \
	Rivet-TVT-Run-I-Jets-3 Rivet-TVT-Run-I-Jets-4 \
	Rivet-TVT-Run-I-Jets-5 Rivet-TVT-Run-I-Jets-6 \
	Rivet-TVT-Run-I-Jets-7 Rivet-TVT-Run-I-Jets-8\
	Rivet-TVT-Run-I-UE\
	Rivet-TVT-630-UE Rivet-TVT-630-Jets-1 \
	Rivet-TVT-630-Jets-2 Rivet-TVT-630-Jets-3
	rivet-rmgaps TVT-Run-I-Jets-4.aida
	rm -rf Rivet-TVT-Jets
	python/merge-TVT-Jets TVT
	rivet-mkhtml -o Rivet-TVT-Jets TVT-Jets.aida:Hw++

	Rivet-LHC-Jets: Rivet-LHC-7-Jets-1 Rivet-LHC-7-Jets-2 \
	Rivet-LHC-7-Jets-3 Rivet-LHC-7-Jets-4 \
	Rivet-LHC-7-Jets-5 Rivet-LHC-7-Jets-6 \
	Rivet-LHC-7-Jets-7 Rivet-LHC-7-Jets-8 \
	Rivet-LHC-7-Jets-9 Rivet-LHC-7-Jets-10 \
	Rivet-LHC-7-Jets-11 Rivet-LHC-7-Jets-12 \
	Rivet-LHC-7-Jets-13 Rivet-LHC-7-UE \
	Rivet-LHC-2360-UE Rivet-LHC-900-UE
	rm -rf Rivet-LHC-Jets
	python/merge-LHC-Jets
	rivet-mkhtml -o Rivet-LHC-Jets LHC-Jets.aida:Hw++

	Rivet-Star: Rivet-Star-UE Rivet-Star-Jets-1 \
	Rivet-Star-Jets-2 Rivet-Star-Jets-3 \
	Rivet-Star-Jets-4
	rm -rf Rivet-Star
	rivet-rmgaps Star-UE.aida
	python/merge-Star Star
	rivet-mkhtml -o Rivet-Star Star.aida

	Rivet-LHC-EW: Rivet-LHC-W-e Rivet-LHC-Powheg-W-e \
	Rivet-LHC-W-mu Rivet-LHC-Powheg-W-mu \
	Rivet-LHC-Z Rivet-LHC-Powheg-Z \
	Rivet-LHC-WW Rivet-LHC-Powheg-WW\
	Rivet-LHC-ZZ Rivet-LHC-Powheg-ZZ
	rm -rf Rivet-LHC-EW;
	python/merge-LHC-EW LHC-{W-e,W-mu,Z,WW,ZZ}.aida -o LHC-EW.aida;
	python/merge-LHC-EW LHC-Powheg-{W-e,W-mu,Z,WW,ZZ}.aida -o LHC-Powheg-EW.aida;
	rivet-mkhtml -o Rivet-LHC-EW LHC-EW.aida:Hw++ LHC-Powheg-EW.aida:Hw++-Powheg;

	Rivet-LHC-Photon: Rivet-LHC-7-PromptPhoton
	rm -rf Rivet-LHC-Photon
	python/merge-aida LHC-7-PromptPhoton.aida -o LHC-Photon.aida
	rivet-mkhtml -o Rivet-LHC-Photon LHC-Photon.aida:Hw++

	test-Gamma-% : Inputs/Gamma-%.in GammaTest.la
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}

	test-DIS-% : Inputs/DIS-%.in DISTest.la
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}

	if WANT_LIBFASTJET
	test-LHC-% : Inputs/LHC-%.in HadronTest.la GammaTest.la HadronJetTest.la
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
	else
	test-LHC-% : Inputs/LHC-%.in HadronTest.la GammaTest.la
	$(HWREAD) $<
	$(HWRUN) $(notdir $(subst .in,.run,$<)) -N $${NUMEVENTS:-10000}
	endif

	clean-local:
	rm -f .out .log .tex .top .run .dump .mult .Bmult *.aida
	diff --git a/src/LHC-TRP.in b/src/LHC-TRP.in
	--- a/src/LHC-TRP.in
	+++ b/src/LHC-TRP.in
	@@ -1,69 +1,69 @@
	##################################################
	# Example generator based on LHC parameters
	# usage: Herwig++ read LHC.in
	##################################################

	cd /Herwig/Particles
	create ThePEG::ParticleData graviton
	setup graviton 39 graviton 0.0 0.0 0.0 0.0 0 0 0 1

	cd /Herwig/EventHandlers

	#set LHCHandler:CascadeHandler NULL
	#set LHCHandler:HadronizationHandler NULL
	#set LHCHandler:DecayHandler NULL

	##################################################
	# Technical parameters for this run
	##################################################
	cd /Herwig


	cd /Herwig/Generators
	set LHCGenerator:NumberOfEvents 1000000
	set LHCGenerator:RandomNumberGenerator:Seed 31122001
	set LHCGenerator:PrintEvent 10
	set LHCGenerator:MaxErrors 10000

	##################################################
	# LHC physics parameters (override defaults here)
	##################################################
	set LHCGenerator:EventHandler:LuminosityFunction:Energy 7000.0

	# Intrinsic pT tune extrapolated to LHC energy
	set /Herwig/Shower/Evolver:IntrinsicPtGaussian 2.2*GeV

	cd /Herwig/Cuts
	#set JetKtCut:MinKT 100.0*GeV
	-set QCDCuts:MHatMin 9000*GeV
	+set QCDCuts:MHatMin 4500*GeV


	##################################################
	# Matrix Elements for hadron-hadron collisions
	##################################################
	cd /Herwig/MatrixElements/
	create Herwig::METRP2to2 METransplanck HwTransplanck.so

	insert SimpleQCD:MatrixElements[0] METransplanck
	#set METransplanck:Process 6
	set METransplanck:NumberExtraDimensions 6
	set METransplanck:PlanckMass 1500

	cd /Herwig/Generators

	##################################################
	# Useful analysis handlers for HepMC related output
	##################################################
	# Schematic overview of an event (requires --with-hepmc to be set at configure time
	# and the graphviz program 'dot' to produce a plot)
	# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/Plot
	# A HepMC dump file (requires --with-hepmc to be set at configure time)
	# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
	# set /Herwig/Analysis/HepMCFile:PrintEvent 100
	# set /Herwig/Analysis/HepMCFile:Format GenEvent
	# set /Herwig/Analysis/HepMCFile:Units GeV_mm

	##################################################
	# Save run for later usage with 'Herwig++ run'
	##################################################
	saverun LHC-TRP LHCGenerator
	diff --git a/src/Leptoquark.model b/src/Leptoquark.model
	--- a/src/Leptoquark.model
	+++ b/src/Leptoquark.model
	@@ -1,238 +1,238 @@
	##################################################
	# Common setup for Leptoquark models
	#
	# See LHC-LQ.in or ILC-LQ.in for example usage
	#
	# This file does not contain anything that
	# users need to touch. User settings are in
	# ???-LQ.in
	#
	###################################################
	#
	# Particle Data object for the Leptoquarks
	#
	###################################################
	cd /Herwig/Particles

	create /ThePEG/ParticleData S0bar
	setup S0bar 9911561 S0bar 400.0 0.0 0.0 0.0 -1 3 1 0
	create /ThePEG/ParticleData S0
	setup S0 -9911561 S0 400.0 0.0 0.0 0.0 1 -3 1 0
	makeanti S0bar S0

	create /ThePEG/ParticleData ~S0bar
	setup ~S0bar 9921551 ~S0bar 400.0 0.0 0.0 0.0 -4 3 1 0
	create /ThePEG/ParticleData ~S0
	setup ~S0 -9921551 ~S0 400.0 0.0 0.0 0.0 4 -3 1 0
	makeanti ~S0bar ~S0

	create /ThePEG/ParticleData S1pbar
	setup S1pbar 9931551 S1pbar 400.0 0.0 0.0 0.0 -4 3 1 0
	create /ThePEG/ParticleData S1p
	setup S1p -9931551 S1p 400.0 0.0 0.0 0.0 4 -3 1 0
	makeanti S1pbar S1p

	create /ThePEG/ParticleData S1mbar
	setup S1mbar 9931661 S1mbar 400.0 0.0 0.0 0.0 2 3 1 0
	create /ThePEG/ParticleData S1m
	setup S1m -9931661 S1m 400.0 0.0 0.0 0.0 -2 -3 1 0
	makeanti S1mbar S1m

	create /ThePEG/ParticleData S1zbar
	setup S1zbar 9931561 S1zbar 400.0 0.0 0.0 0.0 -1 3 1 0
	create /ThePEG/ParticleData S1z
	setup S1z -9931561 S1z 400.0 0.0 0.0 0.0 1 -3 1 0
	makeanti S1zbar S1z

	create /ThePEG/ParticleData S12p
	setup S12p 9941561 S12p 400.0 0.0 0.0 0.0 5 3 1 0
	create /ThePEG/ParticleData S12pbar
	setup S12pbar -9941561 S12pbar 400.0 0.0 0.0 0.0 -5 -3 1 0
	makeanti S12p S12pbar

	create /ThePEG/ParticleData S12m
	setup S12m 9941551 S12m 400.0 0.0 0.0 0.0 2 3 1 0
	create /ThePEG/ParticleData S12mbar
	setup S12mbar -9941551 S12mbar 400.0 0.0 0.0 0.0 -2 -3 1 0
	makeanti S12m S12mbar

	create /ThePEG/ParticleData S12tp
	setup S12tp 9951551 S12tp 400.0 0.0 0.0 0.0 2 3 1 0
	create /ThePEG/ParticleData S12tpbar
	setup S12tpbar -9951551 S12tpbar 400.0 0.0 0.0 0.0 -2 -3 1 0
	makeanti S12tp S12tpbar

	create /ThePEG/ParticleData S12tm
	setup S12tm 9951651 S12tm 400.0 0.0 0.0 0.0 -1 3 1 0
	create /ThePEG/ParticleData S12tmbar
	setup S12tmbar -9951651 S12tmbar 400.0 0.0 0.0 0.0 1 -3 1 0
	makeanti S12tm S12tmbar

	create /ThePEG/ParticleData dS0
	setup dS0 9961551 dS0 400.0 0.0 0.0 0.0 2 3 1 0
	create /ThePEG/ParticleData dS0bar
	setup dS0bar -9961551 dS0bar 400.0 0.0 0.0 0.0 -2 -3 1 0
	makeanti dS0 dS0bar

	create /ThePEG/ParticleData ~dS0
	setup ~dS0 9971561 ~dS0 400.0 0.0 0.0 0.0 5 3 1 0
	create /ThePEG/ParticleData ~dS0bar
	setup ~dS0bar -9971561 ~dS0bar 400.0 0.0 0.0 0.0 -5 -3 1 0
	makeanti ~dS0 ~dS0bar

	create /ThePEG/ParticleData dS1p
	setup dS1p 9981561 dS1p 400.0 0.0 0.0 0.0 5 3 1 0
	create /ThePEG/ParticleData dS1pbar
	setup dS1pbar -9981561 dS1pbar 400.0 0.0 0.0 0.0 -5 -3 1 0
	makeanti dS1p dS1pbar

	create /ThePEG/ParticleData dS1m
	setup dS1m 9981651 dS1m 400.0 0.0 0.0 0.0 -1 3 1 0
	create /ThePEG/ParticleData dS1mbar
	setup dS1mbar -9981651 dS1mbar 400.0 0.0 0.0 0.0 1 -3 1 0
	makeanti dS1m dS1mbar

	create /ThePEG/ParticleData dS1z
	setup dS1z 9981551 dS1z 400.0 0.0 0.0 0.0 2 3 1 0
	create /ThePEG/ParticleData dS1zbar
	setup dS1zbar -9981551 dS1zbar 400.0 0.0 0.0 0.0 -2 -3 1 0
	makeanti dS1z dS1zbar

	create /ThePEG/ParticleData dS12pbar
	setup dS12pbar 9991551 dS12pbar 400.0 0.0 0.0 0.0 -4 3 1 0
	create /ThePEG/ParticleData dS12p
	setup dS12p -9991551 dS12p 400.0 0.0 0.0 0.0 4 -3 1 0
	makeanti dS12pbar dS12p

	create /ThePEG/ParticleData dS12mbar
	setup dS12mbar 9991561 dS12mbar 400.0 0.0 0.0 0.0 -1 3 1 0
	create /ThePEG/ParticleData dS12m
	setup dS12m -9991561 dS12m 400.0 0.0 0.0 0.0 1 -3 1 0
	makeanti dS12mbar dS12m

	create /ThePEG/ParticleData dS12tpbar
	-setup dS12tpbar 9901561 dS12tpbar 400.0 0.0 0.0 0.0 -4 3 1 0
	+setup dS12tpbar 9901561 dS12tpbar 400.0 0.0 0.0 0.0 -1 3 1 0
	create /ThePEG/ParticleData dS12tp
	-setup dS12tp -9901561 dS12tp 400.0 0.0 0.0 0.0 4 -3 1 0
	+setup dS12tp -9901561 dS12tp 400.0 0.0 0.0 0.0 1 -3 1 0
	makeanti dS12tpbar dS12tp

	create /ThePEG/ParticleData dS12tmbar
	-setup dS12tmbar 9901661 dS12tmbar 400.0 0.0 0.0 0.0 -1 3 1 0
	+setup dS12tmbar 9901661 dS12tmbar 400.0 0.0 0.0 0.0 2 3 1 0
	create /ThePEG/ParticleData dS12tm
	-setup dS12tm -9901661 dS12tm 400.0 0.0 0.0 0.0 1 -3 1 0
	+setup dS12tm -9901661 dS12tm 400.0 0.0 0.0 0.0 -2 -3 1 0
	makeanti dS12tmbar dS12tm



	###################################################
	#
	# Main directory and model object
	#
	###################################################
	mkdir /Herwig/NewPhysics/Leptoquark
	cd /Herwig/NewPhysics/Leptoquark
	create Herwig::LeptoquarkModel Model HwLeptoquarkModel.so
	# SM couplings
	set Model:QCD/RunningAlphaS /Herwig/AlphaS
	set Model:EW/RunningAlphaEM /Herwig/AlphaEM
	set Model:EW/CKM /Herwig/CKM
	set Model:RunningMass /Herwig/RunningMass
	set Model:LQCoupling 0.312
	set Model:g_S0_L 1.0
	set Model:g_S0_R 1.0
	set Model:g_S0t_R 1.0
	set Model:g_S1_L 1.0
	set Model:g_S12_L 1.0
	set Model:g_S12_R 1.0
	set Model:g_S12t_L 1.0
	set Model:g_dS0_L 1.0
	set Model:g_dS0_R 1.0
	set Model:g_dS0t_R 1.0
	set Model:g_dS1_L 1.0
	set Model:g_dS12_L 1.0
	set Model:g_dS12_R 1.0
	set Model:g_dS12t_L 1.0
	set Model:derivscale 500.0


	###################################################
	#
	# Vertices
	#
	###################################################
	# create Leptoquark model vertices
	mkdir /Herwig/Vertices/Leptoquark
	cd /Herwig/Vertices/Leptoquark
	library HwLeptoquarkModel.so
	create Herwig::LeptoquarkModelSLQSLQGVertex Leptoquark_SLQSLQGVertex
	create Herwig::LeptoquarkModelSLQSLQGGVertex Leptoquark_SLQSLQGGVertex
	create Herwig::LeptoquarkModelSLQFFVertex Leptoquark_SLQFFVertex


	cd /Herwig/NewPhysics/Leptoquark
	# SM vertices
	set Model:Vertex/FFZ /Herwig/Vertices/FFZVertex
	set Model:Vertex/FFW /Herwig/Vertices/FFWVertex
	set Model:Vertex/FFH /Herwig/Vertices/FFHVertex
	set Model:Vertex/FFG /Herwig/Vertices/FFGVertex
	set Model:Vertex/FFP /Herwig/Vertices/FFPVertex
	set Model:Vertex/GGG /Herwig/Vertices/GGGVertex
	set Model:Vertex/GGGG /Herwig/Vertices/GGGGVertex
	set Model:Vertex/WWH /Herwig/Vertices/WWHVertex
	set Model:Vertex/WWW /Herwig/Vertices/WWWVertex
	set Model:Vertex/WWWW /Herwig/Vertices/WWWWVertex
	set Model:Vertex/HGG /Herwig/Vertices/HGGVertex
	set Model:Vertex/HHH /Herwig/Vertices/HHHVertex
	set Model:Vertex/WWHH /Herwig/Vertices/WWHHVertex

	set Model:Vertex/HHH /Herwig/Vertices/HHHVertex
	set Model:Vertex/HPP /Herwig/Vertices/HPPVertex
	# Leptoquark model vertices
	set Model:Vertex/SLQSLQG /Herwig/Vertices/Leptoquark/Leptoquark_SLQSLQGVertex
	set Model:Vertex/SLQSLQGG /Herwig/Vertices/Leptoquark/Leptoquark_SLQSLQGGVertex
	set Model:Vertex/SLQFF /Herwig/Vertices/Leptoquark/Leptoquark_SLQFFVertex


	###################################################
	#
	# Set up spin correlation Decayers
	#
	###################################################
	cd /Herwig/NewPhysics

	set TwoBodyDC:CreateDecayModes Yes
	set ThreeBodyDC:CreateDecayModes No

	insert NewModel:DecayParticles 0 /Herwig/Particles/S0bar
	insert NewModel:DecayParticles 1 /Herwig/Particles/~S0bar
	insert NewModel:DecayParticles 2 /Herwig/Particles/S1pbar
	insert NewModel:DecayParticles 3 /Herwig/Particles/S1mbar
	insert NewModel:DecayParticles 4 /Herwig/Particles/S1zbar
	insert NewModel:DecayParticles 5 /Herwig/Particles/S12p
	insert NewModel:DecayParticles 6 /Herwig/Particles/S12m
	insert NewModel:DecayParticles 7 /Herwig/Particles/S12tp
	insert NewModel:DecayParticles 8 /Herwig/Particles/S12tm
	insert NewModel:DecayParticles 9 /Herwig/Particles/dS0
	insert NewModel:DecayParticles 10 /Herwig/Particles/~dS0
	insert NewModel:DecayParticles 11 /Herwig/Particles/dS1p
	insert NewModel:DecayParticles 12 /Herwig/Particles/dS1m
	insert NewModel:DecayParticles 13 /Herwig/Particles/dS1z
	insert NewModel:DecayParticles 14 /Herwig/Particles/dS12pbar
	insert NewModel:DecayParticles 15 /Herwig/Particles/dS12mbar
	insert NewModel:DecayParticles 16 /Herwig/Particles/dS12tpbar
	insert NewModel:DecayParticles 17 /Herwig/Particles/dS12tmbar

	###################################################
	# Set up the model framework
	###################################################
	set Leptoquark/Model:ModelGenerator NewModel
	###################################################
	#
	# Choose Leptoquark over SM
	#
	###################################################
	cd /Herwig/Generators
	set LEPGenerator:StandardModelParameters /Herwig/NewPhysics/Leptoquark/Model
	set LHCGenerator:StandardModelParameters /Herwig/NewPhysics/Leptoquark/Model

	#insert /Herwig/Shower/ShowerHandler:DecayInShower 0 9911561

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 9:33 PM (1 d, 1 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4242486
Default Alt Text: (70 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions