No OneTemporary
Actions

Size

103 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc.in b/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc.in
	--- a/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc.in
	+++ b/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc.in
	@@ -1,773 +1,802 @@
	// -- C++ --
	//
	// GoSamAmplitude.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 GoSamAmplitude class.
	//

	#include "GoSamAmplitude.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/Interface/Command.h"

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

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

	#include <boost/progress.hpp>
	#include <boost/filesystem.hpp>

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

	using namespace Herwig;

	namespace bfs = boost::filesystem;

	////////////////////////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////////////////////////


	GoSamAmplitude::GoSamAmplitude(){
	theCodeExists=false;
	isitDR=false;
	theFormOpt=true;
	theNinja=true;
	theHiggsEff=false;}

	GoSamAmplitude::~GoSamAmplitude() {}

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

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


	////////////////////////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////////////////////////


	void GoSamAmplitude::doinit() {
	optionalContractFile() = name() + ".OLPContract.lh";
	MatchboxOLPME::doinit();
	}

	void GoSamAmplitude::doinitrun() {
	optionalContractFile() = name() + ".OLPContract.lh";
	MatchboxOLPME::doinitrun();
	}


	////////////////////////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////////////////////////


	extern "C" void OLP_Start(const char, int i);
	extern "C" void OLP_Polvec(double, double, double*);
	extern "C" void OLP_SetParameter(char, double, double, int);
	extern "C" void OLP_PrintParameter(char*);
	extern "C" void OLP_EvalSubProcess2(int, double, double, double, double*);


	////////////////////////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////////////////////////


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

	char char_cwd[256];
	getcwd(char_cwd, sizeof(char_cwd));
	string cwd = string(char_cwd);

	string folderMatchboxBuild = factory()->buildStorage();
	folderMatchboxBuild.erase(folderMatchboxBuild.begin());

	// set all necessary path and file names
	gosamPath = gosamPathInterface == "" ? cwd + folderMatchboxBuild + "GoSam" : gosamPathInterface;
	// When transitioning to C++ 11 this length()-1 workaround can be replaced by string.back()
	if (gosamPath.at(gosamPath.length()-1) != '/') gosamPath.append("/");
	gosamSourcePath = gosamPath + "source/";
	gosamInstallPath = gosamPath + "build/";

	// create all the directories
	if (!bfs::is_directory(gosamPath)){
	try {
	bfs::create_directory(gosamPath);
	} catch (exception& e) {
	cerr << "--------------------------------------------------------------------------------\n";
	cerr << "The following exception occured:\n\n";
	cerr << " " << e.what() << "\n\n";
	cerr << " -> Please create the parent directory of\n";
	cerr << " " << gosamPath << "\n";
	cerr << " manually!\n";
	cerr << "--------------------------------------------------------------------------------\n";
	abort();
	}
	}
	if (!bfs::is_directory(gosamSourcePath)) bfs::create_directory(gosamSourcePath);
	if (!bfs::is_directory(gosamInstallPath)) bfs::create_directory(gosamInstallPath);

	contractFileTitle = name() + ".OLPContract.lh";
	contractFileName = gosamPath + "/" + contractFileTitle;

	string orderFileName = gosamPath + "/" + name() + ".OLPOrder.lh";

	gosamSetupInFileName = gosamSetupInFileNameInterface == "" ? gosamPath + "/setup.gosam.in" : gosamSetupInFileNameInterface;

	// if it doesn't, then get it from src/defaults
	string cmd = "python @prefix@/bin/GoSamHelper.py ";
	cmd+=" --infile="+gosamSetupInFileName;
	cmd+=" --definfile=@prefix@/share/Herwig++/defaults/setup.gosam.in";
	// set form's tempdir to the current working directory
	cmd+=" --formtempdir="+StringUtils::replace(gosamSourcePath, string("/"), string("\\/")); //@FORMTEMPDIR@
	// set reduction programs to what is specified through the repository
	cmd+=" --reduction="+(theNinja ? string("ninja,golem95") : string("samurai,golem95")); //@REDUCTIONPROGRAMS@
	if (theNinja) generator()->log() << "\nGoSam will use Ninja as reduction program.\n" << flush;
	else if (!theNinja) generator()->log() << "\nGoSam will use Samurai as reduction program.\n" << flush;
	if (theFormOpt) generator()->log() << "Form optimization switched on.\n" << flush;
	else if (!theFormOpt) generator()->log() << "Form optimization switched off.\n" << flush;
	if (theNinja && !theFormOpt) throw Exception() << "Ninja reduction needs form optimization!\n" << Exception::abortnow;
	// use or don't use formopt
	cmd+=" --formopt="+(theFormOpt ? string("") : string(", noformopt")); //@FORMOPT@
	cmd+=" --higgseff="+(theHiggsEff ? string("smehc") : string("smdiag")); //model
	std::system(cmd.c_str());


	// check for old order file and create it if it doesn't already exist
	fillOrderFile(procs, orderFileName);

	ifstream ifile(contractFileName.c_str());
	if(!ifile){
	signOLP(orderFileName, contractFileName);
	}

	if ( !checkOLPContract(contractFileName) ) {
	throw Exception() << "failed to start GoSam" << Exception::abortnow;
	}

	- if (!( DynamicLoader::load(gosamInstallPath+"/lib/libgolem_olp.so")
	- \|\| DynamicLoader::load(gosamInstallPath+"/lib64/libgolem_olp.so"))) buildGoSam();
	+ if (!( DynamicLoader::load(gosamInstallPath+"/lib/libgolem_olp.so")
	+ \|\| DynamicLoader::load(gosamInstallPath+"/lib64/libgolem_olp.so")
	+ \|\| DynamicLoader::load(gosamInstallPath+"/lib/libgolem_olp.dylib")
	+ \|\| DynamicLoader::load(gosamInstallPath+"/lib64/libgolem_olp.dylib"))) buildGoSam();

	int status = -1;
	startOLP(contractFileTitle, status);

	if ( status != 1 ) return false;

	return true;

	}

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

	string tempcontract = contract;

	char char_cwd[256];
	getcwd(char_cwd, sizeof(char_cwd));
	string cwd = string(char_cwd);

	string folderMatchboxBuild = factory()->buildStorage();
	folderMatchboxBuild.erase(folderMatchboxBuild.begin());

	gosamPath = gosamPathInterface == "" ? cwd + folderMatchboxBuild + "GoSam" : gosamPathInterface;
	// When transitioning to C++ 11 this length()-1 workaround can be replaced by string.back()
	if (gosamPath.at(gosamPath.length()-1) != '/') gosamPath.append("/");

	- if (!( DynamicLoader::load(gosamPath+"build/lib/libgolem_olp.so")
	- \|\| DynamicLoader::load(gosamPath+"build/lib64/libgolem_olp.so")))
	+ if (!( DynamicLoader::load(gosamPath+"build/lib/libgolem_olp.so")
	+ \|\| DynamicLoader::load(gosamPath+"build/lib64/libgolem_olp.so")
	+ \|\| DynamicLoader::load(gosamPath+"build/lib/libgolem_olp.dylib")
	+ \|\| DynamicLoader::load(gosamPath+"build/lib64/libgolem_olp.dylib")))
	throw Exception() << "Failed to load GoSam. Please check the log file.\n"
	<< Exception::abortnow;
	tempcontract = gosamPath + tempcontract;

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

	// hand over input parameters for EW scheme considered
	int pStatus = 0;
	double zero = 0.0;
	if ( SM().ewScheme() == 0 \|\| SM().ewScheme() == 6 ) { // EW/Scheme Default and EW/Scheme Independent
	throw Exception() << "`Best value' schemes are not supported by GoSam"
	<< Exception::abortnow;
	} else if ( SM().ewScheme() == 4 ) { // EW/Scheme mW (uses mW,GF,sin2thetaW) seems not to be supported by GoSam
	throw Exception() << "`mW' scheme is not supported by GoSam"
	<< Exception::abortnow;
	} else if ( SM().ewScheme() == 1 ) { // EW/Scheme GMuScheme (uses mW,mZ,GF)
	double in1=getParticleData(ParticleID::Z0)->mass()/GeV;
	double in2=getParticleData(ParticleID::Wplus)->mass()/GeV;
	double in3=SM().fermiConstant()*GeV2;
	OLP_SetParameter((char *)"mass(23)",&in1,&zero,&pStatus);
	OLP_SetParameter((char *)"mass(24)",&in2,&zero,&pStatus);
	OLP_SetParameter((char *)"Gf",&in3,&zero,&pStatus);
	} else if ( SM().ewScheme() == 2 ) { // EW/Scheme alphaMZScheme (uses mW,mZ,alpha(mZ))
	double in1=getParticleData(ParticleID::Z0)->mass()/GeV;
	double in2=getParticleData(ParticleID::Wplus)->mass()/GeV;
	double in3=SM().alphaEMMZ();
	OLP_SetParameter((char *)"mass(23)",&in1,&zero,&pStatus);
	OLP_SetParameter((char *)"mass(24)",&in2,&zero,&pStatus);
	OLP_SetParameter((char *)"alpha",&in3,&zero,&pStatus);
	} else if ( SM().ewScheme() == 3 ) { // EW/Scheme NoMass (uses alpha(mZ),GF,sin2thetaW)
	double in1=SM().fermiConstant()*GeV2;
	double in2=SM().alphaEMMZ();
	double in3=SM().sin2ThetaW();
	OLP_SetParameter((char *)"Gf",&in1,&zero,&pStatus);
	OLP_SetParameter((char *)"alpha",&in2,&zero,&pStatus);
	OLP_SetParameter((char *)"sw2",&in3,&zero,&pStatus);
	} else if ( SM().ewScheme() == 5 ) { // EW/Scheme mZ (uses mZ,alphaEM,sin2thetaW)
	double in1=getParticleData(ParticleID::Z0)->mass()/GeV;
	double in2=SM().alphaEMMZ();
	double in3=SM().sin2ThetaW();
	OLP_SetParameter((char *)"mass(23)",&in1,&zero,&pStatus);
	OLP_SetParameter((char *)"alpha",&in2,&zero,&pStatus);
	OLP_SetParameter((char *)"sw2",&in3,&zero,&pStatus);
	}

	// hand over mass and width of the Higgs
	double wH = getParticleData(25)->width()/GeV;
	double mH = getParticleData(25)->mass()/GeV;
	OLP_SetParameter((char*)"width(25)",&wH,&zero,&pStatus);
	OLP_SetParameter((char*)"mass(25)",&mH,&zero,&pStatus);

	// hand over input parameter for alphaS
	double as = SM().alphaS();
	OLP_SetParameter((char *)"alphaS", &as, &zero, &pStatus);

	// fill massive Particle vector
	if (massiveParticles.empty())
	for (int i=1; i<=6; ++i)
	if (getParticleData(i)->mass()/GeV > 0.0) massiveParticles.push_back(i);
	-
	- // hand over quark masses and widths (iff massive)
	+
	+ // hand over quark masses and widths (iff massive)
	if ( massiveParticles.size() != 0 ) {
	for ( vector<int>::const_iterator mID = massiveParticles.begin(); mID != massiveParticles.end(); ++mID ) {
	- double mass=getParticleData(*mID)->mass()/GeV;
	- double width=getParticleData(*mID)->width()/GeV;
	- string mstr="mass("+static_cast<ostringstream>(&(ostringstream()<<(mID)))->str()+")";
	- string wstr="width("+static_cast<ostringstream>(&(ostringstream()<<(mID)))->str()+")";
	- OLP_SetParameter((char *)&mstr,&mass, &zero, &pStatus );
	- OLP_SetParameter((char *)&wstr,&width, &zero, &pStatus );
	+ string mstr;
	+ string wstr;
	+ int mInt = *mID;
	+ double mass=getParticleData(mInt)->mass()/GeV;
	+ double width=getParticleData(mInt)->width()/GeV;
	+ std::stringstream ss;
	+ ss << mInt;
	+ string str = ss.str();
	+ mstr="mass("+str+")";
	+ wstr="width("+str+")";
	+ char * mchar = new char[mstr.size()+1];
	+ char * wchar = new char[wstr.size()+1];
	+ std::copy(mstr.begin(),mstr.end(),mchar);
	+ std::copy(wstr.begin(),wstr.end(),wchar);
	+ mchar[mstr.size()] = '\0';
	+ wchar[wstr.size()] = '\0';
	+ OLP_SetParameter( mchar, &mass, &zero, &pStatus );
	+ OLP_SetParameter( wchar, &width, &zero, &pStatus );
	+ delete[] mchar;
	+ delete[] wchar;
	+
	+// Nicer but not working properly:
	+
	+// double mass=getParticleData(*mID)->mass()/GeV;
	+// double width=getParticleData(*mID)->width()/GeV;
	+// string mstr="mass("+static_cast<ostringstream>(&(ostringstream()<<(mID)))->str()+")";
	+// string wstr="width("+static_cast<ostringstream>(&(ostringstream()<<(mID)))->str()+")";
	+// cout<<"\n massiv "<<mstr;
	+//
	+// OLP_SetParameter((char *)&mstr,&mass, &zero, &pStatus );
	+// OLP_SetParameter((char *)&wstr,&width, &zero, &pStatus );
	}
	}

	// Note: In the GoSam input file, the standard is to set the parameter
	// 'symmetries' for quark families and lepton generations, which allow
	// for flavour changing only between families/generations. If this pa-
	// rameter is set, GoSam won't allow to set electron and muon mass and
	// width via the interface. Also setting mass and width for the tau is
	// not yet considered.

	// print OLP parameters
	if ( Debug::level > 1 ) {
	string ppstr = factory()->runStorage() + name() + ".OLPParameters.lh";
	OLP_PrintParameter(const_cast<char*>(ppstr.c_str()));
	}

	didStartOLP() = true;

	}


	////////////////////////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////////////////////////


	void GoSamAmplitude::fillOrderFile(const map<pair<Process, int>, int>& procs, string orderFileName) {

	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::loopInducedME2 ) {
	Typestr << "LoopInduced";
	} else if ( (*p).first.second == ProcessType::colourCorrelatedME2 ) {
	Typestr << "ccTree";
	} else if ( (*p).first.second == ProcessType::spinColourCorrelatedME2 ) {
	Typestr << "scTree";
	} else if ( (*p).first.second == ProcessType::oneLoopInterference ) {
	Typestr << "Loop";
	}
	gosamprocinfo pro = gosamprocinfo((*p).second, -1, Processstr.str(), Typestr.str());
	pro.setOAs(p->first.first.orderInAlphaS);
	pro.setOAew(p->first.first.orderInAlphaEW);
	processmap[(*p).second] = pro;
	}

	ifstream oldOrderFileStream(orderFileName.c_str());
	if (oldOrderFileStream){
	oldOrderFileStream.close();
	return;
	}

	ofstream orderFile(orderFileName.c_str());
	int asPower = 100;
	int minlegs = 100;
	int maxlegs = -1;
	int maxasPower = -1;
	int aewPower = 100;
	int maxaewPower = -1;

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

	orderFile << "# OLP order file created by Herwig++/Matchbox for GoSam\n\n";
	orderFile << "InterfaceVersion BLHA2\n";
	orderFile << "MatrixElementSquareType CHsummed\n";
	orderFile << "CorrectionType QCD\n";
	orderFile << "IRregularisation " << (isDR() ? "DRED" : "CDR") << "\n";

	// loop over quarks to check if they have non-zero masses
	for (int i=1; i<=6; ++i) if (getParticleData(i)->mass()/GeV > 0.0) massiveParticles.push_back(i);
	if ( massiveParticles.size() != 0 ) {
	orderFile << "MassiveParticles ";
	for ( vector<int>::const_iterator mID = massiveParticles.begin(); mID != massiveParticles.end(); ++mID ) {
	int mInt = *mID;
	orderFile << mInt << " ";
	}
	orderFile << "\n";
	}
	orderFile << "\n";

	vector < string > types;
	types.push_back("Tree");
	types.push_back("LoopInduced");
	types.push_back("ccTree");
	types.push_back("scTree");
	types.push_back("Loop");

	for ( int i = asPower ; i != maxasPower + 1 ; i++ ) {
	for ( int j = aewPower ; j != maxaewPower + 1 ; j++ ) {
	orderFile << "\nAlphasPower " << i << "\n";
	orderFile << "AlphaPower " << j << "\n";
	for ( vector<string>::iterator it = types.begin() ; it != types.end() ; it++ ) {
	if ( *it == "LoopInduced" ) continue;
	for ( map<int, gosamprocinfo>::iterator p = processmap.begin() ; p != processmap.end() ; ++p )
	if ( (p).second.Tstr() == it && i == (p).second.orderAs() && j == (p).second.orderAew() ) {
	orderFile << "\nAmplitudeType " << *it << "\n";
	break;
	}
	for ( map<int, gosamprocinfo>::iterator p = processmap.begin() ; p != processmap.end() ; ++p )
	if ( (p).second.Tstr() == it && i == (p).second.orderAs() && j == (p).second.orderAew() ) {
	orderFile << (*p).second.Pstr() << "\n";
	}
	}
	}
	}

	// Write out the loop induced processes separately
	int asPowerLI = 100;
	int aewPowerLI = 100;
	for ( map<int, gosamprocinfo>::iterator p = processmap.begin() ; p != processmap.end() ; ++p ) {
	if ( (*p).second.Tstr() != "LoopInduced" ) continue;
	if ( (p).second.orderAs() != asPowerLI \|\| (p).second.orderAew() != aewPowerLI ) {
	asPowerLI = (*p).second.orderAs();
	aewPowerLI = (*p).second.orderAew();
	// At the moment GoSam requires for qcd loop induced processes the as coupling power
	// which would correspond to an associated fictitious Born process
	orderFile << "\nAlphasPower " << (asPowerLI-2) << "\n";
	orderFile << "AlphaPower " << aewPowerLI << "\n";
	orderFile << "\nAmplitudeType " << "LoopInduced" << "\n";
	}
	orderFile << (*p).second.Pstr() << "\n";
	}
	orderFile << flush;

	}

	void GoSamAmplitude::signOLP(const string& order, const string& contract) {
	if(!theCodeExists){
	char char_cwd[256];
	getcwd(char_cwd, sizeof(char_cwd));
	string cwd = string(char_cwd);

	string folderMatchboxBuild = factory()->buildStorage();
	folderMatchboxBuild.erase(folderMatchboxBuild.begin());

	generator()->log() << "\n>>> generating GoSam amplitudes. This may take some time, please be patient.\n"
	<< ">>> see " + cwd + folderMatchboxBuild + "gosam-amplitudes.log for details.\n" << flush;
	string cmd = "@GOSAMPREFIX@/bin/gosam.py --olp --output-file=" + contract + " --config=" +
	gosamSetupInFileName + " --destination=" + gosamSourcePath + " " + order + " > " + cwd + folderMatchboxBuild + "gosam-amplitudes.log 2>&1";
	std::system(cmd.c_str());
	cmd = "python @prefix@/bin/GoSamHelper.py ";
	cmd += " --makelink ";
	cmd += " --makelinkfrom=contract ";
	cmd += " --makelinkto="+factory()->buildStorage() + name() + ".OLPContract.lh";
	std::system(cmd.c_str());
	}
	}

	bool GoSamAmplitude::checkOLPContract(string contractFileName) {

	ifstream infile(contractFileName.c_str());
	string line;
	vector < string > contractfile;

	while (std::getline(infile, line)) contractfile.push_back(line);

	for ( map<int, gosamprocinfo>::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 && (p).second.Pstr().length() == (linex).find("\|") ) {
	string sub = (linex).substr((linex).find("\|") + 1, (linex).find("#") - (linex).find("\|") - 1); // \| 1 23 # buggy??
	if ( sub.find(" 1 ") != 0 )
	throw Exception() << "Failed to check contractfile. Please check the logfile.\n"
	<< Exception::abortnow;
	string subx = sub.substr(3);
	int subint;
	istringstream(subx) >> subint;
	(*p).second.setGID(subint);
	}
	}
	}
	}

	string ids = factory()->buildStorage() + "GoSam.ids.dat";
	ofstream IDS(ids.c_str());

	idpair.clear();
	for ( map<int, gosamprocinfo>::iterator p = processmap.begin() ; p != processmap.end() ; p++ )
	idpair.push_back(-1);
	idpair.push_back(-1);

	for ( map<int, gosamprocinfo>::iterator p = processmap.begin() ; p != processmap.end() ; p++ ) {

	idpair[(p).second.HID()]=(p).second.GID();
	IDS << (p).second.HID() << " " << (p).second.GID() << " " << (*p).second.Tstr() << "\n";
	if ( (*p).second.GID() == -1 ) return 0;
	}
	IDS << flush;

	return 1;

	}

	bool GoSamAmplitude::buildGoSam() {
	if(!theCodeExists){
	generator()->log() << "\n>>> compiling GoSam amplitudes. This may take some time, please be patient.\n"
	<< ">>> see " + gosamSourcePath + "gosam-build.log for details.\n\n" << flush;
	string cmd = "cd " + gosamSourcePath + " && sh autogen.sh FCFLAGS=-g --prefix=" +
	gosamInstallPath + " --disable-static > gosam-build.log 2>&1";
	std::system(cmd.c_str());
	if (!gosamBuildScript.empty()) {
	cmd = "cd " + gosamSourcePath + " && " + gosamBuildScript + " >> gosam-build.log 2>&1";
	std::system(cmd.c_str());
	}
	std::system(cmd.c_str());
	cmd = "cd " + gosamSourcePath + " && make install >> gosam-build.log 2>&1";
	std::system(cmd.c_str());
	}
	theCodeExists=true;
	return 1;
	}

	void GoSamAmplitude::getids() const {
	string line = factory()->buildStorage() + "GoSam.ids.dat";
	ifstream infile(line.c_str());
	int hid;
	int gid;
	string type;
	while (std::getline(infile, line)) {
	idpair.push_back(-1);
	idtypepair.push_back(" ");
	}
	infile.close();
	string line2 = factory()->buildStorage() + "GoSam.ids.dat";
	ifstream infile2(line2.c_str());
	idpair.push_back(-1);
	idtypepair.push_back(" ");
	while (std::getline(infile2, line2)) {
	istringstream(line2) >> hid >> gid >> type;
	idpair[hid]=gid;

	idtypepair[hid]=type;
	}
	infile.close();
	}


	////////////////////////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////////////////////////


	void GoSamAmplitude::evalSubProcess() const {

	useMe();

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

	double scale = sqrt(mu2() / GeV2);

	double out[7] = { };
	if ( idpair.size() == 0 ){ getids(); }
	int id = -99;
	if ( olpId()[ProcessType::loopInducedME2] ) id = olpId()[ProcessType::loopInducedME2];
	else if ( olpId()[ProcessType::oneLoopInterference] ) id = olpId()[ProcessType::oneLoopInterference];
	else id = olpId()[ProcessType::treeME2];

	int callid(idpair[id]); // If id denotes the Herwig ID, this returns the GoSam ID
	string calltype(idtypepair[id]); // If id denotes the Herwig ID, this returns the amplitude type
	double acc;

	OLP_EvalSubProcess2(&(callid), 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[3] * units);
	} else if ( olpId()[ProcessType::loopInducedME2] ) {
	lastTreeME2(out[2] * units);
	}

	}

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

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

	double scale = sqrt(mu2() / GeV2);

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

	if ( idpair.size() == 0 ) getids();

	int callid(idpair[olpId()[ProcessType::colourCorrelatedME2]]);
	double acc;
	OLP_EvalSubProcess2(&(callid), olpMomenta(), &scale, &colourCorrelatorResults[0], &acc);

	cPDVector particles = lastXComb().matrixElement()->mePartonData();

	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 GoSamAmplitude::evalSpinColourCorrelator(pair<int , int > ) const {

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

	double scale = sqrt(mu2() / GeV2);

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

	if ( idpair.size() == 0 ) getids();

	double acc;
	int callid(idpair[olpId()[ProcessType::spinColourCorrelatedME2]]);
	OLP_EvalSubProcess2(&(callid), olpMomenta(), &scale, &spinColourCorrelatorResults[0], &acc);

	for ( int i = 0; i < n; ++i ) {
	for ( int j = 0; j < n; ++j ) {
	Complex scc(spinColourCorrelatorResults[2i+2nj]units, spinColourCorrelatorResults[2i+2nj+1]units);
	lastColourSpinCorrelator(make_pair(i,j),scc);
	}
	}

	}

	LorentzVector<Complex> GoSamAmplitude::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] ={ };
	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;

	}


	////////////////////////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////////////////////////


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


	void GoSamAmplitude::persistentOutput(PersistentOStream & os) const {
	os << idpair << idtypepair << processmap << gosamPathInterface
	<< gosamSetupInFileNameInterface << gosamBuildScript << gosamPath
	<< gosamSourcePath << gosamInstallPath << gosamSetupInFileName
	<< orderFileTitle << contractFileTitle
	<< contractFileName << orderFileName
	<< theCodeExists << theFormOpt << theNinja << isitDR << massiveParticles<<theHiggsEff;
	}

	void GoSamAmplitude::persistentInput(PersistentIStream & is, int) {
	is >> idpair >> idtypepair >> processmap >> gosamPathInterface
	>> gosamSetupInFileNameInterface >> gosamBuildScript >> gosamPath
	>> gosamSourcePath >> gosamInstallPath >> gosamSetupInFileName
	>> orderFileTitle >> contractFileTitle
	>> contractFileName >> orderFileName
	>> theCodeExists >> theFormOpt >> theNinja >> isitDR >> massiveParticles>>theHiggsEff;
	}


	// * 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<GoSamAmplitude, MatchboxOLPME> describeHerwigGoSamAmplitude("Herwig::GoSamAmplitude", "HwMatchboxGoSam.so");

	void GoSamAmplitude::Init() {

	static ClassDocumentation<GoSamAmplitude> documentation("GoSamAmplitude implements an interface to GoSam.", "Matrix elements have been calculated using GoSam");

	static Parameter<GoSamAmplitude,string> interfaceProcessPath
	("ProcessPath",
	"Prefix for the process source code, include files and library produced by GoSam.",
	&GoSamAmplitude::gosamPathInterface, "",
	false, false);

	static Parameter<GoSamAmplitude,string> interfaceSetupInFilename
	("SetupInFilename",
	"File name of the GoSam infile (typically setup.gosam.in) to be used. If left empty a new setup.gosam.in is created in the location specified in Path",
	&GoSamAmplitude::gosamSetupInFileNameInterface, "",
	false, false);

	static Switch<GoSamAmplitude,bool> interfaceCodeExists
	("CodeExists",
	"Switch on or off if Code already exists/not exists.",
	&GoSamAmplitude::theCodeExists, true, false, false);

	static SwitchOption interfaceCodeExistsOn
	(interfaceCodeExists,
	"True",
	"Switch True if Code already exists.",
	true);

	static SwitchOption interfaceCodeExistsOff
	(interfaceCodeExists,
	"False",
	"Switch False if Code has to be build.",
	false);

	static Switch<GoSamAmplitude,bool> interfaceisitDR
	("isDR",
	"Switch on or off DR.",
	&GoSamAmplitude::isitDR, false, false, false);
	static SwitchOption interfaceisitDROn
	(interfaceisitDR,
	"True",
	"Switch True.",
	true);
	static SwitchOption interfaceisitDROff
	(interfaceisitDR,
	"False",
	"Switch False.",
	false);

	static Switch<GoSamAmplitude,bool> interfaceFormOpt
	("FormOpt",
	"Switch On/Off formopt",
	&GoSamAmplitude::theFormOpt, true, false, false);
	static SwitchOption interfaceFormOptOn
	(interfaceFormOpt,
	"On",
	"On",
	true);
	static SwitchOption interfaceFormOptOff
	(interfaceFormOpt,
	"Off",
	"Off",
	false);

	static Switch<GoSamAmplitude,bool> interfaceNinja
	("Ninja",
	"Switch On/Off for reduction with Ninja. If Off then Samurai is used.",
	&GoSamAmplitude::theNinja, true, false, false);
	static SwitchOption interfaceNinjaOn
	(interfaceNinja,
	"On",
	"On",
	true);
	static SwitchOption interfaceNinjaOff
	(interfaceNinja,
	"Off",
	"Off",
	false);

	static Switch<GoSamAmplitude,bool> interfaceHiggsEff
	("HiggsEff",
	"Switch On/Off for effective higgs model.",
	&GoSamAmplitude::theHiggsEff, false, false, false);
	static SwitchOption interfaceHiggsEffOn
	(interfaceHiggsEff,
	"On",
	"On",
	true);
	static SwitchOption interfaceHiggsEffOff
	(interfaceHiggsEff,
	"Off",
	"Off",
	false);

	static Parameter<GoSamAmplitude,string> interfaceBuildScript
	("BuildScript",
	"File name of a custom build script, which is called between 'autogen.sh'"
	"and 'make install'. It can be used for parallelization.",
	&GoSamAmplitude::gosamBuildScript, "",
	false, false);

	}

	diff --git a/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc.in b/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc.in
	--- a/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc.in
	+++ b/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc.in
	@@ -1,700 +1,712 @@
	// -- C++ --
	//
	// MadGraphAmplitude.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 MadGraphAmplitude class.
	//

	#include "MadGraphAmplitude.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 "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "ThePEG/PDT/EnumParticles.h"
	#include <boost/lexical_cast.hpp>
	#include <boost/filesystem.hpp>


	#include <cstdlib>
	#include <dlfcn.h>
	#include <errno.h>
	+#include <sstream>

	using namespace Herwig;


	extern "C" void mginitproc_(char *i,int);
	extern "C" void MG_Calculate_wavefunctions_virt(int* proc,double,double);
	extern "C" void MG_Calculate_wavefunctions_born(int* proc,double, int);
	extern "C" void MG_Jamp (int* proc,int, double);
	extern "C" void MG_LNJamp (int* proc,int, double);
	extern "C" void MG_Virt (int* proc,double*);
	extern "C" void MG_NCol (int* proc,int*);
	extern "C" void MG_vxxxxx (double* p,double* n,int* inc,double* );
	extern "C" void MG_Colour (int* proc,int* i,int* j ,int* color);




	MadGraphAmplitude::MadGraphAmplitude()
	: theMGmodel("loop_sm"),keepinputtopmass(false) {}

	MadGraphAmplitude::~MadGraphAmplitude() {

	}

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

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

	map<string,int>& MadGraphAmplitude::processmap() {
	static map<string,int> theProcessMap;
	return theProcessMap;
	}

	template<class ProcType>
	string processId(const ProcType& p) {
	string orderFileName = "";
	for (typename ProcType::const_iterator it=p.begin();it!=p.end();it++)
	orderFileName +=boost::lexical_cast<string>( (*it)->id());
	orderFileName+=".MG";
	return orderFileName;
	}

	-bool initializedMad=false;
	+bool MadGraphAmplitude::initializedMad=false;


	void MadGraphAmplitude::initProcess(const cPDVector& ) {

	if ( lastMatchboxXComb()->initialized() )
	return;

	string MGProcLibPath = theProcessPath == "" ? factory()->buildStorage()+"MadGraphAmplitudes" : theProcessPath;
	if (MGProcLibPath.at(MGProcLibPath.length()-1) != '/') MGProcLibPath.append("/");

	if ( !DynamicLoader::load(MGProcLibPath+"InterfaceMadGraph.so") )
	throw Exception() << "Failed to load MadGraph amplitudes\n"
	<< DynamicLoader::lastErrorMessage
	<< Exception::abortnow;

	if (!initializedMad){
	string mstr=(factory()->runStorage()+"MadGraphAmplitudes"+"/param_card"+((theMGmodel=="loop_sm")?"":("_"+theMGmodel))+".dat");
	size_t len = mstr.size();
	mginitproc_(const_cast<char*>(mstr.c_str()),len);
	initializedMad=true;
	}
	lastMatchboxXComb()->isInitialized();
	}

	-bool ranMadGraphInitializeExternal = false;
	+

	bool MadGraphAmplitude::initializeExternal() {

	loadProcessMap();

	string MGProcLibPath = theProcessPath == "" ? factory()->buildStorage()+"MadGraphAmplitudes" : theProcessPath;
	if (MGProcLibPath.at(MGProcLibPath.length()-1) != '/') MGProcLibPath.append("/");

	if ( boost::filesystem::exists(MGProcLibPath) ) {
	if ( !boost::filesystem::is_directory(MGProcLibPath) )
	throw Exception() << "MadGraph amplitude storage '"
	<< MGProcLibPath << "' existing but not a directory."
	<< Exception::abortnow;
	} else {
	boost::filesystem::create_directories(MGProcLibPath);
	}

	string born= MGProcLibPath+"BornAmplitudes.dat";
	if ( !boost::filesystem::exists(born) ) {
	ofstream borns(born.c_str());
	for (vector<string>::iterator amps=BornAmplitudes.begin();amps!=BornAmplitudes.end();amps++)
	borns<<*amps<<endl;
	borns.close();
	}

	string virt= MGProcLibPath+"VirtAmplitudes.dat";
	if ( !boost::filesystem::exists(virt) ) {
	ofstream virts(virt.c_str());
	for (vector<string>::iterator amps=VirtAmplitudes.begin();amps!=VirtAmplitudes.end();amps++)
	virts<<*amps<<endl;
	virts.close();
	}

	string runAmplitudes = factory()->runStorage() + "/MadGraphAmplitudes";
	-
	if ( boost::filesystem::exists(runAmplitudes) ) {
	if ( !boost::filesystem::is_directory(runAmplitudes) )
	throw Exception() << "MadGraph amplitude storage '"
	<< runAmplitudes << "' existing but not a directory."
	<< Exception::abortnow;
	} else {
	boost::filesystem::create_directories(runAmplitudes);
	}
	-
	+
	+
	+ //EW-consistency check:
	Energy MW=getParticleData(ParticleID::Wplus)->mass();
	Energy MZ=getParticleData(ParticleID::Z0)->mass();
	-
	if( MW!= sqrt(MZMZ/2.0+sqrt(MZMZMZMZ/4.0-Constants::piSM().alphaEMMZ()MZ*MZ/ sqrt(2.0)/SM().fermiConstant()))){
	cerr<<"\n\n-----!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!-----";
	cerr << "\nYou are using a EW scheme which is inconsistent with the MadGraph parametisation:\n\n"
	<<MW/GeV<< " GeV==MW!= sqrt(MZ^2/2+sqrt(MZ^4/4.0-pialphaEMMZMZ^2/ sqrt(2)/G_f))=="<<
	sqrt(MZMZ/2.0+sqrt(MZMZMZMZ/4.0-Constants::piSM().alphaEMMZ()MZ*MZ/ sqrt(2.0)/SM().fermiConstant()))/GeV
	<<" GeV\n\n-----!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!-----\n";
	}



	+
	+ string para= factory()->runStorage()+"/MadGraphAmplitudes"+"/MG-Parameter.dat";
	+
	+
	+

	-
	- string para= factory()->runStorage()+"/MadGraphAmplitudes"+"/MG-Parameter.dat";
	ofstream params(para.c_str());
	params<<"$WZ$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Z0)->width() /GeV;
	params<<"\n$WW$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Wplus)->width()/GeV;
	params<<"\n$alphas$ " <<std::setiosflags(ios::scientific) <<SM().alphaS();
	params<<"\n$GF$ " <<std::setiosflags(ios::scientific) <<SM().fermiConstant()*GeV2 ;
	params<<"\n$alphaMZ$ " <<std::setiosflags(ios::scientific) <<1/SM().alphaEMMZ();
	params<<"\n$MZ$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Z0)->mass() /GeV<<flush;
	params<<"\n$MW$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Wplus)->mass() /GeV<<flush;
	params<<"\n$sw2$ " <<std::setiosflags(ios::scientific) << SM().sin2ThetaW() <<flush;
	if(theMGmodel=="heft"&&!keepinputtopmass){
	cerr<<"\n---------------------------------------------------------------";
	cerr<<"\n---------------------------------------------------------------";
	cerr<<"\nNote: You are using the Higgs Effective model (heft) in ";
	cerr<<"\n Madgraph. We assume you try to calculate NLO with ";
	cerr<<"\n the GoSam virtual amplitudes. To match the models we ";
	cerr<<"\n therefore set the topmass to 10000000 GeV.";
	cerr<<"\n\n For more information see the \\tau parameter in:";
	cerr<<"\n https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Models/HiggsEffective";
	cerr<<"\n\n The Effective Higgs model in Gosam is using mT=infinity";
	cerr<<"\n\n\n If you want to use the LO matrixelements of MadGraph with finite' topmass you need to add: ";
	cerr<<"\n\n set Madgraph:KeepInputTopMass True";
	cerr<<"\n\n to your input file.";
	cerr<<"\n---------------------------------------------------------------";
	cerr<<"\n---------------------------------------------------------------\n";
	params<<"\n$MT$ 10000000." <<flush;
	}else{
	params<<"\n$MT$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::t)->mass() /GeV <<flush;
	}
	params<<"\n$WT$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::t)->width() /GeV <<flush;
	params<<"\n$MB$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::b)->mass() /GeV <<flush;
	params<<"\n$MH$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::h0)->mass() /GeV <<flush;
	params<<"\n$WH$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::h0)->width() /GeV <<flush;
	params<<"\n$MTA$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::tauplus)->mass() /GeV <<flush;


	string cmd = "python @prefix@/bin/mg2Matchbox.py ";
	cmd +=" --buildpath "+MGProcLibPath;
	cmd +=" --model "+theMGmodel;
	cmd +=" --runpath "+factory()->runStorage()+"/MadGraphAmplitudes " ;
	-
	- cmd +=" --orderas "+static_cast<ostringstream*>(&(ostringstream()<<(factory()->orderInAlphaS())))->str() ;
	-
	- cmd +=" --orderew "+static_cast<ostringstream*>(&(ostringstream()<<(factory()->orderInAlphaEW())))->str();
	+ std::stringstream as,aem;
	+ as << factory()->orderInAlphaS();
	+ cmd +=" --orderas "+as.str() ;
	+ aem <<factory()->orderInAlphaEW();
	+ cmd +=" --orderew "+aem.str();

	// TODO move to boost::system

	if ( !processmap().empty() ){
	//set the parameters

	std::system(cmd.c_str());
	ranMadGraphInitializeExternal = true;
	return true;
	}
	char cwd[1024];
	if ( !getcwd(cwd,sizeof(cwd)) )
	throw Exception() << "failed to determine current working directory\n"
	<< Exception::abortnow;

	cmd +=" --madgraph @MADGRAPHPREFIX@/bin " ;
	cmd +="--build > ";
	cmd += MGProcLibPath+"MG.log 2>&1";


	generator()->log() << "\n\nCompiling MadGraph amplitudes. This may take some time -- please be patient.\n"
	<< "In case of problems see " << MGProcLibPath << "MG.log for details.\n\n"
	<< flush;
	std::system(cmd.c_str());


	cmd = "python @prefix@/bin/mg2Matchbox.py ";
	cmd +=" --buildpath "+MGProcLibPath;
	cmd +=" --model "+theMGmodel;
	cmd +=" --runpath "+factory()->runStorage()+"/MadGraphAmplitudes " ;
	-
	- cmd +=" --orderas "+static_cast<ostringstream*>(&(ostringstream()<<(factory()->orderInAlphaS())))->str() ;
	-
	- cmd +=" --orderew "+static_cast<ostringstream*>(&(ostringstream()<<(factory()->orderInAlphaEW())))->str();
	+ as.clear();
	+ aem.clear();
	+ as << factory()->orderInAlphaS();
	+ cmd +=" --orderas "+as.str() ;
	+ aem <<factory()->orderInAlphaEW();
	+ cmd +=" --orderew "+aem.str();
	+
	+
	+
	+
	std::system(cmd.c_str());






	loadProcessMap();

	ranMadGraphInitializeExternal = true;

	return !processmap().empty();

	}

	int MadGraphAmplitude::externalId(const cPDVector& proc) {
	for (int i=0;i<100;i++){
	colourindex.push_back(-2);
	}
	assert(!BornAmplitudes.empty()\|\|!VirtAmplitudes.empty());

	int res=0;
	string amp="";
	int k=0;
	for (cPDVector::const_iterator it=proc.begin();it!=proc.end();it++,k++){
	amp+=boost::lexical_cast<string>( (*it)->id())+" ";if (k==1)amp+=" > ";
	}


	for (vector<string>::iterator amps=BornAmplitudes.begin();amps!=BornAmplitudes.end();amps++,res++)
	if(amp==*amps)return res+1;
	for (vector<string>::iterator amps=VirtAmplitudes.begin();amps!=VirtAmplitudes.end();amps++,res++)
	if(amp==*amps)return res+1;

	throw Exception() << "One amplitude has no externalId. Please remove the MadGraphAmplitude-folder and rebuild.\n" << Exception::abortnow;



	return res;
	}

	void MadGraphAmplitude::loadProcessMap() {
	if ( !processmap().empty() )
	return;
	string MGProcLibPath = theProcessPath == "" ? factory()->buildStorage()+"MadGraphAmplitudes" : theProcessPath;
	if (MGProcLibPath.at(MGProcLibPath.length()-1) != '/') MGProcLibPath.append("/");
	string line = MGProcLibPath+"processmap.dat";
	ifstream infile(line.c_str());
	if ( !infile )
	return;
	string process;
	int id;
	while (std::getline(infile, line)) {
	std::stringstream ss(line);
	getline( ss, process, ' ' );
	ss >> id;
	processmap().insert(make_pair(process,id));
	}
	}
	+bool MadGraphAmplitude::ranMadGraphInitializeExternal = false;

	void MadGraphAmplitude::doinit() {
	+
	if ( !ranMadGraphInitializeExternal ) {
	initializeExternal();
	}
	MatchboxAmplitude::doinit();
	}

	void MadGraphAmplitude::doinitrun() {
	if ( !ranMadGraphInitializeExternal ) {
	- loadProcessMap();
	- ranMadGraphInitializeExternal = true;
	+ initializeExternal();
	}
	MatchboxAmplitude::doinitrun();
	}

	bool MadGraphAmplitude::canHandle(const PDVector& p,
	Ptr<MatchboxFactory>::tptr factory,
	bool virt) const {
	if ( factory->processData()->diagramMap().find(p) !=
	factory->processData()->diagramMap().end() )
	return true;
	vector<Ptr<Tree2toNDiagram>::ptr> diags =
	factory->diagramGenerator()->generate(p,orderInGs(),orderInGem());
	if ( diags.empty() )
	return false;
	factory->processData()->diagramMap()[p] = diags;
	string amp="";
	int k=0;
	for (PDVector::const_iterator it=p.begin();it!=p.end();it++,k++){
	amp+=boost::lexical_cast<string>( (*it)->id())+" ";if (k==1)amp+=" > ";
	}
	if (virt && factory->highestVirt()>=p.size()){
	VirtAmplitudes.push_back(amp);
	}else{
	BornAmplitudes.push_back(amp);
	}

	return true;
	}

	void MadGraphAmplitude::prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr me) {

	useMe();

	if ( !calculateTreeAmplitudes() ) {
	MatchboxAmplitude::prepareAmplitudes(me);
	return;
	}
	lastMatchboxXComb()->clearheljamp();
	lastMatchboxXComb()->clearhelLNjamp();
	initProcess(mePartonData());
	MatchboxAmplitude::prepareAmplitudes(me);

	}

	Complex MadGraphAmplitude::evaluate(size_t i, const vector<int>& hel, Complex& largeN) {


	//find the colourline:
	int ii = -1;
	int xx=lastMatchboxXComb()->externalId();
	if(colourindex[i]!=-2){

	ii = colourindex[i];
	if (ii==-1) {
	largeN = Complex(0.0);
	return Complex(0.0);
	}
	} else {
	set<vector<size_t> > a = colourOrdering(i);
	int ncol=-1;
	MG_NCol(&xx,&ncol);
	assert(ncol!=-1);
	for( int it = 0; it < ncol; it++ ){
	int n = 0;
	for ( cPDVector::const_iterator nx = mePartonData().begin();
	nx != mePartonData().end(); nx++ )
	if ( (*nx)->coloured() ) n++;
	set<vector<size_t> > tmpset;
	vector<size_t> tmpvek;
	for ( int it2 = 0; it2 < n; it2++ ) {
	int ret=-2;
	MG_Colour(&xx,&it,&it2,&ret);
	assert(ret !=-2);
	if (ret== -1)
	break;
	if ( ret == 0 ) {
	n++;
	tmpset.insert(tmpvek);
	tmpvek.clear();
	} else {
	tmpvek.push_back(ret-1);
	}
	if( it2 == n-1 ) tmpset.insert(tmpvek);
	}
	bool found_all = true;
	for ( set<vector<size_t> >::iterator it3 = a.begin(); it3 != a.end(); it3++ ) {
	bool found_it3=false;
	for ( set<vector<size_t> >::iterator it4 = tmpset.begin(); it4 != tmpset.end(); it4++ ) {
	vector<size_t> it3tmp = gluonsFirst(*it3);
	vector<size_t> it4tmp = (*it4);
	if ( it3tmp.size() != it4tmp.size() ) continue;
	if ( it3tmp == it4tmp ) found_it3 = true;
	}
	found_all = found_all && found_it3;
	}

	if ( found_all ) {
	colourindex[i]=it;


	ii=it;
	}
	}
	}
	if ( ii == -1 ){

	colourindex[i]=ii;
	largeN = Complex(0.0);
	return Complex(0.0);
	}

	const map<vector<int>,vector < complex<double> > >& tmp = lastMatchboxXComb()->heljamp();
	const map<vector<int>,vector < complex<double> > >& tmpLN = lastMatchboxXComb()->helLNjamp();

	if( tmp.find(hel) != tmp.end()) {
	largeN = tmpLN.find(hel)->second[ii];
	return tmp.find(hel)->second[ii];;
	}

	double units = pow(sqrt(lastSHat())/GeV,int(hel.size())-4);



	for(size_t j=0;j<hel.size();j++){
	int cross=crossingMap()[j];
	if( (cross>1&&j<=1)\|\|(cross<=1&&j>1)){
	heltmp[cross]=-1*hel[j];}
	else{heltmp[cross]=hel[j];}
	}


	size_t j=0;
	for (size_t i=0;i<mePartonData().size();i++){
	momenta[j]=abs(meMomenta()[i].e()/GeV)<1.e-13?0.:double(meMomenta()[i].e()/GeV);
	momenta[j+1]=abs(meMomenta()[i].x()/GeV)<1.e-13?0.:double(meMomenta()[i].x()/GeV);
	momenta[j+2]=abs(meMomenta()[i].y()/GeV)<1.e-13?0.:double(meMomenta()[i].y()/GeV);
	momenta[j+3]=abs(meMomenta()[i].z()/GeV)<1.e-13?0.:double(meMomenta()[i].z()/GeV);
	j+=4;
	}

	MG_Calculate_wavefunctions_born(&xx, &momenta[0], &heltmp[0]);


	int ncol=-1;
	MG_NCol(&xx,&ncol);

	Complex res;
	Complex resLN;
	for( int it = 0; it < ncol; it++ ){
	double dd[2];
	MG_Jamp(&xx,&it,&dd[0]);
	Complex d(dd[0],dd[1]);
	if(it==ii)res=d*units;
	lastMatchboxXComb()->pushheljamp(hel,d*units);
	double ddLN[2];
	MG_LNJamp(&xx,&it,&ddLN[0]);
	Complex dLN(ddLN[0],ddLN[1]);
	if(it==ii)resLN=dLN*units;
	lastMatchboxXComb()->pushhelLNjamp(hel,dLN*units);
	}


	largeN = resLN;
	return res;

	}

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




	int tmp=i;
	pg[0]=p.e()/GeV;pg[1]=p.x()/GeV;pg[2]=p.y()/GeV;pg[3]=p.z()/GeV;
	ng[0]=n.e()/GeV;ng[1]=n.x()/GeV;ng[2]=n.y()/GeV;ng[3]=n.z()/GeV;

	MG_vxxxxx(&pg[0],&ng[0],&tmp,&poltmp[0]);

	complex<double> pol[6];
	pol[0]=Complex(poltmp[0],poltmp[1]);
	pol[1]=Complex(poltmp[2],poltmp[3]);
	pol[2]=Complex(poltmp[4],poltmp[5]);
	pol[3]=Complex(poltmp[6],poltmp[7]);

	LorentzVector<Complex> polarization(pol[1],pol[2],pol[3],pol[0]);

	return polarization.conjugate();
	}

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

	vector<size_t> MadGraphAmplitude::gluonsFirst(vector<size_t> vec) {
	vector<size_t> vecout;
	for(vector<size_t>::iterator it= vec.begin();it!= vec.end();++it)
	if ( mePartonData()[crossingMap()[*it]]->id()==21)
	vecout.push_back(crossingMap()[*it]);

	for(vector<size_t>::iterator it= vec.begin();it!= vec.end();++it)
	if ( mePartonData()[crossingMap()[*it]]->id()!=21)
	vecout.push_back(crossingMap()[*it]);

	return vecout;

	}

	double MadGraphAmplitude::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 = crossingMap()[ij.first];
	iCrossed = ij.first;
	for ( unsigned int k = 0; k < crossingMap().size(); ++k )
	if ( crossingMap()[k] == ij.first ) {
	iCrossed = k;
	break;
	}

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

	}




	void MadGraphAmplitude::prepareOneLoopAmplitudes(Ptr<MatchboxMEBase>::tcptr ){
	assert(false);
	}


	double MadGraphAmplitude::oneLoopInterference() const {
	if ( !calculateOneLoopInterference() )
	return lastOneLoopInterference();
	evaloneLoopInterference();
	return lastOneLoopInterference();
	}


	void MadGraphAmplitude::evaloneLoopInterference() const {


	double units = pow(lastSHat()/GeV2,int(mePartonData().size())-4);



	size_t j=0;
	for (size_t i=0;i<mePartonData().size();i++){
	momenta[j]=abs(meMomenta()[i].e()/GeV)<1.e-13?0.:double(meMomenta()[i].e()/GeV);
	momenta[j+1]=abs(meMomenta()[i].x()/GeV)<1.e-13?0.:double(meMomenta()[i].x()/GeV);
	momenta[j+2]=abs(meMomenta()[i].y()/GeV)<1.e-13?0.:double(meMomenta()[i].y()/GeV);
	momenta[j+3]=abs(meMomenta()[i].z()/GeV)<1.e-13?0.:double(meMomenta()[i].z()/GeV);
	j+=4;
	}
	int xx=lastMatchboxXComb()->externalId();

	MG_Calculate_wavefunctions_virt(&xx,&momenta[0],&virt[0]);

	double ifact = 1.;

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

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


	ifact *= lastMatchboxXComb()->matchboxME()->finalStateSymmetry();


	lastOneLoopInterference(virt[1]/ifact*units);
	}

	void MadGraphAmplitude::persistentOutput(PersistentOStream & os) const {
	os << theOrderInGs << theOrderInGem << BornAmplitudes << VirtAmplitudes
	<< colourindex<<crossing << theProcessPath << theMGmodel;
	}

	void MadGraphAmplitude::persistentInput(PersistentIStream & is, int) {
	is >> theOrderInGs >> theOrderInGem >> BornAmplitudes >> VirtAmplitudes
	>> colourindex>>crossing >> theProcessPath >> theMGmodel;
	}

	// * 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<MadGraphAmplitude,MatchboxAmplitude>
	describeHerwigMadGraphAmplitude("Herwig::MadGraphAmplitude", "HwMatchboxMadGraph.so");

	void MadGraphAmplitude::Init() {

	static ClassDocumentation<MadGraphAmplitude> documentation ("MadGraphAmplitude","Matrix elements have been calculated using MadGraph5");

	static Parameter<MadGraphAmplitude,string> interfaceProcessPath
	("ProcessPath",
	"The Process Path.",
	&MadGraphAmplitude::theProcessPath, "",false, false);

	static Parameter<MadGraphAmplitude,string> interfaceModel
	("Model",
	"The MadGraph-Model.",
	&MadGraphAmplitude::theMGmodel, "loop_sm",false, false);

	static Switch<MadGraphAmplitude,bool> interfacekeepinputtopmass
	("KeepInputTopMass",
	"Switch On/Off formopt",
	&MadGraphAmplitude::keepinputtopmass, false, false, false);
	static SwitchOption interfacekeepinputtopmassTrue
	(interfacekeepinputtopmass,
	"On",
	"On",
	true);
	static SwitchOption interfacekeepinputtopmassFalse
	(interfacekeepinputtopmass,
	"Off",
	"Off",
	false);



	}


	diff --git a/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.h b/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.h
	--- a/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.h
	+++ b/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.h
	@@ -1,322 +1,329 @@
	// -- C++ --
	//
	// MadGraphAmplitude.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef Herwig_MadGraphAmplitude_H
	#define Herwig_MadGraphAmplitude_H
	//
	// This is the declaration of the MadGraphAmplitude class.
	//

	#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
	#include "Herwig++/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxCurrents.h"
	#include "ThePEG/Utilities/DynamicLoader.h"

	namespace Herwig {


	using namespace ThePEG;

	/**
	* \ingroup Matchbox
	* \author Johannes Bellm, Simon Platzer
	*
	* \brief MadGraphAmplitude implements an interface to MadGraph
	*/
	class MadGraphAmplitude:
	public MatchboxAmplitude {

	public:

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

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

	public:

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

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

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

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

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

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

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

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

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

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

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

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

	void evaloneLoopInterference() const;


	/**
	* Return true, if one loop corrections are given in the conventions
	* of BDK.
	*/

	virtual bool isCS() const { return false; }
	virtual bool isExpanded() const { return true; }
	virtual bool isBDK() const { return false; }
	virtual bool isDR() const { return false; }
	virtual bool isDRbar() 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 sqr(91.18800000000*GeV);}
	//virtual Energy2 mu2() const { return lastSHat(); }


	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;

	/**
	* Order process in MadGraph conventions
	*/
	vector<size_t> gluonsFirst(vector<size_t> i);

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

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

	/**
	* Initialize this amplitude
	*/
	virtual bool initializeExternal();

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

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

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}


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

	/**
	* The (tree-level) order in \f$g_S\f$ in which this matrix
	* element is given.
	*/
	unsigned int theOrderInGs;

	/**
	* The (tree-level) order in \f$g_{EM}\f$ in which this matrix
	* element is given.
	*/
	unsigned int theOrderInGem;



	/**
	* The path for the process libraries.
	*/
	string theProcessPath;

	/**
	* The path to generate amplitudes in.
	*/
	string theMGmodel;

	bool keepinputtopmass;

	/**
	* The process map.
	*/
	static map<string,int>& processmap();

	/**
	* Load the process map
	*/
	void loadProcessMap();

	/**
	* Initialize the given process
	*/
	void initProcess(const cPDVector&);


	/**
	* Storage for Amplitudes
	*/
	mutable vector<string> BornAmplitudes,VirtAmplitudes;

	/**
	* Helper for color and crossing handling
	*/
	mutable vector<int> colourindex, crossing;

	/**
	* Give this to Madgraph:
	*/

	mutable double virt[20],momenta[50],pg[4],ng[4],poltmp[8];
	mutable int heltmp[10];
	-
	+
	+
	+ /**
	+ * Static Variables to handle initialization.
	+ */
	+
	+ static bool ranMadGraphInitializeExternal;
	+ static bool initializedMad;



	//@}

	private:

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


	};

	}

	#endif /* Herwig_MadGraphAmplitude_H */
	diff --git a/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc.in b/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc.in
	--- a/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc.in
	+++ b/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc.in
	@@ -1,392 +1,393 @@
	// -- 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;

	OpenLoopsAmplitude::OpenLoopsAmplitude() {theCodeExists=false;
	}

	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;

	- bool success = DynamicLoader::load("@OPENLOOPSLIBS@/libopenloops.so");
	+ bool success = DynamicLoader::load("@OPENLOOPSLIBS@/libopenloops.so")
	+ \|\|DynamicLoader::load("@OPENLOOPSLIBS@/libopenloops.dylib") ;

	if ( !success ) {
	- throw Exception() << "Failed to load libopenloops.so\n"
	+ throw Exception() << "Failed to load libopenloops.so/dylib\n"
	<< DynamicLoader::lastErrorMessage
	<< Exception::abortnow;
	}

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

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

	OLP_Start(tempcontract.c_str(), &status);
	int a=0;double null=0.0;double one=1.0;
	int part[8]={1,2,3,4,5,6,23,24};string stri;
	for (int i=0;i<8;i++){
	double mass=getParticleData(part[i])->mass()/GeV;
	double width=getParticleData(part[i])->width()/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);

	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 << "\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());
	double acc ;
	double scale = sqrt(mu2() / GeV2);
	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());
	double acc ;
	double scale = sqrt(mu2() / GeV2);
	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());

	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 <<theCodeExists;
	}

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

	// * 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.");

	static Switch<OpenLoopsAmplitude,bool> interfaceCodeExists
	("CodeExists",
	"Switch on or off if Code already exists/not exists.",
	&OpenLoopsAmplitude::theCodeExists, true, false, false);
	static SwitchOption interfaceCodeExistsOn
	(interfaceCodeExists,
	"True",
	"Switch True if Code already exists.",
	true);
	static SwitchOption interfaceCodeExistsOff
	(interfaceCodeExists,
	"False",
	"Switch False if Code has to be build.",
	false);


	}

	diff --git a/m4/herwig.m4 b/m4/herwig.m4
	--- a/m4/herwig.m4
	+++ b/m4/herwig.m4
	@@ -1,858 +1,899 @@
	# check for gcc bug http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34130
	AC_DEFUN([HERWIG_CHECK_ABS_BUG],
	[
	AC_REQUIRE([HERWIG_COMPILERFLAGS])
	if test "$GCC" = "yes"; then
	AC_MSG_CHECKING([for gcc abs bug])
	AC_RUN_IFELSE([
	AC_LANG_PROGRAM(
	[[ int foo (int i) { return -2 * __builtin_abs(i - 2); } ]],
	[[ if ( foo(1) != -2 \|\| foo(3) != -2 ) return 1; ]]
	)],
	[ AC_MSG_RESULT([not found. Compiler is ok.]) ],
	[
	AC_MSG_RESULT([found. Builtin abs() is buggy.])
	AC_MSG_CHECKING([if -fno-builtin-abs works])
	oldcxxflags=$CXXFLAGS
	CXXFLAGS="$CXXFLAGS -fno-builtin-abs"
	AC_RUN_IFELSE([
	AC_LANG_PROGRAM(
	[[
	#include <cstdlib>
	int foo (int i) { return -2 * std::abs(i - 2); }
	]],
	[[
	if (foo(1) != -2 \|\| foo(3) != -2) return 1;
	]]
	)],
	[
	AC_MSG_RESULT([yes. Setting -fno-builtin-abs.])
	AM_CXXFLAGS="$AM_CXXFLAGS -fno-builtin-abs"
	AM_CFLAGS="$AM_CFLAGS -fno-builtin-abs"
	],
	[
	AC_MSG_RESULT([no. Setting -fno-builtin.])
	AC_MSG_WARN([
	*****************************************************************************
	For this version of gcc, -fno-builtin-abs alone did not work to avoid the
	gcc abs() bug. Instead, all gcc builtin functions are now disabled.
	Update gcc if possible.
	*****************************************************************************])
	AM_CXXFLAGS="$AM_CXXFLAGS -fno-builtin"
	AM_CFLAGS="$AM_CFLAGS -fno-builtin"
	]
	)
	CXXFLAGS=$oldcxxflags
	]
	)
	fi
	])

	dnl ##### THEPEG #####
	AC_DEFUN([HERWIG_CHECK_THEPEG],
	[
	defaultlocation="${prefix}"
	test "x$defaultlocation" = xNONE && defaultlocation="${ac_default_prefix}"
	AC_MSG_CHECKING([for libThePEG in])
	AC_ARG_WITH(thepeg,
	AC_HELP_STRING([--with-thepeg=DIR],[location of ThePEG installation]),
	[],
	[with_thepeg="${defaultlocation}"])
	AC_MSG_RESULT([$with_thepeg])

	if test "x$with_thepeg" = "xno"; then
	AC_MSG_ERROR([Cannot build Herwig++ without ThePEG. Please set --with-thepeg.])
	fi

	THEPEGLDFLAGS="-L${with_thepeg}/lib/ThePEG"

	THEPEGHASLHAPDF="no"
	if test -e ${with_thepeg}/lib/ThePEG/ThePEGLHAPDF.so ; then
	THEPEGHASLHAPDF="yes"
	fi
	if test "${host_cpu}" == "x86_64" -a -e ${with_thepeg}/lib64/ThePEG/libThePEG.so ; then
	THEPEGLDFLAGS="-L${with_thepeg}/lib64/ThePEG"
	if test -e ${with_thepeg}/lib64/ThePEG/ThePEGLHAPDF.so ; then
	THEPEGHASLHAPDF="yes"
	fi
	fi

	if test "x$THEPEGHASLHAPDF" == "xno" ; then
	AC_MSG_ERROR([Herwig++ requires ThePEG to be build with lhapdf.])
	fi

	THEPEGHASFASTJET="no"
	if test -e ${with_thepeg}/lib/ThePEG/FastJetFinder.so ; then
	THEPEGHASFASTJET="yes"
	fi
	if test "${host_cpu}" == "x86_64" -a -e ${with_thepeg}/lib64/ThePEG/libThePEG.so ; then
	THEPEGLDFLAGS="-L${with_thepeg}/lib64/ThePEG"
	if test -e ${with_thepeg}/lib64/ThePEG/FastJetFinder.so ; then
	THEPEGHASFASTJET="yes"
	fi
	fi

	if test "x$THEPEGHASFASTJET" == "xno" ; then
	AC_MSG_ERROR([Herwig++ requires ThePEG to be build with FastJet.])
	fi

	THEPEGPATH="${with_thepeg}"

	oldldflags="$LDFLAGS"
	oldlibs="$LIBS"

	LDFLAGS="$LDFLAGS $THEPEGLDFLAGS"
	AC_CHECK_LIB([ThePEG],[debugThePEG],[],
	[AC_MSG_ERROR([No ThePEG libraries in $THEPEGLDFLAGS. Please set --with-thepeg.])])

	AC_SUBST([THEPEGLIB],[-lThePEG])
	AC_SUBST(THEPEGLDFLAGS)
	AC_SUBST(THEPEGPATH)
	AC_SUBST(THEPEGHASLHAPDF)
	AC_SUBST(THEPEGHASFASTJET)

	LIBS="$oldlibs"
	LDFLAGS="$oldldflags"

	AC_MSG_CHECKING([for ThePEG headers in])
	AC_ARG_WITH([thepeg-headers],
	AC_HELP_STRING([--with-thepeg-headers=DIR],[location of ThePEG include directory]),
	[],
	[with_thepeg_headers="${with_thepeg}/include"])
	AC_MSG_RESULT([$with_thepeg_headers])

	if test "x$with_thepeg_headers" = "xno"; then
	AC_MSG_ERROR([Cannot build Herwig++ without ThePEG headers. Please set --with-thepeg-headers.])
	fi

	THEPEGINCLUDE="-I$with_thepeg_headers"

	oldcppflags="$CPPFLAGS"
	CPPFLAGS="$CPPFLAGS $THEPEGINCLUDE"
	AC_CHECK_HEADER([ThePEG/Config/ThePEG.h],[],
	[AC_MSG_ERROR([No ThePEG headers in $with_thepeg_headers. Please set --with-thepeg-headers.])])
	CPPFLAGS="$oldcppflags"

	AC_SUBST(THEPEGINCLUDE)

	AC_MSG_CHECKING([for HepMCAnalysis.so in ThePEG])


	if test -e "$THEPEGPATH/lib/ThePEG/HepMCAnalysis.so" ; then
	CREATE_HEPMC="create"
	AC_MSG_RESULT([found])
	else
	CREATE_HEPMC="# create"
	AC_MSG_RESULT([not found])
	fi

	AC_SUBST([CREATE_HEPMC])
	])

	dnl ##### LOOPTOOLS #####
	AC_DEFUN([HERWIG_LOOPTOOLS],
	[
	AC_REQUIRE([AC_PROG_FC])
	AC_REQUIRE([AC_FC_LIBRARY_LDFLAGS])
	AC_REQUIRE([AC_PROG_CC])
	AC_REQUIRE([HERWIG_COMPILERFLAGS])

	AC_MSG_CHECKING([if Looptools build works])
	enable_looptools=yes

	if test "x$GCC" = "xyes"; then
	case "${host}" in
	x86_64-\|-darwin1*)
	AM_FCFLAGS="$AM_FCFLAGS -fdefault-integer-8"
	;;
	esac

	AC_LANG_PUSH([Fortran])
	oldFCFLAGS="$FCFLAGS"
	FCFLAGS="$AM_FCFLAGS"
	AC_COMPILE_IFELSE(
	AC_LANG_PROGRAM([],[ print *[,]"Hello"]),
	[],
	[AC_MSG_RESULT([no])
	AC_MSG_ERROR([needs gfortran on 64bit machines])]
	)
	FCFLAGS="$oldFCFLAGS"
	AC_LANG_POP([Fortran])
	fi
	AC_MSG_RESULT([$enable_looptools])

	AC_SUBST([F77],[$FC])
	AC_SUBST([FFLAGS],[$FCFLAGS])
	AC_SUBST([AM_FFLAGS],[$AM_FCFLAGS])
	AC_SUBST([FLIBS],[$FCLIBS])
	])

	dnl ##### VBFNLO #####
	AC_DEFUN([HERWIG_CHECK_VBFNLO],
	[
	AC_MSG_CHECKING([for VBFNLO])

	AC_ARG_WITH([vbfnlo],
	AS_HELP_STRING([--with-vbfnlo=DIR], [Installation path of VBFNLO]),
	[],
	[with_vbfnlo=no]

	)

	AC_MSG_RESULT([$with_vbfnlo])

	AS_IF([test "x$with_vbfnlo" != "xno"],
	[AC_CHECK_FILES(
	${with_vbfnlo}/lib/VBFNLO/libVBFNLO.so,
	[have_vbfnlo=lib], [have_vbfnlo=no])],
	[have_vbfnlo=no])

	AS_IF([test "x$with_vbfnlo" != "xno" -a "x$have_vbfnlo" = "xno" ],
	[AC_CHECK_FILES(
	${with_vbfnlo}/lib64/VBFNLO/libVBFNLO.so,
	[have_vbfnlo=lib64], [have_vbfnlo=no])])

	AS_IF([test "x$have_vbfnlo" = "xlib"],
	[VBFNLOLIBS=${with_vbfnlo}/lib/VBFNLO
	AC_SUBST(VBFNLOLIBS)
	])

	AS_IF([test "x$have_vbfnlo" = "xlib64"],
	[VBFNLOLIBS=${with_vbfnlo}/lib64/VBFNLO
	AC_SUBST(VBFNLOLIBS)
	])

	AS_IF([test "x$with_vbfnlo" != "xno" -a "x$have_vbfnlo" = "xno"],
	[AC_MSG_ERROR([vbfnlo requested but not found])])

	AM_CONDITIONAL(HAVE_VBFNLO,[test "x$have_vbfnlo" = "xlib" -o "x$have_vbfnlo" = "xlib64"])

	if test "x$have_vbfnlo" = "xlib" -o "x$have_vbfnlo" = "xlib64" ; then
	VBFNLOINCLUDE=${with_vbfnlo}/include
	AC_SUBST(VBFNLOINCLUDE)
	VBFNLOLIB=${with_vbfnlo}/${have_vbfnlo}/VBFNLO
	AC_SUBST(VBFNLOLIB)
	LOAD_VBFNLO="library"
	CREATE_VBFNLO="create"
	INSERT_VBFNLO="insert"
	SET_VBFNLO="set"
	MKDIR_VBFNLO="mkdir"
	else
	LOAD_VBFNLO="# library"
	CREATE_VBFNLO="# create"
	INSERT_VBFNLO="# insert"
	SET_VBFNLO="# set"
	MKDIR_VBFNLO="# mkdir"
	fi

	AC_SUBST([LOAD_VBFNLO])
	AC_SUBST([CREATE_VBFNLO])
	AC_SUBST([INSERT_VBFNLO])
	AC_SUBST([SET_VBFNLO])
	AC_SUBST([MKDIR_VBFNLO])

	])

	dnl ##### nlojet #####
	AC_DEFUN([HERWIG_CHECK_NLOJET],
	[
	AC_MSG_CHECKING([for nlojet])

	AC_ARG_WITH([nlojet],
	AS_HELP_STRING([--with-nlojet=DIR], [Installation path of nlojet]),
	[],
	[with_nlojet=no]

	)

	AC_MSG_RESULT([$with_nlojet])

	AS_IF([test "x$with_nlojet" != "xno"],
	[AC_CHECK_FILES(
	${with_nlojet}/lib/libnlo++.so,
	[have_nlojet=lib], [have_nlojet=no])],
	[have_nlojet=no])

	AS_IF([test "x$with_nlojet" != "xno" -a "x$have_nlojet" = "xno" ],
	[AC_CHECK_FILES(
	${with_nlojet}/lib64/libnlo++.so,
	[have_nlojet=lib64], [have_nlojet=no])])

	AS_IF([test "x$have_nlojet" = "xlib"],
	[NLOJETLIBPATH=${with_nlojet}/lib
	AC_SUBST(NLOJETLIBPATH)
	NLOJETINCLUDEPATH=${with_nlojet}/include
	AC_SUBST(NLOJETINCLUDEPATH)
	])

	AS_IF([test "x$have_nlojet" = "xlib64"],
	[NLOJETLIBPATH=${with_nlojet}/lib64
	AC_SUBST(NLOJETLIBPATH)
	NLOJETINCLUDEPATH=${with_nlojet}/include
	AC_SUBST(NLOJETINCLUDEPATH)
	])

	AS_IF([test "x$with_nlojet" != "xno" -a "x$have_nlojet" = "xno"],
	[AC_MSG_ERROR([nlojet requested but not found])])

	AM_CONDITIONAL(HAVE_NLOJET,[test "x$have_nlojet" = "xlib" -o "x$have_nlojet" = "xlib64"])

	if test "x$have_nlojet" = "xlib" -o "x$have_nlojet" = "xlib64" ; then
	LOAD_NLOJET="library"
	CREATE_NLOJET="create"
	INSERT_NLOJET="insert"
	else
	LOAD_NLOJET="# library"
	CREATE_NLOJET="# create"
	INSERT_NLOJET="# insert"
	fi

	AC_SUBST([LOAD_NLOJET])
	AC_SUBST([CREATE_NLOJET])
	AC_SUBST([INSERT_NLOJET])

	])

	dnl ##### njet #####
	AC_DEFUN([HERWIG_CHECK_NJET],
	[
	AC_MSG_CHECKING([for njet])

	AC_ARG_WITH([njet],
	AS_HELP_STRING([--with-njet=DIR], [Installation path of njet]),
	[],
	[with_njet=no]

	)

	AC_MSG_RESULT([$with_njet])

	AS_IF([test "x$with_njet" != "xno"],
	[AC_CHECK_FILES(
	${with_njet}/lib/libnjet2.so,
	[have_njet=lib], [have_njet=no])],
	[have_njet=no])

	AS_IF([test "x$with_njet" != "xno" -a "x$have_njet" = "xno" ],
	[AC_CHECK_FILES(
	${with_njet}/lib64/libnjet2.so,
	[have_njet=lib64], [have_njet=no])])

	AS_IF([test "x$have_njet" = "xlib"],
	[NJETLIBPATH=${with_njet}/lib
	AC_SUBST(NJETLIBPATH)
	NJETINCLUDEPATH=${with_njet}/include
	AC_SUBST(NJETINCLUDEPATH)
	NJETPREFIX=${with_njet}
	AC_SUBST(NJETPREFIX)
	])

	AS_IF([test "x$have_njet" = "xlib64"],
	[NJETLIBPATH=${with_njet}/lib64
	AC_SUBST(NJETLIBPATH)
	NJETINCLUDEPATH=${with_njet}/include
	AC_SUBST(NJETINCLUDEPATH)
	NJETPREFIX=${with_njet}
	AC_SUBST(NJETPREFIX)
	])

	AS_IF([test "x$with_njet" != "xno" -a "x$have_njet" = "xno"],
	[AC_MSG_ERROR([njet requested but not found])])

	AM_CONDITIONAL(HAVE_NJET,[test "x$have_njet" = "xlib" -o "x$have_njet" = "xlib64"])

	if test "x$have_njet" = "xlib" -o "x$have_njet" = "xlib64" ; then
	LOAD_NJET="library"
	CREATE_NJET="create"
	INSERT_NJET="insert"
	else
	LOAD_NJET="# library"
	CREATE_NJET="# create"
	INSERT_NJET="# insert"
	fi

	AC_SUBST([LOAD_NJET])
	AC_SUBST([CREATE_NJET])
	AC_SUBST([INSERT_NJET])

	])



	dnl ##### gosam #####
	AC_DEFUN([HERWIG_CHECK_GOSAM],
	[
	AC_MSG_CHECKING([for gosam])

	AC_ARG_WITH([gosam],
	AS_HELP_STRING([--with-gosam=DIR], [Installation path of gosam]),
	[],
	[with_gosam=no]
	)

	AC_MSG_RESULT([$with_gosam])

	AS_IF([test "x$with_gosam" != "xno"],
	[AC_CHECK_FILES(
	${with_gosam}/bin/gosam.py,
	[have_gosam=lib], [have_gosam=no])],
	[have_gosam=no])

	AS_IF([test "x$have_gosam" = "xlib"],
	[GOSAMPREFIX=${with_gosam}
	AC_SUBST(GOSAMPREFIX)
	])

	AS_IF([test "x$with_gosam" != "xno" -a "x$have_gosam" = "xno"],
	[AC_MSG_ERROR([GoSam requested but not found])])

	AM_CONDITIONAL(HAVE_GOSAM,[test "x$have_gosam" = "xlib" ])

	if test "x$have_gosam" = "xlib" ; then
	LOAD_GOSAM="library"
	CREATE_GOSAM="create"
	INSERT_GOSAM="insert"
	else
	LOAD_GOSAM="# library"
	CREATE_GOSAM="# create"
	INSERT_GOSAM="# insert"
	fi

	AC_SUBST([LOAD_GOSAM])
	AC_SUBST([CREATE_GOSAM])
	AC_SUBST([INSERT_GOSAM])


	])


	dnl ##### gosam-contrib #####
	AC_DEFUN([HERWIG_CHECK_GOSAM_CONTRIB],
	[
	AC_MSG_CHECKING([for gosam-contrib])

	AC_ARG_WITH([gosam-contrib],
	AS_HELP_STRING([--with-gosam-contrib=DIR], [Installation path of gosam-contrib]),
	[],
	[with_gosam_contrib=no]
	)

	AC_MSG_RESULT([$with_gosam_contrib])

	AS_IF([test "x$with_gosam_contrib" = "xno" -a "x$with_gosam" != "xno"],
	[AC_CHECK_FILES(
	${with_gosam}/lib/libsamurai.so,
	[with_gosam_contrib=${with_gosam}], [])
	])

	AS_IF([test "x$with_gosam_contrib" = "xno" -a "x$with_gosam" != "xno"],
	[AC_CHECK_FILES(
	${with_gosam}/lib64/libsamurai.so,
	[with_gosam_contrib=${with_gosam}], [])
	])

	+AS_IF([test "x$with_gosam_contrib" = "xno" -a "x$with_gosam" != "xno"],
	+ [AC_CHECK_FILES(
	+ ${with_gosam}/lib/libsamurai.dylib,
	+ [with_gosam_contrib=${with_gosam}], [])
	+])

	+AS_IF([test "x$with_gosam_contrib" = "xno" -a "x$with_gosam" != "xno"],
	+ [AC_CHECK_FILES(
	+ ${with_gosam}/lib64/libsamurai.dylib,
	+ [with_gosam_contrib=${with_gosam}], [])
	+])

	AS_IF([test "x$with_gosam_contrib" = "xno" -a "x$with_gosam" != "xno"],
	[AC_MSG_ERROR([GoSam requested without requesting GoSam-Contrib])])

	AS_IF([test "x$with_gosam_contrib" != "xno"],
	[AC_CHECK_FILES(
	${with_gosam_contrib}/lib/libsamurai.so,
	[have_gosam_contrib=lib], [have_gosam_contrib=no])],
	[have_gosam_contrib=no])

	AS_IF([test "x$with_gosam_contrib" != "xno" -a "x$have_gosam_contrib" = "xno" ],
	[AC_CHECK_FILES(
	${with_gosam_contrib}/lib64/libsamurai.so,
	[have_gosam_contrib=lib64], [have_gosam_contrib=no])])

	+AS_IF([test "x$with_gosam_contrib" != "xno" -a "x$have_gosam_contrib" = "xno" ],
	+ [AC_CHECK_FILES(
	+ ${with_gosam_contrib}/lib/libsamurai.dylib,
	+ [have_gosam_contrib=lib], [have_gosam_contrib=no])])
	+
	+AS_IF([test "x$with_gosam_contrib" != "xno" -a "x$have_gosam_contrib" = "xno" ],
	+ [AC_CHECK_FILES(
	+ ${with_gosam_contrib}/lib64/libsamurai.dylib,
	+ [have_gosam_contrib=lib64], [have_gosam_contrib=no])])
	+
	+
	+
	+
	+
	+
	+
	AS_IF([test "x$have_gosam_contrib" != "xno"],
	[GOSAMCONTRIBPREFIX=${with_gosam_contrib}
	AC_SUBST(GOSAMCONTRIBPREFIX)
	])

	AS_IF([test "x$have_gosam_contrib" = "xlib"],
	[GOSAMCONTRIBLIBS=${with_gosam_contrib}/lib
	AC_SUBST(GOSAMCONTRIBLIBS)
	])

	AS_IF([test "x$have_gosam_contrib" = "xlib64"],
	[GOSAMCONTRIBLIBS=${with_gosam_contrib}/lib64
	AC_SUBST(GOSAMCONTRIBLIBS)
	])

	AS_IF([test "x$with_gosam_contrib" != "xno" -a "x$have_gosam_contrib" = "xno"],
	[AC_MSG_ERROR([GoSam-Contrib requested but not found])])

	AM_CONDITIONAL(HAVE_GOSAM_CONTRIB,[test "x$have_gosam_contrib" = "xlib" -o "x$have_gosam_contrib" = "xlib64"])

	if test "x$have_gosam_contrib" = "xlib" -o "x$have_gosam_contrib" = "xlib64" ; then
	LOAD_GOSAM_CONTRIB="library"
	CREATE_GOSAM_CONTRIB="create"
	INSERT_GOSAM_CONTRIB="insert"
	else
	LOAD_GOSAM_CONTRIB="# library"
	CREATE_GOSAM_CONTRIB="# create"
	INSERT_GOSAM_CONTRIB="# insert"
	fi

	AC_SUBST([LOAD_GOSAM_CONTRIB])
	AC_SUBST([CREATE_GOSAM_CONTRIB])
	AC_SUBST([INSERT_GOSAM_CONTRIB])


	])


	dnl ##### OpenLoops #####
	AC_DEFUN([HERWIG_CHECK_OPENLOOPS],
	[
	AC_MSG_CHECKING([for OpenLoops])

	AC_ARG_WITH([openloops],
	AS_HELP_STRING([--with-openloops=DIR], [Installation path of OpenLoops]),
	[],
	[with_openloops=no]

	)

	AC_MSG_RESULT([$with_openloops])

	AS_IF([test "x$with_openloops" != "xno"],
	[AC_CHECK_FILES(
	${with_openloops}/lib/libopenloops.so,
	[have_openloops=lib], [have_openloops=no])],
	[have_openloops=no])

	AS_IF([test "x$with_openloops" != "xno" -a "x$have_openloops" = "xno" ],
	[AC_CHECK_FILES(
	+ ${with_openloops}/lib/libopenloops.dylib,
	+ [have_openloops=lib], [have_openloops=no])])
	+
	+AS_IF([test "x$with_openloops" != "xno" -a "x$have_openloops" = "xno" ],
	+ [AC_CHECK_FILES(
	${with_openloops}/lib64/libopenloops.so,
	[have_openloops=lib64], [have_openloops=no])])

	+
	+AS_IF([test "x$with_openloops" != "xno" -a "x$have_openloops" = "xno" ],
	+ [AC_CHECK_FILES(
	+ ${with_openloops}/lib64/libopenloops.dylib,
	+ [have_openloops=lib64], [have_openloops=no])])
	+
	+
	+
	+
	+
	AS_IF([test "x$have_openloops" = "xlib"],
	[OPENLOOPSLIBS=${with_openloops}/lib
	AC_SUBST(OPENLOOPSLIBS)
	])

	AS_IF([test "x$have_openloops" = "xlib64"],
	[OPENLOOPSLIBS=${with_openloops}/lib64
	AC_SUBST(OPENLOOPSLIBS)
	])

	AS_IF([test "x$with_openloops" != "xno" -a "x$have_openloops" = "xno"],
	[AC_MSG_ERROR([OpenLoops requested but not found])])

	AM_CONDITIONAL(HAVE_OPENLOOPS,[test "x$have_openloops" = "xlib" -o "x$have_openloops" = "xlib64"])

	if test "x$have_openloops" = "xlib" -o "x$have_openloops" = "xlib64" ; then
	OPENLOOPSPREFIX=${with_openloops}
	LOAD_OPENLOOPS="library"
	CREATE_OPENLOOPS="create"
	INSERT_OPENLOOPS="insert"
	SET_OPENLOOPS="set"
	MKDIR_OPENLOOPS="mkdir"
	else
	LOAD_OPENLOOPS="# library"
	CREATE_OPENLOOPS="# create"
	INSERT_OPENLOOPS="# insert"
	SET_OPENLOOPS="# set"
	MKDIR_OPENLOOPS="# mkdir"
	fi

	AC_SUBST([OPENLOOPSPREFIX])
	AC_SUBST([LOAD_OPENLOOPS])
	AC_SUBST([CREATE_OPENLOOPS])
	AC_SUBST([INSERT_OPENLOOPS])
	AC_SUBST([SET_OPENLOOPS])
	AC_SUBST([MKDIR_OPENLOOPS])

	])

	#########################################

	dnl ##### madgraph #####
	AC_DEFUN([HERWIG_CHECK_MADGRAPH],
	[
	AC_MSG_CHECKING([for MadGraph])

	AC_ARG_WITH([madgraph],
	AS_HELP_STRING([--with-madgraph=DIR], [Installation path of MadGraph]),
	[],
	[with_madgraph=no]
	)

	AC_MSG_RESULT([$with_madgraph])

	AS_IF([test "x$with_madgraph" != "xno"],
	[AC_CHECK_FILES(
	${with_madgraph}/bin/mg5_aMC,
	[have_madgraph=yes], [have_madgraph=no])],
	[have_madgraph=no])

	AS_IF([test "x$have_madgraph" = "xyes"],
	[MADGRAPHPREFIX=${with_madgraph}
	AC_SUBST(MADGRAPHPREFIX)
	])

	AS_IF([test "x$with_madgraph" != "xno" -a "x$have_madgraph" = "xno"],
	[AC_MSG_ERROR([MadGraph requested but not found])])

	AM_CONDITIONAL(HAVE_MADGRAPH,[test "x$have_madgraph" = "xyes" ])

	if test "x$have_madgraph" = "xyes" ; then
	LOAD_MADGRAPH="library"
	CREATE_MADGRAPH="create"
	INSERT_MADGRAPH="insert"
	SET_MADGRAPH="set"
	else
	LOAD_MADGRAPH="# library"
	CREATE_MADGRAPH="# create"
	INSERT_MADGRAPH="# insert"
	SET_MADGRAPH="# set"
	fi

	AC_SUBST([LOAD_MADGRAPH])
	AC_SUBST([CREATE_MADGRAPH])
	AC_SUBST([INSERT_MADGRAPH])
	AC_SUBST([SET_MADGRAPH])

	])



	dnl ##### hej #####
	AC_DEFUN([HERWIG_CHECK_HEJ],
	[
	AC_MSG_CHECKING([for hej])

	AC_ARG_WITH([hej],
	AS_HELP_STRING([--with-hej=DIR], [Installation path of HEJ]),
	[],
	[with_hej=no]

	)

	AC_MSG_RESULT([$with_hej])

	AS_IF([test "x$with_hej" != "xno"],
	[AC_CHECK_FILES(
	${with_hej}/lib/libJets.so
	${with_hej}/include/HEJ/Jets.h,
	[have_hej=yes], [have_hej=no])],
	[have_hej=no])

	AS_IF([test "x$have_hej" = "xyes"],
	[HEJLIBPATH=${with_hej}/lib
	AC_SUBST(HEJLIBPATH)
	HEJCPPFLAGS=`${with_hej}/bin/hej-config --cppflags`
	AC_SUBST(HEJCPPFLAGS)
	],
	[AS_IF([test "x$with_hej" != "xno"],
	[AC_MSG_ERROR([hej requested but not found])
	])
	])
	AM_CONDITIONAL(HAVE_HEJ,[test "x$have_hej" = "xyes"])
	])

	dnl ##### PDF PATH #####
	AC_DEFUN([HERWIG_PDF_PATH],
	[
	AC_MSG_CHECKING([which Herwig++ PDF data to use])
	AC_ARG_WITH(pdf,
	AC_HELP_STRING([--with-pdf=DIR],[installation path of Herwig++PDF data tarball]),
	[],
	[with_pdf=${prefix}]
	)
	HERWIG_PDF_PREFIX=${with_pdf}/share/Herwig++PDF

	if test -f "${HERWIG_PDF_PREFIX}/mrst/2008/mrstMCal.dat"; then
	AC_MSG_RESULT([$with_pdf])
	localPDFneeded=false
	else
	AC_MSG_RESULT([Using built-in PDF data set. For other data sets, set --with-pdf.])
	HERWIG_PDF_PREFIX=PDF
	localPDFneeded=true
	fi
	HERWIG_PDF_DEFAULT=${HERWIG_PDF_PREFIX}/mrst/2008/mrstMCal.dat
	HERWIG_PDF_NLO=${HERWIG_PDF_PREFIX}/mrst/2002/mrst2002nlo.dat
	HERWIG_PDF_POMERON=${HERWIG_PDF_PREFIX}/diffraction/

	AM_CONDITIONAL(WANT_LOCAL_PDF,[test "x$localPDFneeded" = "xtrue"])
	AC_SUBST(HERWIG_PDF_DEFAULT)
	AC_SUBST(HERWIG_PDF_NLO)
	AC_SUBST(HERWIG_PDF_POMERON)
	])

	dnl ###### GSL ######
	AC_DEFUN([HERWIG_CHECK_GSL],
	[
	AC_MSG_CHECKING([for gsl location])
	GSLINCLUDE=""
	GSLLIBS=""

	AC_ARG_WITH(gsl,
	AC_HELP_STRING([--with-gsl=DIR],[location of gsl installation @<:@default=system libs@:>@]),
	[],
	[with_gsl=system])

	if test "x$with_gsl" = "xno"; then
	AC_MSG_ERROR([libgsl is required. Please install the GNU scientific library and header files.])
	fi

	if test "x$with_gsl" = "xsystem"; then
	AC_MSG_RESULT([in system libraries])
	oldlibs="$LIBS"
	AC_CHECK_LIB(m,main)
	AC_CHECK_LIB(gslcblas,main)
	AC_CHECK_LIB(gsl,main,[],
	[
	AC_MSG_ERROR([Cannot find libgsl. Please install the GNU scientific library and header files or use --with-gsl=.])
	]
	)
	GSLLIBS="$LIBS"
	LIBS=$oldlibs
	else
	if test "`uname -m`" = "x86_64" -a -e "$with_gsl/lib64/libgsl.a" -a -d "$with_gsl/include/gsl"; then
	AC_MSG_RESULT([found in $with_gsl])
	GSLLIBS="-L$with_gsl/lib64 -R$with_gsl/lib64 -lgslcblas -lgsl"
	GSLINCLUDE="-I$with_gsl/include"
	elif test -e "$with_gsl/lib/libgsl.a" -a -d "$with_gsl/include/gsl"; then
	AC_MSG_RESULT([found in $with_gsl])
	GSLLIBS="-L$with_gsl/lib -R$with_gsl/lib -lgslcblas -lgsl"
	GSLINCLUDE="-I$with_gsl/include"
	else
	AC_MSG_RESULT([not found])
	AC_MSG_ERROR([Can't find $with_gsl/lib/libgsl.a or the headers in $with_gsl/include])
	fi
	fi

	AC_SUBST(GSLINCLUDE)
	AC_SUBST(GSLLIBS)
	])

	AC_DEFUN([HERWIG_VERSIONSTRING],
	[
	if test -d $srcdir/.svn; then
	AC_CHECK_PROG(have_svnversion,[svnversion],[yes],[no])
	fi
	AM_CONDITIONAL(USE_SVNVERSION,[test "x$have_svnversion" = "xyes"])
	])

	dnl ##### COMPILERFLAGS #####
	AC_DEFUN([HERWIG_COMPILERFLAGS],
	[
	AC_REQUIRE([HERWIG_CHECK_GSL])
	AC_REQUIRE([HERWIG_CHECK_THEPEG])
	AC_REQUIRE([BOOST_REQUIRE])

	AM_CPPFLAGS="-I\$(top_builddir)/include $THEPEGINCLUDE \$(GSLINCLUDE) \$(BOOST_CPPFLAGS)"

	AC_MSG_CHECKING([for debugging mode])
	AC_ARG_ENABLE(debug,
	AC_HELP_STRING([--enable-debug],[debug mode, use --enable-debug=slow for additional options that slow down the run.]),
	[],
	[enable_debug=no]
	)
	AC_MSG_RESULT([$enable_debug])

	if test "x$enable_debug" = "xno"; then
	debugflags=""
	else
	debugflags="-g"
	fi

	dnl -Wfloat-equal -fvisibility-inlines-hidden -Wctor-dtor-privacy -Weffc++
	if test -n "$GCC"; then
	warnflags="-ansi -pedantic -Wall -Wextra -Wno-overloaded-virtual"

	if test "x$enable_debug" = "xslow"; then
	debugflags="$debugflags -fno-inline"
	AM_CPPFLAGS="$AM_CPPFLAGS -D_GLIBCXX_DEBUG"
	fi
	fi

	dnl do an actual capability check on ld instead of this workaround
	case "${host}" in
	-darwin)
	;;
	*)
	AM_LDFLAGS="-Wl,--enable-new-dtags"
	;;
	esac

	AC_SUBST(AM_CPPFLAGS)
	AC_SUBST(AM_CFLAGS, ["$warnflags $debugflags"])
	AC_SUBST(AM_CXXFLAGS,["$warnflags $debugflags"])
	AC_SUBST(AM_FCFLAGS, ["$debugflags"])
	AC_SUBST(AM_LDFLAGS)
	])

	AC_DEFUN([HERWIG_ENABLE_MODELS],
	[
	AC_MSG_CHECKING([if BSM models should be built])

	AC_ARG_ENABLE(models,
	AC_HELP_STRING([--disable-models],[Turn off compilation of BSM models.]),
	[],
	[enable_models=yes]
	)
	AC_MSG_RESULT([$enable_models])

	LOAD_BSM=""
	if test "$enable_models" = "yes"; then
	LOAD_BSM="read BSMlibs.in"
	fi
	AC_SUBST(LOAD_BSM)

	AM_CONDITIONAL(WANT_BSM,[test "$enable_models" = "yes"])
	])

	AC_DEFUN([HERWIG_OVERVIEW],
	[
	FCSTRING=`$FC --version \| head -1`
	CXXSTRING=`$CXX --version \| head -1`
	CCSTRING=`$CC --version \| head -1`
	if test "x$PYTHON" != "x:"
	then
	python_was_found="yes, using Python $PYTHON_VERSION"
	else
	python_was_found="no, requires Python >= 2.6"
	fi
	cat << _HW_EOF_ > config.herwig
	*****************************************************
	*** $PACKAGE_STRING configuration summary
	*** Please include this information in bug reports!
	***--------------------------------------------------
	*** Prefix: $prefix
	***
	*** BSM models: $enable_models
	*** UFO converter: ${python_was_found}
	***
	*** Herwig debug mode: $enable_debug
	***
	*** ThePEG: $with_thepeg
	*** ThePEG headers: $with_thepeg_headers
	***
	*** VBFNLO: $with_vbfnlo
	*** nlojet: $with_nlojet
	*** njet: $with_njet
	*** GoSam: $with_gosam
	*** GoSam-Contrib: $with_gosam_contrib
	*** OpenLoops: $with_openloops
	*** MadGraph: $with_madgraph
	*** HEJ: $with_hej
	***
	*** GSL: $with_gsl
	*** boost: ${BOOST_CPPFLAGS:-system}
	*** Fastjet: ${fjconfig}
	***
	*** Host: $host
	*** CC: $CCSTRING
	*** CXX: $CXXSTRING
	*** FC: $FCSTRING
	***
	*** CXXFLAGS: $CXXFLAGS
	*****************************************************
	_HW_EOF_
	])

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 6:36 PM (41 m, 19 s)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023834
Default Alt Text: (103 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions