GoSamAmplitude.cc.in
No OneTemporary
Actions

Size

32 KB

Subscribers

None

GoSamAmplitude.cc.in
View Options

	// -- 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;thePrintParameter=true;theFormOpt=true;theNinja=true;
	}

	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() {
	// MatchboxOLPME::doinitrun();
	// }

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

	// string cmd = "";
	// char char_cwd[256];
	// getcwd(char_cwd, sizeof(char_cwd));
	// string cwd = string(char_cwd);
	//
	// // set all necessary path and file names
	// gosamPath = gosamPathInterface == "" ? cwd + "/GoSam" : gosamPathInterface;
	// gosamSourcePath = gosamPath + "/source";
	// gosamInstallPath = gosamPath + "/build";
	//
	// // create all the directories
	// std::system(("mkdir " + gosamPath).c_str());
	// std::system(("mkdir " + gosamSourcePath).c_str());
	// std::system(("mkdir " + gosamInstallPath).c_str());
	//
	// contractFileTitle = name() + ".OLPContract.lh";
	// parametersFileTitle = name() + ".OLPParameters.lh";
	// string orderFileName = gosamPath + "/" + name() + ".OLPOrder.lh";
	// contractFileName = gosamPath + "/" + contractFileTitle;
	// gosamSetupInFileName = gosamSetupInFilenameInterface == "" ? gosamPath + "/setup.gosam.in" : gosamSetupInFilenameInterface;

	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/";

	// cout << "GoSamAmplitude::startOLP(const map<pair<Process, int>, int>& procs): gosamPath = " << gosamPath << endl;

	// 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";

	//parametersFileTitle = name() + ".OLPParameters.lh";

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

	// check if setup.gosam.in exists
	ifstream infile(gosamSetupInFileName.c_str());
	if (!infile) {
	// if it doesn't, then get it from src/defaults
	string cmd = "cp @prefix@/share/Herwig++/defaults/setup.gosam.in " + gosamSetupInFileName;
	std::system(cmd.c_str());
	// set form's tempdir to the current working directory
	cmd = "sed -i 's/@FORMTEMPDIR@/" + StringUtils::replace(gosamSourcePath, string("/"), string("\\/")) + "/g' " + gosamSetupInFileName;
	std::system(cmd.c_str());
	// set reduction programs to what is specified through the repository
	cmd = "sed -i 's/@REDUCTIONPROGRAMS@/" + (theNinja ? string("ninja, golem95") : string("samurai, golem95")) + "/g' " + gosamSetupInFileName;
	std::system(cmd.c_str());
	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 = "sed -i 's/@FORMOPT@/" + (theFormOpt ? string("") : string(", noformopt")) + "/g' " + gosamSetupInFileName;
	std::system(cmd.c_str());
	}
	infile.close();

	// 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 (gosamPathInterface=="") {
	// if ( !( DynamicLoader::load(gosamPath+"/build/lib/libgolem_olp.so") \|\| DynamicLoader::load(gosamPath+"/build/lib64/libgolem_olp.so") ) ) assert(0);
	// buildGoSam();
	// }
	// else {
	// if ( !( DynamicLoader::load(gosamPathInterface+"/build/lib/libgolem_olp.so") \|\| DynamicLoader::load(gosamPathInterface+"/build/lib64/libgolem_olp.so") ) ) assert(0);
	// 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("/");

	// cout << "GoSamAmplitude::startOLP(const string& contract, int& status): gosamPath = " << gosamPath << endl;

	if (gosamPathInterface==""){
	// if ( !( DynamicLoader::load("GoSam/build/lib/libgolem_olp.so") \|\| DynamicLoader::load("GoSam/build/lib64/libgolem_olp.so") ) ) assert(0);
	// tempcontract = "GoSam/" + tempcontract;
	// if ( !( DynamicLoader::load(gosamPath+"/build/lib/libgolem_olp.so") \|\| DynamicLoader::load(gosamPath+"/build/lib64/libgolem_olp.so") ) ) assert(0);
	// tempcontract = gosamPath + "/" + tempcontract;
	if ( !( DynamicLoader::load(gosamPath+"build/lib/libgolem_olp.so") \|\| DynamicLoader::load(gosamPath+"build/lib64/libgolem_olp.so") ) ) assert(0);
	tempcontract = gosamPath + tempcontract;
	}
	else {
	if ( !( DynamicLoader::load(gosamPathInterface+"/build/lib/libgolem_olp.so") \|\| DynamicLoader::load(gosamPathInterface+"/build/lib64/libgolem_olp.so") ) ) assert(0);
	tempcontract = gosamPathInterface + "/" + 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);

	// hand over quark masses and widths (iff massive)
	string mstr;
	string wstr;
	if (massiveParticles.empty())
	for (int i=1; i<=6; ++i)
	if (getParticleData(i)->mass()/GeV > 0.0) massiveParticles.push_back(i);
	if ( massiveParticles.size() != 0 ) {
	for ( vector<int>::const_iterator mID = massiveParticles.begin(); mID != massiveParticles.end(); ++mID ) {
	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;
	}
	}

	// 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 (thePrintParameter) {
	// char char_cwd[256];
	// getcwd(char_cwd, sizeof(char_cwd));
	// gosamPath = gosamPathInterface == "" ? string(char_cwd) + "/GoSam" : gosamPathInterface;
	string ppstr = gosamPath + "/" + name() + ".OLPParameters.lh";
	char * ppchar = new char[ppstr.size()+1];
	std::copy(ppstr.begin(),ppstr.end(),ppchar);
	ppchar[ppstr.size()] = '\0';
	OLP_PrintParameter(ppchar);
	delete[] ppchar;
	}

	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){
	// generator()->log() << "\n>>> generating GoSam amplitudes. This may take some time, please be patient.\n"
	// << ">>> see gosam-amplitudes.log for details.\n" << flush;
	// string cmd = "@GOSAMPREFIX@/bin/gosam.py --olp -o " + contract + " --config=" +
	// gosamSetupInFileName + " --destination=" + gosamSourcePath + " " + order + " > gosam-amplitudes.log 2>&1";
	// std::system(cmd.c_str());

	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 = "ln -s " + contract + " " + 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 ) assert(0);
	string subx = sub.substr(3);
	int subint;
	istringstream(subx) >> subint;
	(*p).second.setGID(subint);
	}
	}
	}
	}

	// string ids = "GoSam_ids.dat";
	// ofstream IDS(ids.c_str());

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

	for ( map<int, gosamprocinfo>::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() << " " << (*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;
	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 = "GoSam_ids.dat";
	// string line = factory()->runStorage() + "GoSam.ids.dat";
	string line = factory()->buildStorage() + "GoSam.ids.dat";
	ifstream infile(line.c_str());
	int hid;
	int gid;
	string type;
	while (std::getline(infile, line)) {
	istringstream(line) >> hid >> gid >> type;
	idpair.insert ( std::pair<int,int>(hid,gid) );
	idtypepair.insert ( std::pair<int,string>(hid,type) );
	}
	}


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


	void GoSamAmplitude::evalSubProcess() const {

	useMe();

	double units = pow(lastSHat() / GeV2, 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.find(id)->second); // If id denotes the Herwig ID, this returns the GoSam ID
	string calltype(idtypepair.find(id)->second); // 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, 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.find(olpId()[ProcessType::colourCorrelatedME2]))).second);

	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, 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.find(olpId()[ProcessType::spinColourCorrelatedME2]))).second);

	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 << gosamPathRelative
	<< gosamSourcePathRelative << gosamInstallPathRelative << gosamSetupInFileName
	<< orderFileTitle << contractFileTitle << parametersFileTitle
	<< contractFileName << orderFileName << parametersFileName
	<< theCodeExists << theFormOpt << theNinja << isitDR << massiveParticles
	<< thePrintParameter;
	}

	void GoSamAmplitude::persistentInput(PersistentIStream & is, int) {
	is >> idpair >> idtypepair >> processmap >> gosamPathInterface
	>> gosamSetupInFilenameInterface >> gosamBuildScript >> gosamPath
	>> gosamSourcePath >> gosamInstallPath >> gosamPathRelative
	>> gosamSourcePathRelative >> gosamInstallPathRelative >> gosamSetupInFileName
	>> orderFileTitle >> contractFileTitle >> parametersFileTitle
	>> contractFileName >> orderFileName >> parametersFileName
	>> theCodeExists >> theFormOpt >> theNinja >> isitDR >> massiveParticles
	>> thePrintParameter;
	}


	// * 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> interfacePrintParameter
	("PrintParameter",
	"Switch On if OLP parameters are to be printed into an output file"
	"'GoSamAmplitude'.OLPParameters.lh, otherwise Off. Default is Off.",
	&GoSamAmplitude::thePrintParameter, false, false, false);
	static SwitchOption interfacePrintParameterOn
	(interfacePrintParameter,
	"On",
	"On",
	true);
	static SwitchOption interfacePrintParameterOff
	(interfacePrintParameter,
	"Off",
	"Off",
	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 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);

	}

File Metadata

Mime Type: text/x-c
Expires: Sat, May 3, 6:55 AM (18 h, 13 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4983195
Default Alt Text: GoSamAmplitude.cc.in (32 KB)

GoSamAmplitude.cc.inNo OneTemporaryActions

GoSamAmplitude.cc.inView Options

File Metadata

Event Timeline

GoSamAmplitude.cc.in
No OneTemporary
Actions

GoSamAmplitude.cc.in
View Options