No OneTemporary
Actions

Size

28 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc b/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc
	--- a/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc
	+++ b/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc
	@@ -1,514 +1,565 @@
	// -- C++ --
	//
	// OpenLoopsAmplitude.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the OpenLoopsAmplitude class.
	//

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

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

	#include "ThePEG/Utilities/DynamicLoader.h"

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

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

	using namespace Herwig;

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

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

	OpenLoopsAmplitude::OpenLoopsAmplitude() :
	- theHiggsEff(false), use_cms(true), psp_tolerance(12){
	+ theHiggsEff(false), use_cms(true), psp_tolerance(12), theDiagonal(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();
	}

	vector<int> OpenLoopsAmplitude::idpair=vector<int>();

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



	string tempcontract=contract;


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

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

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

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

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

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

	didStartOLP() = true;

	}

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

	string orderFileName =
	optionalContractFile().empty() ?
	factory()->buildStorage() + name() + ".OLPContract.lh" :
	optionalContractFile();


	ofstream orderFile(orderFileName.c_str());
	size_t asPower = 100;
	size_t minlegs = 100;
	size_t maxlegs = 0;

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

	orderFile << "# OLP order file created by Herwig/Matchbox for OpenLoops\n\n";
	orderFile << "CorrectionType QCD\n";
	orderFile << "IRregularization " << (isDR() ? "DRED" : "CDR") << "\n";
	orderFile << "extra answerfile " << (factory()->buildStorage() + name() + ".OLPAnswer.lh") << "\n";
	orderFile << "extra psp_tolerance "<<psp_tolerance<<"\n";
	orderFile << "extra use_cms "<<(use_cms?"1":"0")<< "\n";
	+
	+ // To enable non-diagonal CKM processes
	+ if (!theDiagonal) {
	+
	+ orderFile << "model sm_ckm\n";
	+ orderFile << "extra ckmorder 1\n";
	+
	+ std::vector<std::vector<std::complex<double>>> ckmelement = StandardCKM().getUnsquaredMatrix(6);
	+
	+ // write CKM matrix elements to OLP order file
	+ orderFile << "extra vckmdu " << ckmelement[0][0].real() << "\n";
	+ orderFile << "extra vckmsu " << ckmelement[0][1].real() << "\n";
	+ orderFile << "extra vckmbu " << ckmelement[0][2].real() << "\n";
	+
	+ orderFile << "extra vckmidu " << ckmelement[0][0].imag() << "\n";
	+ orderFile << "extra vckmisu " << ckmelement[0][1].imag() << "\n";
	+ orderFile << "extra vckmibu " << ckmelement[0][2].imag() << "\n";
	+
	+ orderFile << "extra vckmdc " << ckmelement[1][0].real() << "\n";
	+ orderFile << "extra vckmsc " << ckmelement[1][1].real() << "\n";
	+ orderFile << "extra vckmbc " << ckmelement[1][2].real() << "\n";
	+
	+ orderFile << "extra vckmidc " << ckmelement[1][0].imag() << "\n";
	+ orderFile << "extra vckmisc " << ckmelement[1][1].imag() << "\n";
	+ orderFile << "extra vckmibc " << ckmelement[1][2].imag() << "\n";
	+
	+ orderFile << "extra vckmdt " << ckmelement[2][0].real() << "\n";
	+ orderFile << "extra vckmst " << ckmelement[2][1].real() << "\n";
	+ orderFile << "extra vckmbt " << ckmelement[2][2].real() << "\n";
	+
	+ orderFile << "extra vckmidt " << ckmelement[2][0].imag() << "\n";
	+ orderFile << "extra vckmist " << ckmelement[2][1].imag() << "\n";
	+ orderFile << "extra vckmibt " << ckmelement[2][2].imag() << "\n";
	+
	+ }
	+
	if (theCollierLib) {
	orderFile << "extra preset 2 "<<"\n";
	if(theHiggsEff){
	orderFile << "extra stability_mode 14\n";
	orderFile << "extra redlib1 1\n";
	}
	}
	if (theHiggsEff){
	orderFile << "model heft\n";
	}
	orderFile << "\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);
	}
	}
	}
	}
	}

	idpair.clear();
	for (size_t i=0;i<processmap.size();i++)idpair.push_back(-1);
	idpair.push_back(-1);
	for ( map<int, OpenLoopsProcInfo>::iterator p = processmap.begin() ; p != processmap.end() ; p++ ) {
	idpair[(p).second.HID()]=(p).second.GID();
	if ( (*p).second.GID() == -1 ) return 0;
	}
	return 1;
	}


	bool OpenLoopsAmplitude::startOLP(const map<pair<Process, int>, int>& procs) {
	string contractFileName = factory()->buildStorage() + name() + ".OLPAnswer.lh";

	string orderFileName =
	optionalContractFile().empty() ?
	factory()->buildStorage() + name() + ".OLPContract.lh" :
	optionalContractFile();

	fillOrderFile(procs);
	int status = -1;
	startOLP(orderFileName, status);
	if ( !checkOLPContract() ) {
	return false;
	}
	if ( status != 1 ) return false;
	return true;
	}


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

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

	assert ( idpair.size() != 0 );

	OLP_EvalSubProcess2(&idpair[id], olpMomenta(), &scale, out,&acc );


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

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

	colourCorrelatorResults.resize(n * (n - 1) / 2);
	assert ( idpair.size() != 0 );


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

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

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


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

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

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

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

	assert ( idpair.size() != 0 ) ;
	int id =idpair[olpId()[ProcessType::spinColourCorrelatedME2]];
	//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;

	}





	string OpenLoopsAmplitude::OpenLoopsLibs_=OPENLOOPSLIBS;
	string OpenLoopsAmplitude::OpenLoopsPrefix_=OPENLOOPSPREFIX;

	void OpenLoopsAmplitude::setOpenLoopsLibs(string p){
	OpenLoopsLibs_=p;
	}
	string OpenLoopsAmplitude::getOpenLoopsLibs() const{
	return OpenLoopsLibs_;
	}


	void OpenLoopsAmplitude::setOpenLoopsPrefix(string p){
	OpenLoopsPrefix_=p;
	}
	string OpenLoopsAmplitude::getOpenLoopsPrefix() const{
	return OpenLoopsPrefix_;
	}


	// 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 << theHiggsEff << use_cms << theCollierLib << OpenLoopsLibs_ << OpenLoopsPrefix_;
	+ os << idpair << theHiggsEff << use_cms << theDiagonal << theCollierLib << OpenLoopsLibs_ << OpenLoopsPrefix_;
	OpenLoopsLibs_.clear();
	OpenLoopsPrefix_.clear();
	}

	void OpenLoopsAmplitude::persistentInput(PersistentIStream & is, int) {
	- is >> idpair >> theHiggsEff >> use_cms >> theCollierLib ;
	+ is >> idpair >> theHiggsEff >> use_cms >> theDiagonal >> theCollierLib ;
	string input=""; is>>input; if (!input.empty())OpenLoopsLibs_=input;
	input=""; is>>input; if (!input.empty())OpenLoopsPrefix_=input;
	}

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

	void OpenLoopsAmplitude::Init() {

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

	static Switch<OpenLoopsAmplitude,bool> interfaceHiggsEff
	("HiggsEff",
	"Switch On/Off for effective higgs model.",
	&OpenLoopsAmplitude::theHiggsEff, false, false, false);
	static SwitchOption interfaceHiggsEffYes
	(interfaceHiggsEff,
	"Yes",
	"Yes",
	true);
	static SwitchOption interfaceHiggsEffNo
	(interfaceHiggsEff,
	"No",
	"No",
	false);

	+ static Switch<OpenLoopsAmplitude,bool> interfaceDiagonal
	+ ("Diagonal",
	+ "Use a diagonal CKM matrix (ignoring the CKM object of the StandardModel).",
	+ &OpenLoopsAmplitude::theDiagonal, false, false, false);
	+ static SwitchOption interfaceDiagonalYes
	+ (interfaceDiagonal,
	+ "Yes",
	+ "Use a diagonal CKM matrix.",
	+ true);
	+ static SwitchOption interfaceDiagonalNo
	+ (interfaceDiagonal,
	+ "No",
	+ "Use the CKM object as used by the StandardModel.",
	+ false);
	+
	static Switch<OpenLoopsAmplitude,bool> interfaceUseComplMass
	("ComplexMassScheme",
	"Switch on or off if Complex Masses.",
	&OpenLoopsAmplitude::use_cms, true, false, false);
	static SwitchOption interfaceUseComplMassYes
	(interfaceUseComplMass,
	"Yes",
	"True for Complex Masses.",
	true);
	static SwitchOption interfaceUseComplMassNo
	(interfaceUseComplMass,
	"No",
	"False for no Complex Masses.",
	false);

	static Switch<OpenLoopsAmplitude,bool> interfaceCollier
	("UseCollier",
	"Switch On/Off for using the Collier Lib (arXiv:1604.06792).",
	&OpenLoopsAmplitude::theCollierLib, true, false, false);
	static SwitchOption interfaceCollierYes
	(interfaceCollier,
	"Yes",
	"Yes",
	true);
	static SwitchOption interfaceCollierNo
	(interfaceCollier,
	"No",
	"No",
	false);


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

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

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

	}

	diff --git a/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.h b/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.h
	--- a/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.h
	+++ b/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.h
	@@ -1,350 +1,358 @@
	// -- C++ --
	//
	// OpenLoopsAmplitude.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef Herwig_OpenLoopsAmplitude_H
	#define Herwig_OpenLoopsAmplitude_H
	//
	// This is the declaration of the OpenLoopsAmplitude class.
	//

	#include "Herwig/MatrixElement/Matchbox/Base/MatchboxOLPME.h"
	+#include "Herwig/Models/StandardModel/StandardCKM.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup Matchbox
	* \author Johannes Bellm, Simon Platzer
	*
	* \brief Process information for OpenLoops
	*/
	class OpenLoopsProcInfo{

	public:

	/**
	* Default constructor
	*/
	OpenLoopsProcInfo() {}

	/**
	* Construct giving data
	*/
	OpenLoopsProcInfo(int HID,int GID, string procstr,string typestr)
	: theHOlpId(HID), theGOlpId(GID), theProcstr(procstr), theTypestr(typestr) {}

	/**
	* Document me
	*/
	int HID() const { return theHOlpId; }

	/**
	* Document me
	*/
	int GID() const { return theGOlpId; }

	/**
	* Document me
	*/
	const string& Pstr() const { return theProcstr; }

	/**
	* Document me
	*/
	const string& Tstr() const { return theTypestr; }

	/**
	* Document me
	*/
	void setGID(int g) { theGOlpId=g; }

	/**
	* Document me
	*/
	void setOAs(int i) { orderAlphas=i; }

	/**
	* Document me
	*/
	int orderAs() { return orderAlphas; }

	private:

	/**
	* Document me
	*/
	int theHOlpId;

	/**
	* Document me
	*/
	int theGOlpId;

	/**
	* Document me
	*/
	string theProcstr;

	/**
	* Document me
	*/
	string theTypestr;

	/**
	* Document me
	*/
	int orderAlphas;

	public:

	/**
	* Write to persistent stream
	*/
	void persistentOutput(PersistentOStream & os) const{
	os << theHOlpId << theGOlpId << theProcstr << theTypestr << orderAlphas;
	}

	/**
	* Read from persistent stream
	*/
	void persistentInput(PersistentIStream &is) {
	is >> theHOlpId >> theGOlpId >> theProcstr >> theTypestr >> orderAlphas;
	}

	};

	/**
	* \ingroup Matchbox
	* \author Johannes Bellm, Simon Platzer
	*
	* \brief OpenLoopsAmplitude implements an interface to OpenLoops
	*/
	class OpenLoopsAmplitude: public MatchboxOLPME {

	public:

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

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

	public:

	virtual void fillOrderFile(const map<pair<Process,int>,int>& procs);


	virtual bool isCS() const { return false; }
	virtual bool isExpanded() const { return true; }
	virtual bool isBDK() const { return false; }
	//virtual bool isDR() const { return true; }
	/**
	* Start the one loop provider, if appropriate, giving order and
	* contract files
	*/

	virtual bool checkOLPContract();

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

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

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


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

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

	/**
	* Call OLP_EvalSubProcess and fill in the results
	*/
	void evalSubProcess() const;

	/**
	* Fill in results for the given colour correlator
	*/
	virtual void evalColourCorrelator(pair<int,int> ij) const;


	/**
	* Fill in results for the given colour/spin correlator
	*/
	virtual void evalSpinColourCorrelator(pair<int,int> ij) const;


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


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

	virtual void doinitrun();


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


	private:

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

	/**
	* Store colour correlator results
	*/
	mutable vector<double> colourCorrelatorResults;

	/**
	* Store spin colour correlator results
	*/
	mutable vector<double> spinColourCorrelatorResults;


	/**
	* first is the olp id from herwig, second the answer from openloops
	*/
	static vector< int > idpair;


	/**
	* Helper map to store information in different procs.
	*/

	map<int , OpenLoopsProcInfo > processmap;


	/**
	* Interface for Higgs Effective
	*/
	bool theHiggsEff;

	/**
	* Complex Mass Scheme.
	*/

	bool use_cms;


	/**
	* Use of Collier Lib (arXiv:1604.06792), available since OpenLoops 1.3.0.
	*/
	bool theCollierLib=true;

	/**
	* parameter to set Phase space tolerance for massiv particles.
	* Should not be used. Better: set Openloops:Massless 11
	*/

	int psp_tolerance;

	/**
	+ * True, if a diagonal CKM matrix should be assumed. This ignores
	+ * the CKM object of the StandardModel.
	+ */
	+ bool theDiagonal;
	+
	+
	+ /**
	* Location of the OpenLoops libraries
	*/
	static string OpenLoopsLibs_;

	/**
	* Location of the OpenLoops
	*/
	static string OpenLoopsPrefix_;


	/**
	* Helper functions to make long strings static
	*/

	void setOpenLoopsLibs(string p);
	string getOpenLoopsLibs() const;

	void setOpenLoopsPrefix(string p);
	string getOpenLoopsPrefix() const;




	};

	}

	#endif /* Herwig_OpenLoopsAmplitude_H */
	diff --git a/src/Matchbox/DiagonalCKM.in b/src/Matchbox/DiagonalCKM.in
	--- a/src/Matchbox/DiagonalCKM.in
	+++ b/src/Matchbox/DiagonalCKM.in
	@@ -1,16 +1,16 @@
	# -- ThePEG-repository --

	set /Herwig/MatrixElements/Matchbox/Factory:QuarkFlavourDiagonal Yes
	set /Herwig/Merging/MergingFactory:QuarkFlavourDiagonal Yes

	set /Herwig/Vertices/FFWMatchboxVertex:Diagonal Yes


	cd /Herwig/MatrixElements/Matchbox/Amplitudes/Builtin
	set Amplitudelnuqqbar:Diagonal Yes
	set Amplitudelnuqqbarg:Diagonal Yes
	set Amplitudelnuqqbargg:Diagonal Yes
	set Amplitudelnuqqbarqqbar:Diagonal Yes


	-
	+set /Herwig/MatrixElements/Matchbox/Amplitudes/OpenLoops:Diagonal Yes
	diff --git a/src/Matchbox/NonDiagonalCKM.in b/src/Matchbox/NonDiagonalCKM.in
	--- a/src/Matchbox/NonDiagonalCKM.in
	+++ b/src/Matchbox/NonDiagonalCKM.in
	@@ -1,7 +1,7 @@
	# -- ThePEG-repository --

	set /Herwig/MatrixElements/Matchbox/Factory:QuarkFlavourDiagonal No
	-set /Herwig/MatrixElements/Matchbox/Factory:QuarkFlavourDiagonal No
	+set /Herwig/Merging/MergingFactory:QuarkFlavourDiagonal No

	set /Herwig/Vertices/FFWMatchboxVertex:Diagonal No

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 6:01 AM (16 h, 6 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4957692
Default Alt Text: (28 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions