No OneTemporary
Actions

Size

65 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc b/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc
	--- a/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc
	+++ b/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc
	@@ -1,546 +1,546 @@
	// -- C++ --
	//
	// MatchboxAmplitude.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 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 MatchboxAmplitude class.
	//

	#include "MatchboxAmplitude.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/StandardModel/StandardModelBase.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/SpinorHelicity.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/SU2Helper.h"
	#include "MatchboxMEBase.h"

	#include <boost/numeric/ublas/io.hpp>
	#include <boost/numeric/ublas/matrix_proxy.hpp>
	#include <iterator>
	using std::ostream_iterator;

	using namespace Herwig;

	MatchboxAmplitude::MatchboxAmplitude()
	: Amplitude() {}

	MatchboxAmplitude::~MatchboxAmplitude() {}

	void MatchboxAmplitude::persistentOutput(PersistentOStream & os) const {
	os << theLastXComb << theColourBasis;
	}

	void MatchboxAmplitude::persistentInput(PersistentIStream & is, int) {
	is >> theLastXComb >> theColourBasis;
	lastMatchboxXComb(theLastXComb);
	}

	void MatchboxAmplitude::cloneDependencies(const std::string&) {}

	Ptr<MatchboxMEBase>::ptr MatchboxAmplitude::makeME(const vector<PDVector>&) const {
	return new_ptr(MatchboxMEBase());
	}

	struct orderPartonData {

	bool operator()(const pair<tcPDPtr,int>& a,
	const pair<tcPDPtr,int>& b) const {

	if ( a.first == b.first )
	return a.second < b.second;

	int acolour = a.first->iColour();
	int bcolour = b.first->iColour();

	if ( abs(acolour) != abs(bcolour) )
	return abs(acolour) < abs(bcolour);

	if ( a.first->iSpin() != b.first->iSpin() )
	return a.first->iSpin() < b.first->iSpin();

	int acharge = a.first->iCharge();
	int bcharge = b.first->iCharge();

	if ( abs(acharge) != abs(bcharge) )
	return abs(acharge) < abs(bcharge);

	if ( abs(a.first->id()) != abs(b.first->id()) )
	return abs(a.first->id()) < abs(b.first->id());

	return a.first->id() > b.first->id();

	}

	};

	void MatchboxAmplitude::setXComb(tStdXCombPtr xc) {
	theLastXComb = xc;
	lastMatchboxXComb(xc);
	fillCrossingMap();
	for ( size_t k = 0 ; k < meMomenta().size(); ++k )
	amplitudeMomenta()[k] = amplitudeMomentum(k);
	}

	void MatchboxAmplitude::fillCrossingMap(size_t shift) {

	if ( !amplitudePartonData().empty() )
	return;

	double csign = 1.;
	set<pair<tcPDPtr,int>,orderPartonData > processLegs;
	for ( unsigned int l = 0; l < mePartonData().size(); ++l ) {
	if ( l > 1 )
	processLegs.insert(make_pair(mePartonData()[l],l));
	else {
	if ( mePartonData()[l]->CC() ) {
	processLegs.insert(make_pair(mePartonData()[l]->CC(),l));
	if ( mePartonData()[l]->iSpin() == PDT::Spin1Half )
	csign *= -1.;
	} else {
	processLegs.insert(make_pair(mePartonData()[l],l));
	}
	}
	}

	crossingSign(csign);

	set<pair<tcPDPtr,int> > amplitudeLegs;
	crossingMap().resize(mePartonData().size());
	amplitudePartonData().resize(mePartonData().size());
	amplitudeMomenta().resize(mePartonData().size());

	int ampCount = 0;

	// process legs are already sorted, we only need to arrange for
	// adjacent particles and anti-particles
	while ( !processLegs.empty() ) {
	set<pair<tcPDPtr,int>,orderPartonData >::iterator next
	= processLegs.begin();
	while ( next->first->id() < 0 ) {
	if ( ++next == processLegs.end() )
	break;
	}
	assert(next != processLegs.end());
	crossingMap()[ampCount] = next->second - shift;
	amplitudeLegs.insert(make_pair(next->first,ampCount));
	tcPDPtr check = next->first;
	processLegs.erase(next);
	++ampCount;
	if ( check->CC() ) {
	set<pair<tcPDPtr,int>,orderPartonData>::iterator checkcc
	= processLegs.end();
	for ( set<pair<tcPDPtr,int>,orderPartonData>::iterator c = processLegs.begin();
	c != processLegs.end(); ++c ) {
	if ( c->first == check->CC() ) {
	checkcc = c; break;
	}
	}
	if ( checkcc == processLegs.end() )
	for ( set<pair<tcPDPtr,int>,orderPartonData>::iterator c = processLegs.begin();
	c != processLegs.end(); ++c ) {
	if ( !SU2Helper::SU2CC(check) )
	continue;
	if ( c->first == SU2Helper::SU2CC(check)->CC() ) {
	checkcc = c; break;
	}
	}
	if ( checkcc == processLegs.end() ) {
	int f = SU2Helper::family(check);
	for ( int i = 1 - f; i < 5 - f; i++ ) {
	bool gotone = false;
	for ( set<pair<tcPDPtr,int>,orderPartonData>::iterator c = processLegs.begin();
	c != processLegs.end(); ++c ) {
	if ( !SU2Helper::SU2CC(check,i) )
	continue;
	if ( c->first == SU2Helper::SU2CC(check,i)->CC() ) {
	checkcc = c; gotone = true; break;
	}
	}
	if ( gotone )
	break;
	}
	}
	// default to just pick the next available anti-particle
	if ( checkcc == processLegs.end() ) {
	checkcc = processLegs.begin();
	while ( checkcc->first->id() > 0 )
	if ( ++checkcc == processLegs.end() )
	break;
	}
	assert(checkcc != processLegs.end());
	crossingMap()[ampCount] = checkcc->second - shift;
	amplitudeLegs.insert(make_pair(checkcc->first,ampCount));
	processLegs.erase(checkcc);
	++ampCount;
	}
	}

	for ( set<pair<tcPDPtr,int> >::const_iterator l = amplitudeLegs.begin();
	l != amplitudeLegs.end(); ++l )
	amplitudePartonData()[l->second] = l->first;

	if ( colourBasis() ) {
	assert(colourBasis()->indexMap().find(mePartonData()) !=
	colourBasis()->indexMap().end());
	const map<size_t,size_t> colourCross =
	colourBasis()->indexMap().find(mePartonData())->second;
	for ( size_t k = 0; k < crossingMap().size(); ++k ) {
	if ( colourCross.find(crossingMap()[k]) !=
	colourCross.end() ) {
	size_t ccross = colourCross.find(crossingMap()[k])->second;
	amplitudeToColourMap()[k] = ccross;
	colourToAmplitudeMap()[ccross] = k;
	}
	}
	}

	}

	const string& MatchboxAmplitude::colourOrderingString(size_t id) const {

	static string empty = "";
	if ( !colourBasis() ) {
	return empty;
	}

	return colourBasis()->orderingString(mePartonData(),colourToAmplitudeMap(),id);

	}

	const vector<vector<size_t> >& MatchboxAmplitude::colourOrdering(size_t id) const {

	static vector<vector<size_t> > empty;
	if ( !colourBasis() ) {
	return empty;
	}

	return colourBasis()->ordering(mePartonData(),colourToAmplitudeMap(),id);

	}

	Lorentz5Momentum MatchboxAmplitude::amplitudeMomentum(int i) const {
	int iCrossed = crossingMap()[i];
	Lorentz5Momentum res = meMomenta()[iCrossed];
	if ( iCrossed < 2 )
	res = -res;
	res.setMass(meMomenta()[iCrossed].mass());
	res.rescaleRho();
	return res;
	}

	set<vector<int> > MatchboxAmplitude::generateHelicities() const {
	set<vector<int> > res;
	vector<int> current(amplitudePartonData().size());
	doGenerateHelicities(res,current,0);
	return res;
	}

	void MatchboxAmplitude::doGenerateHelicities(set<vector<int> >& res,
	vector<int>& current,
	size_t pos) const {

	if ( pos == amplitudePartonData().size() ) {
	res.insert(current);
	return;
	}

	if ( amplitudePartonData()[pos]->iSpin() == PDT::Spin0 \|\|
	( amplitudePartonData()[pos]->iSpin() == PDT::Spin1 &&
	amplitudePartonData()[pos]->mass() != ZERO ) ) {
	current[pos] = 0;
	doGenerateHelicities(res,current,pos+1);
	} else if ( amplitudePartonData()[pos]->iSpin() == PDT::Spin1Half \|\|
	amplitudePartonData()[pos]->iSpin() == PDT::Spin1 ) {
	current[pos] = 1;
	doGenerateHelicities(res,current,pos+1);
	current[pos] = -1;
	doGenerateHelicities(res,current,pos+1);
	}

	}

	void MatchboxAmplitude::prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr) {

	if ( !calculateTreeAmplitudes() )
	return;

	bool initialized = !lastAmplitudes().empty();

	if ( !initialized ) {
	set<vector<int> > helicities = generateHelicities();
	for ( set<vector<int> >::const_iterator h = helicities.begin();
	h != helicities.end(); ++h ) {
	lastAmplitudes().insert(make_pair(*h,CVector(colourBasisDim())));
	lastLargeNAmplitudes().insert(make_pair(*h,CVector(colourBasisDim())));
	}
	}

	AmplitudeIterator amp = lastAmplitudes().begin();
	AmplitudeIterator lamp = lastLargeNAmplitudes().begin();
	for ( ;amp != lastAmplitudes().end(); ++amp, ++lamp ) {
	for ( size_t k = 0; k < colourBasisDim(); ++k )
	amp->second(k) = evaluate(k,amp->first,lamp->second(k));
	}

	if ( !initialized ) {
	map<vector<int>,CVector> clean;
	for ( map<vector<int>,CVector>::const_iterator amp = lastAmplitudes().begin();
	amp != lastAmplitudes().end(); ++amp ) {
	bool nonZero = false;
	for ( size_t k = 0; k < colourBasisDim(); ++k ) {
	if ( amp->second(k) != Complex(0.0) ) {
	nonZero = true;
	break;
	}
	}
	if ( nonZero )
	clean.insert(*amp);
	}
	lastAmplitudes() = clean;
	clean.clear();
	for ( map<vector<int>,CVector>::const_iterator amp = lastLargeNAmplitudes().begin();
	amp != lastLargeNAmplitudes().end(); ++amp ) {
	bool nonZero = false;
	for ( size_t k = 0; k < colourBasisDim(); ++k ) {
	if ( amp->second(k) != Complex(0.0) ) {
	nonZero = true;
	break;
	}
	}
	if ( nonZero )
	clean.insert(*amp);
	}
	lastLargeNAmplitudes() = clean;
	}

	haveTreeAmplitudes();

	}

	void MatchboxAmplitude::prepareOneLoopAmplitudes(Ptr<MatchboxMEBase>::tcptr) {

	if ( !calculateOneLoopAmplitudes() )
	return;

	bool initialized = !lastOneLoopAmplitudes().empty();

	if ( !initialized ) {
	set<vector<int> > helicities = generateHelicities();
	for ( set<vector<int> >::const_iterator h = helicities.begin();
	h != helicities.end(); ++h ) {
	lastOneLoopAmplitudes().insert(make_pair(*h,CVector(colourBasisDim())));
	}
	}

	for ( AmplitudeIterator amp = lastOneLoopAmplitudes().begin();
	amp != lastOneLoopAmplitudes().end(); ++amp ) {
	for ( size_t k = 0; k < colourBasisDim(); ++k )
	amp->second(k) = evaluateOneLoop(k,amp->first);
	}

	if ( !initialized ) {
	map<vector<int>,CVector> clean;
	for ( map<vector<int>,CVector>::const_iterator amp = lastOneLoopAmplitudes().begin();
	amp != lastOneLoopAmplitudes().end(); ++amp ) {
	bool nonZero = false;
	for ( size_t k = 0; k < colourBasisDim(); ++k ) {
	if ( amp->second(k) != Complex(0.0) ) {
	nonZero = true;
	break;
	}
	}
	if ( nonZero )
	clean.insert(*amp);
	}
	lastOneLoopAmplitudes() = clean;
	}

	haveOneLoopAmplitudes();

	}

	Complex MatchboxAmplitude::value(const tcPDVector&,
	const vector<Lorentz5Momentum>&,
	const vector<int>&) {
	assert(false && "ThePEG::Amplitude interface is not sufficient at the moment.");
	throw Exception() << "ThePEG::Amplitude interface is not sufficient at the moment."
	<< Exception::abortnow;
	return 0.;
	}

	double MatchboxAmplitude::me2() const {
	if ( !calculateTreeME2() )
	return lastTreeME2();
	lastTreeME2(crossingSign()*colourBasis()->me2(mePartonData(),lastAmplitudes()));
	return lastTreeME2();
	}

	double MatchboxAmplitude::oneLoopInterference() const {
	if ( !calculateOneLoopInterference() )
	return lastOneLoopInterference();
	lastOneLoopInterference(crossingSign()*
	colourBasis()->interference(mePartonData(),
	lastOneLoopAmplitudes(),lastAmplitudes()));
	return lastOneLoopInterference();
	}

	double MatchboxAmplitude::colourCorrelatedME2(pair<int,int> ij) const {
	if ( !calculateColourCorrelator(ij) )
	return lastColourCorrelator(ij);
	double Nc = generator()->standardModel()->Nc();
	double cfac = 1.;
	if ( mePartonData()[ij.first]->iColour() == PDT::Colour8 ) {
	cfac = Nc;
	} else if ( mePartonData()[ij.first]->iColour() == PDT::Colour3 \|\|
	mePartonData()[ij.first]->iColour() == PDT::Colour3bar ) {
	cfac = (sqr(Nc)-1.)/(2.*Nc);
	} else assert(false);
	double res =
	crossingSign()*colourBasis()->colourCorrelatedME2(ij,mePartonData(),lastAmplitudes());
	lastColourCorrelator(ij,res);
	return res/cfac;
	}

	double MatchboxAmplitude::largeNColourCorrelatedME2(pair<int,int> ij,
	Ptr<ColourBasis>::tptr largeNBasis) const {
	if ( !calculateLargeNColourCorrelator(ij) )
	return lastLargeNColourCorrelator(ij);
	double Nc = generator()->standardModel()->Nc();
	double cfac = 1.;
	if ( mePartonData()[ij.first]->iColour() == PDT::Colour8 ) {
	cfac = Nc;
	} else if ( mePartonData()[ij.first]->iColour() == PDT::Colour3 \|\|
	mePartonData()[ij.first]->iColour() == PDT::Colour3bar ) {
	cfac = Nc/2.;
	} else assert(false);
	double res =
	- crossingSign()*largeNBasis->colourCorrelatedME2(ij,mePartonData(),lastAmplitudes());
	+ crossingSign()*largeNBasis->colourCorrelatedME2(ij,mePartonData(),lastLargeNAmplitudes());
	lastLargeNColourCorrelator(ij,res);
	return res/cfac;
	}

	// compare int vectors modulo certain element
	// which needs to differe between the two
	bool equalsModulo(unsigned int i, const vector<int>& a, const vector<int>& b) {
	assert(a.size()==b.size());
	if ( a[i] == b[i] )
	return false;
	for ( unsigned int k = 0; k < a.size(); ++k ) {
	if ( k == i )
	continue;
	if ( a[k] != b[k] )
	return false;
	}
	return true;
	}

	LorentzVector<Complex> MatchboxAmplitude::plusPolarization(const Lorentz5Momentum& p,
	const Lorentz5Momentum& n) const {

	using namespace SpinorHelicity;

	LorentzVector<complex<Energy> > num =
	PlusSpinorCurrent(PlusConjugateSpinor(n),MinusSpinor(p)).eval();

	complex<Energy> den =
	sqrt(2.)*PlusSpinorProduct(PlusConjugateSpinor(n),PlusSpinor(p)).eval();

	LorentzVector<Complex> polarization(num.x()/den,num.y()/den,num.z()/den,num.t()/den);

	return polarization;

	}

	double MatchboxAmplitude::spinColourCorrelatedME2(pair<int,int> ij,
	const SpinCorrelationTensor& c) const {

	Lorentz5Momentum p = meMomenta()[ij.first];
	Lorentz5Momentum n = meMomenta()[ij.second];

	LorentzVector<Complex> polarization = plusPolarization(p,n);

	Complex pFactor = (polarization*c.momentum())/sqrt(abs(c.scale()));

	double avg =
	colourCorrelatedME2(ij)(-c.diagonal()+ (c.scale() > ZERO ? 1. : -1.)norm(pFactor));

	int iCrossed = -1;
	for ( unsigned int k = 0; k < crossingMap().size(); ++k )
	if ( crossingMap()[k] == ij.first ) {
	iCrossed = k;
	break;
	}
	assert(iCrossed >= 0);

	Complex csCorr = 0.0;

	if ( calculateColourSpinCorrelator(ij) ) {
	set<const CVector*> done;
	for ( AmplitudeConstIterator a = lastAmplitudes().begin();
	a != lastAmplitudes().end(); ++a ) {
	if ( done.find(&(a->second)) != done.end() )
	continue;
	AmplitudeConstIterator b = lastAmplitudes().begin();
	while ( !equalsModulo(iCrossed,a->first,b->first) )
	if ( ++b == lastAmplitudes().end() )
	break;
	if ( b == lastAmplitudes().end() \|\| done.find(&(b->second)) != done.end() )
	continue;
	done.insert(&(a->second)); done.insert(&(b->second));
	if ( a->first[iCrossed] == 1 )
	swap(a,b);
	csCorr += colourBasis()->colourCorrelatedInterference(ij,mePartonData(),a->second,b->second);
	}
	lastColourSpinCorrelator(ij,csCorr);
	} else {
	csCorr = lastColourSpinCorrelator(ij);
	}

	double corr =
	2.real(csCorrsqr(pFactor));

	double Nc = generator()->standardModel()->Nc();
	double cfac = 1.;
	if ( mePartonData()[ij.first]->iColour() == PDT::Colour8 ) {
	cfac = Nc;
	} else if ( mePartonData()[ij.first]->iColour() == PDT::Colour3 \|\|
	mePartonData()[ij.first]->iColour() == PDT::Colour3bar ) {
	cfac = (sqr(Nc)-1.)/(2.*Nc);
	} else assert(false);

	return
	avg + crossingSign()(c.scale() > ZERO ? 1. : -1.)corr/cfac;

	}

	void MatchboxAmplitude::Init() {

	static ClassDocumentation<MatchboxAmplitude> documentation
	("MatchboxAmplitude is the base class for amplitude "
	"implementations inside Matchbox.");

	static Reference<MatchboxAmplitude,ColourBasis> interfaceColourBasis
	("ColourBasis",
	"Set the colour basis implementation.",
	&MatchboxAmplitude::theColourBasis, false, false, true, true, false);

	}

	// * Attention * The following static variable is needed for the type
	// description system in ThePEG. Please check that the template arguments
	// are correct (the class and its base class), and that the constructor
	// arguments are correct (the class name and the name of the dynamically
	// loadable library where the class implementation can be found).
	DescribeAbstractClass<MatchboxAmplitude,Amplitude>
	describeMatchboxAmplitude("Herwig::MatchboxAmplitude", "HwMatchbox.so");
	diff --git a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
	--- a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
	+++ b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
	@@ -1,1200 +1,1201 @@
	// -- C++ --
	//
	// MatchboxMEBase.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 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 MatchboxMEBase class.
	//

	#include "MatchboxMEBase.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDF/PDF.h"
	#include "ThePEG/PDT/PDT.h"
	#include "ThePEG/StandardModel/StandardModelBase.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/Handlers/StdXCombGroup.h"
	#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/DiagramDrawer.h"
	#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"

	#include <iterator>
	using std::ostream_iterator;

	using namespace Herwig;

	MatchboxMEBase::MatchboxMEBase()
	: MEBase(),
	theOneLoop(false),
	theOneLoopNoBorn(false) {}

	MatchboxMEBase::~MatchboxMEBase() {}

	Ptr<MatchboxFactory>::tcptr MatchboxMEBase::factory() const { return theFactory; }

	void MatchboxMEBase::factory(Ptr<MatchboxFactory>::tcptr f) { theFactory = f; }

	Ptr<Tree2toNGenerator>::tptr MatchboxMEBase::diagramGenerator() const { return factory()->diagramGenerator(); }

	Ptr<ProcessData>::tptr MatchboxMEBase::processData() const { return factory()->processData(); }

	unsigned int MatchboxMEBase::getNLight() const { return factory()->nLight(); }

	double MatchboxMEBase::factorizationScaleFactor() const { return factory()->factorizationScaleFactor(); }

	double MatchboxMEBase::renormalizationScaleFactor() const { return factory()->renormalizationScaleFactor(); }

	bool MatchboxMEBase::fixedCouplings() const { return factory()->fixedCouplings(); }

	bool MatchboxMEBase::fixedQEDCouplings() const { return factory()->fixedQEDCouplings(); }

	bool MatchboxMEBase::checkPoles() const { return factory()->checkPoles(); }

	bool MatchboxMEBase::verbose() const { return factory()->verbose(); }

	bool MatchboxMEBase::initVerbose() const { return factory()->initVerbose(); }

	void MatchboxMEBase::getDiagrams() const {

	if ( diagramGenerator() && processData() ) {

	vector<Ptr<Tree2toNDiagram>::ptr> diags;
	for ( vector<PDVector>::const_iterator p = subProcesses().begin();
	p != subProcesses().end(); ++p ) {
	vector<Ptr<Tree2toNDiagram>::ptr>& res =
	processData()->diagramMap()[*p];
	if ( res.empty() ) {
	res = diagramGenerator()->generate(*p,orderInAlphaS(),orderInAlphaEW());
	}
	copy(res.begin(),res.end(),back_inserter(diags));
	processData()->fillMassGenerators(*p);
	}

	if ( diags.empty() )
	return;

	for ( vector<Ptr<Tree2toNDiagram>::ptr>::iterator d = diags.begin();
	d != diags.end(); ++d ) {
	add(*d);
	}

	return;

	}

	throw Exception()
	<< "MatchboxMEBase::getDiagrams() expects a Tree2toNGenerator and ProcessData object.\n"
	<< "Please check your setup." << Exception::abortnow;

	}

	Selector<MEBase::DiagramIndex>
	MatchboxMEBase::diagrams(const DiagramVector & diags) const {

	if ( phasespace() ) {
	return phasespace()->selectDiagrams(diags);
	}

	throw Exception()
	<< "MatchboxMEBase::diagrams() expects a MatchboxPhasespace object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return Selector<MEBase::DiagramIndex>();

	}

	Selector<const ColourLines *>
	MatchboxMEBase::colourGeometries(tcDiagPtr diag) const {

	if ( matchboxAmplitude() ) {
	if ( !matchboxAmplitude()->haveColourFlows() )
	throw Exception() << "A colour flow implementation is not present."
	<< Exception::abortnow;
	if ( matchboxAmplitude()->treeAmplitudes() )
	matchboxAmplitude()->prepareAmplitudes(this);
	return matchboxAmplitude()->colourGeometries(diag);
	}

	throw Exception()
	<< "MatchboxMEBase::colourGeometries() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return Selector<const ColourLines *>();

	}

	unsigned int MatchboxMEBase::orderInAlphaS() const {
	if ( matchboxAmplitude() ) {
	return matchboxAmplitude()->orderInGs();
	}
	throw Exception()
	<< "MatchboxMEBase::orderInAlphaS() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0;
	}

	unsigned int MatchboxMEBase::orderInAlphaEW() const {
	if ( matchboxAmplitude() ) {
	return matchboxAmplitude()->orderInGem();
	}
	throw Exception()
	<< "MatchboxMEBase::orderInAlphaEW() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0;
	}

	void MatchboxMEBase::setXComb(tStdXCombPtr xc) {

	MEBase::setXComb(xc);
	lastMatchboxXComb(xc);
	if ( phasespace() )
	phasespace()->setXComb(xc);
	if ( scaleChoice() )
	scaleChoice()->setXComb(xc);
	if ( matchboxAmplitude() )
	matchboxAmplitude()->setXComb(xc);

	}

	double MatchboxMEBase::generateIncomingPartons(const double* r1, const double* r2) {

	// shamelessly stolen from PartonExtractor.cc

	Energy2 shmax = lastCuts().sHatMax();
	Energy2 shmin = lastCuts().sHatMin();
	Energy2 sh = shminpow(shmax/shmin, r1);
	double ymax = lastCuts().yHatMax();
	double ymin = lastCuts().yHatMin();
	double km = log(shmax/shmin);
	ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/sh)));
	ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/sh)));

	double y = ymin + (r2)(ymax - ymin);
	double x1 = exp(-0.5*log(lastS()/sh) + y);
	double x2 = exp(-0.5*log(lastS()/sh) - y);

	Lorentz5Momentum P1 = lastParticles().first->momentum();
	LorentzMomentum p1 = lightCone((P1.rho() + P1.e())*x1, Energy());
	p1.rotateY(P1.theta());
	p1.rotateZ(P1.phi());
	meMomenta()[0] = p1;

	Lorentz5Momentum P2 = lastParticles().second->momentum();
	LorentzMomentum p2 = lightCone((P2.rho() + P2.e())*x2, Energy());
	p2.rotateY(P2.theta());
	p2.rotateZ(P2.phi());
	meMomenta()[1] = p2;

	lastXCombPtr()->lastX1X2(make_pair(x1,x2));
	lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());

	return km*(ymax - ymin);

	}

	bool MatchboxMEBase::generateKinematics(const double * r) {

	if ( phasespace() ) {

	jacobian(phasespace()->generateKinematics(r,meMomenta()));
	if ( jacobian() == 0.0 )
	return false;

	setScale();
	logGenerateKinematics(r);

	assert(lastMatchboxXComb());

	if ( nDimAmplitude() > 0 ) {
	amplitudeRandomNumbers().resize(nDimAmplitude());
	copy(r + nDimPhasespace(),
	r + nDimPhasespace() + nDimAmplitude(),
	amplitudeRandomNumbers().begin());
	}

	if ( nDimInsertions() > 0 ) {
	insertionRandomNumbers().resize(nDimInsertions());
	copy(r + nDimPhasespace() + nDimAmplitude(),
	r + nDimPhasespace() + nDimAmplitude() + nDimInsertions(),
	insertionRandomNumbers().begin());
	}

	return true;

	}

	throw Exception()
	<< "MatchboxMEBase::generateKinematics() expects a MatchboxPhasespace object.\n"
	<< "Please check your setup." << Exception::abortnow;

	return false;

	}

	int MatchboxMEBase::nDim() const {

	if ( lastMatchboxXComb() )
	return nDimPhasespace() + nDimAmplitude() + nDimInsertions();

	int ampAdd = 0;
	if ( matchboxAmplitude() ) {
	ampAdd = matchboxAmplitude()->nDimAdditional();
	}

	int insertionAdd = 0;
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	insertionAdd = max(insertionAdd,(**v).nDimAdditional());
	}

	return nDimBorn() + ampAdd + insertionAdd;

	}

	int MatchboxMEBase::nDimBorn() const {

	if ( lastMatchboxXComb() )
	return nDimPhasespace();

	if ( phasespace() ) {
	size_t nout = diagrams().front()->partons().size()-2;
	int n = phasespace()->nDim(nout);
	if ( phasespace()->useMassGenerators() ) {
	for ( cPDVector::const_iterator pd =
	diagrams().front()->partons().begin();
	pd != diagrams().front()->partons().end(); ++pd ) {
	if ( processData()->massGenerator(*pd) \|\|
	(**pd).width() != ZERO ) {
	++n;
	}
	}
	}
	return n;
	}

	throw Exception()
	<< "MatchboxMEBase::nDim() expects a MatchboxPhasespace object.\n"
	<< "Please check your setup." << Exception::abortnow;

	return 0;

	}

	void MatchboxMEBase::setScale() const {
	if ( haveX1X2() ) {
	lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
	}
	Energy2 fscale = factorizationScale()*sqr(factorizationScaleFactor());
	Energy2 rscale = renormalizationScale()*sqr(renormalizationScaleFactor());
	Energy2 ewrscale = renormalizationScaleQED();
	lastXCombPtr()->lastScale(fscale);
	if ( !fixedCouplings() ) {
	if ( rscale > lastCuts().scaleMin() )
	lastXCombPtr()->lastAlphaS(SM().alphaS(rscale));
	else
	lastXCombPtr()->lastAlphaS(SM().alphaS(lastCuts().scaleMin()));
	} else {
	lastXCombPtr()->lastAlphaS(SM().alphaS());
	}
	if ( !fixedQEDCouplings() ) {
	lastXCombPtr()->lastAlphaEM(SM().alphaEM(ewrscale));
	} else {
	lastXCombPtr()->lastAlphaEM(SM().alphaEMMZ());
	}
	logSetScale();
	}

	Energy2 MatchboxMEBase::factorizationScale() const {
	if ( scaleChoice() ) {
	return scaleChoice()->factorizationScale();
	}

	throw Exception()
	<< "MatchboxMEBase::factorizationScale() expects a MatchboxScaleChoice object.\n"
	<< "Please check your setup." << Exception::abortnow;

	return ZERO;

	}

	Energy2 MatchboxMEBase::renormalizationScale() const {
	if ( scaleChoice() ) {
	return scaleChoice()->renormalizationScale();
	}

	throw Exception()
	<< "MatchboxMEBase::renormalizationScale() expects a MatchboxScaleChoice object.\n"
	<< "Please check your setup." << Exception::abortnow;

	return ZERO;

	}

	Energy2 MatchboxMEBase::renormalizationScaleQED() const {
	if ( scaleChoice() ) {
	return scaleChoice()->renormalizationScaleQED();
	}
	return renormalizationScale();
	}

	void MatchboxMEBase::setVetoScales(tSubProPtr) const {}

	void MatchboxMEBase::getPDFWeight(Energy2 factorizationScale) const {

	if ( !mePartonData()[0]->coloured() &&
	!mePartonData()[1]->coloured() ) {
	lastMEPDFWeight(1.0);
	logPDFWeight();
	return;
	}

	double w = 1.;

	if ( mePartonData()[0]->coloured() && havePDFWeight1() )
	w *= pdf1(factorizationScale);

	if ( mePartonData()[1]->coloured() && havePDFWeight2() )
	w *= pdf2(factorizationScale);

	lastMEPDFWeight(w);

	logPDFWeight();

	}

	double MatchboxMEBase::pdf1(Energy2 fscale, double xEx) const {

	assert(lastXCombPtr()->partonBins().first->pdf());

	if ( xEx < 1. && lastX1() >= xEx ) {
	return
	( ( 1. - lastX1() ) / ( 1. - xEx ) ) *
	lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
	lastPartons().first->dataPtr(),
	fscale == ZERO ? lastScale() : fscale,
	xEx)/xEx;
	}

	return lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
	lastPartons().first->dataPtr(),
	fscale == ZERO ? lastScale() : fscale,
	lastX1())/lastX1();
	}

	double MatchboxMEBase::pdf2(Energy2 fscale, double xEx) const {

	assert(lastXCombPtr()->partonBins().second->pdf());

	if ( xEx < 1. && lastX2() >= xEx ) {
	return
	( ( 1. - lastX2() ) / ( 1. - xEx ) ) *
	lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
	lastPartons().second->dataPtr(),
	fscale == ZERO ? lastScale() : fscale,
	xEx)/xEx;
	}

	return lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
	lastPartons().second->dataPtr(),
	fscale == ZERO ? lastScale() : fscale,
	lastX2())/lastX2();

	}

	double MatchboxMEBase::me2() const {

	if ( matchboxAmplitude() ) {

	if ( matchboxAmplitude()->treeAmplitudes() )
	matchboxAmplitude()->prepareAmplitudes(this);

	lastME2(matchboxAmplitude()->me2()*
	crossingSign()*
	me2Norm());

	logME2();

	return lastME2();

	}

	throw Exception()
	<< "MatchboxMEBase::me2() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0.;

	}

	double MatchboxMEBase::finalStateSymmetry() const {

	if ( symmetryFactor() > 0.0 )
	return symmetryFactor();

	double sFactor = 1.;

	map<long,int> counts;

	cPDVector checkData;
	copy(mePartonData().begin()+2,mePartonData().end(),back_inserter(checkData));

	cPDVector::iterator p = checkData.begin();
	while ( !checkData.empty() ) {
	if ( counts.find((**p).id()) != counts.end() ) {
	counts[(**p).id()] += 1;
	} else {
	counts[(**p).id()] = 1;
	}
	checkData.erase(p);
	p = checkData.begin();
	continue;
	}

	for ( map<long,int>::const_iterator c = counts.begin();
	c != counts.end(); ++c ) {
	if ( c->second == 1 )
	continue;
	if ( c->second == 2 )
	sFactor /= 2.;
	else if ( c->second == 3 )
	sFactor /= 6.;
	else if ( c->second == 4 )
	sFactor /= 24.;
	}

	symmetryFactor(sFactor);

	return symmetryFactor();

	}

	double MatchboxMEBase::me2Norm(unsigned int addAlphaS) const {

	// assume that we always have incoming
	// spin-1/2 or massless spin-1 particles
	double fac = 1./4.;

	if ( orderInAlphaS() > 0 \|\| addAlphaS != 0 )
	fac *= pow(lastAlphaS()/SM().alphaS(),double(orderInAlphaS()+addAlphaS));
	if ( orderInAlphaEW() > 0 )
	fac *= pow(lastAlphaEM()/SM().alphaEM(),double(orderInAlphaEW()));

	if ( mePartonData()[0]->iColour() == PDT::Colour3 \|\|
	mePartonData()[0]->iColour() == PDT::Colour3bar )
	fac /= SM().Nc();
	else if ( mePartonData()[0]->iColour() == PDT::Colour8 )
	fac /= (SM().Nc()*SM().Nc()-1.);

	if ( mePartonData()[1]->iColour() == PDT::Colour3 \|\|
	mePartonData()[1]->iColour() == PDT::Colour3bar )
	fac /= SM().Nc();
	else if ( mePartonData()[1]->iColour() == PDT::Colour8 )
	fac /= (SM().Nc()*SM().Nc()-1.);

	return finalStateSymmetry()*fac;

	}

	CrossSection MatchboxMEBase::dSigHatDR() const {

	getPDFWeight();

	if ( !lastXCombPtr()->willPassCuts() ) {
	lastME2(0.0);
	lastMECrossSection(ZERO);
	return lastMECrossSection();
	}

	double xme2 = me2();
	lastME2(xme2);

	if ( xme2 == 0. && !oneLoopNoBorn() ) {
	lastME2(0.0);
	lastMECrossSection(ZERO);
	return lastMECrossSection();
	}

	double vme2 = 0.;
	if ( oneLoop() )
	vme2 = oneLoopInterference();

	CrossSection res = ZERO;

	if ( !oneLoopNoBorn() )
	res +=
	(sqr(hbarc)/(2.lastSHat()))
	jacobian()* lastMEPDFWeight() * xme2;

	if ( oneLoop() )
	res +=
	(sqr(hbarc)/(2.lastSHat()))
	jacobian()* lastMEPDFWeight() * vme2;

	if ( !onlyOneLoop() ) {
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	(**v).setXComb(lastXCombPtr());
	res += (**v).dSigHatDR();
	}
	}

	double weight = 0.0;
	bool applied = false;
	for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
	theReweights.begin(); rw != theReweights.end(); ++rw ) {
	(**rw).setXComb(lastXCombPtr());
	if ( !(**rw).apply() )
	continue;
	weight += (**rw).evaluate();
	applied = true;
	}
	if ( applied )
	res *= weight;
	lastMECrossSection(res);

	return lastMECrossSection();

	}

	double MatchboxMEBase::oneLoopInterference() const {

	if ( matchboxAmplitude() ) {

	if ( matchboxAmplitude()->oneLoopAmplitudes() )
	matchboxAmplitude()->prepareOneLoopAmplitudes(this);
	lastME2(matchboxAmplitude()->oneLoopInterference()*
	crossingSign()*
	me2Norm(1));

	logME2();

	return lastME2();

	}

	throw Exception()
	<< "MatchboxMEBase::oneLoopInterference() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0.;

	}

	void MatchboxMEBase::logPoles() const {
	double res2me = oneLoopDoublePole();
	double res1me = oneLoopSinglePole();
	double res2i = 0.;
	double res1i = 0.;
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	res2i += (**v).oneLoopDoublePole();
	res1i += (**v).oneLoopSinglePole();
	}
	double diff2 = abs(res2me) != 0. ? 1.-abs(res2i/res2me) : abs(res2i)-abs(res2me);
	double diff1 = abs(res1me) != 0. ? 1.-abs(res1i/res1me) : abs(res1i)-abs(res1me);
	generator()->log()
	<< "check "
	<< log10(abs(diff2)) << " " << log10(abs(diff1)) << "\n"
	<< flush;
	}

	bool MatchboxMEBase::haveOneLoop() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->haveOneLoop();
	return false;
	}

	bool MatchboxMEBase::onlyOneLoop() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->onlyOneLoop();
	return false;
	}

	bool MatchboxMEBase::isDR() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->isDR();
	return false;
	}

	bool MatchboxMEBase::isCS() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->isCS();
	return false;
	}

	bool MatchboxMEBase::isBDK() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->isBDK();
	return false;
	}

	bool MatchboxMEBase::isExpanded() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->isExpanded();
	return false;
	}

	Energy2 MatchboxMEBase::mu2() const {
	if ( matchboxAmplitude() )
	return matchboxAmplitude()->mu2();
	return 0*GeV2;
	}

	double MatchboxMEBase::oneLoopDoublePole() const {

	if ( matchboxAmplitude() ) {

	return
	matchboxAmplitude()->oneLoopDoublePole()*
	crossingSign()*
	me2Norm(1);

	}

	return 0.;

	}

	double MatchboxMEBase::oneLoopSinglePole() const {

	if ( matchboxAmplitude() ) {

	return
	matchboxAmplitude()->oneLoopSinglePole()*
	crossingSign()*
	me2Norm(1);

	}

	return 0.;

	}

	vector<Ptr<SubtractionDipole>::ptr>
	MatchboxMEBase::getDipoles(const vector<Ptr<SubtractionDipole>::ptr>& dipoles,
	const vector<Ptr<MatchboxMEBase>::ptr>& borns) const {

	vector<Ptr<SubtractionDipole>::ptr> res;

	// keep track of the dipoles we already did set up
	set<pair<pair<pair<int,int>,int>,pair<Ptr<MatchboxMEBase>::tptr,Ptr<SubtractionDipole>::tptr> > > done;

	cPDVector rep = diagrams().front()->partons();
	int nreal = rep.size();

	// now loop over configs
	for ( int emitter = 0; emitter < nreal; ++emitter ) {
	for ( int spectator = 0; spectator < nreal; ++spectator ) {
	if ( emitter == spectator )
	continue;
	for ( int emission = 2; emission < nreal; ++emission ) {
	if ( emission == emitter \|\| emission == spectator )
	continue;
	for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator b =
	borns.begin(); b != borns.end(); ++b ) {
	if ( (**b).onlyOneLoop() )
	continue;
	for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d =
	dipoles.begin(); d != dipoles.end(); ++d ) {
	if ( !rep[emitter]->coloured() \|\|
	!rep[emission]->coloured() \|\|
	!rep[spectator]->coloured() ) {
	continue;
	}
	if ( noDipole(emitter,emission,spectator) ) {
	continue;
	}
	if ( done.find(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(b,d)))
	!= done.end() ) {
	continue;
	}
	if ( !(**d).canHandle(rep,emitter,emission,spectator) ) {
	continue;
	}
	// now get to work
	Ptr<SubtractionDipole>::ptr nDipole = (**d).cloneMe();
	nDipole->realEmitter(emitter);
	nDipole->realEmission(emission);
	nDipole->realSpectator(spectator);
	nDipole->realEmissionME(const_cast<MatchboxMEBase*>(this));
	nDipole->underlyingBornME(*b);
	nDipole->setupBookkeeping();
	if ( !(nDipole->empty()) ) {
	res.push_back(nDipole);
	done.insert(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(b,d)));
	if ( nDipole->isSymmetric() )
	done.insert(make_pair(make_pair(make_pair(emission,emitter),spectator),make_pair(b,d)));
	ostringstream dname;
	dname << fullName() << "." << (**b).name() << "."
	<< (**d).name() << ".[("
	<< emitter << "," << emission << ")," << spectator << "]";
	if ( ! (generator()->preinitRegister(nDipole,dname.str()) ) )
	throw InitException() << "Dipole " << dname.str() << " already existing.";
	nDipole->cloneDependencies(dname.str());
	}
	}
	}
	}
	}
	}

	for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d = res.begin();
	d != res.end(); ++d )
	(**d).partnerDipoles(res);

	return res;

	}

	double MatchboxMEBase::colourCorrelatedME2(pair<int,int> ij) const {

	if ( matchboxAmplitude() ) {

	if ( matchboxAmplitude()->treeAmplitudes() )
	matchboxAmplitude()->prepareAmplitudes(this);
	lastME2(matchboxAmplitude()->colourCorrelatedME2(ij)*
	crossingSign()*
	me2Norm());

	logME2();

	return lastME2();

	}

	throw Exception()
	<< "MatchboxMEBase::colourCorrelatedME2() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0.;

	}

	double MatchboxMEBase::largeNColourCorrelatedME2(pair<int,int> ij,
	Ptr<ColourBasis>::tptr largeNBasis) const {

	if ( matchboxAmplitude() ) {

	if ( matchboxAmplitude()->treeAmplitudes() )
	matchboxAmplitude()->prepareAmplitudes(this);
	+ largeNBasis->prepare(mePartonData(),false);
	lastME2(matchboxAmplitude()->largeNColourCorrelatedME2(ij,largeNBasis)*
	crossingSign()*
	me2Norm());

	logME2();

	return lastME2();

	}

	throw Exception()
	<< "MatchboxMEBase::largeNColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0.;

	}

	double MatchboxMEBase::spinColourCorrelatedME2(pair<int,int> ij,
	const SpinCorrelationTensor& c) const {

	if ( matchboxAmplitude() ) {

	if ( matchboxAmplitude()->treeAmplitudes() )
	matchboxAmplitude()->prepareAmplitudes(this);
	lastME2(matchboxAmplitude()->spinColourCorrelatedME2(ij,c)*
	crossingSign()*
	me2Norm());

	logME2();

	return lastME2();

	}

	throw Exception()
	<< "MatchboxMEBase::spinColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
	<< "Please check your setup." << Exception::abortnow;
	return 0.;

	}

	void MatchboxMEBase::flushCaches() {
	MEBase::flushCaches();
	if ( matchboxAmplitude() )
	matchboxAmplitude()->flushCaches();
	for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r =
	reweights().begin(); r != reweights().end(); ++r ) {
	(**r).flushCaches();
	}
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	(**v).flushCaches();
	}
	}

	void MatchboxMEBase::print(ostream& os) const {

	os << "--- MatchboxMEBase setup -------------------------------------------------------\n";

	os << " '" << name() << "' for subprocesses:\n";
	for ( vector<PDVector>::const_iterator p = subProcesses().begin();
	p != subProcesses().end(); ++p ) {
	os << " ";
	for ( PDVector::const_iterator pp = p->begin();
	pp != p->end(); ++pp ) {
	os << (**pp).PDGName() << " ";
	if ( pp == p->begin() + 1 )
	os << "-> ";
	}
	os << "\n";
	}

	os << " including " << (oneLoop() ? "" : "no ") << "virtual corrections";
	if ( oneLoopNoBorn() )
	os << " without Born contributions";
	os << "\n";

	if ( oneLoop() && !onlyOneLoop() ) {
	os << " using insertion operators\n";
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	os << " '" << (**v).name() << "' with "
	<< ((**v).isDR() ? "" : "C") << "DR/";
	if ( (**v).isCS() )
	os << "CS";
	if ( (**v).isBDK() )
	os << "BDK";
	if ( (**v).isExpanded() )
	os << "expanded";
	os << " conventions\n";
	}
	}

	os << "--------------------------------------------------------------------------------\n";

	os << flush;

	}

	void MatchboxMEBase::printLastEvent(ostream& os) const {

	os << "--- MatchboxMEBase last event information --------------------------------------\n";

	os << " for matrix element '" << name() << "'\n";

	os << " process considered:\n ";

	int in = 0;
	for ( cPDVector::const_iterator p = mePartonData().begin();
	p != mePartonData().end(); ++p ) {
	os << (**p).PDGName() << " ";
	if ( ++in == 2 )
	os << " -> ";
	}

	os << " kinematic environment as set by the XComb " << lastXCombPtr() << ":\n"
	<< " sqrt(shat)/GeV = " << sqrt(lastSHat()/GeV2)
	<< " x1 = " << lastX1() << " x2 = " << lastX2()
	<< " alphaS = " << lastAlphaS() << "\n";

	os << " momenta/GeV generated from random numbers\n ";
	copy(lastXComb().lastRandomNumbers().begin(),
	lastXComb().lastRandomNumbers().end(),ostream_iterator<double>(os," "));
	os << ":\n ";

	for ( vector<Lorentz5Momentum>::const_iterator p = meMomenta().begin();
	p != meMomenta().end(); ++p ) {
	os << (*p/GeV) << "\n ";
	}

	os << "last cross section/nb calculated was:\n "
	<< (lastMECrossSection()/nanobarn) << " (pdf weight " << lastMEPDFWeight() << ")\n";

	os << "--------------------------------------------------------------------------------\n";

	os << flush;

	}

	void MatchboxMEBase::logGenerateKinematics(const double * r) const {

	if ( !verbose() )
	return;

	generator()->log() << "'" << name() << "' generated kinematics\nfrom "
	<< nDim() << " random numbers:\n";
	copy(r,r+nDim(),ostream_iterator<double>(generator()->log()," "));
	generator()->log() << "\n";

	generator()->log() << "storing phase space information in XComb "
	<< lastXCombPtr() << "\n";

	generator()->log() << "generated phase space point (in GeV):\n";

	vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
	cPDVector::const_iterator dit = mePartonData().begin();

	for ( ; pit != meMomenta().end() ; ++pit, ++dit )
	generator()->log() << (**dit).PDGName() << " : "
	<< (*pit/GeV) << "\n";

	generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
	<< "and Jacobian = " << jacobian() << " sHat/GeV2 = "
	<< (lastSHat()/GeV2) << "\n" << flush;

	}

	void MatchboxMEBase::logSetScale() const {

	if ( !verbose() )
	return;

	generator()->log() << "'" << name() << "' set scales using XComb " << lastXCombPtr() << ":\n"
	<< "scale/GeV2 = " << (scale()/GeV2) << " xi_R = "
	<< renormalizationScaleFactor() << " xi_F = "
	<< factorizationScaleFactor() << "\n"
	<< "alpha_s = " << lastAlphaS() << "\n" << flush;

	}

	void MatchboxMEBase::logPDFWeight() const {

	if ( !verbose() )
	return;

	generator()->log() << "'" << name() << "' calculated pdf weight = "
	<< lastMEPDFWeight() << " from XComb "
	<< lastXCombPtr() << "\n"
	<< "x1 = " << lastX1() << " (" << (mePartonData()[0]->coloured() ? "" : "not ") << "used) "
	<< "x2 = " << lastX2() << " (" << (mePartonData()[1]->coloured() ? "" : "not ") << "used)\n"
	<< flush;

	}

	void MatchboxMEBase::logME2() const {

	if ( !verbose() )
	return;

	generator()->log() << "'" << name() << "' evaluated me2 using XComb "
	<< lastXCombPtr() << "\n"
	<< "and phase space point (in GeV):\n";

	vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
	cPDVector::const_iterator dit = mePartonData().begin();

	for ( ; pit != meMomenta().end() ; ++pit, ++dit )
	generator()->log() << (**dit).PDGName() << " : "
	<< (*pit/GeV) << "\n";

	generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
	<< "sHat/GeV2 = " << (lastSHat()/GeV2)
	<< " me2 = " << lastME2() << "\n" << flush;

	}

	void MatchboxMEBase::logDSigHatDR() const {

	if ( !verbose() )
	return;

	generator()->log() << "'" << name() << "' evaluated cross section using XComb "
	<< lastXCombPtr() << "\n"
	<< "Jacobian = " << jacobian() << " sHat/GeV2 = "
	<< (lastSHat()/GeV2) << " dsig/nb = "
	<< (lastMECrossSection()/nanobarn) << "\n" << flush;

	}

	void MatchboxMEBase::cloneDependencies(const std::string& prefix) {

	if ( phasespace() ) {
	Ptr<MatchboxPhasespace>::ptr myPhasespace = phasespace()->cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << myPhasespace->name();
	if ( ! (generator()->preinitRegister(myPhasespace,pname.str()) ) )
	throw InitException() << "Phasespace generator " << pname.str() << " already existing.";
	myPhasespace->cloneDependencies(pname.str());
	phasespace(myPhasespace);
	}

	theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());

	if ( matchboxAmplitude() ) {
	Ptr<MatchboxAmplitude>::ptr myAmplitude = matchboxAmplitude()->cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << myAmplitude->name();
	if ( ! (generator()->preinitRegister(myAmplitude,pname.str()) ) )
	throw InitException() << "Amplitude " << pname.str() << " already existing.";
	myAmplitude->cloneDependencies(pname.str());
	matchboxAmplitude(myAmplitude);
	amplitude(myAmplitude);
	matchboxAmplitude()->orderInGs(orderInAlphaS());
	matchboxAmplitude()->orderInGem(orderInAlphaEW());
	}

	if ( scaleChoice() ) {
	Ptr<MatchboxScaleChoice>::ptr myScaleChoice = scaleChoice()->cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << myScaleChoice->name();
	if ( ! (generator()->preinitRegister(myScaleChoice,pname.str()) ) )
	throw InitException() << "Scale choice " << pname.str() << " already existing.";
	scaleChoice(myScaleChoice);
	}

	for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
	theReweights.begin(); rw != theReweights.end(); ++rw ) {
	Ptr<MatchboxReweightBase>::ptr myReweight = (**rw).cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << (**rw).name();
	if ( ! (generator()->preinitRegister(myReweight,pname.str()) ) )
	throw InitException() << "Reweight " << pname.str() << " already existing.";
	myReweight->cloneDependencies(pname.str());
	*rw = myReweight;
	}

	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	Ptr<MatchboxInsertionOperator>::ptr myIOP = (**v).cloneMe();
	ostringstream pname;
	pname << (prefix == "" ? fullName() : prefix) << "/" << (**v).name();
	if ( ! (generator()->preinitRegister(myIOP,pname.str()) ) )
	throw InitException() << "Insertion operator " << pname.str() << " already existing.";
	*v = myIOP;
	}

	}

	void MatchboxMEBase::prepareXComb(MatchboxXCombData& xc) const {

	if ( phasespace() ) {
	size_t nout = diagrams().front()->partons().size()-2;
	xc.nDimPhasespace(phasespace()->nDim(nout));
	}

	if ( matchboxAmplitude() ) {
	xc.nDimAmplitude(matchboxAmplitude()->nDimAdditional());
	if ( matchboxAmplitude()->colourBasis() ) {
	size_t cdim =
	matchboxAmplitude()->colourBasis()->prepare(diagrams(),matchboxAmplitude()->noCorrelations());
	xc.colourBasisDim(cdim);
	}
	}

	int insertionAdd = 0;
	for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
	virtuals().begin(); v != virtuals().end(); ++v ) {
	insertionAdd = max(insertionAdd,(**v).nDimAdditional());
	}
	xc.nDimInsertions(insertionAdd);

	xc.nLight(getNLight());

	if ( initVerbose() ) {
	string fname = name() + ".diagrams";
	ifstream test(fname.c_str());
	if ( !test ) {
	test.close();
	ofstream out(fname.c_str());
	for ( vector<Ptr<DiagramBase>::ptr>::const_iterator d = diagrams().begin();
	d != diagrams().end(); ++d ) {
	DiagramDrawer::drawDiag(out,dynamic_cast<const Tree2toNDiagram&>(**d));
	out << "\n";
	}
	}
	}

	}

	StdXCombPtr MatchboxMEBase::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
	tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
	tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
	const PBPair & newPartonBins, tCutsPtr newCuts,
	const DiagramVector & newDiagrams, bool mir,
	const PartonPairVec&,
	tStdXCombPtr newHead,
	tMEPtr newME) {

	if ( !newME )
	newME = this;

	Ptr<MatchboxXComb>::ptr xc =
	new_ptr(MatchboxXComb(newMaxEnergy, inc,
	newEventHandler, newSubProcessHandler,
	newExtractor, newCKKW,
	newPartonBins, newCuts, newME,
	newDiagrams, mir,
	newHead));

	prepareXComb(*xc);

	return xc;

	}

	StdXCombPtr MatchboxMEBase::makeXComb(tStdXCombPtr newHead,
	const PBPair & newPartonBins,
	const DiagramVector & newDiagrams,
	tMEPtr newME) {
	if ( !newME )
	newME = this;

	Ptr<MatchboxXComb>::ptr xc =
	new_ptr(MatchboxXComb(newHead, newPartonBins, newME, newDiagrams));

	prepareXComb(*xc);

	return xc;
	}

	void MatchboxMEBase::persistentOutput(PersistentOStream & os) const {
	os << theLastXComb << theFactory << thePhasespace
	<< theAmplitude << theScaleChoice << theVirtuals
	<< theReweights << theSubprocesses << theOneLoop
	<< theOneLoopNoBorn;
	}

	void MatchboxMEBase::persistentInput(PersistentIStream & is, int) {
	is >> theLastXComb >> theFactory >> thePhasespace
	>> theAmplitude >> theScaleChoice >> theVirtuals
	>> theReweights >> theSubprocesses >> theOneLoop
	>> theOneLoopNoBorn;
	lastMatchboxXComb(theLastXComb);
	}

	void MatchboxMEBase::Init() {

	static ClassDocumentation<MatchboxMEBase> documentation
	("MatchboxMEBase is the base class for matrix elements "
	"in the context of the matchbox NLO interface.");

	}

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

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

	void MatchboxMEBase::doinit() {
	MEBase::doinit();
	if ( !theAmplitude )
	theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());
	}

	// * Attention * The following static variable is needed for the type
	// description system in ThePEG. Please check that the template arguments
	// are correct (the class and its base class), and that the constructor
	// arguments are correct (the class name and the name of the dynamically
	// loadable library where the class implementation can be found).
	DescribeClass<MatchboxMEBase,MEBase>
	describeHerwigMatchboxMEBase("Herwig::MatchboxMEBase", "HwMatchbox.so");

	diff --git a/src/defaults/Matchbox.in b/src/defaults/Matchbox.in
	--- a/src/defaults/Matchbox.in
	+++ b/src/defaults/Matchbox.in
	@@ -1,368 +1,372 @@
	################################################################################
	# Input file to setup Matchbox NLO matrix elements
	################################################################################





	################################################################################
	#
	# /!\ Nothing to be done below here, except you exactly know
	# what you're doing.
	#
	# Really.
	#
	################################################################################

	# .....
	# ....
	# ...
	# ..
	# .
	#
















	################################################################################
	# Load libraries
	################################################################################

	library HwExsample2.so
	library HwMatchbox.so
	library JetCuts.so
	library FastJetFinder.so

	mkdir /Herwig/MatrixElements/Matchbox
	cd /Herwig/MatrixElements/Matchbox

	################################################################################
	# For the moment, we can only deal with massless partons
	################################################################################

	set /Herwig/Particles/d:NominalMass 0*GeV
	set /Herwig/Particles/dbar:NominalMass 0*GeV
	set /Herwig/Particles/u:NominalMass 0*GeV
	set /Herwig/Particles/ubar:NominalMass 0*GeV
	set /Herwig/Particles/s:NominalMass 0*GeV
	set /Herwig/Particles/sbar:NominalMass 0*GeV
	set /Herwig/Particles/c:NominalMass 0*GeV
	set /Herwig/Particles/cbar:NominalMass 0*GeV
	set /Herwig/Particles/b:NominalMass 0*GeV
	set /Herwig/Particles/bbar:NominalMass 0*GeV

	################################################################################
	# Cleanup some Herwig++ defaults
	################################################################################

	set /Herwig/Cuts/JetKtCut:MinKT 0.0*GeV
	set /Herwig/Cuts/QCDCuts:MHatMin 0.0*GeV
	set /Herwig/Cuts/NeutralCurrentCut:MinW2 0.0*GeV2
	set /Herwig/Cuts/DISCuts:MHatMin 0.0*GeV

	set /Herwig/Generators/LEPGenerator:EventHandler:CollisionCuts Off
	set /Herwig/Generators/DISGenerator:EventHandler:CollisionCuts Off
	set /Herwig/Generators/LHCGenerator:EventHandler:CollisionCuts Off

	################################################################################
	# Integration/sampling
	################################################################################

	mkdir /Herwig/Samplers
	cd /Herwig/Samplers

	create Herwig::ExSampler ExBinSampler
	create Herwig::GeneralSampler ExSampler
	set ExBinSampler:presampling_points 10000
	set ExSampler:BinSampler ExBinSampler
	set ExSampler:Verbose On

	create Herwig::BinSampler FlatBinSampler
	set FlatBinSampler:InitialPoints 10000

	create Herwig::ProjectingSampler ProjectingBinSampler
	set ProjectingBinSampler:InitialPoints 1250
	set ProjectingBinSampler:NIterations 4
	set ProjectingBinSampler:EnhancementFactor 4.0
	set ProjectingBinSampler:NBins 8

	create Herwig::GeneralSampler Sampler
	set Sampler:BinSampler ExBinSampler
	set Sampler:UpdateAfter 1000

	set /Herwig/EventHandlers/LEPHandler:Sampler Sampler
	set /Herwig/EventHandlers/DISHandler:Sampler Sampler
	set /Herwig/EventHandlers/LHCHandler:Sampler Sampler

	################################################################################
	# Setup the factory object
	################################################################################

	cd /Herwig/MatrixElements/Matchbox

	create Herwig::MatchboxFactory Factory

	set Factory:OrderInAlphaEW 2
	set Factory:OrderInAlphaS 0

	do Factory:StartParticleGroup p
	insert Factory:ParticleGroup 0 /Herwig/Particles/b
	insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/c
	insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/s
	insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/d
	insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/u
	insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
	insert Factory:ParticleGroup 0 /Herwig/Particles/g
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup j
	insert Factory:ParticleGroup 0 /Herwig/Particles/b
	insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/c
	insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/s
	insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/d
	insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
	insert Factory:ParticleGroup 0 /Herwig/Particles/u
	insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
	insert Factory:ParticleGroup 0 /Herwig/Particles/g
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup h0
	insert Factory:ParticleGroup 0 /Herwig/Particles/h0
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup e+
	insert Factory:ParticleGroup 0 /Herwig/Particles/e+
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup e-
	insert Factory:ParticleGroup 0 /Herwig/Particles/e-
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup mu+
	insert Factory:ParticleGroup 0 /Herwig/Particles/mu+
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup mu-
	insert Factory:ParticleGroup 0 /Herwig/Particles/mu-
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup l+
	insert Factory:ParticleGroup 0 /Herwig/Particles/e+
	insert Factory:ParticleGroup 0 /Herwig/Particles/mu+
	do Factory:EndParticleGroup

	do Factory:StartParticleGroup l-
	insert Factory:ParticleGroup 0 /Herwig/Particles/e-
	insert Factory:ParticleGroup 0 /Herwig/Particles/mu-
	do Factory:EndParticleGroup

	################################################################################
	# Setup amplitudes
	################################################################################

	cd /Herwig/MatrixElements/Matchbox
	mkdir Amplitudes
	cd Amplitudes

	create Herwig::SimpleColourBasis SimpleColourBasis

	+create Herwig::SimpleColourBasis LargeNColourBasis
	+set LargeNColourBasis:LargeN On
	+
	create Herwig::MatchboxAmplitudellbarqqbar Amplitudellbarqqbar
	set Amplitudellbarqqbar:ColourBasis SimpleColourBasis

	create Herwig::MatchboxAmplitudellbarqqbarg Amplitudellbarqqbarg
	set Amplitudellbarqqbarg:ColourBasis SimpleColourBasis

	create Herwig::MatchboxAmplitudellbarqqbargg Amplitudellbarqqbargg
	set Amplitudellbarqqbargg:ColourBasis SimpleColourBasis

	create Herwig::MatchboxAmplitudellbarqqbarqqbar Amplitudellbarqqbarqqbar
	set Amplitudellbarqqbarqqbar:ColourBasis SimpleColourBasis

	insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbar
	insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbarg
	insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbargg
	insert /Herwig/MatrixElements/Matchbox/Factory:Amplitudes 0 Amplitudellbarqqbarqqbar

	################################################################################
	# Setup utilities for matching
	################################################################################

	cd /Herwig/MatrixElements/Matchbox

	create Herwig::MEMatching MEMatching
	set MEMatching:RestrictPhasespace Off
	set MEMatching:BornScaleInSubtraction RealScale
	set MEMatching:RealEmissionScaleInSubtraction RealScale
	set MEMatching:EmissionScaleInSubtraction RealScale
	set MEMatching:BornScaleInSplitting ShowerScale
	set MEMatching:RealEmissionScaleInSplitting ShowerScale
	set MEMatching:EmissionScaleInSplitting ShowerScale
	set MEMatching:FFPtCut 1.0*GeV
	set MEMatching:FIPtCut 1.0*GeV
	set MEMatching:IIPtCut 1.0*GeV

	create Herwig::ShowerApproximationGenerator MECorrectionHandler
	set MECorrectionHandler:ShowerApproximation MEMatching
	create Herwig::FlatInvertiblePhasespace InvertiblePhasespace
	set MECorrectionHandler:Phasespace InvertiblePhasespace

	create Herwig::DipoleMatching DipoleMatching
	+set DipoleMatching:LargeNBasis /Herwig/MatrixElements/Matchbox/Amplitudes/LargeNColourBasis
	set DipoleMatching:RestrictPhasespace On
	set DipoleMatching:BornScaleInSubtraction BornScale
	set DipoleMatching:RealEmissionScaleInSubtraction BornScale
	set DipoleMatching:EmissionScaleInSubtraction BornScale
	set DipoleMatching:FFPtCut 1.0*GeV
	set DipoleMatching:FIPtCut 1.0*GeV
	set DipoleMatching:IIPtCut 1.0*GeV

	################################################################################
	# Setup utilities for process generation
	################################################################################

	cd /Herwig/MatrixElements/Matchbox
	mkdir Utility
	cd Utility

	create Herwig::Tree2toNGenerator DiagramGenerator
	insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFGVertex
	insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFPVertex
	insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFZVertex
	insert DiagramGenerator:Vertices 0 /Herwig/Vertices/GGGVertex

	create Herwig::ProcessData ProcessData

	set /Herwig/MatrixElements/Matchbox/Factory:DiagramGenerator DiagramGenerator
	set /Herwig/MatrixElements/Matchbox/Factory:ProcessData ProcessData

	################################################################################
	# Setup phasespace generators
	################################################################################

	cd /Herwig/MatrixElements/Matchbox
	mkdir Phasespace
	cd Phasespace

	create Herwig::MatchboxRambo Rambo

	create Herwig::FlatInvertiblePhasespace DetPhasespace

	create Herwig::TreePhasespaceChannels TreePhasespaceChannels
	create Herwig::TreePhasespace TreePhasespace
	set TreePhasespace:ChannelMap TreePhasespaceChannels
	set TreePhasespace:M0 1*GeV
	set TreePhasespace:MC 0.01*GeV

	set /Herwig/MatrixElements/Matchbox/Factory:Phasespace TreePhasespace

	################################################################################
	# Setup cuts
	################################################################################

	cd /Herwig/Cuts

	create ThePEG::FastJetFinder JetFinder
	set JetFinder:UnresolvedMatcher /Herwig/Matchers/StandardQCDParton
	set JetFinder:Variant AntiKt
	set JetFinder:RecombinationScheme E
	set JetFinder:Mode Inclusive
	set JetFinder:ConeRadius 0.7

	create ThePEG::JetRegion FirstJet
	set FirstJet:PtMin 20.*GeV
	do FirstJet:YRange -5.0 5.0
	insert FirstJet:Accepts[0] 1

	create ThePEG::JetRegion SecondJet
	set SecondJet:PtMin 20.*GeV
	do SecondJet:YRange -5.0 5.0
	insert SecondJet:Accepts[0] 2

	create ThePEG::JetRegion ThirdJet
	set ThirdJet:PtMin 20.*GeV
	do ThirdJet:YRange -5.0 5.0
	insert ThirdJet:Accepts[0] 3

	create ThePEG::NJetsCut NJetsCut
	set NJetsCut:UnresolvedMatcher /Herwig/Matchers/StandardQCDParton
	set NJetsCut:NJetsMin 2

	create ThePEG::JetCuts JetCuts
	set JetCuts:UnresolvedMatcher /Herwig/Matchers/StandardQCDParton
	set JetCuts:Ordering OrderPt

	################################################################################
	# Setup scale choices
	################################################################################

	cd /Herwig/MatrixElements/Matchbox
	mkdir Scales
	cd Scales

	create Herwig::MatchboxScaleChoice SHatScale
	cp SHatScale FixedScale
	set FixedScale:FixedScale 100.*GeV
	create Herwig::MatchboxPtScale MaxPtScale
	set MaxPtScale:JetFinder /Herwig/Cuts/JetFinder
	create Herwig::MatchboxLeptonMassScale LeptonMassScale
	create Herwig::MatchboxLeptonPtScale LeptonPtScale

	set /Herwig/MatrixElements/Matchbox/Factory:ScaleChoice LeptonMassScale

	################################################################################
	# AlphaS
	################################################################################

	cd /Herwig/MatrixElements/Matchbox
	mkdir AlphaS
	cd AlphaS

	library HwDipoleShowerAlphaS.so
	create matchbox::lo_alpha_s LOAlphaS
	set LOAlphaS:min_active_flavours 5

	create matchbox::nlo_alpha_s NLOAlphaS
	set NLOAlphaS:exact_evaluation large_scale
	set NLOAlphaS:min_active_flavours 5

	################################################################################
	# Factories for different colliders
	################################################################################

	cd /Herwig/MatrixElements/Matchbox

	cp Factory EEFactory
	set EEFactory:PartonExtractor /Herwig/Partons/EEExtractor
	set EEFactory:Cuts /Herwig/Cuts/EECuts

	cp Factory DISFactory

	set DISFactory:PartonExtractor /Herwig/Partons/DISExtractor
	set DISFactory:Cuts /Herwig/Cuts/DISCuts

	cp Factory PPFactory

	set PPFactory:PartonExtractor /Herwig/Partons/QCDExtractor
	set PPFactory:Cuts /Herwig/Cuts/QCDCuts

	cd /

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 3:01 PM (1 d, 14 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3804892
Default Alt Text: (65 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions