No OneTemporary
Actions

Size

48 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc b/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc
	--- a/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc
	+++ b/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc
	@@ -1,305 +1,307 @@
	// -- C++ --
	//
	// NJetsAmplitude.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 NJetsAmplitude class.
	//

	#include "NJetsAmplitude.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"

	#include "ThePEG/Utilities/DynamicLoader.h"

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"

	#include "njet.h"

	#include <cstdlib>

	#ifndef NJET_PREFIX
	#error Makefile.am needs to define NJET_PREFIX
	#endif
	#ifndef NJET_LIBS
	#error Makefile.am needs to define NJET_LIBS
	#endif

	using namespace Herwig;

	NJetsAmplitude::NJetsAmplitude() : NJetsPrefix_(NJET_PREFIX),
	NJetsLibs_(NJET_LIBS) {}

	NJetsAmplitude::~NJetsAmplitude() {}

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

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

	void NJetsAmplitude::signOLP(const string& order, const string& contract) {
	string cmd = NJetsPrefix_+"/bin/njet.py -o " + contract + " " + order;
	std::system(cmd.c_str());
	}

	void NJetsAmplitude::startOLP(const string& contract, int& status) {
	NJet::LH_OLP::OLP_Start(contract.c_str(), &status);

	if ( status != 1 )
	return;

	status = 0;

	static double zero = 0.0;
	double param = 0.0;

	param = SM().alphaEMMZ();
	NJet::LH_OLP::OLP_SetParameter("alpha",&param,&zero,&status);
	if ( status != 1 )
	return;

	param = getParticleData(ParticleID::Z0)->hardProcessMass()/GeV;
	NJet::LH_OLP::OLP_SetParameter("mass(23)",&param,&zero,&status);
	if ( status != 1 )
	return;

	param = getParticleData(ParticleID::Wplus)->hardProcessMass()/GeV;
	NJet::LH_OLP::OLP_SetParameter("mass(24)",&param,&zero,&status);
	if ( status != 1 )
	return;

	param = getParticleData(ParticleID::Z0)->hardProcessWidth()/GeV;
	NJet::LH_OLP::OLP_SetParameter("width(23)",&param,&zero,&status);
	if ( status != 1 )
	return;

	param = getParticleData(ParticleID::Wplus)->hardProcessWidth()/GeV;
	NJet::LH_OLP::OLP_SetParameter("width(24)",&param,&zero,&status);
	if ( status != 1 )
	return;

	param = SM().sin2ThetaW();
	NJet::LH_OLP::OLP_SetParameter("sw2",&param,&zero,&status);

	didStartOLP() = true;

	}

	void NJetsAmplitude::loadNJET() {
	if ( ! (DynamicLoader::load(NJetsLibs_+"/libnjet2.so") \|\|
	- DynamicLoader::load("libnjet2.so") ) )
	+ DynamicLoader::load("libnjet2.so") \|\|
	+ DynamicLoader::load(NJetsLibs_+"/libnjet2.dylib") \|\|
	+ DynamicLoader::load("libnjet2.dylib") ) )
	throw Exception() << "NJetsAmplitude: Failed to load libnjet2.so\n"
	<< DynamicLoader::lastErrorMessage
	<< Exception::runerror;
	}

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

	// TODO throw exception on massive leptons in procs

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

	olpOrderFileHeader(orderFile);

	orderFile << "NJetReturnAccuracy yes\n"
	<< "NJetRenormalize yes\n"
	<< "NJetNf " << factory()->nLight() << "\n";

	olpOrderFileProcesses(orderFile,procs);

	orderFile << flush;

	orderFile.close();

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

	signOLP(orderFileName, contractFileName);

	int status = -1;

	startOLP(contractFileName,status);

	if ( status != 1 )
	return false;

	return true;

	}

	LorentzVector<Complex> NJetsAmplitude::plusPolarization(const Lorentz5Momentum& p,
	const Lorentz5Momentum& n,
	int inc) const {

	double pvec[4] = {p.t()/GeV,p.x()/GeV,p.y()/GeV,p.z()/GeV};
	double nvec[4] = {n.t()/GeV,n.x()/GeV,n.y()/GeV,n.z()/GeV};
	double out[8] ={ };
	NJet::LH_OLP::OLP_Polvec(pvec,nvec,out);

	LorentzVector<Complex> res;
	Complex a(out[0],out[1]);
	res.setT(a);
	Complex b(out[2],out[3]);
	res.setX(b);
	Complex c(out[4],out[5]);
	res.setY(c);
	Complex d(out[6],out[7]);
	res.setZ(d);

	if (inc<2)
	return res.conjugate();
	else
	return res;

	}

	void NJetsAmplitude::evalSubProcess() const {

	useMe();

	double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
	fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
	double as;
	if (!hasRunningAlphaS()) as = SM().alphaS();
	else if (hasRunningAlphaS()) as = lastAlphaS();
	double scale = sqrt(mu2()/GeV2);

	double out[7]={};

	int id =
	olpId()[ProcessType::oneLoopInterference] ?
	olpId()[ProcessType::oneLoopInterference] :
	olpId()[ProcessType::treeME2];

	NJet::LH_OLP::OLP_EvalSubProcess(id, olpMomenta(), scale, &as, out);

	if ( olpId()[ProcessType::oneLoopInterference] ) {
	if(calculateTreeME2())lastTreeME2(out[3]*units);
	lastOneLoopInterference(out[2]*units);
	lastOneLoopPoles(pair<double,double>(out[0]units,out[1]units));
	} else if ( olpId()[ProcessType::treeME2] ) {
	lastTreeME2(out[0]*units);
	} else assert(false);

	}

	void NJetsAmplitude::evalColourCorrelator(pair<int,int>) const {

	double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
	fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
	double as;
	if (!hasRunningAlphaS()) as = SM().alphaS();
	else if (hasRunningAlphaS()) as = lastAlphaS();
	double scale = sqrt(mu2()/GeV2);

	int n = lastXComb().meMomenta().size();
	colourCorrelatorResults.resize(n*(n-1)/2);

	NJet::LH_OLP::OLP_EvalSubProcess(olpId()[ProcessType::colourCorrelatedME2],
	olpMomenta(), scale, &as, &colourCorrelatorResults[0]);

	for ( int i = 0; i < n; ++i )
	for ( int j = i+1; j < n; ++j ) {
	lastColourCorrelator(make_pair(i,j),colourCorrelatorResults[i+j(j-1)/2]units);
	}

	}

	void NJetsAmplitude::evalSpinColourCorrelator(pair<int,int>) const {

	double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
	fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
	double as;
	if (!hasRunningAlphaS()) as = SM().alphaS();
	else if (hasRunningAlphaS()) as = lastAlphaS();
	double scale = sqrt(mu2()/GeV2);

	int n = lastXComb().meMomenta().size();
	spinColourCorrelatorResults.resize(2nn);

	NJet::LH_OLP::OLP_EvalSubProcess(olpId()[ProcessType::spinColourCorrelatedME2],
	olpMomenta(), scale, &as, &spinColourCorrelatorResults[0]);

	for ( int i = 0; i < n; ++i )
	for ( int j = 0; j < n; ++j ) {
	if ( i == j \|\| mePartonData()[i]->id() != 21 )
	continue;
	Complex scc(spinColourCorrelatorResults[2i+2nj]units,
	spinColourCorrelatorResults[2i+2nj+1]units);
	lastColourSpinCorrelator(make_pair(i,j),scc);
	}

	}

	void NJetsAmplitude::doinit() {
	loadNJET();
	MatchboxOLPME::doinit();
	}

	void NJetsAmplitude::doinitrun() {
	loadNJET();
	MatchboxOLPME::doinitrun();
	}

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


	void NJetsAmplitude::persistentOutput(PersistentOStream & os) const {
	os << colourCorrelatorResults << spinColourCorrelatorResults
	<< NJetsPrefix_ << NJetsLibs_;
	}

	void NJetsAmplitude::persistentInput(PersistentIStream & is, int) {
	is >> colourCorrelatorResults >> spinColourCorrelatorResults
	>> NJetsPrefix_ >> NJetsLibs_;
	}


	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<NJetsAmplitude,MatchboxOLPME>
	describeHerwigNJetsAmplitude("Herwig::NJetsAmplitude", "HwMatchboxNJet.so");

	void NJetsAmplitude::Init() {

	static ClassDocumentation<NJetsAmplitude> documentation
	("NJetsAmplitude implements an interface to NJets.",
	"Matrix elements have been calculated using NJet \\cite{Badger:2012pg}",
	"%\\cite{Badger:2012pg}\n"
	"\\bibitem{Badger:2012pg}\n"
	"S.~Badger et al.,\n"
	"``Numerical evaluation of virtual corrections to multi-jet production in massless QCD,''\n"
	"arXiv:1209.0100 [hep-ph].\n"
	"%%CITATION = ARXIV:1209.0100;%%");

	static Parameter<NJetsAmplitude,string> interfaceNJetsPrefix
	("NJetsPrefix",
	"The prefix for the location of NJets",
	&NJetsAmplitude::NJetsPrefix_, string(NJET_PREFIX),
	false, false);

	static Parameter<NJetsAmplitude,string> interfaceNJetsLibs
	("NJetsLibs",
	"The location of the NJets library",
	&NJetsAmplitude::NJetsLibs_, string(NJET_LIBS),
	false, false);

	}

	diff --git a/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc b/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc
	--- a/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc
	+++ b/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc
	@@ -1,461 +1,464 @@
	// -- C++ --
	//
	// OpenLoopsAmplitude.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 OpenLoopsAmplitude class.
	//

	#include "OpenLoopsAmplitude.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "ThePEG/Utilities/DynamicLoader.h"

	#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"

	#include <fstream>
	#include <sstream>
	#include <string>
	#include <cstdlib>

	using namespace Herwig;

	#ifndef OPENLOOPSLIBS
	#error Makefile.am needs to define OPENLOOPSLIBS
	#endif

	#ifndef OPENLOOPSPREFIX
	#error Makefile.am needs to define OPENLOOPSPREFIX
	#endif

	OpenLoopsAmplitude::OpenLoopsAmplitude() :
	use_cms(true),psp_tolerance(12),
	OpenLoopsLibs_(OPENLOOPSLIBS), OpenLoopsPrefix_(OPENLOOPSPREFIX) {
	}

	OpenLoopsAmplitude::~OpenLoopsAmplitude() {
	}

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

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

	extern "C" void OLP_Start(const char, int i);
	extern "C" void OLP_SetParameter(const char* ,double* ,double,int);
	extern "C" void ol_setparameter_string(const char, const char);
	extern "C" void OLP_PrintParameter(const char*);
	extern "C" void OLP_EvalSubProcess(int, double, double, double, double*);
	extern "C" void OLP_EvalSubProcess2(int, double, double, double, double*);
	// id ps-point emitter polvec res
	extern "C" void ol_evaluate_sc(int, double, int, double, double*);
	extern "C" void OLP_Polvec(double,double,double*);

	void OpenLoopsAmplitude::doinitrun() {
	MatchboxOLPME::doinitrun();
	}

	void OpenLoopsAmplitude::startOLP(const string& contract, int& status) {
	string tempcontract=contract;

	if ( ! (DynamicLoader::load(OpenLoopsLibs_+"/libopenloops.so") \|\|
	DynamicLoader::load(OpenLoopsPrefix_+"/lib/libopenloops.so") \|\|
	- DynamicLoader::load("libopenloops.so") ) ) {
	+ DynamicLoader::load("libopenloops.so") \|\|
	+ DynamicLoader::load(OpenLoopsLibs_+"/libopenloops.dylib") \|\|
	+ DynamicLoader::load(OpenLoopsPrefix_+"/lib/libopenloops.dylib") \|\|
	+ DynamicLoader::load("libopenloops.dylib") ) ) {
	throw Exception() << "OpenLoopsAmplitude::startOLP(): Failed to load libopenloops.so/dylib\n"
	<< DynamicLoader::lastErrorMessage
	<< Exception::runerror;
	}

	string stabilityPrefix = factory()->runStorage() + "OpenLoops.StabilityLog";
	assert(stabilityPrefix.size() < 256);

	ol_setparameter_string("stability_logdir",stabilityPrefix.c_str());

	ol_setparameter_string("install_path",OpenLoopsPrefix_.c_str());

	int a=0;double null=0.0;double one=1.0;
	int part[10]={1,2,3,4,5,6,15,23,24,25};string stri;
	for (int i=0;i<10;i++){
	map<long,Energy>::const_iterator it=reshuffleMasses().find(part[i]);
	double mass;
	if(it==reshuffleMasses().end())
	mass = getParticleData(part[i])->hardProcessMass()/GeV;
	else
	mass = it->second/GeV;
	double width=getParticleData(part[i])->hardProcessWidth()/GeV;
	std::stringstream ss;
	ss << part[i];
	string str = ss.str();
	stri="mass("+str+")";
	OLP_SetParameter(stri.c_str(),&mass,&null,&a);
	stri="width("+str+")";
	OLP_SetParameter(stri.c_str(),&width,&null,&a);
	}
	stri="alphas";
	one=SM().alphaS();
	OLP_SetParameter( stri.c_str(),&one ,&null,&a);
	stri="alpha";
	one=SM().alphaEMMZ();
	OLP_SetParameter(stri.c_str(),&one ,&null,&a);

	OLP_Start(tempcontract.c_str(), &status);

	didStartOLP() = true;

	}

	void OpenLoopsAmplitude::fillOrderFile(const map<pair<Process, int>, int>& procs) {
	string orderFileName = factory()->buildStorage() + name() + ".OLPContract.lh";
	ofstream orderFile(orderFileName.c_str());
	size_t asPower = 100;
	size_t minlegs = 100;
	size_t maxlegs = 0;

	for ( map<pair<Process, int>, int>::const_iterator t = procs.begin() ; t != procs.end() ; ++t ) {
	asPower = min(asPower, static_cast<size_t>(t->first.first.orderInAlphaS));
	minlegs = min(minlegs, static_cast<size_t>(t->first.first.legs.size()));
	maxlegs = max(maxlegs, static_cast<size_t>(t->first.first.legs.size()));
	}

	orderFile << "# OLP order file created by Herwig/Matchbox for OpenLoops\n\n";
	orderFile << "CorrectionType QCD\n";
	orderFile << "IRregularization " << (isDR() ? "DRED" : "CDR") << "\n";
	orderFile << "extra answerfile " << (factory()->buildStorage() + name() + ".OLPAnswer.lh") << "\n";
	orderFile << "extra psp_tolerance "<<psp_tolerance<<"\n";
	orderFile << "extra use_cms "<<(use_cms?"1":"0")<< "\n";
	orderFile << "\n";

	if (extraOpenLoopsPath!="")
	orderFile << "Extra OpenLoopsPath " << extraOpenLoopsPath << "\n";
	for ( map<pair<Process, int>, int>::const_iterator p = procs.begin() ; p != procs.end() ; ++p ) {
	std::stringstream Processstr;
	std::stringstream Typestr;
	Processstr << (p).first.first.legs[0]->id() << " " << (p).first.first.legs[1]->id() << " -> ";
	for ( PDVector::const_iterator o = (p).first.first.legs.begin() + 2 ; o != (p).first.first.legs.end() ; ++o )
	Processstr << (**o).id() << " ";
	if ( (*p).first.second == ProcessType::treeME2 ) {
	Typestr << "Tree";
	} else if ( (*p).first.second == ProcessType::colourCorrelatedME2 ) {
	Typestr << "ccTree";
	} else if ( (*p).first.second == ProcessType::spinColourCorrelatedME2 ) {
	Typestr << "sctree_polvect";
	} else if ( (*p).first.second == ProcessType::oneLoopInterference ) {
	Typestr << "Loop";
	}
	OpenLoopsProcInfo pro = OpenLoopsProcInfo((*p).second, -1, Processstr.str(), Typestr.str());
	pro.setOAs(p->first.first.orderInAlphaS);
	processmap[(*p).second] = pro;
	}

	vector < string > types;
	types.push_back("Tree");
	types.push_back("ccTree");
	types.push_back("sctree_polvect");
	types.push_back("Loop");
	for ( size_t i = asPower ; i != asPower + maxlegs - minlegs + 1 ; i++ ) {
	orderFile << "\n\nCouplingPower QCD " << i;
	orderFile << "\n\n#AlphasPower " << i;
	for ( vector<string>::iterator it = types.begin() ; it != types.end() ; it++ ) {
	for ( map<int, OpenLoopsProcInfo>::iterator p = processmap.begin() ; p != processmap.end() ; ++p )
	if ( (p).second.Tstr() == it && i == (unsigned int) (*p).second.orderAs() ) {
	orderFile << "\nAmplitudeType " << *it << "\n";
	break;
	}
	for ( map<int, OpenLoopsProcInfo>::iterator p = processmap.begin() ; p != processmap.end() ; ++p )
	if ( (p).second.Tstr() == it && i == (unsigned int) (*p).second.orderAs() ) {
	orderFile << (*p).second.Pstr() << "\n";
	}
	}
	}
	orderFile << flush;
	}

	bool OpenLoopsAmplitude::checkOLPContract() {
	string contractFileName = factory()->buildStorage() + name() + ".OLPAnswer.lh";
	ifstream infile(contractFileName.c_str());
	string line;
	vector < string > contractfile;
	while (std::getline(infile, line)) {
	contractfile.push_back(line);
	}
	for ( map<int, OpenLoopsProcInfo>::iterator p = processmap.begin() ; p != processmap.end() ; p++ ) {
	bool righttype = false;
	for ( vector<string>::iterator linex = contractfile.begin() ; linex != contractfile.end() ; ++linex ) {
	if ( (*linex).find("AmplitudeType ")!= std::string::npos ) {
	if ( (linex).find(" " + (p).second.Tstr() + " ")!= std::string::npos ) {
	righttype = true;
	} else {
	righttype = false;
	}
	}
	if ( righttype ) {
	if ( (linex).find((p).second.Pstr()) != std::string::npos ){
	if( (p).second.Pstr().length() == (linex).find("\|") ) {
	string sub = (linex).substr((linex).find("\|") + 1, (linex).find("#") - (linex).find("\|") - 1); // \| 1 23 # buggy??

	int subint;
	int subint2;
	istringstream(sub) >> subint >> subint2;
	assert(subint==1);
	(*p).second.setGID(subint2);
	}
	}
	}
	}
	}
	string ids = factory()->buildStorage() + "OpenLoops.ids.dat";
	ofstream IDS(ids.c_str());

	for ( map<int, OpenLoopsProcInfo>::iterator p = processmap.begin() ; p != processmap.end() ; p++ ) {
	idpair.insert ( std::pair<int,int>((p).second.HID(),(p).second.GID()) );
	IDS << (p).second.HID() << " " << (p).second.GID() << "\n";
	if ( (*p).second.GID() == -1 ) return 0;
	}
	IDS << flush;
	return 1;
	}

	void OpenLoopsAmplitude::getids() const{
	string line = factory()->buildStorage() + "OpenLoops.ids.dat";
	ifstream infile(line.c_str());
	int hid;
	int gid;
	while (std::getline(infile, line)) {
	istringstream(line) >> hid>>gid;
	idpair.insert ( std::pair<int,int>(hid,gid) );
	}
	}
	bool OpenLoopsAmplitude::startOLP(const map<pair<Process, int>, int>& procs) {
	string contractFileName = factory()->buildStorage() + name() + ".OLPAnswer.lh";
	string orderFileName = factory()->buildStorage() + name() + ".OLPContract.lh";
	fillOrderFile(procs);
	int status = -1;
	startOLP(orderFileName, status);
	if ( !checkOLPContract() ) {
	return false;
	}
	if ( status != 1 ) return false;
	return true;
	}


	void OpenLoopsAmplitude::evalSubProcess() const {
	useMe();
	double units = pow(lastSHat() / GeV2, mePartonData().size() - 4.);
	fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
	double acc ;
	double scale = sqrt(mu2() / GeV2);
	if (hasRunningAlphaS()) {
	int a=0;double null=0.0;double one=1.0;
	string stri="alphas";
	one=lastAlphaS();
	OLP_SetParameter( stri.c_str(),&one ,&null,&a);
	}
	double out[7]={};

	int id = olpId()[ProcessType::oneLoopInterference] ? olpId()[ProcessType::oneLoopInterference] : olpId()[ProcessType::treeME2];

	if ( idpair.size() == 0 ) {
	getids();
	if ( Debug::level > 1 ) {
	string parfile=factory()->runStorage() + name() + ".Parameters.dat";
	OLP_PrintParameter(parfile.c_str());
	}

	}


	OLP_EvalSubProcess2(&((*(idpair.find(id))).second), olpMomenta(), &scale, out,&acc );


	if ( olpId()[ProcessType::oneLoopInterference] ) {
	if(calculateTreeME2())lastTreeME2(out[3] * units);
	lastOneLoopInterference((out[2])* units);
	lastOneLoopPoles(pair<double, double>(out[0] * units, out[1] * units));
	} else if ( olpId()[ProcessType::treeME2] ) {
	lastTreeME2(out[0] * units);
	}
	}

	void OpenLoopsAmplitude::evalColourCorrelator(pair<int, int> ) const {
	double units = pow(lastSHat() / GeV2, mePartonData().size() - 4.);
	fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
	double acc ;
	double scale = sqrt(mu2() / GeV2);
	if (hasRunningAlphaS()) {
	int a=0;double null=0.0;double one=1.0;
	string stri="alphas";
	one=lastAlphaS();
	OLP_SetParameter( stri.c_str(),&one ,&null,&a);
	}
	int n = lastXComb().meMomenta().size();

	colourCorrelatorResults.resize(n * (n - 1) / 2);
	if ( idpair.size() == 0 ) {
	getids();
	if ( Debug::level > 1 ) {
	string parfile=factory()->runStorage() + name() + ".Parameters.dat";
	OLP_PrintParameter(parfile.c_str());
	}

	}


	int id = olpId()[ProcessType::colourCorrelatedME2];

	OLP_EvalSubProcess2(&((*(idpair.find(id))).second), olpMomenta(), &scale, &colourCorrelatorResults[0],&acc );
	for ( int i = 0 ; i < n ; ++i ){
	for ( int j = i + 1 ; j < n ; ++j ) {
	lastColourCorrelator(make_pair(i, j), colourCorrelatorResults[i+j(j-1)/2] units);
	}
	}
	}

	void OpenLoopsAmplitude::evalSpinColourCorrelator(pair<int , int > ) const {
	assert(false);
	}


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

	double units = pow(lastSHat() / GeV2, mePartonData().size() - 4.);
	fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());

	if (hasRunningAlphaS()) {
	int a=0;double null=0.0;double one=1.0;
	string stri="alphas";
	one=lastAlphaS();
	OLP_SetParameter( stri.c_str(),&one ,&null,&a);
	}

	int emitter=ij.first+1;
	int n = lastXComb().meMomenta().size();

	if ( idpair.size() == 0 ) {
	getids();
	if ( Debug::level > 1 ) {
	string parfile=factory()->runStorage() + name() + ".Parameters.dat";
	OLP_PrintParameter(parfile.c_str());
	}

	}
	int id = (*(idpair.find(olpId()[ProcessType::spinColourCorrelatedME2]))).second;
	//double * outx =new double[n];
	spinColourCorrelatorResults.resize(n);
	double polvec[4];

	polvec[0]=c.momentum().e()/GeV;
	polvec[1]=c.momentum().x()/GeV;
	polvec[2]=c.momentum().y()/GeV;
	polvec[3]=c.momentum().z()/GeV;
	double avg= colourCorrelatedME2(ij)*(-c.diagonal());

	ol_evaluate_sc(id, olpMomenta(),emitter,polvec,&spinColourCorrelatorResults[0]);

	double corr =-1.units spinColourCorrelatorResults[ij.second]/c.scale()*c.momentum().dot(c.momentum());


	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 + corr/cfac;

	}





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

	void OpenLoopsAmplitude::persistentOutput(PersistentOStream & os) const {
	os << idpair << OpenLoopsLibs_ << OpenLoopsPrefix_;
	}

	void OpenLoopsAmplitude::persistentInput(PersistentIStream & is, int) {
	is >> idpair >> OpenLoopsLibs_ >> OpenLoopsPrefix_;
	}

	// * 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<OpenLoopsAmplitude, MatchboxOLPME> describeHerwigOpenLoopsAmplitude("Herwig::OpenLoopsAmplitude", "HwMatchboxOpenLoops.so");

	void OpenLoopsAmplitude::Init() {

	static ClassDocumentation<OpenLoopsAmplitude>
	documentation("OpenLoopsAmplitude implements an interface to OpenLoops.",
	"Matrix elements have been calculated using OpenLoops \\cite{Cascioli:2011va}",
	"%\\cite{Cascioli:2011va}\n"
	"\\bibitem{Cascioli:2011va}\n"
	"F.~Cascioli et al.,\n"
	"``Scattering Amplitudes with Open Loops,''\n"
	"arXiv:1111.5206 [hep-ph].\n"
	"%%CITATION = ARXIV:1111.5206;%%");

	static Switch<OpenLoopsAmplitude,bool> interfaceUseComplMass
	("ComplexMassScheme",
	"Switch on or off if Compex Masses.",
	&OpenLoopsAmplitude::use_cms, true, false, false);
	static SwitchOption interfaceUseComplMassOn
	(interfaceUseComplMass,
	"True",
	"True for Complex Masses.",
	true);
	static SwitchOption interfaceUseComplMassOff
	(interfaceUseComplMass,
	"False",
	"False for no Complex Masses.",
	false);

	static Parameter<OpenLoopsAmplitude,int> interfacepsp_tolerance
	("PSP_tolerance",
	"(Debug)Phase Space Tolerance. Better use e.g.: set OpenLoops:Massless 13",
	&OpenLoopsAmplitude::psp_tolerance, 12, 0, 0,
	false, false, Interface::lowerlim);

	static Parameter<OpenLoopsAmplitude,string> interfaceOpenLoopsLibs
	("OpenLoopsLibs",
	"The location of OpenLoops libraries",
	&OpenLoopsAmplitude::OpenLoopsLibs_, string(OPENLOOPSLIBS),
	false, false);

	static Parameter<OpenLoopsAmplitude,string> interfaceOpenLoopsPrefix
	("OpenLoopsPrefix",
	"The location of OpenLoops libraries",
	&OpenLoopsAmplitude::OpenLoopsPrefix_, string(OPENLOOPSPREFIX),
	false, false);

	}

	diff --git a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOAmplitude.cc b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOAmplitude.cc
	--- a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOAmplitude.cc
	+++ b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOAmplitude.cc
	@@ -1,457 +1,459 @@
	// -- C++ --
	//
	// VBFNLOAmplitude.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 VBFNLOAmplitude class.
	//

	#include "VBFNLOAmplitude.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/Switch.h"

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "ThePEG/Utilities/DynamicLoader.h"

	#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"

	#include <cstdlib>

	#include "VBFNLO/utilities/BLHAinterface.h"

	#define DEFSTR(s) CPPSTR(s)
	#define CPPSTR(s) #s

	using namespace Herwig;

	VBFNLOAmplitude::VBFNLOAmplitude()
	: theRanHelSum(false), theAnomCoupl(false), VBFNLOlib_(DEFSTR(VBFNLOLIB))
	{}

	VBFNLOAmplitude::~VBFNLOAmplitude() {}

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

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

	void VBFNLOAmplitude::signOLP(const string& order, const string& contract) {
	int status = 0;
	OLP_Order(const_cast<char*>(order.c_str()),
	const_cast<char*>(contract.c_str()),&status);
	if ( status != 1 )
	throw Exception() << "VBFNLOAmplitude: Failed to sign contract with VBFNLO"
	<< Exception::runerror;
	}

	void VBFNLOAmplitude::setOLPParameter(const string& name, double value) const {

	int pStatus = 0;
	double zero = 0.0;
	OLP_SetParameter(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
	if ( !pStatus )
	throw Exception() << "VBFNLOAmplitude: VBFNLO failed to set parameter '"
	<< name << "' to " << value << "\n"
	<< Exception::runerror;

	}

	void VBFNLOAmplitude::startOLP(const string& contract, int& status) {

	OLP_Start(const_cast<char*>(contract.c_str()), &status);

	map<long,Energy>::const_iterator it=reshuffleMasses().find(ParticleID::b);
	double bmass;
	if(it==reshuffleMasses().end())
	bmass = getParticleData(ParticleID::b)->hardProcessMass()/GeV;
	else
	bmass = it->second/GeV;
	setOLPParameter("mass(5)",bmass);
	setOLPParameter("mass(6)",getParticleData(ParticleID::t)->hardProcessMass()/GeV);

	setOLPParameter("mass(23)",getParticleData(ParticleID::Z0)->hardProcessMass()/GeV);
	setOLPParameter("mass(24)",getParticleData(ParticleID::Wplus)->hardProcessMass()/GeV);
	setOLPParameter("mass(25)",getParticleData(ParticleID::h0)->hardProcessMass()/GeV);

	setOLPParameter("width(23)",getParticleData(ParticleID::Z0)->hardProcessWidth()/GeV);
	setOLPParameter("width(24)",getParticleData(ParticleID::Wplus)->hardProcessWidth()/GeV);
	setOLPParameter("width(25)",getParticleData(ParticleID::h0)->hardProcessWidth()/GeV);

	setOLPParameter("alpha",SM().alphaEMMZ());
	setOLPParameter("sw2",SM().sin2ThetaW());
	setOLPParameter("Gf",SM().fermiConstant()*GeV2);

	setOLPParameter("Nf",factory()->nLight());

	setOLPParameter("alphas",SM().alphaS());

	setOLPParameter("ranhelsum",theRanHelSum);

	setOLPParameter("anomcoupl",theAnomCoupl);

	didStartOLP() = true;

	}

	void VBFNLOAmplitude::loadVBFNLO() {
	if ( ! (DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.so") \|\|
	- DynamicLoader::load("libVBFNLO.so") ) )
	+ DynamicLoader::load("libVBFNLO.so") \|\|
	+ DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.dylib") \|\|
	+ DynamicLoader::load("libVBFNLO.dylib") ) )
	throw Exception() << "VBFNLOAmplitude: failed to load libVBFNLO.so/dylib\n"
	<< DynamicLoader::lastErrorMessage
	<< Exception::runerror;
	}

	bool VBFNLOAmplitude::startOLP(const map<pair<Process,int>,int>& procs) {
	loadVBFNLO();
	string orderFileName = factory()->buildStorage() + name() + ".OLPOrder.lh";
	ofstream orderFile(orderFileName.c_str());

	olpOrderFileHeader(orderFile);

	// add VBFNLO specifics here

	olpOrderFileProcesses(orderFile,procs);

	orderFile << flush;

	orderFile.close();

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

	signOLP(orderFileName, contractFileName);

	int status = -1;

	startOLP(contractFileName,status);

	if ( status != 1 )
	return false;

	return true;


	}

	LorentzVector<Complex> VBFNLOAmplitude::plusPolarization(const Lorentz5Momentum& p,
	const Lorentz5Momentum& n,
	int inc) const {

	// shamelessly stolen from the GoSam interface; mind that we can
	// always cast eq (5.7) in the manual into a form that it only uses
	// <M-\|\|M_+> and then switch bvetween eps_+ for an outgoing and
	// eps_- for an incoming gluon.

	double pvec[4] = {p.t()/GeV,p.x()/GeV,p.y()/GeV,p.z()/GeV};
	double nvec[4] = {n.t()/GeV,n.x()/GeV,n.y()/GeV,n.z()/GeV};
	double out[8] ={ };
	OLP_Polvec(pvec,nvec,out);

	LorentzVector<Complex> res;
	Complex a(out[0],out[1]);
	res.setT(a);
	Complex b(out[2],out[3]);
	res.setX(b);
	Complex c(out[4],out[5]);
	res.setY(c);
	Complex d(out[6],out[7]);
	res.setZ(d);

	if (inc<2)
	return res.conjugate();
	else
	return res;

	}

	void VBFNLOAmplitude::evalSubProcess() const {

	useMe();

	double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
	fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
	double scale = sqrt(mu2()/GeV2);

	if (hasRunningAlphaS()) setOLPParameter("alphas",lastAlphaS());

	double acc = -1.0;
	double out[4]={};

	int id =
	olpId()[ProcessType::oneLoopInterference] ?
	olpId()[ProcessType::oneLoopInterference] :
	olpId()[ProcessType::treeME2];

	if (theRanHelSum) {
	vector<double> helicityrn = amplitudeRandomNumbers();
	if (helicityrn.size()>0) {
	setOLPParameter("HelicityRN",helicityrn[0]);
	}
	}

	OLP_EvalSubProcess2(&id, olpMomenta(), &scale, out, &acc);

	if ( olpId()[ProcessType::oneLoopInterference] ) {
	lastTreeME2(out[3]*units);
	lastOneLoopInterference(out[2]*units);
	lastOneLoopPoles(pair<double,double>(out[0]units,out[1]units));
	} else if ( olpId()[ProcessType::treeME2] ) {
	lastTreeME2(out[0]*units);
	} else assert(false);

	}

	void VBFNLOAmplitude::evalColourCorrelator(pair<int,int>) const {

	double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
	fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
	double scale = sqrt(mu2()/GeV2);

	if (hasRunningAlphaS()) setOLPParameter("alphas",lastAlphaS());

	double acc = -1.0;

	int n = lastXComb().meMomenta().size();
	colourCorrelatorResults.resize(n*(n-1)/2);

	int id = olpId()[ProcessType::colourCorrelatedME2];

	if ( theRanHelSum ) {
	if ( lastHeadMatchboxXComb() ) {
	vector<double> helicityrn = lastHeadMatchboxXComb()->amplitudeRandomNumbers();
	if (helicityrn.size()>0) {
	setOLPParameter("HelicityRN",helicityrn[0]);
	}
	} else if ( amplitudeRandomNumbers().size() > 0 ) {
	vector<double> helicityrn = amplitudeRandomNumbers();
	if (helicityrn.size()>0) {
	setOLPParameter("HelicityRN",helicityrn[0]);
	}
	}
	}

	OLP_EvalSubProcess2(&id, olpMomenta(), &scale, &colourCorrelatorResults[0], &acc);

	for ( int i = 0; i < n; ++i )
	for ( int j = i+1; j < n; ++j ) {
	lastColourCorrelator(make_pair(i,j),colourCorrelatorResults[i+j(j-1)/2]units);
	}

	}

	void VBFNLOAmplitude::evalSpinColourCorrelator(pair<int,int>) const {

	double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
	fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
	double scale = sqrt(mu2()/GeV2);

	if (hasRunningAlphaS()) setOLPParameter("alphas",lastAlphaS());

	double acc = -1.0;

	int n = lastXComb().meMomenta().size();
	spinColourCorrelatorResults.resize(2nn);

	int id = olpId()[ProcessType::spinColourCorrelatedME2];

	if (theRanHelSum && lastHeadMatchboxXComb()) {
	vector<double> helicityrn = lastHeadMatchboxXComb()->amplitudeRandomNumbers();
	if (helicityrn.size()>0) {
	setOLPParameter("HelicityRN",helicityrn[0]);
	}
	}

	OLP_EvalSubProcess2(&id, olpMomenta(), &scale, &spinColourCorrelatorResults[0], &acc);

	for ( int i = 0; i < n; ++i )
	for ( int j = 0; j < n; ++j ) {
	if ( i == j \|\| mePartonData()[i]->id() != 21 )
	continue;
	Complex scc(spinColourCorrelatorResults[2i+2nj]units,
	spinColourCorrelatorResults[2i+2nj+1]units);
	lastColourSpinCorrelator(make_pair(i,j),scc);
	}

	}

	double VBFNLOAmplitude::largeNME2(Ptr<ColourBasis>::tptr cptr) const {
	if ( calculateLargeNME2() )
	evalLargeNSubProcess(cptr);
	return lastLargeNME2();
	}

	void VBFNLOAmplitude::evalLargeNSubProcess(Ptr<ColourBasis>::tptr) const {

	double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
	fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
	double scale = sqrt(mu2()/GeV2);

	if (hasRunningAlphaS()) setOLPParameter("alphas",lastAlphaS());

	double acc = -1.0;
	double out[4]={};

	int id =
	olpId()[ProcessType::oneLoopInterference] ?
	olpId()[ProcessType::oneLoopInterference] :
	olpId()[ProcessType::treeME2];

	if (theRanHelSum) {
	vector<double> helicityrn = amplitudeRandomNumbers();
	if (helicityrn.size()>0) {
	setOLPParameter("HelicityRN",helicityrn[0]);
	}
	}

	setOLPParameter("Nc",-1); // large-N limit
	OLP_EvalSubProcess2(&id, olpMomenta(), &scale, out, &acc);
	setOLPParameter("Nc",generator()->standardModel()->Nc());

	if ( olpId()[ProcessType::oneLoopInterference] ) {
	lastLargeNME2(out[3]*units);
	lastOneLoopInterference(out[2]*units);
	lastOneLoopPoles(pair<double,double>(out[0]units,out[1]units));
	} else if ( olpId()[ProcessType::treeME2] ) {
	lastLargeNME2(out[0]*units);
	} else assert(false);

	}

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

	double cfac = 1.;
	double Nc = generator()->standardModel()->Nc();
	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);
	if ( calculateLargeNColourCorrelator(ij) )
	evalLargeNColourCorrelator(ij,cptr);
	return lastLargeNColourCorrelator(ij)/cfac;

	}

	void VBFNLOAmplitude::evalLargeNColourCorrelator(pair<int,int>,
	Ptr<ColourBasis>::tptr) const {

	double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
	fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
	double scale = sqrt(mu2()/GeV2);

	if (hasRunningAlphaS()) setOLPParameter("alphas",lastAlphaS());

	double acc = -1.0;

	int n = lastXComb().meMomenta().size();
	colourCorrelatorResults.resize(n*(n-1)/2);

	int id = olpId()[ProcessType::colourCorrelatedME2];

	if (theRanHelSum && lastHeadMatchboxXComb()) {
	vector<double> helicityrn = lastHeadMatchboxXComb()->amplitudeRandomNumbers();
	if (helicityrn.size()>0) {
	setOLPParameter("HelicityRN",helicityrn[0]);
	}
	}

	setOLPParameter("Nc",-1); // large-N limit
	OLP_EvalSubProcess2(&id, olpMomenta(), &scale, &colourCorrelatorResults[0], &acc);
	setOLPParameter("Nc",generator()->standardModel()->Nc());

	for ( int i = 0; i < n; ++i )
	for ( int j = i+1; j < n; ++j ) {
	lastLargeNColourCorrelator(make_pair(i,j),colourCorrelatorResults[i+j(j-1)/2]units);
	}

	}

	void VBFNLOAmplitude::doinit() {
	loadVBFNLO();
	MatchboxOLPME::doinit();
	}

	void VBFNLOAmplitude::doinitrun() {
	loadVBFNLO();
	MatchboxOLPME::doinitrun();
	}

	void VBFNLOAmplitude::persistentOutput(PersistentOStream & os) const {
	os << colourCorrelatorResults << spinColourCorrelatorResults << theRanHelSum << theAnomCoupl << VBFNLOlib_;
	}

	void VBFNLOAmplitude::persistentInput(PersistentIStream & is, int) {
	is >> colourCorrelatorResults >> spinColourCorrelatorResults >> theRanHelSum >> theAnomCoupl >> VBFNLOlib_;
	}


	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<VBFNLOAmplitude,MatchboxOLPME>
	describeHerwigVBFNLOAmplitude("Herwig::VBFNLOAmplitude", "HwMatchboxVBFNLO.so");

	void VBFNLOAmplitude::Init() {

	static ClassDocumentation<VBFNLOAmplitude> documentation
	("VBFNLOAmplitude implements an interface to VBFNLO.",
	"Matrix elements have been calculated using VBFNLO "
	"(Ref.~\\cite{VBFNLO} and process-specific references)\n",
	"%\\cite{VBFNLO}\n"
	"\\bibitem{Arnold:2008rz}\n"
	"K.~Arnold, M.~Bahr, G.~Bozzi, F.~Campanario, C.~Englert, T.~Figy, "
	"N.~Greiner and C.~Hackstein {\\it et al.},\n"
	"``VBFNLO: A Parton level Monte Carlo for processes with electroweak bosons,''\n"
	"Comput.\\ Phys.\\ Commun.\\ {\\bf 180} (2009) 1661\n"
	"[arXiv:0811.4559 [hep-ph]];\n"
	"%%CITATION = ARXIV:0811.4559;%%\n"
	"J.~Baglio, J.~Bellm, F.~Campanario, B.~Feigl, J.~Frank, T.~Figy, "
	"M.~Kerner and L.~D.~Ninh {\\it et al.},\n"
	"``Release Note - VBFNLO 2.7.0,''\n"
	"arXiv:1404.3940 [hep-ph].\n"
	"%%CITATION = ARXIV:1404.3940;%%\n");

	static Switch<VBFNLOAmplitude,bool> interfaceRandomHelicitySummation
	("RandomHelicitySummation", "Switch for random helicity summation of leptons and photons",
	&VBFNLOAmplitude::theRanHelSum, false, false, false);
	static SwitchOption interfaceRandomHelicitySummationTrue
	(interfaceRandomHelicitySummation,
	"True",
	"Perform random helicity summation",
	true);
	static SwitchOption interfaceRandomHelicitySummationFalse
	(interfaceRandomHelicitySummation,
	"False",
	"Sum over all helicity combinations",
	false);

	static Switch<VBFNLOAmplitude,bool> interfaceAnomalousCouplings
	("AnomalousCouplings", "Switch for anomalous couplings",
	&VBFNLOAmplitude::theAnomCoupl, false, false, false);
	static SwitchOption interfaceAnomalousCouplingsTrue
	(interfaceAnomalousCouplings,
	"On",
	"Switch anomalous couplings on",
	true);
	static SwitchOption interfaceAnomalousCouplingsFalse
	(interfaceAnomalousCouplings,
	"Off",
	"Switch anomalous couplings off",
	false);

	}

	diff --git a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
	--- a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
	+++ b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
	@@ -1,259 +1,261 @@
	// -- C++ --
	//
	// VBFNLOPhasespace.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 VBFNLOPhasespace class.
	//

	#include "VBFNLOPhasespace.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Reference.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 "Herwig/Utilities/GSLBisection.h"
	#include "ThePEG/Utilities/DynamicLoader.h"

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"

	#include "VBFNLO/utilities/BLHAinterface.h"

	#define DEFSTR(s) CPPSTR(s)
	#define CPPSTR(s) #s

	using namespace Herwig;

	VBFNLOPhasespace::VBFNLOPhasespace() :
	lastSqrtS(0*GeV), needToReshuffle(false), VBFNLOlib_(DEFSTR(VBFNLOLIB))
	{}

	void VBFNLOPhasespace::loadVBFNLO() {
	if ( ! (DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.so") \|\|
	- DynamicLoader::load("libVBFNLO.so") ) )
	+ DynamicLoader::load("libVBFNLO.so") \|\|
	+ DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.dylib") \|\|
	+ DynamicLoader::load("libVBFNLO.dylib") ) )
	throw Exception() << "VBFNLOPhasespace::loadVBFNLO(): Failed to load libVBFNLO.so/dylib\n"
	<< DynamicLoader::lastErrorMessage
	<< Exception::runerror;
	}

	VBFNLOPhasespace::~VBFNLOPhasespace() {}

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

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

	void VBFNLOPhasespace::setXComb(tStdXCombPtr xco) {

	MatchboxPhasespace::setXComb(xco);
	// test for resuffling
	needToReshuffle = false;
	if ( xco ) {
	for ( cPDVector::const_iterator d = mePartonData().begin();
	d != mePartonData().end(); ++d ) {
	// Higgs is massive -> does not need reshuffling
	if ( ( (d).id() != ParticleID::h0 ) && ( (d).hardProcessMass() != ZERO ) ) {
	needToReshuffle = true;
	break;
	}
	}
	}

	// set CMS energy
	int pStatus = 0;
	double zero = 0.0;

	double value = sqrt(lastXCombPtr()->lastS())/GeV;
	if (value && (value != lastSqrtS/GeV)) {
	lastSqrtS = value*GeV;
	string name = "sqrtS";
	OLP_SetParameter(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
	if ( !pStatus )
	throw Exception() << "VBFNLOPhasespace::setXComb(): VBFNLO failed to set parameter '"
	<< name << "' to " << value << "\n"
	<< Exception::runerror;
	}

	}

	double VBFNLOPhasespace::generateTwoToNKinematics(const double* random,
	vector<Lorentz5Momentum>& momenta) {

	double weight;

	int id =
	olpId()[ProcessType::oneLoopInterference] ?
	olpId()[ProcessType::oneLoopInterference] :
	olpId()[ProcessType::treeME2];

	double* p = new double[4*momenta.size()];

	OLP_PhaseSpacePoint(&id, const_cast<double>(random), const_cast<double>(random+1), p, &weight);

	if (weight < 0) {
	throw Exception() << "VBFNLOPhasespace::generateTwoToNKinematics(): Negative weight in VBFNLOPhaseSpace\n"
	<< Exception::runerror;
	}

	if (weight == 0) {
	delete p;
	return 0;
	}

	for ( size_t i = 0; i < momenta.size(); ++i ) {
	momenta[i].setT(p[4i] GeV);
	momenta[i].setX(p[4i+1]GeV);
	momenta[i].setY(p[4i+2]GeV);
	momenta[i].setZ(p[4i+3]GeV);
	momenta[i].rescaleMass();
	}

	delete p;

	Energy beamenergy = sqrt(lastXCombPtr()->lastS())/2.;
	double x1 = momenta[0].e()/beamenergy;
	double x2 = momenta[1].e()/beamenergy;
	Energy2 thisSHat = (momenta[0] + momenta[1]).m2();

	// reshuffle so that particles have correct mass
	if ( needToReshuffle ) {

	// boost final-state into partonic CMS
	Boost toCMS = (momenta[0]+momenta[1]).findBoostToCM();
	for ( size_t i = 2; i < momenta.size(); ++i ) {
	momenta[i].boost(toCMS);
	}

	// copied from MatchboxRambo phasespace
	double xi;

	ReshuffleEquation solve(sqrt(thisSHat),mePartonData().begin()+2,mePartonData().end(),
	momenta.begin()+2,momenta.end());

	GSLBisection solver(1e-10,1e-8,10000);

	try {
	xi = solver.value(solve,0.0,1.1);
	} catch (GSLBisection::GSLerror) {
	return 0.;
	} catch (GSLBisection::IntervalError) {
	return 0.;
	}

	weight = pow(xi,3.(momenta.size()-3.));

	Energy num = ZERO;
	Energy den = ZERO;

	cPDVector::const_iterator d = mePartonData().begin()+2;
	for ( vector<Lorentz5Momentum>::iterator k = momenta.begin()+2;
	k != momenta.end(); ++k, ++d ) {
	num += (k).vect().mag2()/(k).t();
	Energy q = (*k).t();
	(k).setT(sqrt(sqr((d).hardProcessMass())+xixisqr((k).t())));
	(k).setVect(xi(*k).vect());
	weight = q/(k).t();
	den += (k).vect().mag2()/(k).t();
	(k).setMass((*d).hardProcessMass());
	}

	// unboost
	for ( size_t i = 2; i < momenta.size(); ++i ) {
	momenta[i].boost(-toCMS);
	}
	}

	if ( !matchConstraints(momenta) )
	return 0.;

	lastXCombPtr()->lastX1X2(make_pair(x1,x2));
	lastXCombPtr()->lastSHat(thisSHat);

	weight /= pow(thisSHat/GeV2,momenta.size()-4);
	weight /= x1*x2;

	fillDiagramWeights();

	return weight;

	}

	int VBFNLOPhasespace::nDimPhasespace(int nFinal) const {
	return 3*nFinal;

	//get this from within VBFNLO
	int pStatus = 0;
	double value, zero;
	string name = "PSdimension";
	OLP_GetParameter(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
	if ( pStatus != 1) {
	throw Exception() << "VBFNLOPhasespace::nDimPhasespace(): Cannot get phasespace dimension in VBFNLOPhaseSpace\n"
	<< "error code: " << pStatus << "\n"
	<< Exception::runerror;
	}
	// one additional number (first) needed for channel selection
	// one additional number (last) needed for global phi integration
	return value+2;
	}

	Energy VBFNLOPhasespace::ReshuffleEquation::operator() (double xi) const {
	cPDVector::const_iterator d = dataBegin;
	vector<Lorentz5Momentum>::const_iterator p = momentaBegin;
	Energy res = -w;
	for ( ; d != dataEnd; ++d, ++p ) {
	res += sqrt(sqr((**d).hardProcessMass()) +
	xixisqr(p->t()));
	}
	return res;
	}

	void VBFNLOPhasespace::doinit() {
	loadVBFNLO();
	MatchboxPhasespace::doinit();
	}

	void VBFNLOPhasespace::doinitrun() {
	loadVBFNLO();
	MatchboxPhasespace::doinitrun();
	}

	void VBFNLOPhasespace::persistentOutput(PersistentOStream & os) const {
	os << needToReshuffle << theLastXComb;
	}

	void VBFNLOPhasespace::persistentInput(PersistentIStream & is, int) {
	is >> needToReshuffle >> theLastXComb;
	}


	// * 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<VBFNLOPhasespace,MatchboxPhasespace>
	describeHerwigVBFNLOPhasespace("Herwig::VBFNLOPhasespace", "HwMatchboxVBFNLO.so");

	void VBFNLOPhasespace::Init() {

	static ClassDocumentation<VBFNLOPhasespace> documentation
	("VBFNLOPhasespace is an interface to the internal phasespace generator "
	"of VBFNLO. It uses the information passed via the BLHA interface to "
	"obtain information on the required channels.");

	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 6:26 PM (7 h, 56 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023797
Default Alt Text: (48 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions