No OneTemporary
Actions

Size

49 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,956 +1,959 @@
	// -- C++ --
	//
	// MatchboxAmplitude.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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/Interface/Parameter.h"
	#include "ThePEG/Interface/Command.h"
	#include "ThePEG/Interface/Switch.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 "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
	#include "ThePEG/Utilities/StringUtils.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(), theCleanupAfter(20),
	treeLevelHelicityPoints(0),
	oneLoopHelicityPoints(0),
	theTrivialColourLegs(false) {
	}

	MatchboxAmplitude::~MatchboxAmplitude() {}

	void MatchboxAmplitude::persistentOutput(PersistentOStream & os) const {
	os << theLastXComb << theColourBasis << theFactory
	<< theCleanupAfter << treeLevelHelicityPoints << oneLoopHelicityPoints
	<< theTrivialColourLegs << theReshuffleMasses.size();
	if ( !theReshuffleMasses.empty() ) {
	for (auto const & r : theReshuffleMasses )
	os << r.first << ounit(r.second,GeV);
	}
	}

	void MatchboxAmplitude::persistentInput(PersistentIStream & is, int) {
	size_t reshuffleSize;
	is >> theLastXComb >> theColourBasis >> theFactory
	>> theCleanupAfter >> treeLevelHelicityPoints >> oneLoopHelicityPoints
	>> theTrivialColourLegs >> reshuffleSize;
	theReshuffleMasses.clear();
	while ( reshuffleSize > 0 ) {
	long id; Energy m;
	is >> id >> iunit(m,GeV);
	theReshuffleMasses[id] = m;
	--reshuffleSize;
	}
	lastMatchboxXComb(theLastXComb);
	}

	Ptr<MatchboxFactory>::tptr MatchboxAmplitude::factory() const {
	return theFactory;
	}

	void MatchboxAmplitude::factory(Ptr<MatchboxFactory>::tptr f) {
	theFactory = f;
	}

	void MatchboxAmplitude::doinit() {
	Amplitude::doinit();
	if ( colourBasis() ) {
	colourBasis()->factory(factory());
	colourBasis()->init();
	}
	}

	void MatchboxAmplitude::doinitrun() {
	Amplitude::doinitrun();
	if ( colourBasis())
	colourBasis()->initrun();
	}

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

	MatchboxMEBasePtr MatchboxAmplitude::makeME(const PDVector&) const {
	return new_ptr(MatchboxMEBase());
	}

	Selector<const ColourLines *> MatchboxAmplitude::colourGeometries(tcDiagPtr d) const {
	if ( haveColourFlows() )
	return colourBasis()->colourGeometries(d,lastLargeNAmplitudes());
	return Selector<const ColourLines *>();
	}

	void MatchboxAmplitude::olpOrderFileHeader(ostream& os) const {

	os << "# OLP order file created by Herwig/Matchbox\n\n";

	os << "InterfaceVersion BLHA2\n\n";

	os << "Model SM\n"
	<< "CorrectionType QCD\n"
	<< "IRregularisation " << (isDR() ? "DRED" : "CDR") << "\n"
	<< "Extra HelAvgInitial no\n"
	<< "Extra ColAvgInitial no\n"
	<< "Extra MCSymmetrizeFinal no\n";

	os << "\n";

	}

	void MatchboxAmplitude::olpOrderFileProcesses(ostream& os,
	const map<pair<Process,int>,int>& proc) const {

	map<int,pair<Process,int> > sorted;

	for (auto const & p : proc ) {
	sorted[p.second] = p.first;
	}

	unsigned int currentOrderInAlphaS = sorted.begin()->second.first.orderInAlphaS;
	unsigned int currentOrderInAlphaEW = sorted.begin()->second.first.orderInAlphaEW;
	int currentType = sorted.begin()->second.second;

	os << "AlphasPower " << currentOrderInAlphaS << "\n"
	<< "AlphaPower " << currentOrderInAlphaEW << "\n"
	<< "AmplitudeType ";
	if ( currentType == ProcessType::treeME2 ) {
	os << "tree\n";
	} else if ( currentType == ProcessType::oneLoopInterference ) {
	os << "loop\n";
	} else if ( currentType == ProcessType::colourCorrelatedME2 ) {
	os << "cctree\n";
	} else if ( currentType == ProcessType::spinColourCorrelatedME2 ) {
	os << "sctree\n";
	} else if ( currentType == ProcessType::loopInducedME2 ) {
	os << "loopinduced\n";
	} else if ( currentType == ProcessType::spinCorrelatedME2 ) {
	os << "stree\n";
	} else assert(false);


	for ( map<int,pair<Process,int> >::const_iterator p = sorted.begin();
	p != sorted.end(); ++p ) {

	if ( currentOrderInAlphaS != p->second.first.orderInAlphaS ) {
	currentOrderInAlphaS = p->second.first.orderInAlphaS;
	os << "AlphasPower " << currentOrderInAlphaS << "\n";
	}

	if ( currentOrderInAlphaEW != p->second.first.orderInAlphaEW ) {
	currentOrderInAlphaEW = p->second.first.orderInAlphaEW;
	os << "AlphaPower " << currentOrderInAlphaEW << "\n";
	}

	if ( currentType != p->second.second ) {
	currentType = p->second.second;
	os << "AmplitudeType ";
	if ( currentType == ProcessType::treeME2 ) {
	os << "tree\n";
	} else if ( currentType == ProcessType::oneLoopInterference ) {
	os << "loop\n";
	} else if ( currentType == ProcessType::colourCorrelatedME2 ) {
	os << "cctree\n";
	} else if ( currentType == ProcessType::spinColourCorrelatedME2 ) {
	os << "sctree\n";
	} else if ( currentType == ProcessType::spinCorrelatedME2 ) {
	os << "stree\n";
	} else assert(false);
	}

	os << p->second.first.legs[0]->id() << " "
	<< p->second.first.legs[1]->id() << " -> ";
	for ( PDVector::const_iterator o = p->second.first.legs.begin() + 2;
	o != p->second.first.legs.end(); ++o ) {
	os << (**o).id() << " ";
	}
	os << "\n";

	}

	}

	bool MatchboxAmplitude::startOLP(const map<pair<Process,int>,int>& procs) {

	string orderFileName = factory()->buildStorage() + name() + ".OLPOrder.lh";
	ofstream orderFile(orderFileName.c_str());

	olpOrderFileHeader(orderFile);
	olpOrderFileProcesses(orderFile,procs);

	string contractFileName = factory()->buildStorage() + name() + ".OLPContract.lh";

	signOLP(orderFileName, contractFileName);

	// TODO check the contract file

	int status = 0;
	startOLP(contractFileName, status);

	if ( status != 1 )
	return false;

	return true;

	}


	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();
	- if ( treeAmplitudes() \|\| oneLoopAmplitudes() )
	+ if ( ( treeAmplitudes() \|\| oneLoopAmplitudes() ) &&
	+ hasAmplitudeMomenta() )
	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());
	+
	+ if ( hasAmplitudeMomenta() )
	+ 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;
	}


	//This happens for e.g. p p-> W- gamma & p p->W- W- j j
	//Still working for pp->W-H-W- e+ nue jj ???
	if(next == processLegs.end()){
	next = processLegs.begin();
	for (;next!=processLegs.end();next++){
	assert(next->first->id() < 0 );
	crossingMap()[ampCount] = next->second - shift;
	amplitudeLegs.insert(make_pair(next->first,ampCount));
	++ampCount;
	processLegs.erase(next);
	}
	break;
	}

	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 ( processLegs.empty() ) break;
	if ( checkcc == processLegs.end() ) {
	checkcc = processLegs.begin();
	while ( checkcc->first->id() > 0 )
	if ( ++checkcc == processLegs.end() )
	break;
	}
	// if still not there, use whatever is available at the end
	if ( checkcc == processLegs.end() )
	checkcc = processLegs.begin();
	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() && !colourBasis()->indexMap().empty()) {
	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 set<vector<size_t> >& MatchboxAmplitude::colourOrdering(size_t id) const {

	static set<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());
	Energy2 rho = res.t()*res.t() - res.mass2();
	res.setRho(sqrt(abs(rho)));
	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 ) {
	current[pos] = 0;
	doGenerateHelicities(res,current,pos+1);
	} else if ( amplitudePartonData()[pos]->iSpin() == PDT::Spin1Half ) {
	current[pos] = 1;
	doGenerateHelicities(res,current,pos+1);
	current[pos] = -1;
	doGenerateHelicities(res,current,pos+1);
	}else if (amplitudePartonData()[pos]->iSpin() == PDT::Spin1 ) {
	if (amplitudePartonData()[pos]->hardProcessMass() != ZERO){
	current[pos] = 0;
	doGenerateHelicities(res,current,pos+1);
	}
	current[pos] = 1;
	doGenerateHelicities(res,current,pos+1);
	current[pos] = -1;
	doGenerateHelicities(res,current,pos+1);
	}
	}

	vector<unsigned int> MatchboxAmplitude::physicalHelicities(const vector<int>&) const {
	throw Exception()
	<< "MatchboxAmplitude::physicalHelicities(): The amplitude '" << name() << "' does not support the spin correlation algorithm"
	<< Exception::runerror;
	static vector<unsigned int> dummy;
	return dummy;
	}

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

	if ( !calculateTreeAmplitudes() )
	return;

	bool initialized =
	!lastAmplitudes().empty() && treeLevelHelicityPoints > theCleanupAfter;

	if ( !initialized ) {
	treeLevelHelicityPoints++;
	map<vector<int>,CVector> all;
	map<vector<int>,CVector> allLargeN;
	set<vector<int> > helicities = generateHelicities();
	for ( set<vector<int> >::const_iterator h = helicities.begin();
	h != helicities.end(); ++h ) {
	all.insert(make_pair(*h,CVector(max(colourBasisDim(),1))));
	allLargeN.insert(make_pair(*h,CVector(max(colourBasisDim(),1))));
	}
	AmplitudeIterator amp = all.begin();
	AmplitudeIterator lamp = allLargeN.begin();
	for ( ; amp != all.end(); ++amp, ++lamp ) {
	for ( size_t k = 0; k < max(colourBasisDim(),1); ++k ){
	amp->second(k) = evaluate(k,amp->first,lamp->second(k));
	if ( amp->second(k) != Complex(0.0) ) {
	if ( lastAmplitudes().find(amp->first)!=lastAmplitudes().end() ) {
	lastAmplitudes().find(amp->first)->second = amp->second;
	lastLargeNAmplitudes().find(lamp->first)->second = lamp->second;
	} else {
	lastAmplitudes().insert(*amp);
	lastLargeNAmplitudes().insert(*lamp);
	}
	} else if ( lastAmplitudes().find(amp->first)!=lastAmplitudes().end() ){
	lastAmplitudes().find(amp->first)->second = amp->second;
	lastLargeNAmplitudes().find(lamp->first)->second = lamp->second;
	}
	}
	}
	} else {
	AmplitudeIterator amp = lastAmplitudes().begin();
	AmplitudeIterator lamp = lastLargeNAmplitudes().begin();
	for ( ;amp != lastAmplitudes().end(); ++amp, ++lamp ) {
	for ( size_t k = 0; k < max(colourBasisDim(),1); ++k ){
	amp->second(k) = evaluate(k,amp->first,lamp->second(k));
	}
	}
	}

	haveTreeAmplitudes();

	}

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

	if ( !calculateOneLoopAmplitudes() )
	return;

	bool initialized =
	!lastOneLoopAmplitudes().empty() && oneLoopHelicityPoints > theCleanupAfter;

	if ( !initialized ) {
	oneLoopHelicityPoints++;
	map<vector<int>,CVector> all;
	set<vector<int> > helicities = generateHelicities();
	for ( set<vector<int> >::const_iterator h = helicities.begin();
	h != helicities.end(); ++h ) {
	all.insert(make_pair(*h,CVector(max(colourBasisDim(),1))));
	}
	AmplitudeIterator amp = all.begin();
	for ( ; amp != all.end(); ++amp ) {
	for ( size_t k = 0; k < max(colourBasisDim(),1); ++k ){
	amp->second(k) = evaluateOneLoop(k,amp->first);
	if ( amp->second(k) != Complex(0.0) ) {
	if ( lastOneLoopAmplitudes().find(amp->first)!=lastOneLoopAmplitudes().end() ) {
	lastOneLoopAmplitudes().find(amp->first)->second = amp->second;
	} else{
	lastOneLoopAmplitudes().insert(*amp);
	}
	} else if ( lastOneLoopAmplitudes().find(amp->first)!=lastOneLoopAmplitudes().end() ){
	lastOneLoopAmplitudes().find(amp->first)->second = amp->second;
	}
	}
	}
	} else {
	AmplitudeIterator amp = lastOneLoopAmplitudes().begin();
	for ( ;amp != lastOneLoopAmplitudes().end(); ++amp ) {
	for ( size_t k = 0; k < max(colourBasisDim(),1); ++k ){
	amp->second(k) = evaluateOneLoop(k,amp->first);
	}
	}
	}

	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() << "MatchboxAmplitude::value(): ThePEG::Amplitude interface is not sufficient at the moment."
	<< Exception::runerror;
	return 0.;
	}

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

	double MatchboxAmplitude::largeNME2(Ptr<ColourBasis>::tptr largeNBasis) const {
	if ( !calculateLargeNME2() )
	return lastLargeNME2();
	double res = 0.;
	if ( !trivialColourLegs() )
	res = largeNBasis->me2(mePartonData(),lastLargeNAmplitudes());
	else
	res = me2();
	lastLargeNME2(res);
	return res;
	}


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

	double MatchboxAmplitude::colourCharge(tcPDPtr data) const {
	double cfac = 1.;
	double Nc = generator()->standardModel()->Nc();
	if ( data->iColour() == PDT::Colour8 ) {
	cfac = Nc;
	} else if ( data->iColour() == PDT::Colour3 \|\|
	data->iColour() == PDT::Colour3bar ) {
	cfac = (sqr(Nc)-1.)/(2.*Nc);
	} else assert(false);
	return cfac;
	}

	double MatchboxAmplitude::largeNColourCharge(tcPDPtr data) const {
	double cfac = 1.;
	double Nc = generator()->standardModel()->Nc();
	if ( data->iColour() == PDT::Colour8 ) {
	cfac = Nc;
	} else if ( data->iColour() == PDT::Colour3 \|\|
	data->iColour() == PDT::Colour3bar ) {
	cfac = 0.5*Nc;
	} else assert(false);
	return cfac;
	}

	double MatchboxAmplitude::colourCorrelatedME2(pair<int,int> ij) const {
	double cfac = colourCharge(mePartonData()[ij.first]);
	if ( !calculateColourCorrelator(ij) )
	return lastColourCorrelator(ij)/cfac;
	double res = 0.;
	if ( !trivialColourLegs() )
	res = colourBasis()->colourCorrelatedME2(ij,mePartonData(),lastAmplitudes());
	else {
	// two partons TiTj = (-Ti^2-Tj^2)/2
	// three partons TiTj = (-Ti^2-Tj^2+Tk^2)/2
	int n = mePartonData().size();
	for ( int i = 0; i < n; ++i ) {
	if ( !mePartonData()[i]->coloured() )
	continue;
	if ( i == ij.first \|\| i == ij.second )
	res -= colourCharge(mePartonData()[i])*me2();
	else
	res += colourCharge(mePartonData()[i])*me2();
	}
	res *= 0.5;
	}
	lastColourCorrelator(ij,res);
	return res/cfac;
	}

	double MatchboxAmplitude::largeNColourCorrelatedME2(pair<int,int> ij,
	Ptr<ColourBasis>::tptr largeNBasis) const {
	double cfac = largeNColourCharge(mePartonData()[ij.first]);
	if ( !calculateLargeNColourCorrelator(ij) )
	return lastLargeNColourCorrelator(ij)/cfac;
	double res = 0.;
	if ( !trivialColourLegs() )
	res = largeNBasis->colourCorrelatedME2(ij,mePartonData(),lastLargeNAmplitudes());
	else {
	// two partons TiTj = (-Ti^2-Tj^2)/2
	// three partons TiTj = (-Ti^2-Tj^2+Tk^2)/2
	int n = mePartonData().size();
	for ( int i = 0; i < n; ++i ) {
	if ( !mePartonData()[i]->coloured() )
	continue;
	if ( i == ij.first \|\| i == ij.second )
	res -= largeNColourCharge(mePartonData()[i])*largeNME2(largeNBasis);
	else
	res += largeNColourCharge(mePartonData()[i])*largeNME2(largeNBasis);
	}
	res *= 0.5;
	}
	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,
	int) 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,ij.first);

	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 + (c.scale() > ZERO ? 1. : -1.)*corr/cfac;

	}

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

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

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

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

	double avg =
	me2()(-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 ( calculateSpinCorrelator(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()->interference(mePartonData(),a->second,b->second);
	}
	lastSpinCorrelator(ij,csCorr);
	} else {
	csCorr = lastSpinCorrelator(ij);
	}

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

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

	}

	void MatchboxAmplitude::checkReshuffling(Ptr<MatchboxPhasespace>::tptr ps) {
	set<long> noReshuffle;
	for ( map<long,Energy>::const_iterator m = reshuffleMasses().begin();
	m != reshuffleMasses().end(); ++m ) {
	tcPDPtr data = getParticleData(m->first);
	assert(data);
	bool needReshuffle = m->second != data->hardProcessMass();
	needReshuffle \|=
	(data->hardProcessWidth() != ZERO \|\| data->massGenerator()) &&
	ps->useMassGenerators();
	if ( !needReshuffle )
	noReshuffle.insert(m->first);
	}
	for ( set<long>::const_iterator rm = noReshuffle.begin();
	rm != noReshuffle.end(); ++rm )
	theReshuffleMasses.erase(*rm);
	}

	string MatchboxAmplitude::doReshuffle(string in) {
	in = StringUtils::stripws(in);
	if ( in.empty() )
	throw Exception() << "MatchboxAmplitude::doReshuffle(): Expecting PDG id and mass value"
	<< Exception::runerror;
	istringstream ins(in);
	long id;
	ins >> id;
	if ( ins.eof() )
	throw Exception() << "MatchboxAmplitude::doReshuffle(): expecting PDG id and mass value"
	<< Exception::runerror;
	Energy m;
	ins >> iunit(m,GeV);
	theReshuffleMasses[id] = m;
	return "";
	}

	string MatchboxAmplitude::doMassless(string in) {
	in = StringUtils::stripws(in);
	if ( in.empty() )
	throw Exception() << "MatchboxAmplitude::doMassless(): Expecting PDG id"
	<< Exception::runerror;
	istringstream ins(in);
	long id;
	ins >> id;
	theReshuffleMasses[id] = ZERO;
	return "";
	}

	string MatchboxAmplitude::doOnShell(string in) {
	in = StringUtils::stripws(in);
	if ( in.empty() )
	throw Exception() << "MatchboxAmplitude::doOnShell(): Expecting PDG id"
	<< Exception::runerror;
	istringstream ins(in);
	long id;
	ins >> id;
	tcPDPtr data = getParticleData(id);
	assert(data);
	theReshuffleMasses[id] = data->hardProcessMass();
	return "";
	}

	string MatchboxAmplitude::doClearReshuffling(string) {
	theReshuffleMasses.clear();
	return "";
	}

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

	static Parameter<MatchboxAmplitude,int> interfaceCleanupAfter
	("CleanupAfter",
	"The number of points after which helicity combinations are cleaned up.",
	&MatchboxAmplitude::theCleanupAfter, 20, 1, 0,
	false, false, Interface::lowerlim);

	static Command<MatchboxAmplitude> interfaceReshuffle
	("Reshuffle",
	"Reshuffle the mass for the given PDG id to a different mass shell for amplitude evaluation.",
	&MatchboxAmplitude::doReshuffle, false);

	static Command<MatchboxAmplitude> interfaceMassless
	("Massless",
	"Reshuffle the mass for the given PDG id to be massless for amplitude evaluation.",
	&MatchboxAmplitude::doMassless, false);

	static Command<MatchboxAmplitude> interfaceOnShell
	("OnShell",
	"Reshuffle the mass for the given PDG id to be the on-shell mass for amplitude evaluation.",
	&MatchboxAmplitude::doOnShell, false);

	static Command<MatchboxAmplitude> interfaceClearReshuffling
	("ClearReshuffling",
	"Do not perform any reshuffling.",
	&MatchboxAmplitude::doClearReshuffling, false);

	static Switch<MatchboxAmplitude,bool> interfaceTrivialColourLegs
	("TrivialColourLegs",
	"Assume the process considered has trivial colour correllations.",
	&MatchboxAmplitude::theTrivialColourLegs, false, false, false);
	static SwitchOption interfaceTrivialColourLegsYes
	(interfaceTrivialColourLegs,
	"Yes",
	"",
	true);
	static SwitchOption interfaceTrivialColourLegsNo
	(interfaceTrivialColourLegs,
	"No",
	"",
	false);
	interfaceTrivialColourLegs.rank(-1);

	}

	// * 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", "Herwig.so");
	diff --git a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
	--- a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
	+++ b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
	@@ -1,719 +1,724 @@
	// -- C++ --
	//
	// MatchboxAmplitude.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_MatchboxAmplitude_H
	#define HERWIG_MatchboxAmplitude_H
	//
	// This is the declaration of the MatchboxAmplitude class.
	//

	#include "ThePEG/MatrixElement/Amplitude.h"
	#include "ThePEG/Handlers/LastXCombInfo.h"
	#include "Herwig/MatrixElement/Matchbox/Utility/ColourBasis.h"
	#include "Herwig/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
	#include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
	#include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXComb.h"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
	#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
	#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief Process information with coupling order
	*/
	struct Process {

	PDVector legs;
	unsigned int orderInAlphaS;
	unsigned int orderInAlphaEW;

	Process()
	: orderInAlphaS(0), orderInAlphaEW(0) {}

	Process(const PDVector& p,
	unsigned int oas,
	unsigned int oae)
	: legs(p), orderInAlphaS(oas), orderInAlphaEW(oae) {}

	bool operator==(const Process& other) const {
	return
	legs == other.legs &&
	orderInAlphaS == other.orderInAlphaS &&
	orderInAlphaEW == other.orderInAlphaEW;
	}

	bool operator<(const Process& other) const {
	if ( orderInAlphaS != other.orderInAlphaS )
	return orderInAlphaS < other.orderInAlphaS;
	if ( orderInAlphaEW != other.orderInAlphaEW )
	return orderInAlphaEW < other.orderInAlphaEW;
	return legs < other.legs;
	}

	void persistentOutput(PersistentOStream & os) const {
	os << legs << orderInAlphaS << orderInAlphaEW;
	}

	void persistentInput(PersistentIStream & is) {
	is >> legs >> orderInAlphaS >> orderInAlphaEW;
	}

	};

	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief Enumerate the type of calculation required
	*/
	namespace ProcessType {

	enum Types {

	treeME2 = 0,
	colourCorrelatedME2,
	spinColourCorrelatedME2,
	oneLoopInterference,
	loopInducedME2,
	spinCorrelatedME2

	};

	}


	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief MatchboxAmplitude is the base class for amplitude
	* implementations inside Matchbox.
	*
	* @see \ref MatchboxAmplitudeInterfaces "The interfaces"
	* defined for MatchboxAmplitude.
	*/
	class MatchboxAmplitude:
	public Amplitude,
	public LastXCombInfo<StandardXComb>,
	public LastMatchboxXCombInfo {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	MatchboxAmplitude();

	/**
	* The destructor.
	*/
	virtual ~MatchboxAmplitude();
	//@}

	public:

	/**
	* Return the amplitude. Needs to be implemented from
	* ThePEG::Amplitude but is actually ill-defined, as colours of the
	* external particles are not specified. To this extent, this
	* implementation just asserts.
	*/
	virtual Complex value(const tcPDVector & particles,
	const vector<Lorentz5Momentum> & momenta,
	const vector<int> & helicities);

	/**
	* Return the factory which produced this matrix element
	*/
	Ptr<MatchboxFactory>::tptr factory() const;

	/**
	* Set the factory which produced this matrix element
	*/
	virtual void factory(Ptr<MatchboxFactory>::tptr f);

	/** @name Subprocess information */
	//@{

	/**
	* Return true, if this amplitude can handle the given process.
	*/
	virtual bool canHandle(const PDVector& p,
	Ptr<MatchboxFactory>::tptr,
	bool) const { return canHandle(p); }

	/**
	* Return true, if this amplitude can handle the given process.
	*/
	virtual bool canHandle(const PDVector&) const { return false; }

	/**
	* Return the number of random numbers required to evaluate this
	* amplitude at a fixed phase space point.
	*/
	virtual int nDimAdditional() const { return 0; }

	/**
	* Return a ME instance appropriate for this amplitude and the given
	* subprocesses
	*/
	virtual Ptr<MatchboxMEBase>::ptr makeME(const PDVector&) const;

	/**
	* Set the (tree-level) order in \f$g_S\f$ in which this matrix
	* element should be evaluated.
	*/
	virtual void orderInGs(unsigned int) {}

	/**
	* Return the (tree-level) order in \f$g_S\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInGs() const = 0;

	/**
	* Set the (tree-level) order in \f$g_{EM}\f$ in which this matrix
	* element should be evaluated.
	*/
	virtual void orderInGem(unsigned int) {}

	/**
	* Return the (tree-level) order in \f$g_{EM}\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInGem() const = 0;

	/**
	* Return the Herwig StandardModel object
	*/
	Ptr<StandardModel>::tcptr standardModel() {
	if ( !hwStandardModel() )
	hwStandardModel(dynamic_ptr_cast<Ptr<StandardModel>::tcptr>(HandlerBase::standardModel()));
	return hwStandardModel();
	}

	/**
	* Return true, if this amplitude already includes averaging over
	* incoming parton's quantum numbers.
	*/
	virtual bool hasInitialAverage() const { return false; }

	/**
	* Return true, if this amplitude already includes symmetry factors
	* for identical outgoing particles.
	*/
	virtual bool hasFinalStateSymmetry() const { return false; }

	/**
	* Return true, if this amplitude is handled by a BLHA one-loop provider
	*/
	virtual bool isOLPTree() const { return false; }

	/**
	* Return true, if this amplitude is handled by a BLHA one-loop provider
	*/
	virtual bool isOLPLoop() const { return false; }

	/**
	* Return true, if colour and spin correlated matrix elements should
	* be ordered from the OLP
	*/
	virtual bool needsOLPCorrelators() const { return true; }

	/**
	* Write the order file header
	*/
	virtual void olpOrderFileHeader(ostream&) const;

	/**
	* Write the order file process list
	*/
	virtual void olpOrderFileProcesses(ostream&,
	const map<pair<Process,int>,int>& procs) const;

	/**
	* Start the one loop provider, if appropriate, giving order and
	* contract files
	*/
	virtual void signOLP(const string&, const string&) { }

	/**
	* Start the one loop provider, if appropriate
	*/
	virtual void startOLP(const string&, int& status) { status = -1; }

	/**
	* Start the one loop provider, if appropriate. This default
	* implementation writes an BLHA 2.0 order file and starts the OLP
	*/
	virtual bool startOLP(const map<pair<Process,int>,int>& procs);

	/**
	* Return true, if this amplitude needs to initialize an external
	* code.
	*/
	virtual bool isExternal() const { return false; }

	/**
	* Initialize this amplitude
	*/
	virtual bool initializeExternal() { return false; }

	/**
	* Return a generic process id for the given process
	*/
	virtual int externalId(const cPDVector&) { return 0; }

	/**
	* Return the map with masses to be used for amplitude evaluation
	*/
	const map<long,Energy>& reshuffleMasses() const { return theReshuffleMasses; }

	/**
	* Check if reshuffling is needed at all
	*/
	void checkReshuffling(Ptr<MatchboxPhasespace>::tptr);

	+ /**
	+ * Return true, if this amplitude makes use of amplitudeMomenta
	+ */
	+ virtual bool hasAmplitudeMomenta() const { return false; }
	+
	//@}

	/** @name Colour basis. */
	//@{

	/**
	* Return the colour basis.
	*/
	virtual Ptr<ColourBasis>::tptr colourBasis() const { return theColourBasis; }

	/**
	* Return true, if the colour basis is capable of assigning colour
	* flows.
	*/
	virtual bool haveColourFlows() const {
	return colourBasis() ? colourBasis()->haveColourFlows() : false;
	}

	/**
	* Return a Selector with possible colour geometries for the selected
	* diagram weighted by their relative probabilities.
	*/
	virtual Selector<const ColourLines *> colourGeometries(tcDiagPtr diag) const;

	/**
	* Return an ordering identifier for the current subprocess and
	* colour absis tensor index.
	*/
	const string& colourOrderingString(size_t id) const;

	/**
	* Return an ordering identifier for the current subprocess and
	* colour absis tensor index.
	*/
	const set<vector<size_t> >& colourOrdering(size_t id) const;

	//@}

	/** @name Phasespace point, crossing and helicities */
	//@{

	/**
	* Set the xcomb object.
	*/
	virtual void setXComb(tStdXCombPtr xc);

	/**
	* Return the momentum as crossed appropriate for this amplitude.
	*/
	Lorentz5Momentum amplitudeMomentum(int) const;

	/**
	* Perform a normal ordering of external legs and fill the
	* crossing information as. This default implementation sorts
	* lexicographically in (abs(colour)/spin/abs(charge)), putting pairs
	* of particles/anti-particles where possible.
	*/
	virtual void fillCrossingMap(size_t shift = 0);

	/**
	* Generate the helicity combinations.
	*/
	virtual set<vector<int> > generateHelicities() const;

	/**
	* Return the helicity combination of the physical process in the
	* conventions used by the spin correlation algorithm.
	*/
	virtual vector<unsigned int> physicalHelicities(const vector<int>&) const;

	//@}

	/** @name Tree-level amplitudes */
	//@{

	/**
	* Calculate the tree level amplitudes for the phasespace point
	* stored in lastXComb.
	*/
	virtual void prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr);

	/**
	* Return the matrix element squared.
	*/
	virtual double me2() const;

	/**
	* Return the colour charge of a given leg
	*/
	double colourCharge(tcPDPtr) const;

	/**
	* Return the large-N charge of a given leg
	*/
	double largeNColourCharge(tcPDPtr) const;

	/**
	* Return the largeN matrix element squared.
	*/
	virtual double largeNME2(Ptr<ColourBasis>::tptr largeNBasis) const;

	/**
	* Return the colour correlated matrix element.
	*/
	virtual double colourCorrelatedME2(pair<int,int> ij) const;

	/**
	* Return the large-N colour correlated matrix element.
	*/
	virtual double largeNColourCorrelatedME2(pair<int,int> ij,
	Ptr<ColourBasis>::tptr largeNBasis) const;

	/**
	* Return true if trivial colour configuration.
	*/
	bool trivialColourLegs() const { return theTrivialColourLegs; }

	/**
	* Return true, if this amplitude is capable of consistently filling
	* the rho matrices for the spin correllations
	*/
	virtual bool canFillRhoMatrix() const { return false; }

	/**
	* Return a positive helicity polarization vector for a gluon of
	* momentum p (with reference vector n) to be used when evaluating
	* spin correlations.
	*/
	virtual LorentzVector<Complex> plusPolarization(const Lorentz5Momentum& p,
	const Lorentz5Momentum& n,
	int id = -1) const;

	/**
	* Return the colour and spin correlated matrix element.
	*/
	virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
	const SpinCorrelationTensor& c) const;

	/**
	* Return the spin correlated matrix element.
	*/
	virtual double spinCorrelatedME2(pair<int,int> emitterSpectator,
	const SpinCorrelationTensor& c) const;


	/**
	* Return true, if tree-level contributions will be evaluated at amplitude level.
	*/
	virtual bool treeAmplitudes() const { return true; }

	/**
	* Evaluate the amplitude for the given colour tensor id and
	* helicity assignment
	*/
	virtual Complex evaluate(size_t, const vector<int>&, Complex&) { return 0.; }

	//@}

	/** @name One-loop amplitudes */
	//@{

	/**
	* Return the one-loop amplitude, if applicable.
	*/
	virtual Ptr<MatchboxAmplitude>::tptr oneLoopAmplitude() const {
	return Ptr<MatchboxAmplitude>::tptr();
	}

	/**
	* Diasble one-loop functionality if not needed.
	*/
	virtual void disableOneLoop() {}

	/**
	* Return true, if this amplitude is capable of calculating one-loop
	* (QCD) corrections.
	*/
	virtual bool haveOneLoop() const { return false; }

	/**
	* Return true, if this amplitude only provides
	* one-loop (QCD) corrections.
	*/
	virtual bool onlyOneLoop() const { return false; }

	/**
	* Return true, if one-loop contributions will be evaluated at amplitude level.
	*/
	virtual bool oneLoopAmplitudes() const { return true; }

	/**
	* Return true, if one loop corrections have been calculated in
	* dimensional reduction. Otherwise conventional dimensional
	* regularization is assumed. Note that renormalization is always
	* assumed to be MSbar.
	*/
	virtual bool isDR() const { return false; }

	/**
	* Return true, if the amplitude is DRbar renormalized, otherwise
	* MSbar is assumed.
	*/
	virtual bool isDRbar() const { return true; }

	/**
	* Return true, if one loop corrections are given in the conventions
	* of the integrated dipoles.
	*/
	virtual bool isCS() const { return false; }

	/**
	* Return true, if one loop corrections are given in the conventions
	* of BDK.
	*/
	virtual bool isBDK() const { return false; }

	/**
	* Return true, if one loop corrections are given in the conventions
	* of everything expanded.
	*/
	virtual bool isExpanded() const { return false; }

	/**
	* Return the value of the dimensional regularization
	* parameter. Note that renormalization scale dependence is fully
	* restored in DipoleIOperator.
	*/
	virtual Energy2 mu2() const { return 0.*GeV2; }

	/**
	* Indicate that this amplitude is running alphas by itself.
	*/
	virtual bool hasRunningAlphaS() const { return false; }

	/**
	* Indicate that this amplitude is running alphaew by itself.
	*/
	virtual bool hasRunningAlphaEW() const { return false; }

	/**
	* If defined, return the coefficient of the pole in epsilon^2
	*/
	virtual double oneLoopDoublePole() const { return 0.; }

	/**
	* If defined, return the coefficient of the pole in epsilon
	*/
	virtual double oneLoopSinglePole() const { return 0.; }

	/**
	* Calculate the one-loop amplitudes for the phasespace point
	* stored in lastXComb, if provided.
	*/
	virtual void prepareOneLoopAmplitudes(Ptr<MatchboxMEBase>::tcptr);

	/**
	* Return the one-loop/tree interference.
	*/
	virtual double oneLoopInterference() const;

	/**
	* Evaluate the amplitude for the given colour tensor id and
	* helicity assignment
	*/
	virtual Complex evaluateOneLoop(size_t, const vector<int>&) { return 0.; }

	//@}

	/** @name Caching and helpers to setup amplitude objects. */
	//@{

	/**
	* Flush all cashes.
	*/
	virtual void flushCaches() {}

	/**
	* Clone this amplitude.
	*/
	Ptr<MatchboxAmplitude>::ptr cloneMe() const {
	return dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(clone());
	}

	/**
	* Clone the dependencies, using a given prefix.
	*/
	virtual void cloneDependencies(const std::string& prefix="" , bool slim=false);

	//@}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	// If needed, insert declarations of virtual function defined in the
	// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).

	protected:

	/** @name Standard Interfaced functions. */
	//@{

	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit();

	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();
	//@}

	private:

	/**
	* The factory which produced this matrix element
	*/
	Ptr<MatchboxFactory>::tptr theFactory;

	/**
	* Recursively generate helicities
	*/
	void doGenerateHelicities(set<vector<int> >& res,
	vector<int>& current,
	size_t pos) const;

	/**
	* The colour basis implementation to be used.
	*/
	Ptr<ColourBasis>::ptr theColourBasis;

	/**
	* The number of points after which helicity combinations wil be
	* cleaned up
	*/
	int theCleanupAfter;

	/**
	* The number of points that are calculated before a certain
	* helicity is excluded. Needed in pp->V
	*/
	int treeLevelHelicityPoints;

	/**
	* The number of points that are calculated before a certain
	* helicity is excluded. Needed in pp->V
	*/
	int oneLoopHelicityPoints;

	/**
	* The map with masses to be used for amplitude evaluation
	*/
	map<long,Energy> theReshuffleMasses;

	/**
	* True if trivial colour configuration.
	*/
	bool theTrivialColourLegs;

	/**
	* A command to fill the reshuffle mass map
	*/
	string doReshuffle(string);

	/**
	* A command to fill the reshuffle mass map
	*/
	string doMassless(string);

	/**
	* A command to fill the reshuffle mass map
	*/
	string doOnShell(string);

	/**
	* Clear the reshuffling map
	*/
	string doClearReshuffling(string);

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

	};

	inline PersistentOStream& operator<<(PersistentOStream& os,
	const Process& h) {
	h.persistentOutput(os);
	return os;
	}

	inline PersistentIStream& operator>>(PersistentIStream& is,
	Process& h) {
	h.persistentInput(is);
	return is;
	}

	}

	#endif /* HERWIG_MatchboxAmplitude_H */

File Metadata

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

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions