Page MenuHomeHEPForge
Files F10881392

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Matchbox/External/Makefile.am b/MatrixElement/Matchbox/External/Makefile.am
--- a/MatrixElement/Matchbox/External/Makefile.am
+++ b/MatrixElement/Matchbox/External/Makefile.am
@@ -1,96 +1,100 @@
SUBDIRS = BLHAGeneric VBFNLO NJet GoSam OpenLoops MadGraph
pkglib_LTLIBRARIES =
##############
if HAVE_GOSAM
pkglib_LTLIBRARIES += HwMatchboxGoSam.la
endif
HwMatchboxGoSam_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 15:0:0
HwMatchboxGoSam_la_CPPFLAGS = $(AM_CPPFLAGS) \
-DHERWIG_BINDIR="\"$(bindir)\"" \
-DHERWIG_PKGDATADIR="\"$(pkgdatadir)\"" \
-DGOSAM_PREFIX="\"$(GOSAMPREFIX)\""
HwMatchboxGoSam_la_SOURCES = \
GoSam/GoSamAmplitude.cc
###############
if HAVE_VBFNLO
pkglib_LTLIBRARIES += HwMatchboxVBFNLO.la
endif
HwMatchboxVBFNLO_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 6:0:0
HwMatchboxVBFNLO_la_CPPFLAGS = $(AM_CPPFLAGS)
HwMatchboxVBFNLO_la_CPPFLAGS += -I$(VBFNLOINCLUDE)
HwMatchboxVBFNLO_la_CPPFLAGS += -DVBFNLOLIB=$(VBFNLOLIB)
+if HAVE_VBFNLO3
+HwMatchboxVBFNLO_la_CPPFLAGS += -DVBFNLO3
+endif
+
HwMatchboxVBFNLO_la_SOURCES = \
VBFNLO/VBFNLOAmplitude.cc \
VBFNLO/VBFNLOPhasespace.cc
###############
if HAVE_OPENLOOPS
pkglib_LTLIBRARIES += HwMatchboxOpenLoops.la
endif
HwMatchboxOpenLoops_la_SOURCES = \
OpenLoops/OpenLoopsAmplitude.cc
HwMatchboxOpenLoops_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 15:0:0
HwMatchboxOpenLoops_la_CPPFLAGS = $(AM_CPPFLAGS) \
-DOPENLOOPSLIBS="\"$(OPENLOOPSLIBS)\"" \
-DOPENLOOPSPREFIX="\"$(OPENLOOPSPREFIX)\""
##############
if HAVE_NJET
pkglib_LTLIBRARIES += HwMatchboxNJet.la
endif
HwMatchboxNJet_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 15:0:0
HwMatchboxNJet_la_CPPFLAGS = $(AM_CPPFLAGS) -I$(NJETINCLUDEPATH) \
-DNJET_PREFIX="\"$(NJETPREFIX)\"" \
-DNJET_LIBS="\"$(NJETLIBPATH)\"" \
-DNJET_VERSION="$(NJET_VERSION)"
HwMatchboxNJet_la_SOURCES = \
NJet/NJetsAmplitude.cc
##############
if HAVE_MADGRAPH
pkglib_LTLIBRARIES += HwMatchboxMadGraph.la
endif
HwMatchboxMadGraph_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 6:0:0
HwMatchboxMadGraph_la_SOURCES = \
MadGraph/MadGraphAmplitude.cc
HwMatchboxMadGraph_la_CPPFLAGS = $(AM_CPPFLAGS) \
-DHERWIG_BINDIR="\"$(bindir)\"" \
-DHERWIG_INCLUDEDIR="\"$(includedir)\"" \
-DHERWIG_PKGDATADIR="\"$(pkgdatadir)\"" \
-DMADGRAPH_PREFIX="\"$(MADGRAPHPREFIX)\""
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/MatrixElement/Matchbox/External/VBFNLO/VBFNLOAmplitude.cc b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOAmplitude.cc
--- a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOAmplitude.cc
+++ b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOAmplitude.cc
@@ -1,495 +1,536 @@
// -*- C++ -*-
//
// VBFNLOAmplitude.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 VBFNLOAmplitude class.
//
#include "VBFNLOAmplitude.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/DynamicLoader.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include <cstdlib>
#include "VBFNLO/utilities/BLHAinterface.h"
#define DEFSTR(s) CPPSTR(s)
#define CPPSTR(s) #s
using namespace Herwig;
VBFNLOAmplitude::VBFNLOAmplitude()
: theRanHelSum(false), theAnomCoupl(false), VBFNLOlib_(DEFSTR(VBFNLOLIB))
{}
VBFNLOAmplitude::~VBFNLOAmplitude() {}
IBPtr VBFNLOAmplitude::clone() const {
return new_ptr(*this);
}
IBPtr VBFNLOAmplitude::fullclone() const {
return new_ptr(*this);
}
void VBFNLOAmplitude::signOLP(const string& order, const string& contract) {
int status = 0;
+ #ifdef VBFNLO3
+ OLP_Order_VBFNLO(const_cast<char*>(order.c_str()),
+ const_cast<char*>(contract.c_str()),&status);
+ #else
OLP_Order(const_cast<char*>(order.c_str()),
const_cast<char*>(contract.c_str()),&status);
+ #endif
if ( status != 1 )
throw Exception() << "VBFNLOAmplitude: Failed to sign contract with VBFNLO.\n"
<< "The BLHA contract file " << contract << "\n"
<< "may contain further details about the error."
<< Exception::runerror;
}
void VBFNLOAmplitude::setOLPParameter(const string& name, double value) const {
int pStatus = 0;
double zero = 0.0;
+ #ifdef VBFNLO3
+ OLP_SetParameter_VBFNLO(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
+ #else
OLP_SetParameter(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
+ #endif
if ( !pStatus )
throw Exception() << "VBFNLOAmplitude: VBFNLO failed to set parameter '"
<< name << "' to " << value << "\n"
<< Exception::runerror;
}
void VBFNLOAmplitude::startOLP(const string& contract, int& status) {
+ #ifdef VBFNLO3
+ OLP_Start_VBFNLO(const_cast<char*>(contract.c_str()), &status);
+ #else
OLP_Start(const_cast<char*>(contract.c_str()), &status);
+ #endif
map<long,Energy>::const_iterator it=reshuffleMasses().find(ParticleID::b);
double bmass;
if(it==reshuffleMasses().end())
bmass = getParticleData(ParticleID::b)->hardProcessMass()/GeV;
else
bmass = it->second/GeV;
+
setOLPParameter("mass(5)",bmass);
setOLPParameter("mass(6)",getParticleData(ParticleID::t)->hardProcessMass()/GeV);
setOLPParameter("mass(23)",getParticleData(ParticleID::Z0)->hardProcessMass()/GeV);
setOLPParameter("mass(24)",getParticleData(ParticleID::Wplus)->hardProcessMass()/GeV);
setOLPParameter("mass(25)",getParticleData(ParticleID::h0)->hardProcessMass()/GeV);
setOLPParameter("width(23)",getParticleData(ParticleID::Z0)->hardProcessWidth()/GeV);
setOLPParameter("width(24)",getParticleData(ParticleID::Wplus)->hardProcessWidth()/GeV);
setOLPParameter("width(25)",getParticleData(ParticleID::h0)->hardProcessWidth()/GeV);
setOLPParameter("alpha",SM().alphaEMMZ());
setOLPParameter("sw2",SM().sin2ThetaW());
setOLPParameter("Gf",SM().fermiConstant()*GeV2);
setOLPParameter("Nf",factory()->nLight());
setOLPParameter("alphas",SM().alphaS());
setOLPParameter("ranhelsum",theRanHelSum);
setOLPParameter("anomcoupl",theAnomCoupl);
didStartOLP() = true;
}
void VBFNLOAmplitude::loadVBFNLO() {
if ( ! DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.so") ) {
string error1 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.dylib") ) {
string error2 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load("libVBFNLO.so") ) {
string error3 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load("libVBFNLO.dylib") ) {
string error4 = DynamicLoader::lastErrorMessage;
throw Exception() << "VBFNLOAmplitude: failed to load libVBFNLO.so/dylib\n"
<< "Error messages are:\n\n"
<< "* " << VBFNLOlib_ << "/libVBFNLO.so:\n"
<< error1 << "\n"
<< "* " << VBFNLOlib_ << "/libVBFNLO.dylib:\n"
<< error2 << "\n"
<< "* libVBFNLO.so:\n"
<< error3 << "\n"
<< "* libVBFNLO.dylib:\n"
<< error4 << "\n"
<< Exception::runerror;
}
}
}
}
}
bool VBFNLOAmplitude::startOLP(const map<pair<Process,int>,int>& procs) {
loadVBFNLO();
string orderFileName = factory()->buildStorage() + name() + ".OLPOrder.lh";
ofstream orderFile(orderFileName.c_str());
olpOrderFileHeader(orderFile);
// add VBFNLO specifics here
olpOrderFileProcesses(orderFile,procs);
orderFile << flush;
orderFile.close();
string contractFileName = factory()->buildStorage() + name() + ".OLPContract.lh";
signOLP(orderFileName, contractFileName);
int status = -1;
-
- startOLP(contractFileName,status);
+
+ startOLP(contractFileName,status);
if ( status != 1 )
return false;
return true;
}
LorentzVector<Complex> VBFNLOAmplitude::plusPolarization(const Lorentz5Momentum& p,
const Lorentz5Momentum& n,
int inc) const {
// shamelessly stolen from the GoSam interface; mind that we can
// always cast eq (5.7) in the manual into a form that it only uses
// <M-||M_+> and then switch bvetween eps_+ for an outgoing and
// eps_- for an incoming gluon.
double pvec[4] = {p.t()/GeV,p.x()/GeV,p.y()/GeV,p.z()/GeV};
double nvec[4] = {n.t()/GeV,n.x()/GeV,n.y()/GeV,n.z()/GeV};
double out[8] ={ };
+ #ifdef VBFNLO3
+ OLP_Polvec_VBFNLO(pvec,nvec,out);
+ #else
OLP_Polvec(pvec,nvec,out);
+ #endif
LorentzVector<Complex> res;
Complex a(out[0],out[1]);
res.setT(a);
Complex b(out[2],out[3]);
res.setX(b);
Complex c(out[4],out[5]);
res.setY(c);
Complex d(out[6],out[7]);
res.setZ(d);
if (inc<2)
return res.conjugate();
else
return res;
}
void VBFNLOAmplitude::evalSubProcess() const {
useMe();
double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double scale = sqrt(mu2()/GeV2);
- if (hasRunningAlphaS()) setOLPParameter("alphas",lastAlphaS());
-
+ if (hasRunningAlphaS())
+ setOLPParameter("alphas",lastAlphaS());
+
double acc = -1.0;
double out[4]={};
int id =
olpId()[ProcessType::oneLoopInterference] ?
olpId()[ProcessType::oneLoopInterference] :
olpId()[ProcessType::treeME2];
if (theRanHelSum) {
vector<double> helicityrn;
if ( lastHeadMatchboxXComb() ) {
helicityrn = lastHeadMatchboxXComb()->amplitudeRandomNumbers();
} else {
helicityrn = amplitudeRandomNumbers();
}
if (helicityrn.size()>0) {
setOLPParameter("HelicityRN",helicityrn[0]);
}
}
+ #ifdef VBFNLO3
+ OLP_EvalSubProcess2_VBFNLO(&id, olpMomenta(), &scale, out, &acc);
+ #else
OLP_EvalSubProcess2(&id, olpMomenta(), &scale, out, &acc);
+ #endif
if ( olpId()[ProcessType::oneLoopInterference] ) {
lastTreeME2(out[3]*units);
lastOneLoopInterference(out[2]*units);
lastOneLoopPoles(pair<double,double>(out[0]*units,out[1]*units));
} else if ( olpId()[ProcessType::treeME2] ) {
lastTreeME2(out[0]*units);
} else assert(false);
}
void VBFNLOAmplitude::evalColourCorrelator(pair<int,int>) const {
double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double scale = sqrt(mu2()/GeV2);
- if (hasRunningAlphaS()) setOLPParameter("alphas",lastAlphaS());
+ if (hasRunningAlphaS())
+ setOLPParameter("alphas",lastAlphaS());
double acc = -1.0;
int n = lastXComb().meMomenta().size();
colourCorrelatorResults.resize(n*(n-1)/2);
int id = olpId()[ProcessType::colourCorrelatedME2];
if (theRanHelSum) {
vector<double> helicityrn;
if ( lastHeadMatchboxXComb() ) {
helicityrn = lastHeadMatchboxXComb()->amplitudeRandomNumbers();
} else {
helicityrn = amplitudeRandomNumbers();
}
if (helicityrn.size()>0) {
setOLPParameter("HelicityRN",helicityrn[0]);
}
}
+ #ifdef VBFNLO3
+ OLP_EvalSubProcess2_VBFNLO(&id, olpMomenta(), &scale, &colourCorrelatorResults[0], &acc);
+ #else
OLP_EvalSubProcess2(&id, olpMomenta(), &scale, &colourCorrelatorResults[0], &acc);
+ #endif
for ( int i = 0; i < n; ++i )
for ( int j = i+1; j < n; ++j ) {
lastColourCorrelator(make_pair(i,j),colourCorrelatorResults[i+j*(j-1)/2]*units);
}
}
void VBFNLOAmplitude::evalSpinColourCorrelator(pair<int,int>) const {
double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double scale = sqrt(mu2()/GeV2);
- if (hasRunningAlphaS()) setOLPParameter("alphas",lastAlphaS());
+ if (hasRunningAlphaS())
+ setOLPParameter("alphas",lastAlphaS());
double acc = -1.0;
int n = lastXComb().meMomenta().size();
spinColourCorrelatorResults.resize(2*n*n);
int id = olpId()[ProcessType::spinColourCorrelatedME2];
if (theRanHelSum) {
vector<double> helicityrn;
if ( lastHeadMatchboxXComb() ) {
helicityrn = lastHeadMatchboxXComb()->amplitudeRandomNumbers();
} else {
helicityrn = amplitudeRandomNumbers();
}
if (helicityrn.size()>0) {
setOLPParameter("HelicityRN",helicityrn[0]);
}
}
+ #ifdef VBFNLO3
+ OLP_EvalSubProcess2_VBFNLO(&id, olpMomenta(), &scale, &spinColourCorrelatorResults[0], &acc);
+ #else
OLP_EvalSubProcess2(&id, olpMomenta(), &scale, &spinColourCorrelatorResults[0], &acc);
+ #endif
for ( int i = 0; i < n; ++i )
for ( int j = 0; j < n; ++j ) {
if ( i == j || mePartonData()[i]->id() != 21 )
continue;
Complex scc(spinColourCorrelatorResults[2*i+2*n*j]*units,
spinColourCorrelatorResults[2*i+2*n*j+1]*units);
lastColourSpinCorrelator(make_pair(i,j),scc);
}
}
double VBFNLOAmplitude::largeNME2(Ptr<ColourBasis>::tptr cptr) const {
if ( calculateLargeNME2() )
evalLargeNSubProcess(cptr);
return lastLargeNME2();
}
void VBFNLOAmplitude::evalLargeNSubProcess(Ptr<ColourBasis>::tptr) const {
double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double scale = sqrt(mu2()/GeV2);
if (hasRunningAlphaS()) setOLPParameter("alphas",lastAlphaS());
double acc = -1.0;
double out[4]={};
int id =
olpId()[ProcessType::oneLoopInterference] ?
olpId()[ProcessType::oneLoopInterference] :
olpId()[ProcessType::treeME2];
if (theRanHelSum) {
vector<double> helicityrn;
if ( lastHeadMatchboxXComb() ) {
helicityrn = lastHeadMatchboxXComb()->amplitudeRandomNumbers();
} else {
helicityrn = amplitudeRandomNumbers();
}
if (helicityrn.size()>0) {
setOLPParameter("HelicityRN",helicityrn[0]);
}
}
setOLPParameter("Nc",-1); // large-N limit
+ #ifdef VBFNLO3
+ OLP_EvalSubProcess2_VBFNLO(&id, olpMomenta(), &scale, out, &acc);
+ #else
OLP_EvalSubProcess2(&id, olpMomenta(), &scale, out, &acc);
+ #endif
setOLPParameter("Nc",generator()->standardModel()->Nc());
if ( olpId()[ProcessType::oneLoopInterference] ) {
lastLargeNME2(out[3]*units);
lastOneLoopInterference(out[2]*units);
lastOneLoopPoles(pair<double,double>(out[0]*units,out[1]*units));
} else if ( olpId()[ProcessType::treeME2] ) {
lastLargeNME2(out[0]*units);
} else assert(false);
}
double VBFNLOAmplitude::largeNColourCorrelatedME2(pair<int,int> ij,
Ptr<ColourBasis>::tptr cptr) const {
double cfac = 1.;
double Nc = generator()->standardModel()->Nc();
if ( mePartonData()[ij.first]->iColour() == PDT::Colour8 ) {
cfac = Nc;
} else if ( mePartonData()[ij.first]->iColour() == PDT::Colour3 ||
mePartonData()[ij.first]->iColour() == PDT::Colour3bar ) {
cfac = Nc/2.;
} else assert(false);
if ( calculateLargeNColourCorrelator(ij) )
evalLargeNColourCorrelator(ij,cptr);
return lastLargeNColourCorrelator(ij)/cfac;
}
void VBFNLOAmplitude::evalLargeNColourCorrelator(pair<int,int>,
Ptr<ColourBasis>::tptr) const {
double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double scale = sqrt(mu2()/GeV2);
if (hasRunningAlphaS()) setOLPParameter("alphas",lastAlphaS());
double acc = -1.0;
int n = lastXComb().meMomenta().size();
colourCorrelatorResults.resize(n*(n-1)/2);
int id = olpId()[ProcessType::colourCorrelatedME2];
if (theRanHelSum) {
vector<double> helicityrn;
if ( lastHeadMatchboxXComb() ) {
helicityrn = lastHeadMatchboxXComb()->amplitudeRandomNumbers();
} else {
helicityrn = amplitudeRandomNumbers();
}
if (helicityrn.size()>0) {
setOLPParameter("HelicityRN",helicityrn[0]);
}
}
setOLPParameter("Nc",-1); // large-N limit
+ #ifdef VBFNLO3
+ OLP_EvalSubProcess2_VBFNLO(&id, olpMomenta(), &scale, &colourCorrelatorResults[0], &acc);
+ #else
OLP_EvalSubProcess2(&id, olpMomenta(), &scale, &colourCorrelatorResults[0], &acc);
+ #endif
setOLPParameter("Nc",generator()->standardModel()->Nc());
for ( int i = 0; i < n; ++i )
for ( int j = i+1; j < n; ++j ) {
lastLargeNColourCorrelator(make_pair(i,j),colourCorrelatorResults[i+j*(j-1)/2]*units);
}
}
void VBFNLOAmplitude::doinit() {
loadVBFNLO();
MatchboxOLPME::doinit();
}
void VBFNLOAmplitude::doinitrun() {
loadVBFNLO();
MatchboxOLPME::doinitrun();
}
void VBFNLOAmplitude::persistentOutput(PersistentOStream & os) const {
os << colourCorrelatorResults << spinColourCorrelatorResults << theRanHelSum << theAnomCoupl << VBFNLOlib_;
}
void VBFNLOAmplitude::persistentInput(PersistentIStream & is, int) {
is >> colourCorrelatorResults >> spinColourCorrelatorResults >> theRanHelSum >> theAnomCoupl >> VBFNLOlib_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<VBFNLOAmplitude,MatchboxOLPME>
describeHerwigVBFNLOAmplitude("Herwig::VBFNLOAmplitude", "HwMatchboxVBFNLO.so");
void VBFNLOAmplitude::Init() {
static ClassDocumentation<VBFNLOAmplitude> documentation
("VBFNLOAmplitude implements an interface to VBFNLO.",
"Matrix elements have been calculated using VBFNLO "
"(Ref.~\\cite{VBFNLO} and process-specific references)\n",
"%\\cite{VBFNLO}\n"
"\\bibitem{Arnold:2008rz}\n"
"K.~Arnold, M.~Bahr, G.~Bozzi, F.~Campanario, C.~Englert, T.~Figy, "
"N.~Greiner and C.~Hackstein {\\it et al.},\n"
"``VBFNLO: A Parton level Monte Carlo for processes with electroweak bosons,''\n"
"Comput.\\ Phys.\\ Commun.\\ {\\bf 180} (2009) 1661\n"
"[arXiv:0811.4559 [hep-ph]];\n"
"%%CITATION = ARXIV:0811.4559;%%\n"
"J.~Baglio, J.~Bellm, F.~Campanario, B.~Feigl, J.~Frank, T.~Figy, "
"M.~Kerner and L.~D.~Ninh {\\it et al.},\n"
"``Release Note - VBFNLO 2.7.0,''\n"
"arXiv:1404.3940 [hep-ph].\n"
"%%CITATION = ARXIV:1404.3940;%%\n");
static Switch<VBFNLOAmplitude,bool> interfaceRandomHelicitySummation
("RandomHelicitySummation", "Switch for random helicity summation of leptons and photons",
&VBFNLOAmplitude::theRanHelSum, false, false, false);
static SwitchOption interfaceRandomHelicitySummationYes
(interfaceRandomHelicitySummation,
"Yes",
"Perform random helicity summation",
true);
static SwitchOption interfaceRandomHelicitySummationNo
(interfaceRandomHelicitySummation,
"No",
"Sum over all helicity combinations",
false);
static Switch<VBFNLOAmplitude,bool> interfaceAnomalousCouplings
("AnomalousCouplings", "Switch for anomalous couplings",
&VBFNLOAmplitude::theAnomCoupl, false, false, false);
static SwitchOption interfaceAnomalousCouplingsYes
(interfaceAnomalousCouplings,
"Yes",
"Switch anomalous couplings on",
true);
static SwitchOption interfaceAnomalousCouplingsNo
(interfaceAnomalousCouplings,
"No",
"Switch anomalous couplings off",
false);
}
diff --git a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
--- a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
+++ b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
@@ -1,277 +1,289 @@
// -*- C++ -*-
//
// VBFNLOPhasespace.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 VBFNLOPhasespace class.
//
#include "VBFNLOPhasespace.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/Utilities/GSLBisection.h"
#include "ThePEG/Utilities/DynamicLoader.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
#include "VBFNLO/utilities/BLHAinterface.h"
#define DEFSTR(s) CPPSTR(s)
#define CPPSTR(s) #s
using namespace Herwig;
VBFNLOPhasespace::VBFNLOPhasespace() :
lastSqrtS(0*GeV), needToReshuffle(false), VBFNLOlib_(DEFSTR(VBFNLOLIB))
{}
void VBFNLOPhasespace::loadVBFNLO() {
if ( ! DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.so") ) {
string error1 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.dylib") ) {
string error2 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load("libVBFNLO.so") ) {
string error3 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load("libVBFNLO.dylib") ) {
string error4 = DynamicLoader::lastErrorMessage;
throw Exception() << "VBFNLOPhasespace: failed to load libVBFNLO.so/dylib\n"
<< "Error messages are:\n\n"
<< "* " << VBFNLOlib_ << "/libVBFNLO.so:\n"
<< error1 << "\n"
<< "* " << VBFNLOlib_ << "/libVBFNLO.dylib:\n"
<< error2 << "\n"
<< "* libVBFNLO.so:\n"
<< error3 << "\n"
<< "* libVBFNLO.dylib:\n"
<< error4 << "\n"
<< Exception::runerror;
}
}
}
}
}
VBFNLOPhasespace::~VBFNLOPhasespace() {}
IBPtr VBFNLOPhasespace::clone() const {
return new_ptr(*this);
}
IBPtr VBFNLOPhasespace::fullclone() const {
return new_ptr(*this);
}
void VBFNLOPhasespace::setXComb(tStdXCombPtr xco) {
MatchboxPhasespace::setXComb(xco);
// test for resuffling
needToReshuffle = false;
if ( xco ) {
for ( cPDVector::const_iterator d = mePartonData().begin();
d != mePartonData().end(); ++d ) {
// Higgs is massive -> does not need reshuffling
if ( ( (**d).id() != ParticleID::h0 ) && ( (**d).hardProcessMass() != ZERO ) ) {
needToReshuffle = true;
break;
}
}
}
// set CMS energy
int pStatus = 0;
double zero = 0.0;
double value = sqrt(lastXCombPtr()->lastS())/GeV;
if (value && (value != lastSqrtS/GeV)) {
lastSqrtS = value*GeV;
string name = "sqrtS";
+ #ifdef VBFNLO3
+ OLP_SetParameter_VBFNLO(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
+ #else
OLP_SetParameter(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
+ #endif
if ( !pStatus )
throw Exception() << "VBFNLOPhasespace::setXComb(): VBFNLO failed to set parameter '"
<< name << "' to " << value << "\n"
<< Exception::runerror;
}
}
double VBFNLOPhasespace::generateTwoToNKinematics(const double* random,
vector<Lorentz5Momentum>& momenta) {
double weight;
int id =
olpId()[ProcessType::oneLoopInterference] ?
olpId()[ProcessType::oneLoopInterference] :
olpId()[ProcessType::treeME2];
double* p = new double[4*momenta.size()];
+ #ifdef VBFNLO3
+ OLP_PhaseSpacePoint_VBFNLO(&id, const_cast<double*>(random), const_cast<double*>(random+1), p, &weight);
+ #else
OLP_PhaseSpacePoint(&id, const_cast<double*>(random), const_cast<double*>(random+1), p, &weight);
+ #endif
if (weight < 0) {
throw Exception() << "VBFNLOPhasespace::generateTwoToNKinematics(): Negative weight in VBFNLOPhaseSpace\n"
<< Exception::runerror;
}
if (weight == 0) {
delete[] p;
return 0;
}
for ( size_t i = 0; i < momenta.size(); ++i ) {
momenta[i].setT(p[4*i] *GeV);
momenta[i].setX(p[4*i+1]*GeV);
momenta[i].setY(p[4*i+2]*GeV);
momenta[i].setZ(p[4*i+3]*GeV);
momenta[i].rescaleMass();
}
delete[] p;
Energy beamenergy = sqrt(lastXCombPtr()->lastS())/2.;
double x1 = momenta[0].e()/beamenergy;
double x2 = momenta[1].e()/beamenergy;
Energy2 thisSHat = (momenta[0] + momenta[1]).m2();
// reshuffle so that particles have correct mass
if ( needToReshuffle ) {
// boost final-state into partonic CMS
Boost toCMS = (momenta[0]+momenta[1]).findBoostToCM();
for ( size_t i = 2; i < momenta.size(); ++i ) {
momenta[i].boost(toCMS);
}
// copied from MatchboxRambo phasespace
double xi;
ReshuffleEquation solve(sqrt(thisSHat),mePartonData().begin()+2,mePartonData().end(),
momenta.begin()+2,momenta.end());
GSLBisection solver(1e-10,1e-8,10000);
try {
xi = solver.value(solve,0.0,1.1);
} catch (GSLBisection::GSLerror) {
return 0.;
} catch (GSLBisection::IntervalError) {
return 0.;
}
weight *= pow(xi,3.*(momenta.size()-3.));
Energy num = ZERO;
Energy den = ZERO;
cPDVector::const_iterator d = mePartonData().begin()+2;
for ( vector<Lorentz5Momentum>::iterator k = momenta.begin()+2;
k != momenta.end(); ++k, ++d ) {
num += (*k).vect().mag2()/(*k).t();
Energy q = (*k).t();
(*k).setT(sqrt(sqr((**d).hardProcessMass())+xi*xi*sqr((*k).t())));
(*k).setVect(xi*(*k).vect());
weight *= q/(*k).t();
den += (*k).vect().mag2()/(*k).t();
(*k).setMass((**d).hardProcessMass());
}
// unboost
for ( size_t i = 2; i < momenta.size(); ++i ) {
momenta[i].boost(-toCMS);
}
}
if ( !matchConstraints(momenta) )
return 0.;
lastXCombPtr()->lastX1X2(make_pair(x1,x2));
lastXCombPtr()->lastSHat(thisSHat);
weight /= pow(thisSHat/GeV2,momenta.size()-4);
weight /= x1*x2;
fillDiagramWeights();
return weight;
}
int VBFNLOPhasespace::nDimPhasespace(int nFinal) const {
return 3*nFinal;
//get this from within VBFNLO
int pStatus = 0;
double value, zero;
string name = "PSdimension";
+ #ifdef VBFNLO3
+ OLP_GetParameter_VBFNLO(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
+ #else
OLP_GetParameter(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
+ #endif
if ( pStatus != 1) {
throw Exception() << "VBFNLOPhasespace::nDimPhasespace(): Cannot get phasespace dimension in VBFNLOPhaseSpace\n"
<< "error code: " << pStatus << "\n"
<< Exception::runerror;
}
// one additional number (first) needed for channel selection
// one additional number (last) needed for global phi integration
return value+2;
}
Energy VBFNLOPhasespace::ReshuffleEquation::operator() (double xi) const {
cPDVector::const_iterator d = dataBegin;
vector<Lorentz5Momentum>::const_iterator p = momentaBegin;
Energy res = -w;
for ( ; d != dataEnd; ++d, ++p ) {
res += sqrt(sqr((**d).hardProcessMass()) +
xi*xi*sqr(p->t()));
}
return res;
}
void VBFNLOPhasespace::doinit() {
loadVBFNLO();
MatchboxPhasespace::doinit();
}
void VBFNLOPhasespace::doinitrun() {
loadVBFNLO();
MatchboxPhasespace::doinitrun();
}
void VBFNLOPhasespace::persistentOutput(PersistentOStream & os) const {
os << needToReshuffle << theLastXComb;
}
void VBFNLOPhasespace::persistentInput(PersistentIStream & is, int) {
is >> needToReshuffle >> theLastXComb;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<VBFNLOPhasespace,MatchboxPhasespace>
describeHerwigVBFNLOPhasespace("Herwig::VBFNLOPhasespace", "HwMatchboxVBFNLO.so");
void VBFNLOPhasespace::Init() {
static ClassDocumentation<VBFNLOPhasespace> documentation
("VBFNLOPhasespace is an interface to the internal phasespace generator "
"of VBFNLO. It uses the information passed via the BLHA interface to "
"obtain information on the required channels.");
}
diff --git a/Shower/QTilde/Couplings/ShowerAlphaQCDNP.cc b/Shower/QTilde/Couplings/ShowerAlphaQCDNP.cc
new file mode 100644
--- /dev/null
+++ b/Shower/QTilde/Couplings/ShowerAlphaQCDNP.cc
@@ -0,0 +1,419 @@
+// -*- C++ -*-
+//
+// ShowerAlphaQCDNP.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 ShowerAlphaQCDNP class.
+//
+#include "ShowerAlphaQCDNP.h"
+#include "ThePEG/PDT/ParticleData.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/Interface/Parameter.h"
+#include "ThePEG/Interface/ParVector.h"
+#include "ThePEG/Interface/Command.h"
+#include "ThePEG/Interface/Deleted.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/Utilities/Throw.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include "ThePEG/Config/Constants.h"
+#include "Herwig/Utilities/AlphaS.h"
+#include "gsl/gsl_sf_lambert.h"
+
+using namespace Herwig;
+using Herwig::Math::alphaS;
+using Herwig::Math::derivativeAlphaS;
+
+DescribeClass<ShowerAlphaQCDNP,ShowerAlpha>
+describeShowerAlphaQCDNP("Herwig::ShowerAlphaQCDNP","HwShower.so");
+
+IBPtr ShowerAlphaQCDNP::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr ShowerAlphaQCDNP::fullclone() const {
+ return new_ptr(*this);
+}
+
+void ShowerAlphaQCDNP::persistentOutput(PersistentOStream & os) const {
+ os << ounit(_qmin,GeV) << _asMaxNP << _nfMaxNP
+ << ounit(_thresholds,GeV) << ounit(_lambda,GeV)
+ << _nloop << _thresopt
+ << _alphain << _tolerance << _maxtry
+ << _npPower << ounit(_optInputScale,GeV) << _quarkBranching
+ << ounit(_quarkMasses,GeV);
+}
+
+void ShowerAlphaQCDNP::persistentInput(PersistentIStream & is, int) {
+ is >> iunit(_qmin,GeV) >> _asMaxNP >> _nfMaxNP
+ >> iunit(_thresholds,GeV) >> iunit(_lambda,GeV)
+ >> _nloop >> _thresopt
+ >> _alphain >> _tolerance >> _maxtry
+ >> _npPower >> iunit(_optInputScale,GeV) >> _quarkBranching
+ >> iunit(_quarkMasses,GeV);
+}
+
+void ShowerAlphaQCDNP::Init() {
+
+ static ClassDocumentation<ShowerAlphaQCDNP> documentation
+ ("This (concrete) class describes the QCD alpha running.");
+
+ // default such that as(qmin) = 1 in the current parametrization.
+ static Parameter<ShowerAlphaQCDNP,Energy> intQmin
+ ("Qmin", "Q < Qmin is treated with NP parametrization as of (unit [GeV])",
+ &ShowerAlphaQCDNP::_qmin, GeV, 1*GeV, 1*GeV,
+ 100.0*GeV,false,false,false);
+
+ static Parameter<ShowerAlphaQCDNP,unsigned int> interfaceNumberOfLoops
+ ("NumberOfLoops",
+ "The number of loops to use in the alpha_S calculation",
+ &ShowerAlphaQCDNP::_nloop, 3, 1, 3,
+ false, false, Interface::limited);
+
+ static Parameter<ShowerAlphaQCDNP,double> interfaceAlphaMZ
+ ("AlphaIn",
+ "The input value of the strong coupling at the chosen InputScale (default: MZ)",
+ &ShowerAlphaQCDNP::_alphain, 0.118, 0.1, 0.2,
+ false, false, Interface::limited);
+
+ static Parameter<ShowerAlphaQCDNP,double> interfaceTolerance
+ ("Tolerance",
+ "The tolerance for discontinuities in alphaS at thresholds.",
+ &ShowerAlphaQCDNP::_tolerance, 1e-10, 1e-20, 1e-4,
+ false, false, Interface::limited);
+
+ static Parameter<ShowerAlphaQCDNP,unsigned int> interfaceMaximumIterations
+ ("MaximumIterations",
+ "The maximum number of iterations for the Newton-Raphson method to converge.",
+ &ShowerAlphaQCDNP::_maxtry, 100, 10, 1000,
+ false, false, Interface::limited);
+
+ static Switch<ShowerAlphaQCDNP,bool> interfaceThresholdOption
+ ("ThresholdOption",
+ "Whether to use the consistuent or normal masses for the thresholds",
+ &ShowerAlphaQCDNP::_thresopt, true, false, false);
+ static SwitchOption interfaceThresholdOptionCurrent
+ (interfaceThresholdOption,
+ "Current",
+ "Use the current masses",
+ true);
+ static SwitchOption interfaceThresholdOptionConstituent
+ (interfaceThresholdOption,
+ "Constituent",
+ "Use the constitent masses.",
+ false);
+
+ static Command<ShowerAlphaQCDNP> interfaceValue
+ ("Value",
+ "",
+ &ShowerAlphaQCDNP::value, false);
+
+ static Command<ShowerAlphaQCDNP> interfacecheck
+ ("check",
+ "check",
+ &ShowerAlphaQCDNP::check, false);
+
+ static Command<ShowerAlphaQCDNP> interfacecheckNf
+ ("checkNf",
+ "checkNf",
+ &ShowerAlphaQCDNP::checkNf, false);
+
+ static Command<ShowerAlphaQCDNP> interfacecheckscaleNP
+ ("checkscaleNP",
+ "checkscaleNP",
+ &ShowerAlphaQCDNP::checkscaleNP, false);
+
+ static Parameter<ShowerAlphaQCDNP,Energy> interfaceInputScale
+ ("InputScale",
+ "An optional input scale. MZ will be used if not set.",
+ &ShowerAlphaQCDNP::_optInputScale, GeV, 91.1876_GeV, ZERO, 0*GeV,
+ false, false, Interface::lowerlim);
+
+ static ParVector<ShowerAlphaQCDNP,Energy> interfaceQuarkMasses
+ ("QuarkMasses",
+ "The quark masses to be used instead of the masses set in the particle data.",
+ &ShowerAlphaQCDNP::_quarkMasses, GeV, -1, 0.0*GeV, 0.0*GeV, 0*GeV,
+ false, false, Interface::lowerlim);
+
+
+ static Switch<ShowerAlphaQCDNP,bool> interfaceQuarkBranching
+ ("QuarkBranching",
+ "True, if this coupling is used in a gluon to qqbar branching.",
+ &ShowerAlphaQCDNP::_quarkBranching, false, false, false);
+ static SwitchOption interfaceQuarkBranchingYes
+ (interfaceQuarkBranching,
+ "Yes",
+ "Use in gluon to qqbar branching.",
+ true);
+ static SwitchOption interfaceQuarkBranchingNo
+ (interfaceQuarkBranching,
+ "No",
+ "Use in gluon emission.",
+ false);
+
+ static Parameter<ShowerAlphaQCDNP,double> interfaceNPPower
+ ("NPPower",
+ "The non-perturbative power law",
+ &ShowerAlphaQCDNP::_npPower, 2.0, 1.0, 0,
+ false, false, Interface::lowerlim);
+
+}
+
+void ShowerAlphaQCDNP::doinit() {
+ ShowerAlpha::doinit();
+ // calculate the value of 5-flavour lambda
+ // evaluate the initial
+ // value of Lambda from alphas if needed using Newton-Raphson
+ _lambda[2]=computeLambda(_optInputScale,_alphain,5);
+
+ // compute the threshold matching
+ // top threshold
+ for(int ix=1;ix<4;++ix) {
+ if ( _quarkMasses.empty() ) {
+ if(_thresopt)
+ _thresholds[ix]=getParticleData(ix+3)->mass();
+ else
+ _thresholds[ix]=getParticleData(ix+3)->constituentMass();
+ } else {
+ // starting at zero rather than one, cf the other alphas's
+ _thresholds[ix] = _quarkMasses[ix+2];
+ }
+ }
+ // compute 6 flavour lambda by matching at top mass using Newton Raphson
+ _lambda[3]=computeLambda(_thresholds[3],alphaS(_thresholds[3],_lambda[2],5,_nloop),6);
+ // bottom threshold
+ // compute 4 flavour lambda by matching at bottom mass using Newton Raphson
+ _lambda[1]=computeLambda(_thresholds[2],alphaS(_thresholds[2],_lambda[2],5,_nloop),4);
+ // charm threshold
+ // compute 3 flavour lambda by matching at charm mass using Newton Raphson
+ _lambda[0]=computeLambda(_thresholds[1],alphaS(_thresholds[1],_lambda[1],4,_nloop),3);
+
+ //Energy q = scaleNPMin();
+ //auto nflam = getLamNf(q);
+ //_asMaxNP = alphaS(q, nflam.second, nflam.first, _nloop);
+
+ Energy q = scaleThatmMinimizesScaleNP() ;
+ _asMaxNP = value(q*q)*1.05;
+
+ //_absoluteCutoff=_qmin*pow(100.*log(10.),-1./_npPower);
+ _absoluteCutoff=2.0*GeV;
+
+ _nfMaxNP = nfNP(sqr(_absoluteCutoff));
+
+}
+
+string ShowerAlphaQCDNP::check(string args) {
+
+ doinit();
+
+ istringstream argin(args);
+
+ double Q_low, Q_high;
+ long n_steps;
+
+ argin >> Q_low >> Q_high >> n_steps;
+
+ string fname;
+ argin >> fname;
+
+ Repository::clog() << "checking alpha_s in range [" << Q_low << "," << Q_high << "] GeV in "
+ << n_steps << " steps.\nResults are written to " << fname << "\n";
+
+ double step_width = (Q_high-Q_low)/n_steps;
+
+ ofstream out (fname.c_str());
+
+ for (long k = 0; k <= n_steps; ++k) {
+
+ Energy Q = Q_low*GeV + k*step_width*GeV;
+
+ //out << (Q/GeV) << " " << value(Q*Q) << " " << _asMaxNP << "\n";
+ out << (Q/GeV) << " " << value(Q*Q) << "\n";
+
+
+ }
+
+ return "alpha_s check finished";
+
+}
+
+string ShowerAlphaQCDNP::checkNf(string args) {
+
+ doinit();
+
+ istringstream argin(args);
+
+ double Q_low, Q_high;
+ long n_steps;
+
+ argin >> Q_low >> Q_high >> n_steps;
+
+ string fname;
+ argin >> fname;
+
+ Repository::clog() << "checking nfNP in range [" << Q_low << "," << Q_high << "] GeV in "
+ << n_steps << " steps.\nResults are written to " << fname << "\n";
+
+ double step_width = (Q_high-Q_low)/n_steps;
+
+ ofstream out (fname.c_str());
+
+ for (long k = 0; k <= n_steps; ++k) {
+
+ Energy Q = Q_low*GeV + k*step_width*GeV;
+
+ //out << (Q/GeV) << " " << nfNP(Q*Q) << " " << _nfMaxNP << "\n";
+ out << (Q/GeV) << " " << nfNP(Q*Q) << "\n";
+
+ }
+
+ return "nfNP check finished";
+
+}
+
+
+
+string ShowerAlphaQCDNP::checkscaleNP(string args) {
+
+ doinit();
+
+ istringstream argin(args);
+
+ double Q_low, Q_high;
+ long n_steps;
+
+ argin >> Q_low >> Q_high >> n_steps;
+
+ string fname;
+ argin >> fname;
+
+ Repository::clog() << "checking scaleNP (in GeV) in range [" << Q_low << "," << Q_high << "] GeV in "
+ << n_steps << " steps.\nResults are written to " << fname << "\n";
+
+ double step_width = (Q_high-Q_low)/n_steps;
+
+ ofstream out (fname.c_str());
+
+ for (long k = 0; k <= n_steps; ++k) {
+
+ Energy Q = Q_low*GeV + k*step_width*GeV;
+
+ out << (Q/GeV) << " " << scaleNP(scaleFactor()*Q)/(1.*GeV) << "\n";
+
+ }
+
+ return "scaleNP check finished";
+
+}
+
+
+double ShowerAlphaQCDNP::value(const Energy2 scale) const {
+
+ Energy q = scaleFactor()*sqrt(scale);
+
+ if ( q > _absoluteCutoff ) {
+ auto nflam = getLamNf(q);
+ return alphaS(scaleNP(q), nflam.second, nflam.first, _nloop);
+ }
+
+
+ return 0.;
+
+}
+
+
+string ShowerAlphaQCDNP::value(string scale) {
+ istringstream readscale(scale);
+ double inScale; readscale >> inScale;
+ Energy theScale = inScale * GeV;
+ initialize();
+ ostringstream showvalue ("");
+ showvalue << "alpha_s (" << theScale/GeV << " GeV) = "
+ << value (sqr(theScale));
+ return showvalue.str();
+}
+
+pair<short, Energy> ShowerAlphaQCDNP::getLamNf(Energy q) const {
+ short nf = 6;
+ // get lambda and nf according to the thresholds
+ if (q < _thresholds[1]) nf = 3;
+ else if (q < _thresholds[2]) nf = 4;
+ else if (q < _thresholds[3]) nf = 5;
+ return pair<short,Energy>(nf, _lambda[nf-3]);
+}
+
+Energy ShowerAlphaQCDNP::computeLambda(Energy match,
+ double alpha,
+ unsigned int nflav) const {
+
+ Energy lamtest=200.0*MeV;
+ double xtest;
+ unsigned int ntry=0;
+ do {
+ ++ntry;
+ xtest = log(sqr(match/lamtest));
+ xtest += (alpha-alphaS(match,lamtest,nflav,_nloop))/derivativeAlphaS(match,lamtest,nflav,_nloop);
+ Energy newLambda = match/exp(0.5 *xtest);
+ lamtest = newLambda<match ? newLambda : 0.5*(lamtest+match);
+ }
+ while(abs(alpha-alphaS(match,lamtest,nflav,_nloop)) > _tolerance && ntry < _maxtry);
+ return lamtest;
+
+}
+
+// --------------------------------------------------------------------------------
+// main modification for nonperturbative running in terms of changed argument and
+// nonperturbative number of flavours
+// --------------------------------------------------------------------------------
+
+double ShowerAlphaQCDNP::nfNP(const Energy2 q2) const {
+ Energy q=sqrt(q2);
+
+ if (q<_absoluteCutoff){
+ return 0.;
+ }
+ else{
+
+ //get p-nf
+ double nf=getLamNf(q).first;
+
+ //beta function coefficients; b=sum_i (alphaS/(4Pi))^i * sum_j b_ij *nf^j
+ using Constants::pi;
+ double b00 = 11.;
+ double b01 = - 2./3.;
+ double b10 = 51.;
+ double b11 = - 19./3.;
+ double c = (b00+b10*value(q2)/(4.*pi))/(b01+b11*value(q2)/(4.*pi));
+
+ //absolute number of np-flavors
+ double nfNPabs=(c+nf)*(q*scaleNPDerivative(q)/scaleNP(q))-c;
+
+ //return ratio of np-flavors over p-flavors
+ return nfNPabs/nf;
+ }
+}
+
+Energy ShowerAlphaQCDNP::scaleNP(const Energy q) const {
+ return q*(1.+(_qmin/q)*(exp(pow(_qmin/q,_npPower))-1.));
+}
+
+double ShowerAlphaQCDNP::scaleNPDerivative(const Energy q) const {
+ return 1.-exp(pow(_qmin/q,_npPower))*_npPower*pow(_qmin/q,1.+_npPower);
+}
+
+Energy ShowerAlphaQCDNP::scaleThatmMinimizesScaleNP() const {
+ double x=(1./_npPower)*gsl_sf_lambert_W0(pow(_npPower,1./(1.+_npPower))/(1.+_npPower));
+ return _qmin*exp(x)*pow(_npPower,1./(1.+_npPower));
+
+}
+
+Energy ShowerAlphaQCDNP::scaleNPMin() const {
+ double x=((1.+_npPower)/_npPower)*gsl_sf_lambert_W0(pow(_npPower,1./(1.+_npPower))/(1.+_npPower));
+ return _qmin*(-1.+exp(x)+pow(x,-1./_npPower));
+}
diff --git a/Shower/QTilde/Couplings/ShowerAlphaQCDNP.h b/Shower/QTilde/Couplings/ShowerAlphaQCDNP.h
new file mode 100644
--- /dev/null
+++ b/Shower/QTilde/Couplings/ShowerAlphaQCDNP.h
@@ -0,0 +1,382 @@
+// -*- C++ -*-
+//
+// ShowerAlphaQCDNP.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_ShowerAlphaQCDNP_H
+#define HERWIG_ShowerAlphaQCDNP_H
+//
+// This is the declaration of the ShowerAlphaQCDNP class.
+//
+
+#include "Herwig/Shower/ShowerAlpha.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/** \ingroup Shower
+ *
+ * This concrete class provides the definition of the
+ * pure virtual function value() and overestimateValue() for the
+ * strong coupling.
+ *
+ * A number of different options for the running of the coupling
+ * and its initial definition are supported.
+ *
+ * @see \ref ShowerAlphaQCDNPInterfaces "The interfaces"
+ * defined for ShowerAlphaQCDNP.
+ */
+class ShowerAlphaQCDNP: public ShowerAlpha {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ ShowerAlphaQCDNP() :
+ ShowerAlpha(),
+ _qmin(1*GeV), _asMaxNP(1.0), _nfMaxNP(15.0),
+ _thresholds(4), _lambda(4),
+ _nloop(2), _thresopt(false),
+ _alphain(0.118), _tolerance(1e-10),
+ _maxtry(100),
+ _npPower(2.), _optInputScale(91.1876_GeV),
+ _quarkBranching(false) {}
+
+public:
+
+ /**
+ * Methods to return the coupling
+ */
+ //@{
+ /**
+ * It returns the running coupling value evaluated at the input scale
+ * multiplied by the scale factor scaleFactor().
+ * @param scale The scale
+ * @return The coupling
+ */
+ virtual double value(const Energy2 scale) const;
+
+ /**
+ * Virtual method that is supposed to return the
+ * running alpha value evaluated at the input scale.
+ * @param scale The scale
+ * @return The coupling
+ */
+ virtual double showerValue(const Energy2 scale) const {
+ if ( !_quarkBranching ) {
+ return value(scale);
+ }
+ return nfNP(scale)*value(scale);
+ }
+
+ /**
+ * It returns the running coupling value evaluated at the input scale
+ * multiplied by the scale factor scaleFactor().
+ */
+ virtual double overestimateValue() const {
+ return _asMaxNP;
+ }
+
+ /**
+ * Virtual method, which
+ * should be overridden in a derived class to provide an
+ * overestimate approximation of the alpha value.
+ */
+ virtual double showerOverestimateValue() const {
+ if ( !_quarkBranching ) {
+ return overestimateValue();
+ }
+ return overestimateNfNP()*overestimateValue();
+ }
+
+ /**
+ * Return the ratio of the coupling at the scale to the overestimated value
+ */
+ virtual double ratio(const Energy2 scale,double factor =1.) const {
+ Energy2 q2 = sqr(factor)*scale;
+ return value(q2)/overestimateValue();
+ }
+
+ /**
+ * Virtual method which returns the ratio of the running alpha
+ * value at the input scale to the overestimated value.
+ * @param scale The scale
+ * @return The ratio
+ */
+ virtual double showerRatio(const Energy2 scale,double factor=1.) const {
+
+ if ( !_quarkBranching ) {
+ return ratio(scale,factor);
+ }
+ return ratioNfNP(scale,factor)*ratio(scale,factor);
+ }
+
+ /**
+ * Initialize this coupling.
+ */
+ virtual void initialize() { doinit(); }
+
+ /**
+ * A command to initialize the coupling and write
+ * its value at the scale given by the argument (in GeV)
+ */
+ string value(string);
+
+ /**
+ * Match thresholds and write alpha_s
+ * specified file; arguments are
+ * Q_low/GeV Q_high/GeV n_steps filename
+ */
+ string check(string args);
+
+ /**
+ * Match thresholds and write alpha_s
+ * specified file; arguments are
+ * Q_low/GeV Q_high/GeV n_steps filename
+ */
+ string checkNf(string args);
+
+ string checkscaleNP(string args);
+
+ //@}
+
+ /**
+ * Get the value of \f$\Lambda_{\rm QCd}\f$
+ * @param nf number of flavours
+ */
+ Energy lambdaQCD(unsigned int nf) {
+ if (nf <= 3) return _lambda[0];
+ else if (nf==4 || nf==5) return _lambda[nf-3];
+ else return _lambda[3];
+ }
+
+ /**
+ * Return the quark masses to be used; if not empty these masses
+ * should be considered instead of the ones set in the particle data
+ * objects.
+ */
+ const vector<Energy>& quarkMasses() const { return _quarkMasses; }
+
+ /**
+ * Return the non-perturbative running flavour number, normalized to
+ * the number of active flavours.
+ */
+ double nfNP(const Energy2 scale) const;
+
+ /**
+ * Return the modification of the scale argument
+ */
+ Energy scaleNP(const Energy scale) const;
+
+ /**
+ * Return the the derivative of scaleNP
+ */
+ double scaleNPDerivative(const Energy scale) const;
+
+ /**
+ * Return the minimum of the modification of the scale argument
+ */
+ Energy scaleNPMin() const;
+
+
+ /**
+ * Return the scale q that minimizes scaleNP (i.e. scaleNP(q)=scaleNPMin() )
+ */
+ Energy scaleThatmMinimizesScaleNP() const;
+
+ /**
+ * Return the maximum of the non-perturbative running flavour number
+ * normalized to the number of perturbatively active flavours.
+ */
+ double overestimateNfNP() const {
+ return _nfMaxNP;
+ }
+
+ /**
+ * Return ration to the maximum of the non-perturbative running
+ * flavour number normalized to the number of perturbatively active
+ * flavours.
+ */
+ double ratioNfNP(const Energy2 scale,double factor=1.) const {
+ Energy2 q2 = sqr(factor)*scale;
+ return nfNP(q2)/overestimateNfNP();
+ }
+
+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;
+ //@}
+
+
+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();
+ //@}
+
+private:
+
+ /**
+ * Member functions which calculate the coupling
+ */
+ //@{
+ /**
+ * Compute the value of \f$Lambda\f$ needed to get the input value of
+ * the strong coupling at the scale given for the given number of flavours
+ * using the Newton-Raphson method
+ * @param match The scale for the coupling
+ * @param alpha The input coupling
+ * @param nflav The number of flavours
+ */
+ Energy computeLambda(Energy match, double alpha, unsigned int nflav) const;
+
+ /**
+ * Return the value of \f$\Lambda\f$ and the number of flavours at the scale.
+ * @param q The scale
+ * @return The number of flavours at the scale and \f$\Lambda\f$.
+ */
+ pair<short, Energy> getLamNf(Energy q) const;
+ //@}
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ ShowerAlphaQCDNP & operator=(const ShowerAlphaQCDNP &) = delete;
+
+private:
+
+ /**
+ * Transition value
+ */
+ Energy _qmin;
+
+ /**
+ * Maximum value of alphas
+ */
+ double _asMaxNP;
+
+ /**
+ * Maximum of non-perturbative flavour number
+ */
+ double _nfMaxNP;
+
+ /**
+ * Thresholds for the different number of flavours
+ */
+ vector<Energy> _thresholds;
+
+ /**
+ * \f$\Lambda\f$ for the different number of flavours
+ */
+ vector<Energy> _lambda;
+
+ /**
+ * Option for the number of loops
+ */
+ unsigned int _nloop;
+
+ /**
+ * Option for the threshold masses
+ */
+ bool _thresopt;
+
+ /**
+ * Input value of \f$alpha_S(M_Z)\f$
+ */
+ double _alphain;
+
+ /**
+ * Tolerance for discontinuities at the thresholds
+ */
+ double _tolerance;
+
+ /**
+ * Maximum number of iterations for the Newton-Raphson method to converge
+ */
+ unsigned int _maxtry;
+
+
+ /**
+ * Power of non-perturbative modification
+ */
+ double _npPower;
+
+ /**
+ * Lowest cutoff scale
+ */
+ Energy _absoluteCutoff;
+
+ /**
+ * An optional input scale to be used for the input alphas; if zero MZ will
+ * be used out of the particle data object.
+ */
+ Energy _optInputScale;
+
+ /**
+ * Flag if factor for quark branching should be included
+ */
+ bool _quarkBranching;
+
+ /**
+ * The quark masses to be used; if not empty these masses should be
+ * considered instead of the ones set in the particle data objects.
+ */
+ vector<Energy> _quarkMasses;
+
+};
+
+}
+
+#endif /* HERWIG_ShowerAlphaQCDNP_H */
diff --git a/Shower/QTilde/Makefile.am b/Shower/QTilde/Makefile.am
--- a/Shower/QTilde/Makefile.am
+++ b/Shower/QTilde/Makefile.am
@@ -1,57 +1,58 @@
SUBDIRS = Matching
pkglib_LTLIBRARIES = HwShower.la
HwShower_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 30:0:0
HwShower_la_SOURCES = \
Couplings/ShowerAlphaQCD.h Couplings/ShowerAlphaQCD.cc \
Couplings/ShowerAlphaQED.h Couplings/ShowerAlphaQED.cc\
+Couplings/ShowerAlphaQCDNP.h Couplings/ShowerAlphaQCDNP.cc\
QTildeShowerHandler.h QTildeShowerHandler.fh QTildeShowerHandler.cc \
SplittingFunctions/HalfHalfOneSplitFn.h SplittingFunctions/HalfHalfOneSplitFn.cc\
SplittingFunctions/HalfHalfOneEWSplitFn.h SplittingFunctions/HalfHalfOneEWSplitFn.cc\
SplittingFunctions/HalfHalfZeroEWSplitFn.h SplittingFunctions/HalfHalfZeroEWSplitFn.cc\
SplittingFunctions/OneOneOneEWSplitFn.h SplittingFunctions/OneOneOneEWSplitFn.cc\
SplittingFunctions/OneOneOneQEDSplitFn.h SplittingFunctions/OneOneOneQEDSplitFn.cc\
SplittingFunctions/OneOneZeroEWSplitFn.h SplittingFunctions/OneOneZeroEWSplitFn.cc\
SplittingFunctions/OneOneOneSplitFn.h SplittingFunctions/OneOneOneSplitFn.cc\
SplittingFunctions/OneOneOneMassiveSplitFn.h SplittingFunctions/OneOneOneMassiveSplitFn.cc\
SplittingFunctions/ZeroZeroOneSplitFn.h SplittingFunctions/ZeroZeroOneSplitFn.cc\
SplittingFunctions/OneHalfHalfSplitFn.h SplittingFunctions/OneHalfHalfSplitFn.cc\
SplittingFunctions/HalfOneHalfSplitFn.h SplittingFunctions/HalfOneHalfSplitFn.cc\
SplittingFunctions/CMWOneOneOneSplitFn.h SplittingFunctions/CMWOneOneOneSplitFn.cc\
SplittingFunctions/CMWHalfHalfOneSplitFn.h SplittingFunctions/CMWHalfHalfOneSplitFn.cc\
Dark/EnumParticles.h Dark/HiddenValleyAlpha.cc Dark/HiddenValleyAlpha.h \
Dark/HiddenValleyModel.cc Dark/HiddenValleyModel.h \
Dark/HiddenValleyFFZPrimeVertex.cc Dark/HiddenValleyFFZPrimeVertex.h \
Dark/OneHalfHalfDarkSplitFn.cc Dark/OneHalfHalfDarkSplitFn.h \
Dark/OneOneOneDarkSplitFn.cc Dark/OneOneOneDarkSplitFn.h \
Dark/HalfHalfOneDarkSplitFn.cc Dark/HalfHalfOneDarkSplitFn.h \
Kinematics/Decay_QTildeShowerKinematics1to2.cc \
Kinematics/Decay_QTildeShowerKinematics1to2.h \
Kinematics/IS_QTildeShowerKinematics1to2.cc Kinematics/IS_QTildeShowerKinematics1to2.h \
Kinematics/FS_QTildeShowerKinematics1to2.cc Kinematics/FS_QTildeShowerKinematics1to2.h \
Kinematics/KinematicHelpers.h \
Kinematics/KinematicsReconstructor.cc \
Kinematics/KinematicsReconstructor.tcc \
Kinematics/KinematicsReconstructor.h \
Kinematics/KinematicsReconstructor.fh \
Base/HardTree.cc Base/HardTree.h Base/HardTree.fh \
Base/HardBranching.h Base/HardBranching.fh Base/HardBranching.cc\
Base/PartnerFinder.h Base/PartnerFinder.fh Base/PartnerFinder.cc \
Base/ShowerVeto.h Base/ShowerVeto.fh Base/ShowerVeto.cc \
Base/FullShowerVeto.h Base/FullShowerVeto.fh Base/FullShowerVeto.cc \
SplittingFunctions/SplittingGenerator.cc SplittingFunctions/SplittingGenerator.h\
SplittingFunctions/SplittingGenerator.fh \
Base/ShowerTree.h Base/ShowerTree.fh Base/ShowerTree.cc \
ShowerConfig.h ShowerConfig.cc \
Base/Branching.h \
Base/ShowerParticle.cc Base/ShowerParticle.fh Base/ShowerParticle.h \
Kinematics/ShowerKinematics.fh Kinematics/ShowerKinematics.h Kinematics/ShowerKinematics.cc \
Kinematics/ShowerBasis.fh Kinematics/ShowerBasis.h Kinematics/ShowerBasis.cc \
Base/ShowerProgenitor.fh Base/ShowerProgenitor.h \
SplittingFunctions/SudakovFormFactor.cc SplittingFunctions/SudakovFormFactor.h SplittingFunctions/SudakovFormFactor.fh \
SplittingFunctions/SudakovCutOff.cc SplittingFunctions/SudakovCutOff.h SplittingFunctions/SudakovCutOff.fh \
SplittingFunctions/PTCutOff.cc SplittingFunctions/PTCutOff.h \
SplittingFunctions/MassCutOff.cc SplittingFunctions/MassCutOff.h \
SplittingFunctions/VariableMassCutOff.cc SplittingFunctions/VariableMassCutOff.h \
SplittingFunctions/SplittingFunction.h SplittingFunctions/SplittingFunction.fh \
SplittingFunctions/SplittingFunction.cc \
Base/ShowerVertex.cc Base/ShowerVertex.fh Base/ShowerVertex.h
diff --git a/m4/herwig.m4 b/m4/herwig.m4
--- a/m4/herwig.m4
+++ b/m4/herwig.m4
@@ -1,1009 +1,1026 @@
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])
THEPEGHASHEPMC="no"
if test -e ${with_thepeg}/lib/ThePEG/HepMCAnalysis.so ; then
THEPEGHASHEPMC="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/HepMCAnalysis.so ; then
THEPEGHASHEPMC="yes"
fi
fi
if test "x$THEPEGHASHEPMC" == "xno" ; then
CREATE_HEPMC="# create"
AC_MSG_RESULT([not found])
else
CREATE_HEPMC="create"
AC_MSG_RESULT([found])
fi
AC_SUBST([CREATE_HEPMC])
AC_MSG_CHECKING([for RivetAnalysis.so in ThePEG])
THEPEGHASRIVET="no"
if test -e ${with_thepeg}/lib/ThePEG/RivetAnalysis.so ; then
THEPEGHASRIVET="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/RivetAnalysis.so ; then
THEPEGHASRIVET="yes"
fi
fi
if test "x$THEPEGHASRIVET" == "xno" ; then
CREATE_RIVET="# create"
AC_MSG_RESULT([not found])
else
CREATE_RIVET="create"
AC_MSG_RESULT([found])
fi
AM_CONDITIONAL(HAVE_RIVET,[test x$THEPEGHASRIVET = xyes])
AC_SUBST([CREATE_RIVET])
])
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_LANG_PUSH([Fortran])
AC_MSG_CHECKING([checking if fortran compiler compiles argument mismatches])
AC_COMPILE_IFELSE(AC_LANG_SOURCE([[
program temp
call a(1.0D0)
end program
subroutine a(b)
double complex b
end]]),
[AC_MSG_RESULT([yes])],
[oldFCFLAGS="$FCFLAGS"
FCFLAGS="-std=legacy"
AC_MSG_CHECKING([checking if fortran compiler compiles argument mismatches with -std=legacy])
AC_COMPILE_IFELSE(AC_LANG_SOURCE([[
program temp
double precision b
double complex c
b = 1.0D0
c =1.0D0
call a(b)
call a(c)
end program
subroutine a(b)
double complex b
end]]),
[AC_MSG_RESULT([yes])
AM_FCFLAGS="$AM_FCFLAGS -std=legacy"],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([fortran compiler won't compile LoopTools])])
FCFLAGS="$oldFCFLAGS -std=legacy"]
)
AC_MSG_CHECKING([checking if fortran compiler compiles long lines])
AC_COMPILE_IFELSE(AC_LANG_SOURCE([[
program temp
write (*,*) 'aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa'
end program]]),
[AC_MSG_RESULT([yes])],
[AC_MSG_RESULT([no])
oldFCFLAGS="$FCFLAGS"
FCFLAGS="-ffixed-line-length-none"
AC_MSG_CHECKING([checking if fortran compiler compiles long lines with -ffixed-line-length-none])
AC_COMPILE_IFELSE(AC_LANG_SOURCE([[
program temp
write (*,*) 'aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa'
end program]]),
[AC_MSG_RESULT([yes])
AM_FCFLAGS="$AM_FCFLAGS -ffixed-line-length-none"
FCFLAGS="$oldFCFLAGS -ffixed-line-length-none"],
[FCFLAGS="-132"
AC_MSG_CHECKING([checking if fortran compiler compiles long lines with -132])
AC_COMPILE_IFELSE(AC_LANG_SOURCE([[
program temp
write (*,*) 'aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa'
end program]]),
[AC_MSG_RESULT([yes])
AM_FCFLAGS="$AM_FCFLAGS -132"
FCFLAGS="$oldFCFLAGS -132"],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([fortran compiler won't compile LoopTools])])])
]
)
AC_LANG_POP([Fortran])
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$with_vbfnlo" != "xno" -a "x$have_vbfnlo" = "xno" ],
[AC_CHECK_FILES(
${with_vbfnlo}/lib/VBFNLO/libVBFNLO.dylib,
[have_vbfnlo=lib], [have_vbfnlo=no])])
AS_IF([test "x$with_vbfnlo" != "xno" -a "x$have_vbfnlo" = "xno" ],
[AC_CHECK_FILES(
${with_vbfnlo}/lib64/VBFNLO/libVBFNLO.dylib,
[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
- AC_REQUIRE([AC_PROG_SED])
- VBFNLOINCLUDE=${with_vbfnlo}/include
- AC_SUBST(VBFNLOINCLUDE)
- VBFNLOLIB=$(echo ${with_vbfnlo}/${have_vbfnlo}/VBFNLO | $SED -e 's%/\+%/%g')
- AC_SUBST(VBFNLOLIB)
- LOAD_VBFNLO="library"
- CREATE_VBFNLO="create"
- INSERT_VBFNLO="insert"
- SET_VBFNLO="set"
- DO_VBFNLO="do"
- MKDIR_VBFNLO="mkdir"
+# Check if VBFNLO is being used
+AS_IF([test "x$with_vbfnlo" != "xno"],
+[
+ # Check for VBFNLO version
+ AC_MSG_CHECKING([for VBFNLO version >= 3.0.0])
+ tmp_vbfnloversion=$(${with_vbfnlo}/bin/vbfnlo --version | awk '{print $2}')
+ AX_COMPARE_VERSION([${tmp_vbfnloversion}], [lt], [3.0.0],
+ [have_vbfnlo3=no],
+ [have_vbfnlo3=yes])
+ AC_MSG_RESULT([VBFNLO 3.0.0 or higher: $have_vbfnlo3])
+ have_vbfnlo3=no
+],[have_vbfnlo3=no])
+
+# Define conditional based on VBFNLO version
+AM_CONDITIONAL(HAVE_VBFNLO3, [test "x$have_vbfnlo3" = "xyes"])
+
+# Configure library paths and options only if VBFNLO library is specified
+if test "x$have_vbfnlo" = "xlib" -o "x$have_vbfnlo" = "xlib64"; then
+ AC_REQUIRE([AC_PROG_SED])
+ VBFNLOINCLUDE=${with_vbfnlo}/include
+ AC_SUBST(VBFNLOINCLUDE)
+ VBFNLOLIB=$(echo ${with_vbfnlo}/${have_vbfnlo}/VBFNLO | $SED -e 's%/\+%/%g')
+ AC_SUBST(VBFNLOLIB)
+ LOAD_VBFNLO="library"
+ CREATE_VBFNLO="create"
+ INSERT_VBFNLO="insert"
+ SET_VBFNLO="set"
+ DO_VBFNLO="do"
+ MKDIR_VBFNLO="mkdir"
else
- LOAD_VBFNLO="# library"
- CREATE_VBFNLO="# create"
- INSERT_VBFNLO="# insert"
- SET_VBFNLO="# set"
- DO_VBFNLO="# do"
- MKDIR_VBFNLO="# mkdir"
+ LOAD_VBFNLO="# library"
+ CREATE_VBFNLO="# create"
+ INSERT_VBFNLO="# insert"
+ SET_VBFNLO="# set"
+ DO_VBFNLO="# do"
+ MKDIR_VBFNLO="# mkdir"
fi
AC_SUBST([LOAD_VBFNLO])
AC_SUBST([CREATE_VBFNLO])
AC_SUBST([INSERT_VBFNLO])
AC_SUBST([SET_VBFNLO])
AC_SUBST([DO_VBFNLO])
AC_SUBST([MKDIR_VBFNLO])
])
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$with_njet" != "xno" -a "x$have_njet" = "xno" ],
[AC_CHECK_FILES(
${with_njet}/lib/libnjet2.dylib,
[have_njet=lib], [have_njet=no])])
AS_IF([test "x$with_njet" != "xno" -a "x$have_njet" = "xno" ],
[AC_CHECK_FILES(
${with_njet}/lib/libnjet3.so,
[have_njet=lib], [have_njet=no])])
AS_IF([test "x$with_njet" != "xno" -a "x$have_njet" = "xno" ],
[AC_CHECK_FILES(
${with_njet}/lib64/libnjet3.so,
[have_njet=lib], [have_njet=no])])
AS_IF([test "x$with_njet" != "xno" -a "x$have_njet" = "xno" ],
[AC_CHECK_FILES(
${with_njet}/lib/libnjet3.dylib,
[have_njet=lib], [have_njet=no])])
AS_IF([test "x$with_njet" != "xno" ],
[AC_CHECK_FILES(
${with_njet}/include/njet.h,
[njet_include=include], [njet_include=include/njet])])
AS_IF([test "x$have_njet" = "xlib"],
[NJETLIBPATH=${with_njet}/lib
AC_SUBST(NJETLIBPATH)
NJETINCLUDEPATH=${with_njet}/$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}/$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"])
AS_IF([test "x$with_njet" != "xno"],[AC_MSG_CHECKING([for Njet version])
tmp_njetversion=[$(${with_njet}/bin/njet.py --version 2>&1 | grep -oE '[0-9]{4}$')]
AS_IF([test -z "$tmp_njetversion"],
[tmp_njetversion=1023])
AC_MSG_RESULT([${tmp_njetversion}])
NJET_VERSION=${tmp_njetversion}
AC_SUBST(NJET_VERSION)
])
if test "x$have_njet" = "xlib" -o "x$have_njet" = "xlib64" ; then
LOAD_NJET="library"
CREATE_NJET="create"
INSERT_NJET="insert"
DO_NJET="do"
else
LOAD_NJET="# library"
CREATE_NJET="# create"
INSERT_NJET="# insert"
DO_NJET="# do"
fi
AC_SUBST([LOAD_NJET])
AC_SUBST([CREATE_NJET])
AC_SUBST([INSERT_NJET])
AC_SUBST([DO_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])])
AS_IF([test "x$with_gosam" != "xno"],
[AC_MSG_CHECKING([for GoSam version >= 2.0.4])
tmp_gosamversion=[$(${with_gosam}/bin/gosam.py --version | grep 'GoSam.*rev' | cut -d' ' -f2)]
AX_COMPARE_VERSION([${tmp_gosamversion}],[lt],[2.0.4],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([Herwig requires GoSam 2.0.4 or later, found ${tmp_gosamversion}])],
[AC_MSG_RESULT([yes])])])
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"
DO_GOSAM="do"
else
LOAD_GOSAM="# library"
CREATE_GOSAM="# create"
INSERT_GOSAM="# insert"
DO_GOSAM="# do"
fi
AC_SUBST([LOAD_GOSAM])
AC_SUBST([CREATE_GOSAM])
AC_SUBST([INSERT_GOSAM])
AC_SUBST([DO_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"
DO_OPENLOOPS="do"
MKDIR_OPENLOOPS="mkdir"
else
LOAD_OPENLOOPS="# library"
CREATE_OPENLOOPS="# create"
INSERT_OPENLOOPS="# insert"
SET_OPENLOOPS="# set"
DO_OPENLOOPS="# do"
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([DO_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"
DO_MADGRAPH="do"
else
LOAD_MADGRAPH="# library"
CREATE_MADGRAPH="# create"
INSERT_MADGRAPH="# insert"
SET_MADGRAPH="# set"
DO_MADGRAPH="# do"
fi
AC_SUBST([LOAD_MADGRAPH])
AC_SUBST([CREATE_MADGRAPH])
AC_SUBST([INSERT_MADGRAPH])
AC_SUBST([SET_MADGRAPH])
AC_SUBST([DO_MADGRAPH])
])
AC_DEFUN([HERWIG_CHECK_PYTHIA],
[
dnl check if a directory is specified for Pythia
AC_ARG_WITH(pythia,
[AC_HELP_STRING([--with-pythia=dir],
[Assume the given directory for Pythia])])
dnl search for the pythia-config script
if test "$with_pythia" = ""; then
AC_PATH_PROG(pythiaconfig, pythia8-config, no)
else
AC_PATH_PROG(pythiaconfig, pythia8-config, no, ${with_pythia}/bin)
fi
if test "${pythiaconfig}" = "no"; then
AC_MSG_CHECKING(Pythia)
AC_MSG_RESULT(no);
PYTHIA8LIB=
# $2
else
PYTHIA8DATA=`${pythiaconfig} --datadir`/xmldoc
PYTHIA8LIB="-L`${pythiaconfig} --libdir` -lpythia8"
fi
AC_SUBST(PYTHIA8DATA)
AC_SUBST(PYTHIA8LIB)
])
dnl CHECK PYTHIA END
dnl ##### EvtGen #####
AC_DEFUN([HERWIG_CHECK_EVTGEN],
[
AC_MSG_CHECKING([for evtgen])
AC_ARG_WITH([evtgen],
AS_HELP_STRING([--with-evtgen=DIR], [Installation path of EvtGen]),
[],
[with_evtgen=no]
)
AC_MSG_RESULT([$with_evtgen])
AS_IF([test "x$with_evtgen" != "xno"],
[AC_CHECK_FILES(
${with_evtgen}/lib/libEvtGenExternal.so,
[have_evtgen=lib], [have_evtgen=no])],
[have_evtgen=no])
AS_IF([test "x$with_evtgen" != "xno" -a "x$have_evtgen" = "xno"],
[AC_CHECK_FILES(
${with_evtgen}/lib64/libEvtGenExternal.so,
[have_evtgen=lib64], [have_evtgen=no])])
AS_IF([test "x$with_evtgen" != "xno" -a "x$have_evtgen" = "xno" ],
[AC_CHECK_FILES(
${with_evtgen}/lib/libEvtGenExternal.dylib,
[have_evtgen=lib], [have_evtgen=no])])
AS_IF([test "x$have_evtgen" = "xlib" -o "x$have_evtgen" = "xlib64" ],
[EVTGENPREFIX=${with_evtgen}
AC_SUBST(EVTGENPREFIX)
])
AS_IF([test "x$have_evtgen" = "xlib" -o "x$have_evtgen" = "xlib64" ],
[AC_CHECK_FILES(
${with_evtgen}/share/evt.pdl,
[EVTGENSHARE=${with_evtgen}/share], [EVTGENSHARE=${with_evtgen}/share/EvtGen])
AC_SUBST(EVTGENSHARE)])
AS_IF([test "x$with_evtgen" != "xno" -a "x$have_evtgen" = "xno"],
[AC_MSG_ERROR([EvtGen requested but not found])])
AC_SUBST([EVTGENINCLUDE],[-I$EVTGENPREFIX/include])
AM_CONDITIONAL(HAVE_EVTGEN,[test "x$have_evtgen" = "xlib" -o "x$have_evtgen" = "xlib64"])
if test "x$have_evtgen" = "xlib" ; then
LOAD_EVTGEN_DECAYS="read EvtGenBDecays.in"
LOAD_EVTGEN_DECAYER="read EvtGenDecayer.in"
EVTGENLIBS="-L$with_evtgen/lib -lEvtGen -lEvtGenExternal"
elif test "x$have_evtgen" = "xlib64" ; then
LOAD_EVTGEN_DECAYS="read EvtGenBDecays.in"
LOAD_EVTGEN_DECAYER="read EvtGenDecayer.in"
EVTGENLIBS="-L$with_evtgen/lib64 -lEvtGen -lEvtGenExternal"
else
LOAD_EVTGEN_DECAYS="read HerwigBDecays.in"
LOAD_EVTGEN_DECAYER="#read EvtGenDecayer.in"
EVTGENLIBS=""
fi
AC_SUBST([LOAD_EVTGEN_DECAYS])
AC_SUBST([LOAD_EVTGEN_DECAYER])
AC_SUBST([EVTGENLIBS])
])
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"
GSLPATH="$with_gsl"
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"
GSLPATH="$with_gsl"
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"
GSLPATH="$with_gsl"
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_SUBST(GSLPATH)
])
dnl ##### COMPILERFLAGS #####
AC_DEFUN([HERWIG_COMPILERFLAGS],
[
AC_REQUIRE([HERWIG_CHECK_GSL])
AC_REQUIRE([HERWIG_CHECK_THEPEG])
AC_REQUIRE([BOOST_REQUIRE])
AC_REQUIRE([AX_COMPILER_VENDOR])
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="-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
case "${ax_cv_cxx_compiler_vendor}" in
gnu)
AM_CXXFLAGS="-pedantic -Wall -W -Wno-use-after-free"
;;
clang)
AM_CXXFLAGS="-pedantic -Wall -Wno-overloaded-virtual -Wno-unused-function -Wno-unused-parameter"
dnl -Wno-unneeded-internal-declaration
;;
intel)
AM_CXXFLAGS="-strict-ansi -Wall -wd13000,1418,981,444,383,1599,1572,2259,980"
;;
esac
AC_SUBST(AM_CPPFLAGS)
AC_SUBST(AM_CFLAGS, ["$warnflags $debugflags"])
AC_SUBST(AM_CXXFLAGS,["$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
***
*** GoSam: $with_gosam
*** GoSam-Contrib: $with_gosam_contrib
*** MadGraph: $with_madgraph
*** njet: $with_njet ${NJET_VERSION}
*** OpenLoops: $with_openloops
*** VBFNLO: $with_vbfnlo
***
*** EvtGen: $with_evtgen
*** GSL: $with_gsl
*** boost: ${BOOST_CPPFLAGS:-system}
*** Fastjet: ${fjconfig}
***
*** Host: $host
*** CC: $CCSTRING
*** CXX: $CXXSTRING
*** FC: $FCSTRING
***
*** CXXFLAGS: $CXXFLAGS
*** FCFLAGS: $FCFLAGS
*****************************************************
_HW_EOF_
])
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:18 AM (1 d, 13 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4982979
Default Alt Text
(113 KB)

Event Timeline