Page MenuHomeHEPForge
Files F8308910

No OneTemporary
Actions

View Options

diff --git a/.hgignore b/.hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -1,84 +1,80 @@
Makefile$
Makefile\.in$
\.deps$
\.libs$
\.l[ao]$
\.so$
\.so\.
\.o$
~$
\.pyc$
\.codereplace$
\.orig$
\.tar\.(gz|bz2)$
^autom4te.cache$
^config.herwig$
^config.log$
^config.status$
^configure$
^include/Herwig\+\+$
^Config/config.guess$
^Config/config.h$
^Config/config.h.in$
^Config/config.sub$
^Config/depcomp$
^Config/compile$
^Config/install-sh$
^Config/missing$
^Config/stamp-h1$
^Config/ar-lib$
^Contrib/make_makefiles.sh$
^Doc/HerwigDefaults.in$
^Doc/refman-html$
^Doc/fixinterfaces.pl$
^Doc/refman.conf$
^Doc/refman.h$
^Doc/AllInterfaces.h$
^Doc/HerwigDefaults.rpo$
^Doc/Herwig\+\+-refman.tag$
^Doc/tagfileThePEG.tag$
^Doc/.*\.log$
^(src|Utilities)/version.tmp$
^(src|Utilities)/version.tmp.new$
^(src|Utilities)/versionstring.h$
^INSTALL$
^aclocal.m4$
^confdefs.h$
^conftest.c$
^conftest.err$
^include/done-all-links$
^libtool$
\.dirstamp$
^src/herwigopts.h$
^src/herwigopts.c$
^src/defaults/Analysis.in$
^src/herwig-config$
^src/.*\.(run|tex|out|log|rpo|spc|top|dump|dot|aux|pdf|ps|png|svg|hepmc|dvi)$
^src/Makefile-UserModules$
^lib/done-all-links$
^lib/apple-fixes$
^src/Herwig\+\+
^src/defaults/PDF.in$
^src/defaults/done-all-links$
^src/tune$
^src/Herwig$
^src/tests/.*\.(time|mult|Bmult|chisq)$
^Tests/.*/.*\.(top|ps|pyc|info|dat|pdf|png)$
^Tests/.*\.(top|ps|pyc|info|dat|pdf|png|log|out)$
^Tests/.*\.(top|run|tex|mult|Bmult|aida|yoda)$
^Tests/Rivet-.*$
^Tests/plots$
^Tests/.*index.html$
^Herwig\+\+\-
^Models/Feynrules/python/Makefile-FR$
-^MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc$
-^MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc$
^MatrixElement/Matchbox/External/MadGraph/mg2Matchbox.py$
-^MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc$
-^MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc$
^MatrixElement/Matchbox/Scales/MatchboxScale.cc$
^Utilities/Statistics/combineDistributions$
^Utilities/Statistics/combineRuns$
^Utilities/Statistics/makeDistributions$
^src/defaults/MatchboxDefaults.in$
^src/defaults/setup.gosam.in$
diff --git a/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc.in b/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc
rename from MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc.in
rename to MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc
--- a/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc.in
+++ b/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc
@@ -1,1013 +1,1035 @@
// -*- C++ -*-
//
// GoSamAmplitude.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the GoSamAmplitude class.
//
#include "GoSamAmplitude.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
#include <boost/progress.hpp>
#include <boost/filesystem.hpp>
#include <fstream>
#include <sstream>
#include <string>
#include <cstdlib>
#include <exception>
using namespace Herwig;
namespace bfs = boost::filesystem;
+#ifndef HERWIG_BINDIR
+#error Makefile.am needs to define HERWIG_BINDIR
+#endif
+#ifndef HERWIG_PKGDATADIR
+#error Makefile.am needs to define HERWIG_PKGDATADIR
+#endif
+#ifndef GOSAM_PREFIX
+#error Makefile.am needs to define GOSAM_PREFIX
+#endif
+
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
-GoSamAmplitude::GoSamAmplitude(){
- theCodeExists=false;
- isitDR=false;
- theFormOpt=true;
- theNinja=true;
- theHiggsEff=false;
- theAccuracyTarget=6;
- theMassiveLeptons=false;
- theLoopInducedOption=0;
- doneGoSamInit=false;
- doneGoSamInitRun=false;}
+GoSamAmplitude::GoSamAmplitude() :
+ theAccuracyTarget(6),theCodeExists(false),theFormOpt(true),theNinja(true),
+ theHiggsEff(false),theMassiveLeptons(false),theLoopInducedOption(0),
+ isitDR(false),doneGoSamInit(false),doneGoSamInitRun(false),
+ bindir_(HERWIG_BINDIR), pkgdatadir_(HERWIG_PKGDATADIR), GoSamPrefix_(GOSAM_PREFIX)
+{}
GoSamAmplitude::~GoSamAmplitude() {}
IBPtr GoSamAmplitude::clone() const {
return new_ptr(*this);
}
IBPtr GoSamAmplitude::fullclone() const {
return new_ptr(*this);
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void GoSamAmplitude::doinit() {
optionalContractFile() = name() + ".OLPContract.lh";
MatchboxOLPME::doinit();
doneGoSamInit = true;
}
void GoSamAmplitude::doinitrun() {
optionalContractFile() = name() + ".OLPContract.lh";
MatchboxOLPME::doinitrun();
doneGoSamInitRun = true;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
extern "C" void OLP_Start(const char*, int* i);
extern "C" void OLP_Polvec(double*, double*, double*);
extern "C" void OLP_SetParameter(char*, double*, double*, int*);
extern "C" void OLP_PrintParameter(char*);
extern "C" void OLP_EvalSubProcess2(int*, double*, double*, double*, double*);
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
bool GoSamAmplitude::startOLP(const map<pair<Process, int>, int>& procs) {
char char_cwd[256];
getcwd(char_cwd, sizeof(char_cwd));
string cwd = string(char_cwd);
string folderMatchboxBuild = factory()->buildStorage();
folderMatchboxBuild.erase(folderMatchboxBuild.begin());
// set all necessary path and file names
gosamPath = gosamPathInterface == "" ? cwd + folderMatchboxBuild + "GoSam" : gosamPathInterface;
// When transitioning to C++ 11 this length()-1 workaround can be replaced by string.back()
if (gosamPath.at(gosamPath.length()-1) != '/') gosamPath.append("/");
gosamSourcePath = gosamPath + "source/";
gosamInstallPath = gosamPath + "build/";
// create all the directories
if (!bfs::is_directory(gosamPath)){
try {
bfs::create_directory(gosamPath);
} catch (exception& e) {
cerr << "--------------------------------------------------------------------------------\n";
cerr << "The following exception occured:\n\n";
cerr << " " << e.what() << "\n\n";
cerr << " -> Please create the parent directory of\n";
cerr << " " << gosamPath << "\n";
cerr << " manually!\n";
cerr << "--------------------------------------------------------------------------------\n";
abort();
}
}
if (!bfs::is_directory(gosamSourcePath)) bfs::create_directory(gosamSourcePath);
if (!bfs::is_directory(gosamInstallPath)) bfs::create_directory(gosamInstallPath);
contractFileTitle = name() + ".OLPContract.lh";
contractFileName = gosamPath + "/" + contractFileTitle;
string orderFileName = gosamPath + "/" + name() + ".OLPOrder.lh";
// Set the path variable (plus file name) where to find the GoSam specific input file
gosamSetupInFileName = gosamSetupInFileNameInterface == "" ? gosamPath + "/setup.gosam.in" : gosamSetupInFileNameInterface;
// Use the python script GoSamHelper.py to make replacements in the GoSam
// specific input file at gosamSetupInFileName. If the GoSam input file
// does not exist yet at gosamSetupInFileName the python script will get
// it from src/defaults/ before making the replacements.
- string cmd = "python @prefix@/bin/GoSamHelper.py ";
+ string cmd = "python "+bindir_+"/GoSamHelper.py ";
cmd+=" --usrinfile="+gosamSetupInFileNameInterface;
cmd+=" --infile="+gosamSetupInFileName+".tbu";
- cmd+=" --definfile=@prefix@/share/Herwig++/defaults/setup.gosam.in";
+ cmd+=" --definfile="+pkgdatadir_+"/defaults/setup.gosam.in";
cmd+=" --formtempdir="+StringUtils::replace(gosamSourcePath, string("/"), string("\\/")); //@FORMTEMPDIR@
cmd+=" --reduction="+(theNinja ? string("ninja,golem95") : string("samurai,golem95")); //@REDUCTIONPROGRAMS@
cmd+=" --formopt="+(theFormOpt ? string("") : string(", noformopt")); //@FORMOPT@
cmd+=" --higgseff="+(theHiggsEff ? string("smehc") : string("smdiag")); //@MODEL@
std::system(cmd.c_str());
generator()->log() << "\n\n>>> NOTE: According to the repository settings for the GoSam interface:\n" << flush;
if (theHiggsEff) generator()->log() << "\n -- GoSam will use a model with an effective ggH coupling (model=smehc).\n" << flush;
else if (!theHiggsEff) generator()->log() << "\n -- GoSam will use its default model (model=smdiag).\n" << flush;
if (theNinja) generator()->log() << " -- GoSam will use Ninja as reduction program (reduction_programs=ninja,golem95).\n" << flush;
else if (!theNinja) generator()->log() << " -- GoSam will use Samurai as reduction program (reduction_programs=samurai,golem95).\n" << flush;
if (theFormOpt) generator()->log() << " -- Form optimization switched on (extensions=autotools).\n" << flush;
else if (!theFormOpt) generator()->log() << " -- Form optimization switched off (extensions=autotools, noformopt).\n" << flush;
if (theNinja && !theFormOpt) throw Exception() << "\n\n>>> NOTE: Ninja reduction needs form optimization!\n" << Exception::abortnow;
if (gosamSetupInFileNameInterface == "") {
generator()->log() << "\n Please be aware that you are using a copy of the default GoSam input file!\n"
<< " Please note that if you need special options to be considered for the specific\n"
<< " process you are looking at (diagram filtering, etc.) these are not automatically\n"
<< " set for you. In that case please consider to specify your own GoSam input file\n"
<< " via 'set " << name() << ":SetupInFilename' in the input file.\n\n" << flush;
}
// If one uses a custom GoSam input file at gosamSetupInFileName = gosamSetupInFileNameInterface
// then please note that not all options in there might match the corresponding Herwig repository
// options
if (gosamSetupInFileNameInterface != "") {
generator()->log() << "\n Please be aware that you are using a custom GoSam input file!\n"
<< " Please note that if you have set the options for model, reduction_programs,\n"
<< " extensions and/or form.tempdir manually these will of course not be replaced\n"
<< " by the corresponding repository settings mentioned above.\n\n" << flush;
}
generator()->log() << "\n>>> NOTE: GoSam may return the set of used parameters for this process via the OLP_PrintParameter() function:\n\n"
<< " -- If Debug::level > 1, the OLP parameters are being written to file: at " << factory()->runStorage() + name() + ".OLPParameters.lh.\n\n" << flush;
double accuracyTarget = 1.0/pow(10.0,accuracyTargetNegExp());
time_t rawtime;
time (&rawtime);
accuracyFileTitle = name() + ".OLPAccuracy.lh";
accuracyFile = factory()->buildStorage() + accuracyFileTitle;
ofstream accuracyFileStream;
if ( Debug::level > 1 ) {
accuracyFileStream.open(accuracyFile.c_str()); // Opening accuracyFile once here removes all previous content before the read step
accuracyFileStream << "\nFile to contain those PSPs for which GoSam evaluated one-loop interference terms or loop induced ME2s\n"
<< "with acc > target accuracy = " << accuracyTarget << ". Date/Time: " << ctime(&rawtime) << endl;
}
generator()->log() << "\n>>> NOTE: GoSam will return the accuracy of one-loop interference terms or loop induced ME2s\n"
<< " at every PSP via the BLHA2 acc parameter:\n\n"
<< " -- In cases where acc > 10^-AccuracyTarget = " << accuracyTarget << " the corresponding PSPs are being dis-\n"
<< " carded.\n"
<< " -- The default value for AccuracyTarget is 6, but you may consider setting it otherwise\n"
<< " via 'set " << name() << ":AccuracyTarget' in the input file.\n"
<< " -- Currently the value for AccuracyTarget is set to " << accuracyTargetNegExp() << ".\n"
<< " -- If Debug::level > 1, the discarded PSPs are being written to file: at " + accuracyFile << ".\n"
<< " -- If the amount of PSPs with acc > " << accuracyTarget << " is significant, please consider to re-evaluate\n"
<< " your process setup (accuracy target, masses, cuts, etc.)!\n\n\n" << flush;
// check for old order file and create it if it doesn't already exist
fillOrderFile(procs, orderFileName);
ifstream ifile(contractFileName.c_str());
if(!ifile){
signOLP(orderFileName, contractFileName);
}
if ( !checkOLPContract(contractFileName) ) {
throw Exception() << "failed to start GoSam" << Exception::abortnow;
}
if (!( DynamicLoader::load(gosamInstallPath+"/lib/libgolem_olp.so")
|| DynamicLoader::load(gosamInstallPath+"/lib64/libgolem_olp.so")
|| DynamicLoader::load(gosamInstallPath+"/lib/libgolem_olp.dylib")
|| DynamicLoader::load(gosamInstallPath+"/lib64/libgolem_olp.dylib"))) buildGoSam();
int status = -1;
startOLP(contractFileTitle, status);
if ( status != 1 ) return false;
return true;
}
void GoSamAmplitude::startOLP(const string& contract, int& status) {
string tempcontract = contract;
char char_cwd[256];
getcwd(char_cwd, sizeof(char_cwd));
string cwd = string(char_cwd);
string folderMatchboxBuild = factory()->buildStorage();
folderMatchboxBuild.erase(folderMatchboxBuild.begin());
gosamPath = gosamPathInterface == "" ? cwd + folderMatchboxBuild + "GoSam" : gosamPathInterface;
// When transitioning to C++ 11 this length()-1 workaround can be replaced by string.back()
if (gosamPath.at(gosamPath.length()-1) != '/') gosamPath.append("/");
if (!( DynamicLoader::load(gosamPath+"build/lib/libgolem_olp.so")
|| DynamicLoader::load(gosamPath+"build/lib64/libgolem_olp.so")
|| DynamicLoader::load(gosamPath+"build/lib/libgolem_olp.dylib")
|| DynamicLoader::load(gosamPath+"build/lib64/libgolem_olp.dylib")))
throw Exception() << "Failed to load GoSam. Please check the log file.\n"
<< Exception::abortnow;
tempcontract = gosamPath + tempcontract;
OLP_Start(tempcontract.c_str(), &status);
// hand over input parameters for EW scheme considered
int pStatus = 0;
double zero = 0.0;
if ( SM().ewScheme() == 0 || SM().ewScheme() == 6 ) { // EW/Scheme Default and EW/Scheme Independent
throw Exception() << "`Best value' schemes are not supported by GoSam"
<< Exception::abortnow;
} else if ( SM().ewScheme() == 4 ) { // EW/Scheme mW (uses mW,GF,sin2thetaW) seems not to be supported by GoSam
throw Exception() << "`mW' scheme is not supported by GoSam"
<< Exception::abortnow;
} else if ( SM().ewScheme() == 1 ) { // EW/Scheme GMuScheme (uses mW,mZ,GF)
double in1=getParticleData(ParticleID::Z0)->mass()/GeV;
double in2=getParticleData(ParticleID::Wplus)->mass()/GeV;
double in3=SM().fermiConstant()*GeV2;
OLP_SetParameter((char *)"mass(23)",&in1,&zero,&pStatus);
OLP_SetParameter((char *)"mass(24)",&in2,&zero,&pStatus);
OLP_SetParameter((char *)"Gf",&in3,&zero,&pStatus);
} else if ( SM().ewScheme() == 2 ) { // EW/Scheme alphaMZScheme (uses mW,mZ,alpha(mZ))
double in1=getParticleData(ParticleID::Z0)->mass()/GeV;
double in2=getParticleData(ParticleID::Wplus)->mass()/GeV;
double in3=SM().alphaEMMZ();
OLP_SetParameter((char *)"mass(23)",&in1,&zero,&pStatus);
OLP_SetParameter((char *)"mass(24)",&in2,&zero,&pStatus);
OLP_SetParameter((char *)"alpha",&in3,&zero,&pStatus);
} else if ( SM().ewScheme() == 3 ) { // EW/Scheme NoMass (uses alpha(mZ),GF,sin2thetaW)
double in1=SM().fermiConstant()*GeV2;
double in2=SM().alphaEMMZ();
double in3=SM().sin2ThetaW();
OLP_SetParameter((char *)"Gf",&in1,&zero,&pStatus);
OLP_SetParameter((char *)"alpha",&in2,&zero,&pStatus);
OLP_SetParameter((char *)"sw2",&in3,&zero,&pStatus);
} else if ( SM().ewScheme() == 5 ) { // EW/Scheme mZ (uses mZ,alphaEM,sin2thetaW)
double in1=getParticleData(ParticleID::Z0)->mass()/GeV;
double in2=SM().alphaEMMZ();
double in3=SM().sin2ThetaW();
OLP_SetParameter((char *)"mass(23)",&in1,&zero,&pStatus);
OLP_SetParameter((char *)"alpha",&in2,&zero,&pStatus);
OLP_SetParameter((char *)"sw2",&in3,&zero,&pStatus);
}
// hand over mass and width of the Higgs
double wH = getParticleData(25)->width()/GeV;
double mH = getParticleData(25)->mass()/GeV;
OLP_SetParameter((char*)"width(25)",&wH,&zero,&pStatus);
OLP_SetParameter((char*)"mass(25)",&mH,&zero,&pStatus);
// hand over initial input parameter for alphaS
double as = SM().alphaS();
OLP_SetParameter((char *)"alphaS", &as, &zero, &pStatus);
// fill massive Particle vector
if (massiveParticles.empty()) {
// with quark masses
for (int i=1; i<=6; ++i)
if (getParticleData(i)->mass()/GeV > 0.0) massiveParticles.push_back(i);
// with lepton masses
if (theMassiveLeptons && getParticleData(11)->mass()/GeV > 0.0) massiveParticles.push_back(11);
if (theMassiveLeptons && getParticleData(13)->mass()/GeV > 0.0) massiveParticles.push_back(13);
if (theMassiveLeptons && getParticleData(15)->mass()/GeV > 0.0) massiveParticles.push_back(15);
}
// hand over quark (and possibly lepton) masses and widths (iff massive)
if ( massiveParticles.size() != 0 ) {
for ( vector<int>::const_iterator mID = massiveParticles.begin(); mID != massiveParticles.end(); ++mID ) {
string mstr;
string wstr;
int mInt = *mID;
double mass=getParticleData(mInt)->mass()/GeV;
double width=getParticleData(mInt)->width()/GeV;
std::stringstream ss;
ss << mInt;
string str = ss.str();
mstr="mass("+str+")";
wstr="width("+str+")";
char * mchar = new char[mstr.size()+1];
char * wchar = new char[wstr.size()+1];
std::copy(mstr.begin(),mstr.end(),mchar);
std::copy(wstr.begin(),wstr.end(),wchar);
mchar[mstr.size()] = '\0';
wchar[wstr.size()] = '\0';
OLP_SetParameter( mchar, &mass, &zero, &pStatus );
OLP_SetParameter( wchar, &width, &zero, &pStatus );
delete[] mchar;
delete[] wchar;
// Nicer but not working properly:
// double mass=getParticleData(*mID)->mass()/GeV;
// double width=getParticleData(*mID)->width()/GeV;
// string mstr="mass("+static_cast<ostringstream*>(&(ostringstream()<<(*mID)))->str()+")";
// string wstr="width("+static_cast<ostringstream*>(&(ostringstream()<<(*mID)))->str()+")";
// cout<<"\n massiv "<<mstr;
//
// OLP_SetParameter((char *)&mstr,&mass, &zero, &pStatus );
// OLP_SetParameter((char *)&wstr,&width, &zero, &pStatus );
}
}
// Note: In the GoSam input file, the standard is to set the parameter
// 'symmetries' for quark families and lepton generations, which allow
// for flavour changing only between families/generations. If this pa-
// rameter is set, GoSam won't allow to set electron and muon mass and
// width via the interface. Also setting mass and width for the tau is
// not yet considered.
// print OLP parameters
if ( Debug::level > 1 ) {
string ppstr = factory()->runStorage() + name() + ".OLPParameters.lh";
OLP_PrintParameter(const_cast<char*>(ppstr.c_str()));
}
didStartOLP() = true;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void GoSamAmplitude::fillOrderFile(const map<pair<Process, int>, int>& procs, string orderFileName) {
for ( map<pair<Process, int>, int>::const_iterator p = procs.begin() ; p != procs.end() ; ++p ) {
std::stringstream Processstr;
std::stringstream Typestr;
Processstr << (*p).first.first.legs[0]->id() << " " << (*p).first.first.legs[1]->id() << " -> ";
for ( PDVector::const_iterator o = (*p).first.first.legs.begin() + 2 ; o != (*p).first.first.legs.end() ; ++o )
Processstr << (**o).id() << " ";
if ( (*p).first.second == ProcessType::treeME2 ) {
Typestr << "Tree";
} else if ( (*p).first.second == ProcessType::loopInducedME2 ) {
Typestr << "LoopInduced";
} else if ( (*p).first.second == ProcessType::colourCorrelatedME2 ) {
Typestr << "ccTree";
} else if ( (*p).first.second == ProcessType::spinColourCorrelatedME2 ) {
Typestr << "scTree";
} else if ( (*p).first.second == ProcessType::oneLoopInterference ) {
Typestr << "Loop";
}
gosamprocinfo pro = gosamprocinfo((*p).second, -1, Processstr.str(), Typestr.str());
pro.setOAs(p->first.first.orderInAlphaS);
pro.setOAew(p->first.first.orderInAlphaEW);
processmap[(*p).second] = pro;
}
ifstream oldOrderFileStream(orderFileName.c_str());
if (oldOrderFileStream){
oldOrderFileStream.close();
return;
}
ofstream orderFile(orderFileName.c_str());
int asPower = 100;
int minlegs = 100;
int maxlegs = -1;
int maxasPower = -1;
int aewPower = 100;
int maxaewPower = -1;
for ( map<pair<Process, int>, int>::const_iterator t = procs.begin() ; t != procs.end() ; ++t ) {
asPower = min(asPower, static_cast<int>(t->first.first.orderInAlphaS));
minlegs = min(minlegs, static_cast<int>(t->first.first.legs.size()));
maxlegs = max(maxlegs, static_cast<int>(t->first.first.legs.size()));
maxasPower = max(maxasPower, static_cast<int>(t->first.first.orderInAlphaS));
aewPower = min(aewPower, static_cast<int>(t->first.first.orderInAlphaEW));
maxaewPower = max(maxaewPower, static_cast<int>(t->first.first.orderInAlphaEW));
}
orderFile << "# OLP order file created by Herwig++/Matchbox for GoSam\n\n";
orderFile << "InterfaceVersion BLHA2\n";
orderFile << "MatrixElementSquareType CHsummed\n";
orderFile << "CorrectionType QCD\n";
orderFile << "IRregularisation " << (isDR() ? "DRED" : "CDR") << "\n";
// loop over quarks to check if they have non-zero masses
for (int i=1; i<=6; ++i) if (getParticleData(i)->mass()/GeV > 0.0) massiveParticles.push_back(i);
// check if leptons have non-zero masses (iff theMassiveLeptons==true)
if (theMassiveLeptons && getParticleData(11)->mass()/GeV > 0.0) massiveParticles.push_back(11);
if (theMassiveLeptons && getParticleData(13)->mass()/GeV > 0.0) massiveParticles.push_back(13);
if (theMassiveLeptons && getParticleData(15)->mass()/GeV > 0.0) massiveParticles.push_back(15);
if ( massiveParticles.size() != 0 ) {
orderFile << "MassiveParticles ";
for ( vector<int>::const_iterator mID = massiveParticles.begin(); mID != massiveParticles.end(); ++mID ) {
int mInt = *mID;
orderFile << mInt << " ";
}
orderFile << "\n";
}
orderFile << "\n";
vector < string > types;
types.push_back("Tree");
types.push_back("LoopInduced");
types.push_back("ccTree");
types.push_back("scTree");
types.push_back("Loop");
for ( int i = asPower ; i != maxasPower + 1 ; i++ ) {
for ( int j = aewPower ; j != maxaewPower + 1 ; j++ ) {
orderFile << "\nAlphasPower " << i << "\n";
orderFile << "AlphaPower " << j << "\n";
for ( vector<string>::iterator it = types.begin() ; it != types.end() ; it++ ) {
if ( *it == "LoopInduced" ) continue;
for ( map<int, gosamprocinfo>::iterator p = processmap.begin() ; p != processmap.end() ; ++p )
if ( (*p).second.Tstr() == *it && i == (*p).second.orderAs() && j == (*p).second.orderAew() ) {
orderFile << "\nAmplitudeType " << *it << "\n";
break;
}
for ( map<int, gosamprocinfo>::iterator p = processmap.begin() ; p != processmap.end() ; ++p )
if ( (*p).second.Tstr() == *it && i == (*p).second.orderAs() && j == (*p).second.orderAew() ) {
orderFile << (*p).second.Pstr() << "\n";
}
}
}
}
// Write out the loop induced processes separately
int asPowerLI = 100;
int aewPowerLI = 100;
for ( map<int, gosamprocinfo>::iterator p = processmap.begin() ; p != processmap.end() ; ++p ) {
if ( (*p).second.Tstr() != "LoopInduced" ) continue;
if ( (*p).second.orderAs() != asPowerLI || (*p).second.orderAew() != aewPowerLI ) {
asPowerLI = (*p).second.orderAs();
aewPowerLI = (*p).second.orderAew();
// At the moment GoSam requires for qcd loop induced processes the as coupling power
// which would correspond to an associated fictitious Born process
orderFile << "\nAlphasPower " << (asPowerLI-2) << "\n";
orderFile << "AlphaPower " << aewPowerLI << "\n";
orderFile << "\nAmplitudeType " << "LoopInduced" << "\n";
}
orderFile << (*p).second.Pstr() << "\n";
}
orderFile << flush;
}
void GoSamAmplitude::signOLP(const string& order, const string& contract) {
if(!theCodeExists){
char char_cwd[256];
getcwd(char_cwd, sizeof(char_cwd));
string cwd = string(char_cwd);
string folderMatchboxBuild = factory()->buildStorage();
folderMatchboxBuild.erase(folderMatchboxBuild.begin());
generator()->log() << "\n>>> generating GoSam amplitudes. This may take some time, please be patient.\n"
<< ">>> see " + cwd + folderMatchboxBuild + "gosam-amplitudes.log for details.\n" << flush;
- string cmd = "@GOSAMPREFIX@/bin/gosam.py --olp --output-file=" + contract + " --config=" +
+ string cmd = GoSamPrefix_+"/bin/gosam.py --olp --output-file=" + contract + " --config=" +
gosamSetupInFileName+".tbu" + " --destination=" + gosamSourcePath + " " + order + " > " + cwd + folderMatchboxBuild + "gosam-amplitudes.log 2>&1";
std::system(cmd.c_str());
- cmd = "python @prefix@/bin/GoSamHelper.py ";
+ cmd = "python "+bindir_+"/GoSamHelper.py ";
cmd += " --makelink ";
cmd += " --makelinkfrom=contract ";
cmd += " --makelinkto="+factory()->buildStorage() + name() + ".OLPContract.lh";
std::system(cmd.c_str());
}
}
bool GoSamAmplitude::checkOLPContract(string contractFileName) {
ifstream infile(contractFileName.c_str());
string line;
vector < string > contractfile;
while (std::getline(infile, line)) contractfile.push_back(line);
for ( map<int, gosamprocinfo>::iterator p = processmap.begin() ; p != processmap.end() ; p++ ) {
bool righttype = false;
for ( vector<string>::iterator linex = contractfile.begin() ; linex != contractfile.end() ; ++linex ) {
if ( (*linex).find("AmplitudeType ")!= std::string::npos ) {
if ( (*linex).find(" " + (*p).second.Tstr() + " ")!= std::string::npos ) {
righttype = true;
} else {
righttype = false;
}
}
if ( righttype ) {
if ( (*linex).find((*p).second.Pstr()) != std::string::npos && (*p).second.Pstr().length() == (*linex).find("|") ) {
string sub = (*linex).substr((*linex).find("|") + 1, (*linex).find("#") - (*linex).find("|") - 1); // | 1 23 # buggy??
if ( sub.find(" 1 ") != 0 )
throw Exception() << "Failed to check contractfile. Please check the logfile.\n"
<< Exception::abortnow;
string subx = sub.substr(3);
int subint;
istringstream(subx) >> subint;
(*p).second.setGID(subint);
}
}
}
}
string ids = factory()->buildStorage() + "GoSam.ids.dat";
ofstream IDS(ids.c_str());
idpair.clear();
for ( map<int, gosamprocinfo>::iterator p = processmap.begin() ; p != processmap.end() ; p++ )
idpair.push_back(-1);
idpair.push_back(-1);
for ( map<int, gosamprocinfo>::iterator p = processmap.begin() ; p != processmap.end() ; p++ ) {
idpair[(*p).second.HID()]=(*p).second.GID();
IDS << (*p).second.HID() << " " << (*p).second.GID() << " " << (*p).second.Tstr() << "\n";
if ( (*p).second.GID() == -1 ) return 0;
}
IDS << flush;
return 1;
}
bool GoSamAmplitude::buildGoSam() {
if(!theCodeExists){
generator()->log() << "\n>>> compiling GoSam amplitudes. This may take some time, please be patient.\n"
<< ">>> see " + gosamSourcePath + "gosam-build.log for details.\n\n" << flush;
string cmd = "cd " + gosamSourcePath + " && sh autogen.sh FCFLAGS=-g --prefix=" +
gosamInstallPath + " --disable-static > gosam-build.log 2>&1";
std::system(cmd.c_str());
if (!gosamBuildScript.empty()) {
cmd = "cd " + gosamSourcePath + " && " + gosamBuildScript + " >> gosam-build.log 2>&1";
std::system(cmd.c_str());
}
std::system(cmd.c_str());
cmd = "cd " + gosamSourcePath + " && make install >> gosam-build.log 2>&1";
std::system(cmd.c_str());
}
theCodeExists=true;
return 1;
}
void GoSamAmplitude::getids() const {
string line = factory()->buildStorage() + "GoSam.ids.dat";
ifstream infile(line.c_str());
int hid;
int gid;
string type;
while (std::getline(infile, line)) {
idpair.push_back(-1);
idtypepair.push_back(" ");
}
infile.close();
string line2 = factory()->buildStorage() + "GoSam.ids.dat";
ifstream infile2(line2.c_str());
idpair.push_back(-1);
idtypepair.push_back(" ");
while (std::getline(infile2, line2)) {
istringstream(line2) >> hid >> gid >> type;
idpair[hid]=gid;
idtypepair[hid]=type;
}
infile.close();
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void GoSamAmplitude::evalSubProcess() const {
useMe();
double units = pow(lastSHat() / GeV2, int(mePartonData().size()) - 4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double scale = sqrt(mu2() / GeV2);
if (hasRunningAlphaS()) {
int pStatus = 0;
double zero = 0.0;
double as;
as = lastAlphaS();
OLP_SetParameter((char *)"alphaS", &as, &zero, &pStatus);
}
double out[7] = { };
double acc;
if ( idpair.size() == 0 ){ getids(); }
int id = -99;
if ( olpId()[ProcessType::loopInducedME2] ) id = olpId()[ProcessType::loopInducedME2];
else if ( olpId()[ProcessType::oneLoopInterference] ) id = olpId()[ProcessType::oneLoopInterference];
else id = olpId()[ProcessType::treeME2];
int callid(idpair[id]); // If id denotes the Herwig ID, this returns the GoSam ID
string calltype(idtypepair[id]); // If id denotes the Herwig ID, this returns the amplitude type
OLP_EvalSubProcess2(&(callid), olpMomenta(), &scale, out, &acc);
double accuracyTarget = 1.0/pow(10.0,accuracyTargetNegExp());
accuracyFileTitle = name() + ".OLPAccuracy.lh";
accuracyFile = factory()->buildStorage() + accuracyFileTitle;
ofstream accuracyFileStream;
if ( Debug::level > 1 ) {
accuracyFileStream.open(accuracyFile.c_str(),ios::app);
}
if ( (olpId()[ProcessType::oneLoopInterference]||olpId()[ProcessType::loopInducedME2]) && acc > accuracyTarget ) {
if ( Debug::level > 1 ) {
vector<Lorentz5Momentum> currentpsp = lastXComb().meMomenta();
time_t rawtime;
time (&rawtime);
if (doneGoSamInit) accuracyFileStream << "READ phase: ";
else if (doneGoSamInitRun) accuracyFileStream << "RUN phase: ";
accuracyFileStream << "Sub-process with Herwig ID = " << id << " and GoSam ID = " << callid << ", " << ctime(&rawtime);
accuracyFileStream << "GoSam evaluated one-loop interference or loop induced ME2 with acc = " << acc
<< " > target accuracy = " << accuracyTarget << ", at PSP [in units of GeV]:" << endl;
for (size_t i=0; i!=currentpsp.size(); ++i) {
accuracyFileStream << "(t,x,y,z,mass;m)[" << i << "]=("
<< currentpsp[i].t()/GeV << ","
<< currentpsp[i].x()/GeV << ","
<< currentpsp[i].y()/GeV << ","
<< currentpsp[i].z()/GeV << ","
<< currentpsp[i].mass()/GeV << ";"
<< currentpsp[i].m()/GeV << ")"
<< endl;
}
accuracyFileStream << endl;
}
throw Veto(); // Dispose of PSP
}
if ( olpId()[ProcessType::oneLoopInterference] ) {
if (calculateTreeME2()) lastTreeME2(out[3] * units);
lastOneLoopInterference((out[2])* units);
lastOneLoopPoles(pair<double, double>(out[0] * units, out[1] * units));
} else if ( olpId()[ProcessType::treeME2] ) {
lastTreeME2(out[3] * units);
} else if ( olpId()[ProcessType::loopInducedME2] ) {
lastTreeME2(out[2] * units);
}
}
void GoSamAmplitude::evalColourCorrelator(pair<int, int> ) const {
double units = pow(lastSHat() / GeV2, int(mePartonData().size()) - 4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double scale = sqrt(mu2() / GeV2);
if (hasRunningAlphaS()) {
int pStatus = 0;
double zero = 0.0;
double as;
as = lastAlphaS();
OLP_SetParameter((char *)"alphaS", &as, &zero, &pStatus);
}
int n = lastXComb().meMomenta().size();
colourCorrelatorResults.resize(n * (n - 1) / 2);
if ( idpair.size() == 0 ) getids();
int callid(idpair[olpId()[ProcessType::colourCorrelatedME2]]);
double acc;
OLP_EvalSubProcess2(&(callid), olpMomenta(), &scale, &colourCorrelatorResults[0], &acc);
cPDVector particles = lastXComb().matrixElement()->mePartonData();
for ( int i = 0 ; i < n ; ++i ) {
for ( int j = i + 1 ; j < n ; ++j ) {
lastColourCorrelator(make_pair(i, j), colourCorrelatorResults[i+j*(j-1)/2] * units);
}
}
}
void GoSamAmplitude::evalSpinColourCorrelator(pair<int , int > ) const {
double units = pow(lastSHat() / GeV2, int(mePartonData().size()) - 4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double scale = sqrt(mu2() / GeV2);
if (hasRunningAlphaS()) {
int pStatus = 0;
double zero = 0.0;
double as;
as = lastAlphaS();
OLP_SetParameter((char *)"alphaS", &as, &zero, &pStatus);
}
int n = lastXComb().meMomenta().size();
spinColourCorrelatorResults.resize(2*n*n);
if ( idpair.size() == 0 ) getids();
double acc;
int callid(idpair[olpId()[ProcessType::spinColourCorrelatedME2]]);
OLP_EvalSubProcess2(&(callid), olpMomenta(), &scale, &spinColourCorrelatorResults[0], &acc);
for ( int i = 0; i < n; ++i ) {
for ( int j = 0; j < n; ++j ) {
Complex scc(spinColourCorrelatorResults[2*i+2*n*j]*units, spinColourCorrelatorResults[2*i+2*n*j+1]*units);
lastColourSpinCorrelator(make_pair(i,j),scc);
}
}
}
LorentzVector<Complex> GoSamAmplitude::plusPolarization(const Lorentz5Momentum& p, const Lorentz5Momentum& n, int inc) const {
double pvec[4] = {p.t()/GeV,p.x()/GeV,p.y()/GeV,p.z()/GeV};
double nvec[4] = {n.t()/GeV,n.x()/GeV,n.y()/GeV,n.z()/GeV};
double out[8] ={ };
OLP_Polvec(pvec,nvec,out);
LorentzVector<Complex> res;
Complex a(out[0],out[1]);
res.setT(a);
Complex b(out[2],out[3]);
res.setX(b);
Complex c(out[4],out[5]);
res.setY(c);
Complex d(out[6],out[7]);
res.setZ(d);
if (inc<2)
return res.conjugate();
else
return res;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void GoSamAmplitude::persistentOutput(PersistentOStream & os) const {
os << idpair << idtypepair << processmap << gosamPathInterface
<< gosamSetupInFileNameInterface << gosamBuildScript << gosamPath
<< gosamSourcePath << gosamInstallPath << gosamSetupInFileName
<< orderFileTitle << contractFileTitle
<< contractFileName << orderFileName
<< theCodeExists << theFormOpt << theNinja << isitDR << massiveParticles << theHiggsEff
<< theAccuracyTarget << theMassiveLeptons << theLoopInducedOption
<< doneGoSamInit << doneGoSamInitRun;
}
void GoSamAmplitude::persistentInput(PersistentIStream & is, int) {
is >> idpair >> idtypepair >> processmap >> gosamPathInterface
>> gosamSetupInFileNameInterface >> gosamBuildScript >> gosamPath
>> gosamSourcePath >> gosamInstallPath >> gosamSetupInFileName
>> orderFileTitle >> contractFileTitle
>> contractFileName >> orderFileName
>> theCodeExists >> theFormOpt >> theNinja >> isitDR >> massiveParticles >> theHiggsEff
>> theAccuracyTarget >> theMassiveLeptons >> theLoopInducedOption
>> doneGoSamInit >> doneGoSamInitRun;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<GoSamAmplitude, MatchboxOLPME> describeHerwigGoSamAmplitude("Herwig::GoSamAmplitude", "HwMatchboxGoSam.so");
void GoSamAmplitude::Init() {
static ClassDocumentation<GoSamAmplitude> documentation("GoSamAmplitude implements an interface to GoSam.", "Matrix elements have been calculated using GoSam");
static Parameter<GoSamAmplitude,string> interfaceProcessPath
("ProcessPath",
"Prefix for the process source code, include files and library produced by GoSam.",
&GoSamAmplitude::gosamPathInterface, "",
false, false);
static Parameter<GoSamAmplitude,string> interfaceSetupInFilename
("SetupInFilename",
"File name of the GoSam infile (typically setup.gosam.in) to be used. If left empty a new setup.gosam.in is created in the location specified in Path",
&GoSamAmplitude::gosamSetupInFileNameInterface, "",
false, false);
static Switch<GoSamAmplitude,bool> interfaceCodeExists
("CodeExists",
"Switch on or off if Code already exists/not exists.",
&GoSamAmplitude::theCodeExists, true, false, false);
static SwitchOption interfaceCodeExistsOn
(interfaceCodeExists,
"True",
"Switch True if Code already exists.",
true);
static SwitchOption interfaceCodeExistsOff
(interfaceCodeExists,
"False",
"Switch False if Code has to be build.",
false);
static Switch<GoSamAmplitude,bool> interfaceisitDR
("isDR",
"Switch on or off DR.",
&GoSamAmplitude::isitDR, false, false, false);
static SwitchOption interfaceisitDROn
(interfaceisitDR,
"True",
"Switch True.",
true);
static SwitchOption interfaceisitDROff
(interfaceisitDR,
"False",
"Switch False.",
false);
static Switch<GoSamAmplitude,bool> interfaceFormOpt
("FormOpt",
"Switch On/Off formopt",
&GoSamAmplitude::theFormOpt, true, false, false);
static SwitchOption interfaceFormOptOn
(interfaceFormOpt,
"On",
"On",
true);
static SwitchOption interfaceFormOptOff
(interfaceFormOpt,
"Off",
"Off",
false);
static Switch<GoSamAmplitude,bool> interfaceNinja
("Ninja",
"Switch On/Off for reduction with Ninja. If Off then Samurai is used.",
&GoSamAmplitude::theNinja, true, false, false);
static SwitchOption interfaceNinjaOn
(interfaceNinja,
"On",
"On",
true);
static SwitchOption interfaceNinjaOff
(interfaceNinja,
"Off",
"Off",
false);
static Switch<GoSamAmplitude,bool> interfaceHiggsEff
("HiggsEff",
"Switch On/Off for effective higgs model.",
&GoSamAmplitude::theHiggsEff, false, false, false);
static SwitchOption interfaceHiggsEffOn
(interfaceHiggsEff,
"On",
"On",
true);
static SwitchOption interfaceHiggsEffOff
(interfaceHiggsEff,
"Off",
"Off",
false);
static Parameter<GoSamAmplitude,string> interfaceBuildScript
("BuildScript",
"File name of a custom build script, which is called between 'autogen.sh'"
"and 'make install'. It can be used for parallelization.",
&GoSamAmplitude::gosamBuildScript, "",
false, false);
static Parameter<GoSamAmplitude,int> interfaceAccuracyTarget
("AccuracyTarget",
"Integer to parametrize the threshold value for the BLHA2 acc parameter, returned by GoSam in the case of "
"sub-processes with one-loop intereference terms or loop induced sub-processes."
"If acc > 10^-AccuracyTarget the corresponding PSP is being discarded. Discarded PSPs are written to file "
"if Debug::level > 1.",
&GoSamAmplitude::theAccuracyTarget, 6, 0, 0,
false, false, Interface::lowerlim);
static Switch<GoSamAmplitude,bool> interfaceMassiveLeptons
("MassiveLeptons",
"If set to Yes, then pass on the light lepton masses - as well as the tau mass - to GoSam."
"Otherwise GoSam will use light leptons of zero mass as default, as well as its own default tau mass.",
&GoSamAmplitude::theMassiveLeptons, false, false, false);
static SwitchOption interfaceMassiveLeptonsNo
(interfaceMassiveLeptons,
"No",
"No",
false);
static SwitchOption interfaceMassiveLeptonsYes
(interfaceMassiveLeptons,
"Yes",
"Yes",
true);
static Switch<GoSamAmplitude,int> interfaceLoopInducedOption
("LoopInducedOption",
"Options for the GoSam interface, in the case that a loop induced process is being considered. The default "
"option is 0, for which only the squared one-loop amplitude in the Standard Model is being considered. All "
"other options consider additional contributions from a model with an effective interaction, which lead to "
"the same final state, such as the squared effective amplitude, or the interference term between the one- "
"loop amplitude in the Standard Model and the effective amplitude, or any additive combinations therefrom. "
"In order to use those options an appropriate model has to be used.",
&GoSamAmplitude::theLoopInducedOption, 0, false, false);
static SwitchOption interfaceLoopInducedOptionLI2
(interfaceLoopInducedOption,
"LI2",
"Only consider the squared one-loop amplitude in the Standard Model.",
0);
static SwitchOption interfaceLoopInducedOptionEff2
(interfaceLoopInducedOption,
"Eff2",
"Only consider the squared effective amplitude.",
1);
static SwitchOption interfaceLoopInducedOptionLIEffInterference
(interfaceLoopInducedOption,
"LIEffInterference",
"Only consider the interference term between the one-loop amplitude "
"in the Standard Model and the effective amplitude.",
2);
static SwitchOption interfaceLoopInducedOptionLI2plusEff2
(interfaceLoopInducedOption,
"LI2plusEff2",
"Consider the sum of the squared one-loop amplitude in the Standard "
"Model plus the squared effective amplitude.",
3);
static SwitchOption interfaceLoopInducedOptionLI2plusLIEffInterference
(interfaceLoopInducedOption,
"LI2plusEffInterference",
"Consider the sum of the squared one-loop amplitude in the Standard "
"Model plus the interference term between the one-loop amplitude in "
"the Standard Model and the effective amplitude.",
4);
static SwitchOption interfaceLoopInducedOptionEff2plusLIEffInterference
(interfaceLoopInducedOption,
"Eff2plusEffInterference",
"Consider the sum of the squared effective amplitude plus the inter- "
"ference term between the one-loop amplitude in the Standard Model "
"and the effective amplitude.",
5);
static SwitchOption interfaceLoopInducedOptionAllAdditions
(interfaceLoopInducedOption,
"AllAdditions",
"Consider the sum of the squared one-loop amplitude in the Standard "
"Model plus all other contributions, which come with the effective "
"Model.",
6);
+
+ static Parameter<GoSamAmplitude,string> interfaceBinDir
+ ("BinDir",
+ "The location for the installed executable",
+ &GoSamAmplitude::bindir_, string(HERWIG_BINDIR),
+ false, false);
+ static Parameter<GoSamAmplitude,string> interfacePKGDATADIR
+ ("DataDir",
+ "The location for the installed Herwig++ data files",
+ &GoSamAmplitude::pkgdatadir_, string(HERWIG_PKGDATADIR),
+ false, false);
+
+ static Parameter<GoSamAmplitude,string> interfaceGoSamPrefix
+ ("GoSamPrefix",
+ "The prefix for the location of GoSam",
+ &GoSamAmplitude::GoSamPrefix_, string(GOSAM_PREFIX),
+ false, false);
}
diff --git a/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.h b/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.h
--- a/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.h
+++ b/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.h
@@ -1,312 +1,329 @@
// -*- C++ -*-
//
// GoSamAmplitude.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_GoSamAmplitude_H
#define Herwig_GoSamAmplitude_H
//
// This is the declaration of the GoSamAmplitude class.
//
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxOLPME.h"
#include "ThePEG/Utilities/DynamicLoader.h"
namespace Herwig {
using namespace ThePEG;
class gosamprocinfo{
public:
gosamprocinfo(){};
gosamprocinfo(int HID,int GID, string procstr,string typestr):
theHOlpId(HID),theGOlpId(GID),theProcstr(procstr),theTypestr(typestr){
}
~gosamprocinfo(){}
int HID() const {return theHOlpId;}
int GID() const {return theGOlpId;}
string Pstr() const {return theProcstr;}
string Tstr() const {return theTypestr;}
void setGID(int g){theGOlpId=g;}
void setOAs(int i){ orderAlphas=i;}
int orderAs(){return orderAlphas;}
void setOAew(int j){ orderAlphaew=j;}
int orderAew(){return orderAlphaew;}
private:
int theHOlpId;
int theGOlpId;
string theProcstr;
string theTypestr;
int orderAlphas;
int orderAlphaew;
public:
void persistentOutput(PersistentOStream & os) const{os<<theHOlpId<<theGOlpId<<theProcstr<<theTypestr<<orderAlphas<<orderAlphaew;}
void persistentInput(PersistentIStream &is) {is>>theHOlpId>>theGOlpId>>theProcstr>>theTypestr>>orderAlphas>>orderAlphaew;}
};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief GoSamAmplitude implements an interface to GoSam
*/
class GoSamAmplitude: public MatchboxOLPME {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
GoSamAmplitude();
/**
* The destructor.
*/
virtual ~GoSamAmplitude();
//@}
public:
virtual void fillOrderFile(const map<pair<Process,int>,int>& procs, string OrderFileName);
virtual bool isCS() const { return false; }
virtual bool isExpanded() const { return true; }
virtual bool isBDK() const { return false; }
virtual bool isDR() const { return isitDR; }
virtual bool isDRbar() const {return false;}
/**
* Start the one loop provider, if appropriate, giving order and
* contract files
*/
virtual void signOLP(const string&, const string&);
virtual bool checkOLPContract(string contractFileName);
/**
* Return true, if this amplitude already includes symmetry factors
* for identical outgoing particles.
*/
// virtual bool hasFinalStateSymmetry() const { return true; }
virtual bool hasFinalStateSymmetry() const { return false; }
virtual bool buildGoSam();
/**
* 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);
/**
* 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;
/**
* Return a positive helicity polarization vector for a gluon of
* momentum p (with reference vector n) to be used when evaluating
* spin correlations.
*/
virtual LorentzVector<Complex> plusPolarization(const Lorentz5Momentum& p,
const Lorentz5Momentum& n, int inc) const;
/**
* Fill in results for the given colour/spin correlator
*/
virtual void evalSpinColourCorrelator(pair<int,int> ij) const;
void getids() const;
int accuracyTargetNegExp() const { return theAccuracyTarget; }
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 Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
GoSamAmplitude & operator=(const GoSamAmplitude &);
/**
* 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 gosam
*/
mutable vector<int> idpair;
/**
* first is the olp id from herwig, second the amplitude type
*/
mutable vector<string> idtypepair;
/**
* Map to store all processes handled by this Amplitude
*/
map<int , gosamprocinfo > processmap;
mutable string gosamPathInterface;
mutable string gosamSetupInFileNameInterface;
mutable string gosamBuildScript;
mutable string gosamPath;
mutable string gosamSourcePath;
mutable string gosamInstallPath;
mutable string gosamSetupInFileName;
mutable string orderFileTitle;
mutable string contractFileTitle;
mutable string contractFileName;
mutable string orderFileName;
mutable string accuracyFileTitle;
mutable string accuracyFile;
int theAccuracyTarget;
bool theCodeExists;
bool theFormOpt;
bool theNinja;
bool theHiggsEff;
bool theMassiveLeptons;
int theLoopInducedOption;
bool isitDR;
mutable bool doneGoSamInit;
mutable bool doneGoSamInitRun;
/**
* The PDG codes of those quarks with mass
*/
vector<int> massiveParticles; //theMassiveParticles;
/**
* Method to create the setup.in file for GoSam
*/
void setupGoSamIn(string setupGoSamInFile);
+protected:
+
+ /**
+ * Location of the installed executables
+ */
+ string bindir_;
+
+ /**
+ * Location of the data files
+ */
+ string pkgdatadir_;
+
+ /**
+ * Location of GOSAM
+ */
+ string GoSamPrefix_;
+
}; // end "class GoSamAmplitude: public MatchboxOLPME"
inline PersistentOStream& operator<<(PersistentOStream& os, const gosamprocinfo& p) {
p.persistentOutput(os); return os;
}
inline PersistentIStream& operator>>(PersistentIStream& is, gosamprocinfo& p) {
p.persistentInput(is); return is;
}
} // end "namespace Herwig"
#endif /* Herwig_GoSamAmplitude_H */
diff --git a/MatrixElement/Matchbox/External/GoSam/Makefile.am b/MatrixElement/Matchbox/External/GoSam/Makefile.am
--- a/MatrixElement/Matchbox/External/GoSam/Makefile.am
+++ b/MatrixElement/Matchbox/External/GoSam/Makefile.am
@@ -1,13 +1,18 @@
pkglib_LTLIBRARIES = HwMatchboxGoSam.la
HwMatchboxGoSam_la_LDFLAGS = -module -version-info 11:1:0
HwMatchboxGoSam_la_SOURCES = \
GoSamAmplitude.h GoSamAmplitude.cc
+HwMatchboxGoSam_la_CPPFLAGS = $(AM_CPPFLAGS) \
+-DHERWIG_BINDIR="\"$(bindir)\"" \
+-DHERWIG_PKGDATADIR="\"$(pkgdatadir)\"" \
+-DGOSAM_PREFIX="\"$(GOSAMPREFIX)\""
+
install-exec-local:
$(install_sh_SCRIPT) $(srcdir)/GoSamHelper.py $(DESTDIR)$(bindir)/GoSamHelper.py
uninstall-local:
rm -f $(DESTDIR)$(bindir)/GoSamHelper.py
diff --git a/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc.in b/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc
rename from MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc.in
rename to MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc
--- a/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc.in
+++ b/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc
@@ -1,800 +1,811 @@
// -*- C++ -*-
//
// MadGraphAmplitude.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MadGraphAmplitude class.
//
#include "MadGraphAmplitude.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
-
#include "ThePEG/PDT/EnumParticles.h"
#include <boost/lexical_cast.hpp>
#include <boost/filesystem.hpp>
-
-
#include <cstdlib>
#include <dlfcn.h>
#include <errno.h>
#include <sstream>
using namespace Herwig;
-
+#ifndef HERWIG_BINDIR
+#error Makefile.am needs to define HERWIG_BINDIR
+#endif
+#ifndef HERWIG_PKGDATADIR
+#error Makefile.am needs to define HERWIG_PKGDATADIR
+#endif
+#ifndef MADGRAPH_PREFIX
+#error Makefile.am needs to define MADGRAPH_PREFIX
+#endif
extern "C" void mginitproc_(char *i,int);
extern "C" void MG_Calculate_wavefunctions_virt(int* proc,double*,double*);
extern "C" void MG_Calculate_wavefunctions_born(int* proc,double*, int*);
extern "C" void MG_Jamp (int* proc,int*, double*);
extern "C" void MG_LNJamp (int* proc,int*, double*);
extern "C" void MG_Virt (int* proc,double*);
extern "C" void MG_NCol (int* proc,int*);
extern "C" void MG_vxxxxx (double* p,double* n,int* inc,double* );
extern "C" void MG_Colour (int* proc,int* i,int* j ,int* color);
-
-
-
MadGraphAmplitude::MadGraphAmplitude()
: theMGmodel("loop_sm"),keepinputtopmass(false),
- prefix_("@prefix@"), madgraphPrefix_("@MADGRAPHPREFIX@")
- {}
+ bindir_(HERWIG_BINDIR), pkgdatadir_(HERWIG_PKGDATADIR), madgraphPrefix_(MADGRAPH_PREFIX)
+{}
MadGraphAmplitude::~MadGraphAmplitude() {
}
IBPtr MadGraphAmplitude::clone() const {
return new_ptr(*this);
}
IBPtr MadGraphAmplitude::fullclone() const {
return new_ptr(*this);
}
bool MadGraphAmplitude::initializedMad=false;
vector<string> MadGraphAmplitude::BornAmplitudes=vector<string>();
vector<string> MadGraphAmplitude::VirtAmplitudes=vector<string>();
void MadGraphAmplitude::initProcess(const cPDVector& ) {
if ( lastMatchboxXComb()->initialized() )
return;
if ( !DynamicLoader::load(mgProcLibPath()+"InterfaceMadGraph.so") )
throw Exception() << "Failed to load MadGraph amplitudes\n"
<< DynamicLoader::lastErrorMessage
<< Exception::abortnow;
if (!initializedMad){
string mstr=(factory()->runStorage()+"MadGraphAmplitudes"+"/param_card"+((theMGmodel=="loop_sm")?"":("_"+theMGmodel))+".dat");
size_t len = mstr.size();
mginitproc_(const_cast<char*>(mstr.c_str()),len);
initializedMad=true;
}
lastMatchboxXComb()->isInitialized();
}
bool MadGraphAmplitude::writeAmplitudesDat(){
bool res=false;
string born= mgProcLibPath()+"BornAmplitudes.dat";
if ( !boost::filesystem::exists(born) ) {
ofstream borns(born.c_str());
for (vector<string>::iterator amps=BornAmplitudes.begin();amps!=BornAmplitudes.end();amps++)
borns<<*amps<<endl;
borns.close();
res=true;
}
string virt= mgProcLibPath()+"VirtAmplitudes.dat";
if ( !boost::filesystem::exists(virt) ) {
ofstream virts(virt.c_str());
for (vector<string>::iterator amps=VirtAmplitudes.begin();amps!=VirtAmplitudes.end();amps++)
virts<<*amps<<endl;
virts.flush();
virts.close();
res=true;
}
return res;
}
bool MadGraphAmplitude::checkAmplitudes(){
string born= mgProcLibPath()+"BornAmplitudes.dat";
string virt= mgProcLibPath()+"VirtAmplitudes.dat";
assert ( boost::filesystem::exists(born)|| boost::filesystem::exists(virt));
bool foundallborns=true;
for (vector<string>::iterator amps=BornAmplitudes.begin();amps!=BornAmplitudes.end();amps++){
ifstream borns(born.c_str());
string line;
bool foundthisborn=false;
while (std::getline(borns, line)) {
if(line==*amps)foundthisborn=true;
}
foundallborns&=foundthisborn;
}
bool foundallvirts=true;
for (vector<string>::iterator amps=VirtAmplitudes.begin();amps!=VirtAmplitudes.end();amps++){
ifstream virts(virt.c_str());
string line;
bool foundthisvirt=false;
while (std::getline(virts, line)) {
if(line==*amps)foundthisvirt=true;
}
foundallvirts&=foundthisvirt;
}
if (!foundallborns||!foundallvirts)
throw Exception() << "One amplitude has no externalId. Please remove the MadGraphAmplitude-folder and rebuild.\n" << Exception::abortnow;
return foundallborns && foundallvirts;
}
string MadGraphAmplitude::mgProcLibPath(){
string res=theProcessPath == "" ? factory()->buildStorage()+"MadGraphAmplitudes" : theProcessPath;
if (res.at(res.length()-1) != '/') res.append("/");
return res;
}
bool MadGraphAmplitude::initializeExternal() {
if ( boost::filesystem::exists(mgProcLibPath()) ) {
if ( !boost::filesystem::is_directory(mgProcLibPath()) )
throw Exception() << "MadGraph amplitude storage '"
<< mgProcLibPath() << "' existing but not a directory."
<< Exception::abortnow;
} else {
boost::filesystem::create_directories(mgProcLibPath());
}
string runAmplitudes = factory()->runStorage() + "/MadGraphAmplitudes";
if ( boost::filesystem::exists(runAmplitudes) ) {
if ( !boost::filesystem::is_directory(runAmplitudes) )
throw Exception() << "MadGraph amplitude storage '"
<< runAmplitudes << "' existing but not a directory."
<< Exception::abortnow;
} else {
boost::filesystem::create_directories(runAmplitudes);
}
//EW-consistency check:
Energy MW=getParticleData(ParticleID::Wplus)->mass();
Energy MZ=getParticleData(ParticleID::Z0)->mass();
if( MW!= sqrt(MZ*MZ/2.0+sqrt(MZ*MZ*MZ*MZ/4.0-Constants::pi*SM().alphaEMMZ()*MZ*MZ/ sqrt(2.0)/SM().fermiConstant()))){
cerr<<"\n\n-----!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!-----";
cerr << "\nYou are using a EW scheme which is inconsistent with the MadGraph parametisation:\n\n"
<<MW/GeV<< " GeV==MW!= sqrt(MZ^2/2+sqrt(MZ^4/4.0-pi*alphaEMMZ*MZ^2/ sqrt(2)/G_f))=="<<
sqrt(MZ*MZ/2.0+sqrt(MZ*MZ*MZ*MZ/4.0-Constants::pi*SM().alphaEMMZ()*MZ*MZ/ sqrt(2.0)/SM().fermiConstant()))/GeV
<<" GeV\n\n-----!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!-----\n";
}
string para= factory()->runStorage()+"/MadGraphAmplitudes"+"/MG-Parameter.dat";
ofstream params(para.c_str());
params<<"$WZ$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Z0)->width() /GeV;
params<<"\n$WW$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Wplus)->width()/GeV;
params<<"\n$alphas$ " <<std::setiosflags(ios::scientific) <<SM().alphaS();
params<<"\n$GF$ " <<std::setiosflags(ios::scientific) <<SM().fermiConstant()*GeV2 ;
params<<"\n$alphaMZ$ " <<std::setiosflags(ios::scientific) <<1/SM().alphaEMMZ();
params<<"\n$MZ$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Z0)->mass() /GeV<<flush;
params<<"\n$MW$ " <<std::setiosflags(ios::scientific) <<getParticleData(ParticleID::Wplus)->mass() /GeV<<flush;
params<<"\n$sw2$ " <<std::setiosflags(ios::scientific) << SM().sin2ThetaW() <<flush;
if(theMGmodel=="heft"&&!keepinputtopmass){
cerr<<"\n---------------------------------------------------------------";
cerr<<"\n---------------------------------------------------------------";
cerr<<"\nNote: You are using the Higgs Effective model (heft) in ";
cerr<<"\n Madgraph. We assume you try to calculate NLO with ";
cerr<<"\n the GoSam virtual amplitudes. To match the models we ";
cerr<<"\n therefore set the topmass to 10000000 GeV.";
cerr<<"\n\n For more information see the \\tau parameter in:";
cerr<<"\n https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Models/HiggsEffective";
cerr<<"\n\n The Effective Higgs model in Gosam is using mT=infinity";
cerr<<"\n\n\n If you want to use the LO matrixelements of MadGraph with finite' topmass you need to add: ";
cerr<<"\n\n set Madgraph:KeepInputTopMass True";
cerr<<"\n\n to your input file.";
cerr<<"\n---------------------------------------------------------------";
cerr<<"\n---------------------------------------------------------------\n";
params<<"\n$MT$ 10000000." <<flush;
}else{
params<<"\n$MT$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::t)->mass() /GeV <<flush;
}
params<<"\n$WT$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::t)->width() /GeV <<flush;
params<<"\n$MB$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::b)->mass() /GeV <<flush;
params<<"\n$MH$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::h0)->mass() /GeV <<flush;
params<<"\n$WH$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::h0)->width() /GeV <<flush;
params<<"\n$MTA$ " <<std::setiosflags(ios::scientific) << getParticleData(ParticleID::tauplus)->mass() /GeV <<flush;
- string cmd = "python " + prefix_ + "/bin/mg2Matchbox.py ";
+ string cmd = "python " + bindir_ + "/mg2Matchbox.py ";
cmd +=" --buildpath "+mgProcLibPath();
cmd +=" --model "+theMGmodel;
cmd +=" --runpath "+factory()->runStorage()+"/MadGraphAmplitudes " ;
+ cmd +=" --datadir "+pkgdatadir_;
std::stringstream as,aem;
as << factory()->orderInAlphaS();
cmd +=" --orderas "+as.str() ;
aem <<factory()->orderInAlphaEW();
cmd +=" --orderew "+aem.str();
// TODO move to boost::system
writeAmplitudesDat();
if (boost::filesystem::exists(mgProcLibPath()+"InterfaceMadGraph.so") ){
//set the parameters
checkAmplitudes();
std::system(cmd.c_str());
ranMadGraphInitializeExternal = true;
return true;
}
char cwd[1024];
if ( !getcwd(cwd,sizeof(cwd)) )
throw Exception() << "failed to determine current working directory\n"
<< Exception::abortnow;
cmd +=" --madgraph " + madgraphPrefix_ + "/bin " ;
cmd +="--build > ";
cmd += mgProcLibPath()+"MG.log 2>&1";
generator()->log() << "\n\nCompiling MadGraph amplitudes. This may take some time -- please be patient.\n"
<< "In case of problems see " << mgProcLibPath() << "MG.log for details.\n\n"
<< flush;
std::system(cmd.c_str());
- cmd = "python " + prefix_ + "/bin/mg2Matchbox.py ";
+ cmd = "python " + bindir_ + "/mg2Matchbox.py ";
cmd +=" --buildpath "+mgProcLibPath();
cmd +=" --model "+theMGmodel;
cmd +=" --runpath "+factory()->runStorage()+"/MadGraphAmplitudes " ;
as.clear();
aem.clear();
as << factory()->orderInAlphaS();
cmd +=" --orderas "+as.str() ;
aem <<factory()->orderInAlphaEW();
cmd +=" --orderew "+aem.str();
std::system(cmd.c_str());
ranMadGraphInitializeExternal = true;
return boost::filesystem::exists(mgProcLibPath()+"InterfaceMadGraph.so");
}
int MadGraphAmplitude::externalId(const cPDVector& proc) {
for (int i=0;i<100;i++){
colourindex.push_back(-2);
}
assert(!BornAmplitudes.empty()||!VirtAmplitudes.empty());
writeAmplitudesDat();
int res=0;
string amp="";
int k=0;
for (cPDVector::const_iterator it=proc.begin();it!=proc.end();it++,k++){
amp+=boost::lexical_cast<string>( (*it)->id())+" ";if (k==1)amp+=" > ";
}
string born= mgProcLibPath()+"BornAmplitudes.dat";
string virt= mgProcLibPath()+"VirtAmplitudes.dat";
assert ( boost::filesystem::exists(born)|| boost::filesystem::exists(virt));
ifstream borns(born.c_str());
string line;
while (std::getline(borns, line)) {
res+=1;
if(line==amp)return res;
}
ifstream virts(virt.c_str());
while (std::getline(virts, line)) {
res+=1;
if(line==amp)return res;
}
throw Exception() << "One amplitude has no externalId. Please remove the MadGraphAmplitude-folder and rebuild.\n" << Exception::abortnow;
return res;
}
bool MadGraphAmplitude::ranMadGraphInitializeExternal = false;
void MadGraphAmplitude::doinit() {
if ( !ranMadGraphInitializeExternal ) {
initializeExternal();
}
MatchboxAmplitude::doinit();
}
void MadGraphAmplitude::doinitrun() {
if ( !ranMadGraphInitializeExternal ) {
initializeExternal();
}
MatchboxAmplitude::doinitrun();
}
bool MadGraphAmplitude::canHandle(const PDVector& p,
Ptr<MatchboxFactory>::tptr factory,
bool virt) const {
if ( factory->processData()->diagramMap().find(p) !=
factory->processData()->diagramMap().end() )
return true;
vector<Ptr<Tree2toNDiagram>::ptr> diags =
factory->diagramGenerator()->generate(p,orderInGs(),orderInGem());
if ( diags.empty() )
return false;
factory->processData()->diagramMap()[p] = diags;
string amp="";
int k=0;
for (PDVector::const_iterator it=p.begin();it!=p.end();it++,k++){
amp+=boost::lexical_cast<string>( (*it)->id())+" ";if (k==1)amp+=" > ";
}
if (virt && factory->highestVirt()>=p.size()){
VirtAmplitudes.push_back(amp);
}else{
BornAmplitudes.push_back(amp);
}
return true;
}
void MadGraphAmplitude::prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr me) {
useMe();
if ( !calculateTreeAmplitudes() ) {
MatchboxAmplitude::prepareAmplitudes(me);
return;
}
if (colourindex.empty()) {
for (int i=0;i<100;i++){
colourindex.push_back(-2);
}
}
lastMatchboxXComb()->clearheljamp();
lastMatchboxXComb()->clearhelLNjamp();
initProcess(mePartonData());
MatchboxAmplitude::prepareAmplitudes(me);
}
Complex MadGraphAmplitude::evaluate(size_t i, const vector<int>& hel, Complex& largeN) {
//find the colourline:
int ii = -1;
int xx=lastMatchboxXComb()->externalId();
if(colourindex[i]!=-2){
ii = colourindex[i];
if (ii==-1) {
largeN = Complex(0.0);
return Complex(0.0);
}
} else {
set<vector<size_t> > a = colourOrdering(i);
int ncol=-1;
MG_NCol(&xx,&ncol);
assert(ncol!=-1);
for( int it = 0; it < ncol; it++ ){
int n = 0;
for ( cPDVector::const_iterator nx = mePartonData().begin();
nx != mePartonData().end(); nx++ )
if ( (*nx)->coloured() ) n++;
set<vector<size_t> > tmpset;
vector<size_t> tmpvek;
for ( int it2 = 0; it2 < n; it2++ ) {
int ret=-2;
MG_Colour(&xx,&it,&it2,&ret);
assert(ret !=-2);
if (ret== -1)
break;
if ( ret == 0 ) {
n++;
tmpset.insert(tmpvek);
tmpvek.clear();
} else {
tmpvek.push_back(ret-1);
}
if( it2 == n-1 ) tmpset.insert(tmpvek);
}
bool found_all = true;
for ( set<vector<size_t> >::iterator it3 = a.begin(); it3 != a.end(); it3++ ) {
bool found_it3=false;
for ( set<vector<size_t> >::iterator it4 = tmpset.begin(); it4 != tmpset.end(); it4++ ) {
vector<size_t> it3tmp = gluonsFirst(*it3);
vector<size_t> it4tmp = (*it4);
if ( it3tmp.size() != it4tmp.size() ) continue;
if ( it3tmp == it4tmp ) found_it3 = true;
}
found_all = found_all && found_it3;
}
if ( found_all ) {
colourindex[i]=it;
ii=it;
}
}
}
if ( ii == -1 ){
colourindex[i]=ii;
largeN = Complex(0.0);
return Complex(0.0);
}
const map<vector<int>,vector < complex<double> > >& tmp = lastMatchboxXComb()->heljamp();
const map<vector<int>,vector < complex<double> > >& tmpLN = lastMatchboxXComb()->helLNjamp();
if( tmp.find(hel) != tmp.end()) {
largeN = tmpLN.find(hel)->second[ii];
return tmp.find(hel)->second[ii];;
}
double units = pow(sqrt(lastSHat())/GeV,int(hel.size())-4);
int heltmp[10];
for(size_t j=0;j<hel.size();j++){
int cross=crossingMap()[j];
if( (cross>1&&j<=1)||(cross<=1&&j>1)){
heltmp[cross]=-1*hel[j];}
else{heltmp[cross]=hel[j];}
}
vector<Lorentz5Momentum> reshuffled = meMomenta();
if ( !reshuffleMasses().empty() && reshuffled.size() > 3 ) {
const cPDVector& pdata = mePartonData();
const map<long,Energy>& masses = reshuffleMasses();
lastMatchboxXComb()->reshuffle(reshuffled,pdata,masses);
}
double momenta[50];
size_t j=0;
for (size_t i=0;i<mePartonData().size();i++){
momenta[j]=abs(reshuffled[i].e()/GeV)<1.e-13?0.:double(reshuffled[i].e()/GeV);
momenta[j+1]=abs(reshuffled[i].x()/GeV)<1.e-13?0.:double(reshuffled[i].x()/GeV);
momenta[j+2]=abs(reshuffled[i].y()/GeV)<1.e-13?0.:double(reshuffled[i].y()/GeV);
momenta[j+3]=abs(reshuffled[i].z()/GeV)<1.e-13?0.:double(reshuffled[i].z()/GeV);
j+=4;
}
MG_Calculate_wavefunctions_born(&xx, &momenta[0], &heltmp[0]);
int ncol=-1;
MG_NCol(&xx,&ncol);
Complex res;
Complex resLN;
for( int it = 0; it < ncol; it++ ){
double dd[2];
MG_Jamp(&xx,&it,&dd[0]);
Complex d(dd[0],dd[1]);
if(it==ii)res=d*units;
lastMatchboxXComb()->pushheljamp(hel,d*units);
double ddLN[2];
MG_LNJamp(&xx,&it,&ddLN[0]);
Complex dLN(ddLN[0],ddLN[1]);
if(it==ii)resLN=dLN*units;
lastMatchboxXComb()->pushhelLNjamp(hel,dLN*units);
}
largeN = resLN;
return res;
}
vector<unsigned int> MadGraphAmplitude::physicalHelicities(const vector<int>& hel) const {
vector<unsigned int> res(hel.size(),0);
for ( size_t j = 0; j < hel.size(); ++j ) {
int cross = crossingMap()[j];
int xhel = 0;
if ( (cross > 1 && j <= 1) || (cross <= 1 && j > 1) )
xhel = -1*hel[j];
else
xhel = hel[j];
if ( mePartonData()[cross]->iSpin() == PDT::Spin1Half )
res[cross] = (xhel == -1 ? 0 : 1);
else if ( mePartonData()[cross]->iSpin() == PDT::Spin1 )
res[cross] = (unsigned int)(xhel + 1);
else assert(false);
}
return res;
}
LorentzVector<Complex> MadGraphAmplitude::plusPolarization(const Lorentz5Momentum& p,
const Lorentz5Momentum& n,
int i) const {
int tmp=i;
double pg[4],ng[4],poltmp[8];
pg[0]=p.e()/GeV;pg[1]=p.x()/GeV;pg[2]=p.y()/GeV;pg[3]=p.z()/GeV;
ng[0]=n.e()/GeV;ng[1]=n.x()/GeV;ng[2]=n.y()/GeV;ng[3]=n.z()/GeV;
MG_vxxxxx(&pg[0],&ng[0],&tmp,&poltmp[0]);
complex<double> pol[6];
pol[0]=Complex(poltmp[0],poltmp[1]);
pol[1]=Complex(poltmp[2],poltmp[3]);
pol[2]=Complex(poltmp[4],poltmp[5]);
pol[3]=Complex(poltmp[6],poltmp[7]);
LorentzVector<Complex> polarization(pol[1],pol[2],pol[3],pol[0]);
return polarization.conjugate();
}
bool equalsModulo(unsigned int i, const vector<int>& a, const vector<int>& b) {
assert(a.size()==b.size());
if ( a[i] == b[i] )
return false;
for ( unsigned int k = 0; k < a.size(); ++k ) {
if ( k == i )
continue;
if ( a[k] != b[k] )
return false;
}
return true;
}
vector<size_t> MadGraphAmplitude::gluonsFirst(vector<size_t> vec) {
vector<size_t> vecout;
for(vector<size_t>::iterator it= vec.begin();it!= vec.end();++it)
if ( mePartonData()[crossingMap()[*it]]->id()==21)
vecout.push_back(crossingMap()[*it]);
for(vector<size_t>::iterator it= vec.begin();it!= vec.end();++it)
if ( mePartonData()[crossingMap()[*it]]->id()!=21)
vecout.push_back(crossingMap()[*it]);
return vecout;
}
double MadGraphAmplitude::spinColourCorrelatedME2(pair<int,int> ij,
const SpinCorrelationTensor& c) const {
vector<Lorentz5Momentum> reshuffled = meMomenta();
if ( !reshuffleMasses().empty() && reshuffled.size() > 3 ) {
const cPDVector& pdata = mePartonData();
const map<long,Energy>& masses = reshuffleMasses();
lastMatchboxXComb()->reshuffle(reshuffled,pdata,masses);
}
Lorentz5Momentum p = reshuffled[ij.first];
Lorentz5Momentum n = reshuffled[ij.second];
LorentzVector<Complex> polarization = plusPolarization(p,n,ij.first);
Complex pFactor = (polarization*c.momentum())/sqrt(abs(c.scale()));
double avg =
colourCorrelatedME2(ij)*(-c.diagonal()+ (c.scale() > ZERO ? 1. : -1.)*norm(pFactor));
int iCrossed = crossingMap()[ij.first];
iCrossed = ij.first;
for ( unsigned int k = 0; k < crossingMap().size(); ++k )
if ( crossingMap()[k] == ij.first ) {
iCrossed = k;
break;
}
Complex csCorr = 0.0;
if ( calculateColourSpinCorrelator(ij) ) {
set<const CVector*> done;
for ( AmplitudeConstIterator a = lastAmplitudes().begin();
a != lastAmplitudes().end(); ++a ) {
if ( done.find(&(a->second)) != done.end() )
continue;
AmplitudeConstIterator b = lastAmplitudes().begin();
while ( !equalsModulo(iCrossed,a->first,b->first) )
if ( ++b == lastAmplitudes().end() )
break;
if ( b == lastAmplitudes().end() || done.find(&(b->second)) != done.end() )
continue;
done.insert(&(a->second)); done.insert(&(b->second));
if ( a->first[iCrossed] == 1 )
swap(a,b);
csCorr -= colourBasis()->colourCorrelatedInterference(ij,mePartonData(),a->second,b->second);
}
lastColourSpinCorrelator(ij,csCorr);
} else {
csCorr = lastColourSpinCorrelator(ij);
}
double corr =
2.*real(csCorr*sqr(pFactor));
double Nc = generator()->standardModel()->Nc();
double cfac = 1.;
if ( mePartonData()[ij.first]->iColour() == PDT::Colour8 ) {
cfac = Nc;
} else if ( mePartonData()[ij.first]->iColour() == PDT::Colour3 ||
mePartonData()[ij.first]->iColour() == PDT::Colour3bar ) {
cfac = (sqr(Nc)-1.)/(2.*Nc);
} else assert(false);
return
( avg +(c.scale() > ZERO ? 1. : -1.)*corr/cfac);
}
void MadGraphAmplitude::prepareOneLoopAmplitudes(Ptr<MatchboxMEBase>::tcptr ){
assert(false);
}
double MadGraphAmplitude::oneLoopInterference() const {
if ( !calculateOneLoopInterference() )
return lastOneLoopInterference();
evaloneLoopInterference();
return lastOneLoopInterference();
}
void MadGraphAmplitude::evaloneLoopInterference() const {
double units = pow(lastSHat()/GeV2,int(mePartonData().size())-4);
vector<Lorentz5Momentum> reshuffled = meMomenta();
if ( !reshuffleMasses().empty() && reshuffled.size() > 3 ) {
const cPDVector& pdata = mePartonData();
const map<long,Energy>& masses = reshuffleMasses();
lastMatchboxXComb()->reshuffle(reshuffled,pdata,masses);
}
double virt[20];
double momenta[50];
size_t j=0;
for (size_t i=0;i<mePartonData().size();i++){
momenta[j]=abs(reshuffled[i].e()/GeV)<1.e-13?0.:double(reshuffled[i].e()/GeV);
momenta[j+1]=abs(reshuffled[i].x()/GeV)<1.e-13?0.:double(reshuffled[i].x()/GeV);
momenta[j+2]=abs(reshuffled[i].y()/GeV)<1.e-13?0.:double(reshuffled[i].y()/GeV);
momenta[j+3]=abs(reshuffled[i].z()/GeV)<1.e-13?0.:double(reshuffled[i].z()/GeV);
j+=4;
}
int xx=lastMatchboxXComb()->externalId();
MG_Calculate_wavefunctions_virt(&xx,&momenta[0],&virt[0]);
double ifact = 1.;
ifact = 1./4.;
if (lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour3 ||
lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour3bar )
ifact /= SM().Nc();
else if ( lastMatchboxXComb()->matchboxME()->mePartonData()[0]->iColour() == PDT::Colour8 )
ifact /= (SM().Nc()*SM().Nc()-1.);
if ( lastMatchboxXComb()->matchboxME()->mePartonData()[1]->iColour() == PDT::Colour3 ||
lastMatchboxXComb()->matchboxME()->mePartonData()[1]->iColour() == PDT::Colour3bar )
ifact /= SM().Nc();
else if ( mePartonData()[1]->iColour() == PDT::Colour8 )
ifact /= (SM().Nc()*SM().Nc()-1.);
ifact *= lastMatchboxXComb()->matchboxME()->finalStateSymmetry();
lastOneLoopInterference(virt[1]/ifact*units);
lastOneLoopPoles(pair<double, double>(virt[2]/ifact*units,virt[3]/ifact*units));
}
void MadGraphAmplitude::persistentOutput(PersistentOStream & os) const {
os << theOrderInGs << theOrderInGem << BornAmplitudes << VirtAmplitudes
- << colourindex<<crossing << theProcessPath << theMGmodel << prefix_ << madgraphPrefix_;
+ << colourindex<<crossing << theProcessPath << theMGmodel << bindir_
+ << pkgdatadir_ << madgraphPrefix_;
}
void MadGraphAmplitude::persistentInput(PersistentIStream & is, int) {
is >> theOrderInGs >> theOrderInGem >> BornAmplitudes >> VirtAmplitudes
- >> colourindex>>crossing >> theProcessPath >> theMGmodel >> prefix_ >> madgraphPrefix_;
+ >> colourindex>>crossing >> theProcessPath >> theMGmodel >> bindir_
+ >> pkgdatadir_ >> madgraphPrefix_;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<MadGraphAmplitude,MatchboxAmplitude>
describeHerwigMadGraphAmplitude("Herwig::MadGraphAmplitude", "HwMatchboxMadGraph.so");
void MadGraphAmplitude::Init() {
static ClassDocumentation<MadGraphAmplitude> documentation ("MadGraphAmplitude","Matrix elements have been calculated using MadGraph5");
static Parameter<MadGraphAmplitude,string> interfaceProcessPath
("ProcessPath",
"The Process Path.",
&MadGraphAmplitude::theProcessPath, "",false, false);
static Parameter<MadGraphAmplitude,string> interfaceModel
("Model",
"The MadGraph-Model.",
&MadGraphAmplitude::theMGmodel, "loop_sm",false, false);
static Switch<MadGraphAmplitude,bool> interfacekeepinputtopmass
("KeepInputTopMass",
"Switch On/Off formopt",
&MadGraphAmplitude::keepinputtopmass, false, false, false);
static SwitchOption interfacekeepinputtopmassTrue
(interfacekeepinputtopmass,
"On",
"On",
true);
static SwitchOption interfacekeepinputtopmassFalse
(interfacekeepinputtopmass,
"Off",
"Off",
false);
- static Parameter<MadGraphAmplitude,string> interfacePrefix
- ("Prefix",
- "The prefix for the installed location of the code",
- &MadGraphAmplitude::prefix_, "@prefix@",
+ static Parameter<MadGraphAmplitude,string> interfaceBinDir
+ ("BinDir",
+ "The location for the installed executable",
+ &MadGraphAmplitude::bindir_, string(HERWIG_BINDIR),
+ false, false);
+
+ static Parameter<MadGraphAmplitude,string> interfacePKGDATADIR
+ ("DataDir",
+ "The location for the installed Herwig++ data files",
+ &MadGraphAmplitude::pkgdatadir_, string(HERWIG_PKGDATADIR),
false, false);
static Parameter<MadGraphAmplitude,string> interfaceMadgraphPrefix
("MadgraphPrefix",
"The prefix for the location of MadGraph",
- &MadGraphAmplitude::madgraphPrefix_, "@MADGRAPHPREFIX@",
+ &MadGraphAmplitude::madgraphPrefix_, string(MADGRAPH_PREFIX),
false, false);
}
diff --git a/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.h b/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.h
--- a/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.h
+++ b/MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.h
@@ -1,352 +1,357 @@
// -*- C++ -*-
//
// MadGraphAmplitude.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_MadGraphAmplitude_H
#define Herwig_MadGraphAmplitude_H
//
// This is the declaration of the MadGraphAmplitude class.
//
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
#include "Herwig++/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxCurrents.h"
#include "ThePEG/Utilities/DynamicLoader.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Johannes Bellm, Simon Platzer
*
* \brief MadGraphAmplitude implements an interface to MadGraph
*/
class MadGraphAmplitude:
public MatchboxAmplitude {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MadGraphAmplitude();
/**
* The destructor.
*/
virtual ~MadGraphAmplitude();
//@}
public:
/**
* Return true, if this amplitude can handle the given process.
*/
virtual bool canHandle(const PDVector&,
Ptr<MatchboxFactory>::tptr,
bool) const;
/**
* Return true, if this amplitude already includes symmetry factors
* for identical outgoing particles.
*/
virtual bool hasFinalStateSymmetry() const { return false; }
/**
* Set the (tree-level) order in \f$g_S\f$ in which this matrix
* element should be evaluated.
*/
virtual void orderInGs(unsigned int ogs) { theOrderInGs = ogs; }
/**
* Return the (tree-level) order in \f$g_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInGs() const { return theOrderInGs; }
/**
* Set the (tree-level) order in \f$g_{EM}\f$ in which this matrix
* element should be evaluated.
*/
virtual void orderInGem(unsigned int oge) { theOrderInGem = oge; }
/**
* Return the (tree-level) order in \f$g_{EM}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInGem() const { return theOrderInGem; }
/**
* Return true, if this amplitude is capable of calculating one-loop
* (QCD) corrections.
*/
virtual bool haveOneLoop() const { return true; }
/**
* Calculate the tree level amplitudes for the phasespace point
* stored in lastXComb.
*/
virtual void prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr);
/**
* Calculate the one-loop amplitudes for the phasespace point
* stored in lastXComb, if provided.
*/
virtual void prepareOneLoopAmplitudes(Ptr<MatchboxMEBase>::tcptr);
/**
* Evaluate the amplitude for the given colour tensor id and
* helicity assignment
*/
virtual Complex evaluate(size_t, const vector<int>&, Complex&);
/**
* Return true, if one-loop contributions will be evaluated at amplitude level.
*/
virtual bool oneLoopAmplitudes() const { return false; }
/**
* Return the one-loop/tree interference.
*/
virtual double oneLoopInterference() const;
void evaloneLoopInterference() const;
/**
* Return true, if one loop corrections are given in the conventions
* of BDK.
*/
virtual bool isCS() const { return false; }
virtual bool isExpanded() const { return true; }
virtual bool isBDK() const { return false; }
virtual bool isDR() const { return false; }
virtual bool isDRbar() const {return false;}
/**
* Return the value of the dimensional regularization
* parameter. Note that renormalization scale dependence is fully
* restored in DipoleIOperator.
*/
virtual Energy2 mu2() const { return sqr(91.18800000000*GeV);}
//virtual Energy2 mu2() const { return lastSHat(); }
virtual LorentzVector<Complex> plusPolarization(const Lorentz5Momentum& p,
const Lorentz5Momentum& n,
int id = -1) const;
/**
* Return the colour and spin correlated matrix element.
*/
virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
const SpinCorrelationTensor& c) const;
/**
* Order process in MadGraph conventions
*/
vector<size_t> gluonsFirst(vector<size_t> i);
/**
* Flush all cashes.
*/
virtual void flushCaches() {
MatchboxAmplitude::flushCaches();
}
/**
* Return true, if this amplitude needs to initialize an external
* code.
*/
virtual bool isExternal() const { return true; }
/**
* Initialize this amplitude
*/
virtual bool initializeExternal();
/**
* Return a generic process id for the given process
*/
virtual int externalId(const cPDVector&);
/**
* Write out BornAmplitudes.dat and VirtualAmplitudes.dat,
* return true if new files are written.
*/
virtual bool writeAmplitudesDat();
/**
* CHeck if all amplitudes are in the Libary.
*/
virtual bool checkAmplitudes();
string mgProcLibPath();
/**
* Return true, if this amplitude is capable of consistently filling
* the rho matrices for the spin correllations
*/
virtual bool canFillRhoMatrix() const { return true; }
/**
* Return the helicity combination of the physical process in the
* conventions used by the spin correlation algorithm.
*/
virtual vector<unsigned int> physicalHelicities(const vector<int>&) 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;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* The (tree-level) order in \f$g_S\f$ in which this matrix
* element is given.
*/
unsigned int theOrderInGs;
/**
* The (tree-level) order in \f$g_{EM}\f$ in which this matrix
* element is given.
*/
unsigned int theOrderInGem;
/**
* The path for the process libraries.
*/
string theProcessPath;
/**
* The path to generate amplitudes in.
*/
string theMGmodel;
bool keepinputtopmass;
/**
* Initialize the given process
*/
void initProcess(const cPDVector&);
/**
* Storage for Amplitudes
*/
static vector<string> BornAmplitudes,VirtAmplitudes;
/**
* Helper for color and crossing handling
*/
mutable vector<int> colourindex, crossing;
/**
* Static Variables to handle initialization.
*/
static bool ranMadGraphInitializeExternal;
static bool initializedMad;
//@}
protected:
/**
* Location of the installed executables
*/
- string prefix_;
+ string bindir_;
+
+ /**
+ * Location of the data files
+ */
+ string pkgdatadir_;
/**
* Location of MADGRAPH
*/
string madgraphPrefix_;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MadGraphAmplitude & operator=(const MadGraphAmplitude &);
};
}
#endif /* Herwig_MadGraphAmplitude_H */
diff --git a/MatrixElement/Matchbox/External/MadGraph/Makefile.am b/MatrixElement/Matchbox/External/MadGraph/Makefile.am
--- a/MatrixElement/Matchbox/External/MadGraph/Makefile.am
+++ b/MatrixElement/Matchbox/External/MadGraph/Makefile.am
@@ -1,21 +1,26 @@
pkglib_LTLIBRARIES = HwMatchboxMadGraph.la
HwMatchboxMadGraph_la_LDFLAGS = -module -version-info 1:0:0
HwMatchboxMadGraph_la_SOURCES = \
MadGraphAmplitude.h MadGraphAmplitude.cc
madgraphdir = ${pkgdatadir}/MadGraphInterface
INPUTFILES = \
InterfaceMadGraph.f.in \
param_card.dat.in \
param_card_heft.dat.in
dist_madgraph_DATA = $(INPUTFILES)
+HwMatchboxMadGraph_la_CPPFLAGS = $(AM_CPPFLAGS) \
+-DHERWIG_BINDIR="\"$(bindir)\"" \
+-DHERWIG_PKGDATADIR="\"$(pkgdatadir)\"" \
+-DMADGRAPH_PREFIX="\"$(MADGRAPHPREFIX)\""
+
install-exec-local:
$(install_sh_SCRIPT) $(builddir)/mg2Matchbox.py $(DESTDIR)$(bindir)/mg2Matchbox.py
uninstall-local:
rm -f $(DESTDIR)$(bindir)/mg2Matchbox.py
diff --git a/MatrixElement/Matchbox/External/MadGraph/mg2Matchbox.py.in b/MatrixElement/Matchbox/External/MadGraph/mg2Matchbox.py.in
--- a/MatrixElement/Matchbox/External/MadGraph/mg2Matchbox.py.in
+++ b/MatrixElement/Matchbox/External/MadGraph/mg2Matchbox.py.in
@@ -1,406 +1,405 @@
#! /usr/bin/env python
import os,sys,glob,errno,shutil,argparse,time
# helper to replace all sourceText in fileName with replaceText
def replacetext(fileName, sourceText, replaceText):
file = open(fileName, "r")
text = file.read()
file.close()
file = open(fileName, "w")
file.write(text.replace(sourceText, replaceText))
file.close()
# helper to build recursivly path
def mkdir_p(path):
try:
os.makedirs(path)
except OSError as exc: # Python >2.5
if exc.errno == errno.EEXIST and os.path.isdir(path):
pass
else: raise
# helper to find all files of with name in path
def find(name, path):
for root, dirs, files in os.walk(path):
if name in files:
return os.path.join(root, name)
# fill the proc.dat file from BornAmplitudes.dat and VirtAmplitudes.dat.
def fillprocs(model,oras,orew):
bornlist=[]
virtlist=[]
fileproc=open("proc.dat","w")
fileproc.write("import model "+model+"\n")
borns="BornAmplitudes.dat"
virts="VirtAmplitudes.dat"
first=True
procnr=0
virtlines=""
bornlines=""
minlegs=100
legs=0
for i in [borns, virts]:
file = open(i, "r")
for line in file:
if (len(line.split(" "))<minlegs):
minlegs=len(line.split(" "))
for it in line.split(" "):
if it.replace("-","").isdigit():
legs+=1
file.close()
#conversion for heft model to go from (2QCD+1QED)->1HIG for each FS HIGGS.
HIG=0
if (model=="heft"):
HIG=(int(oras)+int(orew)-legs+2)/2
if (int(oras)+int(orew)-legs+2)%2!=0:
print "Warning: No possible coupling power:(int(oras)+int(orew)-legs+2)%2!=0 "
exit()
return
file = open(borns, "r")
for line in file:
#this assumes extra QCD emmissions
addalphas=len(line.split(" "))-minlegs
linetmp=line.rstrip()
procnr+=1
bornlist+=[str(procnr)]
if first:
if HIG ==0 :
bornlines+="generate "+linetmp+" QCD="+str(int(oras)+addalphas)+" QED="+str(orew)+" @"+str(procnr)+"\n"
else:
bornlines+="generate "+linetmp+" HIG="+str(HIG)+" QCD="+str(int(oras)+addalphas-2*HIG)+" QED="+str(int(orew)-HIG)+" @"+str(procnr)+"\n"
first=False
else:
if HIG ==0 :
bornlines+="add process "+linetmp+" QCD="+str(int(oras)+addalphas)+" QED="+str(orew)+" @"+str(procnr)+"\n"
else:
bornlines+="add process "+linetmp+" HIG="+str(HIG)+" QCD="+str(int(oras)+addalphas-2*HIG)+" QED="+str(int(orew)-HIG)+" @"+str(procnr)+"\n"
file.close()
first=True
file = open(virts, "r")
for line in file:
addalphas=len(line.split(" "))-minlegs
linetmp=line.rstrip()+" QED="+str(int(orew))+" [ virt=QCD ]"
procnr+=1
virtlist+=[str(procnr)]
if first:
virtlines+="generate "+linetmp+" @"+str(procnr)+"\n"
first=False
else:
virtlines+="add process "+linetmp+" @"+str(procnr)+"\n"
file.close()
fileproc.write(bornlines)
if virtlines!="" and bornlines!="":
fileproc.write("output matchbox MG5 --postpone_model\n")
fileproc.write(virtlines)
fileproc.write("output matchbox MG5 -f\n")
fileproc.close()
return bornlist,virtlist
def build_matchbox_tmp(pwd,buildpath):
cwd=os.getcwd()
os.chdir(pwd)
mkdir_p(pwd+"/Herwig/MG_tmp/")
if not buildpath.startswith("/"):
buildpath=pwd+"/"+buildpath.lstrip("./")
resources=glob.glob(buildpath +"/MG5/SubProcesses/MadLoop5_resources/*")
resources+=glob.glob(buildpath +"/MG5/Cards/*")
resources+=glob.glob(buildpath +"/MG5/Cards/SubProcesses/*")
for i in resources:
if not os.path.isfile( pwd+"/Herwig/MG_tmp/"+os.path.basename(i)) \
and not os.path.islink( pwd+"/Herwig/MG_tmp/"+os.path.basename(i)):
os.symlink(i, pwd+"/Herwig/MG_tmp/"+os.path.basename(i))
os.chdir(cwd)
parser = argparse.ArgumentParser()
parser.add_argument('--buildpath', help='installpath')
parser.add_argument('--build', help='build', action="store_true")
parser.add_argument('--madgraph', help='madgraph_installpath')
parser.add_argument('--runpath', help='runpath')
parser.add_argument('--model', help='model')
parser.add_argument('--orderas', help='orderas')
parser.add_argument('--orderew', help='orderew')
+parser.add_argument('--datadir', help='datadir')
args = parser.parse_args()
pwd=os.getcwd()
param_card=""
mkdir_p(pwd+"/Herwig/MG_tmp/")
if args.model=="loop_sm" or args.model=="heft":
if args.model=="loop_sm":
param_card="param_card.dat"
else:
param_card="param_card_"+args.model+".dat"
- file = open("@prefix@/share/Herwig++/MadGraphInterface/"+param_card+".in", "r")
+ file = open("%s/MadGraphInterface/%s.in" % (args.datadir,param_card) , "r")
paramcard = file.read()
file.close()
file = open(args.runpath+"/"+param_card, "w")
params=open(args.runpath+"/MG-Parameter.dat", "r")
for line in params:
a=line.split()
paramcard=paramcard.replace(a[0],a[1])
params.close()
file.write(paramcard)
file.close()
else:
print "---------------------------------------------------------------"
print "---------------------------------------------------------------"
print "Warning: The model set for the MadGraph Interface "
print " needs a parameter setting by hand."
print " Please fill the param_card_"+args.model+".dat"
print " with your favourite assumptions."
print " And make sure Herwig uses the same parameters."
print "---------------------------------------------------------------"
print "---------------------------------------------------------------"
if os.path.isfile(args.buildpath +"/MG5/Cards/param_card.dat") and not os.path.isfile(args.runpath+"/"+"param_card_"+args.model+".dat"):
shutil.copyfile(args.buildpath +"/MG5/Cards/param_card.dat", args.runpath+"/"+"param_card_"+args.model+".dat")
time.sleep(1)
if not os.path.isdir(args.buildpath):
print "The MadGraph Install path was not existend. It has been created for you."
print "Just start Herwig++ read again.."
mkdir_p(args.buildpath)
exit()
os.chdir(args.buildpath)
if os.path.isfile("InterfaceMadGraph.so"):
build_matchbox_tmp(pwd,args.buildpath)
exit()
Bornlist,Virtlist=fillprocs(args.model,args.orderas,args.orderew)
if not args.madgraph and not os.path.isfile("InterfaceMadGraph.so"):
print "*** warning *** MadGraph build failed, check logfile for details"
exit()
os.system("python "+args.madgraph+"/mg5_aMC proc.dat")
routines=[["","BORN(momenta,hel)"],
["","SLOOPMATRIX(momenta,virt)"],
["","GET_JAMP(color,Jamp)"],
["","GET_LNJAMP(color,Jamp)"],
["","GET_NCOL(color)"],
["","GET_NCOLOR(i,j,color)"]]
for routine in routines:
for i in Bornlist + list(set(Virtlist) - set(Bornlist)):
if routine[1]=="Virt(amp)" or routine[1]=="SLOOPMATRIX(momenta,virt)" and i not in Virtlist:
continue
if routine[0]=="":
routine[0]+=" IF (proc .EQ. "+i+") THEN\n CALL "
routine[0]+= "MG5_"+i+"_"+routine[1]+"\n"
else:
routine[0]+=" ELSE IF (proc .EQ. "+i+") THEN\n"\
" CALL "
routine[0]+= "MG5_"+i+"_"+routine[1]+"\n"
if routine[0]!="":
routine[0]+=" ELSE\n"
routine[0]+=" WRITE(*,*) '##W02A WARNING No id found '\n"
routine[0]+=" ENDIF \n"
-
-
-
-shutil.copyfile("@prefix@/share/Herwig++/MadGraphInterface/InterfaceMadGraph.f.in", "InterfaceMadGraph.f")
+
+shutil.copyfile("%s/MadGraphInterface/InterfaceMadGraph.f.in" % args.datadir, "InterfaceMadGraph.f")
replacetext("InterfaceMadGraph.f","MG_CalculateBORNtxt",routines[0][0])
replacetext("InterfaceMadGraph.f","MG_CalculateVIRTtxt",routines[1][0])
replacetext("InterfaceMadGraph.f","MG_Jamptxt", routines[2][0])
replacetext("InterfaceMadGraph.f","MG_LNJamptxt", routines[3][0])
replacetext("InterfaceMadGraph.f","MG_NColtxt", routines[4][0])
replacetext("InterfaceMadGraph.f","MG_ColourMattxt",routines[5][0])
MG_vxxxxxtxt=""
if routines[1][0]!="":
MG_vxxxxxtxt=""" subroutine MG_vxxxxx(p, n,inc,VC)
$ bind(c, name='MG_vxxxxx')
IMPLICIT NONE
double precision p(0:3)
double precision n(0:3)
INTEGER inc
double precision VC(0:7)
double complex VCtmp(0:4)
double complex Ninplus(8)
double complex Noutminus(8)
double complex Pinplus(8)
double complex Poutplus(8)
double complex denom
double complex IMAG1
PARAMETER (IMAG1=(0D0,1D0))
CALL IXXXXX(n, 0.0d0, +1, +1, Ninplus); ! |n+>
CALL OXXXXX(p, 0.0d0, +1, +1, Poutplus); ! <p+|
CALL OXXXXX(n, 0.0d0, -1, +1, Noutminus);! <n-|
CALL IXXXXX(p, 0.0d0, +1, +1, Pinplus); ! |p+>
!<p+| gamma_mu |n+>
VCtmp(0)=Ninplus(5)*Poutplus(7)+Ninplus(6)*Poutplus(8)+
$ Ninplus(7)*Poutplus(5)+Ninplus(8)*Poutplus(6)
VCtmp(1)=(Ninplus(7)*Poutplus(6)+Ninplus(8)*Poutplus(5)-
$ Ninplus(5)*Poutplus(8)-Ninplus(6)*Poutplus(7))
VCtmp(2)=(-IMAG1*(Ninplus(5)*Poutplus(8)+Ninplus(8)
$ * Poutplus(5))+IMAG1*(Ninplus(6)*Poutplus(7)
$ + Ninplus(7)*Poutplus(6)));
VCtmp(3)=(Ninplus(6)*Poutplus(8)+Ninplus(7)*Poutplus(5)-
$ Ninplus(5)*Poutplus(7)-Ninplus(8)*Poutplus(6))
denom= (Pinplus(5)*Noutminus(5)+Pinplus(6)
$ *Noutminus(6)+Pinplus(7)*Noutminus(7)
$ +Pinplus(8)*Noutminus(8))*sqrt(2.d0)
if (inc.gt.1 )THEN
denom=denom*-1.0d0
if (abs(VCtmp(3)).ne.0.d0)THEN
denom=denom* (VCtmp(3))/abs(VCtmp(3))
endif
if (abs(VCtmp(1)).ne.0.d0)THEN
denom=denom* (VCtmp(1))/abs(VCtmp(1))
endif
endif
VCtmp(0)=-1.0d0*VCtmp(0)/denom
VCtmp(1)=-1.0d0*VCtmp(1)/denom
VCtmp(2)=-1.0d0*VCtmp(2)/denom
VCtmp(3)=-1.0d0*VCtmp(3)/denom
VC(0)= real(VCtmp(0))
VC(1)=aimag(VCtmp(0))
VC(2)= real(VCtmp(1))
VC(3)=aimag(VCtmp(1))
VC(4)= real(VCtmp(2))
VC(5)=aimag(VCtmp(2))
VC(6)= real(VCtmp(3))
VC(7)=aimag(VCtmp(3))
END"""
else:
MG_vxxxxxtxt=""" subroutine MG_vxxxxx(p, n,inc,VC)
$ bind(c, name='MG_vxxxxx')
IMPLICIT NONE
double precision p(0:3)
double precision n(0:3)
INTEGER inc
double precision VC(0:7)
double complex VCtmp(0:4)
double complex Ninplus(6)
double complex Noutminus(6)
double complex Pinplus(6)
double complex Poutplus(6)
double complex denom
double complex IMAG1
PARAMETER (IMAG1=(0D0,1D0))
CALL IXXXXX(n, 0.0d0, +1, +1, Ninplus); ! |n+>
CALL OXXXXX(p, 0.0d0, +1, +1, Poutplus); ! <p+|
CALL OXXXXX(n, 0.0d0, -1, +1, Noutminus);! <n-|
CALL IXXXXX(p, 0.0d0, +1, +1, Pinplus); ! |p+>
!<p+| gamma_mu |n+>
VCtmp(0)=Ninplus(3)*Poutplus(5)+Ninplus(4)*Poutplus(6)+
$ Ninplus(5)*Poutplus(3)+Ninplus(6)*Poutplus(4)
VCtmp(1)=(Ninplus(5)*Poutplus(4)+Ninplus(6)*Poutplus(3)-
$ Ninplus(3)*Poutplus(6)-Ninplus(4)*Poutplus(5))
VCtmp(2)=(-IMAG1*(Ninplus(3)*Poutplus(6)+Ninplus(6)
$ * Poutplus(3))+IMAG1*(Ninplus(4)*Poutplus(5)
$ + Ninplus(5)*Poutplus(4)));
VCtmp(3)=(Ninplus(4)*Poutplus(6)+Ninplus(5)*Poutplus(3)-
$ Ninplus(3)*Poutplus(5)-Ninplus(6)*Poutplus(4))
denom= (Pinplus(3)*Noutminus(3)+Pinplus(4)
$ *Noutminus(4)+Pinplus(5)*Noutminus(5)
$ +Pinplus(6)*Noutminus(6))*sqrt(2.d0)
if (inc.gt.1 )THEN
denom=denom*-1.0d0
if (abs(VCtmp(3)).ne.0.d0)THEN
denom=denom* (VCtmp(3))/abs(VCtmp(3))
endif
if (abs(VCtmp(1)).ne.0.d0)THEN
denom=denom* (VCtmp(1))/abs(VCtmp(1))
endif
endif
VCtmp(0)=-1.0d0*VCtmp(0)/denom
VCtmp(1)=-1.0d0*VCtmp(1)/denom
VCtmp(2)=-1.0d0*VCtmp(2)/denom
VCtmp(3)=-1.0d0*VCtmp(3)/denom
VC(0)= real(VCtmp(0))
VC(1)=aimag(VCtmp(0))
VC(2)= real(VCtmp(1))
VC(3)=aimag(VCtmp(1))
VC(4)= real(VCtmp(2))
VC(5)=aimag(VCtmp(2))
VC(6)= real(VCtmp(3))
VC(7)=aimag(VCtmp(3))
END"""
replacetext("InterfaceMadGraph.f","MG_vxxxxxtxt",MG_vxxxxxtxt)
make=" "
fortanfiles=glob.glob('*/*/*.f')+glob.glob('*/*/*/*.f')
for i in fortanfiles:
if "check_sa" not in i:
if not os.path.islink(i):
make += " "+i+"\\\n "
incfiles=glob.glob('*/*/*.inc')+glob.glob('*/*/*/*.inc')
coefdir=""
for i in incfiles:
if "nexternal.inc" in i:
coefdir+=" -I"+i.replace("nexternal.inc"," ")
file=open("makefile","w")
file.write("include MG5/Source/make_opts ")
if Virtlist!=[]:
file.write("\nLIBDIR = MG5/lib\nLINKLIBS = -L$(LIBDIR) -lcts -liregi -L$(LIBDIR)/golem95_lib -lgolem")
file.write("\nLIBS = $(LIBDIR)libcts.$(libext) $(LIBDIR)libgolem.$(libext) $(LIBDIR)libiregi.$(libext)")
file.write("\nPROCESS= InterfaceMadGraph.f "+make+"\n\nall: \n\t @FC@ @FFLAGS@ -w -fbounds-check -ffixed-line-length-132 -fPIC -fno-f2c -shared -s -o InterfaceMadGraph.so -IMG5/SubProcesses/" )
if Virtlist!=[]:
file.write(" -IMG5/lib/golem95_include ")
if coefdir != "":
file.write(coefdir)
file.write(" $(PROCESS) $(LINKLIBS) ")
file.close()
os.chdir(pwd)
os.chdir(args.buildpath)
replacetext("MG5/Source/MODEL/lha_read.f", "ident_card.dat","Herwig/MG_tmp/ident_card.dat")
replacetext("MG5/Source/MODEL/lha_read.f", "param.log","Herwig/MG_tmp/param.log")
if Virtlist!=[]:
replacetext("MG5/SubProcesses/MadLoopCommons.f", "PREFIX='./'","PREFIX='./Herwig/MG_tmp/'")
os.system("make")
build_matchbox_tmp(pwd,args.buildpath)
diff --git a/MatrixElement/Matchbox/External/NJet/Makefile.am b/MatrixElement/Matchbox/External/NJet/Makefile.am
--- a/MatrixElement/Matchbox/External/NJet/Makefile.am
+++ b/MatrixElement/Matchbox/External/NJet/Makefile.am
@@ -1,9 +1,8 @@
pkglib_LTLIBRARIES = HwMatchboxNJet.la
HwMatchboxNJet_la_LDFLAGS = -module -version-info 11:1:0
-HwMatchboxNJet_la_CPPFLAGS = $(AM_CPPFLAGS)
-HwMatchboxNJet_la_CPPFLAGS += -I$(NJETINCLUDEPATH)
+HwMatchboxNJet_la_CPPFLAGS = $(AM_CPPFLAGS) -I$(NJETINCLUDEPATH) \
+-DNJET_PREFIX="\"$(NJETPREFIX)\""
HwMatchboxNJet_la_SOURCES = \
NJetsAmplitude.h NJetsAmplitude.cc
-
diff --git a/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc.in b/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc
rename from MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc.in
rename to MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc
--- a/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc.in
+++ b/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc
@@ -1,280 +1,290 @@
// -*- C++ -*-
//
// NJetsAmplitude.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the NJetsAmplitude class.
//
#include "NJetsAmplitude.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Utilities/DynamicLoader.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
#include "njet.h"
#include <cstdlib>
+#ifndef NJET_PREFIX
+#error Makefile.am needs to define NJET_PREFIX
+#endif
+
using namespace Herwig;
-NJetsAmplitude::NJetsAmplitude() {}
+NJetsAmplitude::NJetsAmplitude() : NJetsPrefix_(NJET_PREFIX) {}
NJetsAmplitude::~NJetsAmplitude() {}
IBPtr NJetsAmplitude::clone() const {
return new_ptr(*this);
}
IBPtr NJetsAmplitude::fullclone() const {
return new_ptr(*this);
}
void NJetsAmplitude::signOLP(const string& order, const string& contract) {
- string cmd = "@NJETPREFIX@/bin/njet.py -o " + contract + " " + order;
+ string cmd = NJetsPrefix_+"/bin/njet.py -o " + contract + " " + order;
std::system(cmd.c_str());
}
void NJetsAmplitude::startOLP(const string& contract, int& status) {
NJet::LH_OLP::OLP_Start(contract.c_str(), &status);
if ( status != 1 )
return;
status = 0;
static double zero = 0.0;
double param = 0.0;
param = SM().alphaEMMZ();
NJet::LH_OLP::OLP_SetParameter("alpha",&param,&zero,&status);
if ( status != 1 )
return;
param = getParticleData(ParticleID::Z0)->mass()/GeV;
NJet::LH_OLP::OLP_SetParameter("mass(23)",&param,&zero,&status);
if ( status != 1 )
return;
param = getParticleData(ParticleID::Wplus)->mass()/GeV;
NJet::LH_OLP::OLP_SetParameter("mass(24)",&param,&zero,&status);
if ( status != 1 )
return;
param = getParticleData(ParticleID::Z0)->width()/GeV;
NJet::LH_OLP::OLP_SetParameter("width(23)",&param,&zero,&status);
if ( status != 1 )
return;
param = getParticleData(ParticleID::Wplus)->width()/GeV;
NJet::LH_OLP::OLP_SetParameter("width(24)",&param,&zero,&status);
if ( status != 1 )
return;
param = SM().sin2ThetaW();
NJet::LH_OLP::OLP_SetParameter("sw2",&param,&zero,&status);
didStartOLP() = true;
}
bool NJetsAmplitude::startOLP(const map<pair<Process,int>,int>& procs) {
if ( !DynamicLoader::load("libnjet2.so") )
throw Exception() << "failed to load libnjet2.so\n"
<< DynamicLoader::lastErrorMessage
<< Exception::abortnow;
// TODO throw exception on massive leptons in procs
string orderFileName = factory()->buildStorage() + name() + ".OLPOrder.lh";
ofstream orderFile(orderFileName.c_str());
olpOrderFileHeader(orderFile);
orderFile << "NJetReturnAccuracy yes\n"
<< "NJetRenormalize yes\n"
<< "NJetNf " << factory()->nLight() << "\n";
olpOrderFileProcesses(orderFile,procs);
orderFile << flush;
orderFile.close();
string contractFileName = factory()->buildStorage() + name() + ".OLPContract.lh";
signOLP(orderFileName, contractFileName);
int status = -1;
startOLP(contractFileName,status);
if ( status != 1 )
return false;
return true;
}
LorentzVector<Complex> NJetsAmplitude::plusPolarization(const Lorentz5Momentum& p,
const Lorentz5Momentum& n,
int inc) const {
double pvec[4] = {p.t()/GeV,p.x()/GeV,p.y()/GeV,p.z()/GeV};
double nvec[4] = {n.t()/GeV,n.x()/GeV,n.y()/GeV,n.z()/GeV};
double out[8] ={ };
NJet::LH_OLP::OLP_Polvec(pvec,nvec,out);
LorentzVector<Complex> res;
Complex a(out[0],out[1]);
res.setT(a);
Complex b(out[2],out[3]);
res.setX(b);
Complex c(out[4],out[5]);
res.setY(c);
Complex d(out[6],out[7]);
res.setZ(d);
if (inc<2)
return res.conjugate();
else
return res;
}
void NJetsAmplitude::evalSubProcess() const {
useMe();
double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double as;
if (!hasRunningAlphaS()) as = SM().alphaS();
else if (hasRunningAlphaS()) as = lastAlphaS();
double scale = sqrt(mu2()/GeV2);
double out[7]={};
int id =
olpId()[ProcessType::oneLoopInterference] ?
olpId()[ProcessType::oneLoopInterference] :
olpId()[ProcessType::treeME2];
NJet::LH_OLP::OLP_EvalSubProcess(id, olpMomenta(), scale, &as, out);
if ( olpId()[ProcessType::oneLoopInterference] ) {
if(calculateTreeME2())lastTreeME2(out[3]*units);
lastOneLoopInterference(out[2]*units);
lastOneLoopPoles(pair<double,double>(out[0]*units,out[1]*units));
} else if ( olpId()[ProcessType::treeME2] ) {
lastTreeME2(out[0]*units);
} else assert(false);
}
void NJetsAmplitude::evalColourCorrelator(pair<int,int>) const {
double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double as;
if (!hasRunningAlphaS()) as = SM().alphaS();
else if (hasRunningAlphaS()) as = lastAlphaS();
double scale = sqrt(mu2()/GeV2);
int n = lastXComb().meMomenta().size();
colourCorrelatorResults.resize(n*(n-1)/2);
NJet::LH_OLP::OLP_EvalSubProcess(olpId()[ProcessType::colourCorrelatedME2],
olpMomenta(), scale, &as, &colourCorrelatorResults[0]);
for ( int i = 0; i < n; ++i )
for ( int j = i+1; j < n; ++j ) {
lastColourCorrelator(make_pair(i,j),colourCorrelatorResults[i+j*(j-1)/2]*units);
}
}
void NJetsAmplitude::evalSpinColourCorrelator(pair<int,int>) const {
double units = pow(lastSHat()/GeV2,mePartonData().size()-4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double as;
if (!hasRunningAlphaS()) as = SM().alphaS();
else if (hasRunningAlphaS()) as = lastAlphaS();
double scale = sqrt(mu2()/GeV2);
int n = lastXComb().meMomenta().size();
spinColourCorrelatorResults.resize(2*n*n);
NJet::LH_OLP::OLP_EvalSubProcess(olpId()[ProcessType::spinColourCorrelatedME2],
olpMomenta(), scale, &as, &spinColourCorrelatorResults[0]);
for ( int i = 0; i < n; ++i )
for ( int j = 0; j < n; ++j ) {
if ( i == j || mePartonData()[i]->id() != 21 )
continue;
Complex scc(spinColourCorrelatorResults[2*i+2*n*j]*units,
spinColourCorrelatorResults[2*i+2*n*j+1]*units);
lastColourSpinCorrelator(make_pair(i,j),scc);
}
}
void NJetsAmplitude::doinit() {
if ( !DynamicLoader::load("libnjet2.so") )
throw Exception() << "failed to load libnjet2.so\n"
<< DynamicLoader::lastErrorMessage
<< Exception::abortnow;
MatchboxOLPME::doinit();
}
void NJetsAmplitude::doinitrun() {
if ( !DynamicLoader::load("libnjet2.so") )
throw Exception() << "failed to load libnjet2.so\n"
<< DynamicLoader::lastErrorMessage
<< Exception::abortnow;
MatchboxOLPME::doinitrun();
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void NJetsAmplitude::persistentOutput(PersistentOStream & os) const {
os << colourCorrelatorResults << spinColourCorrelatorResults;
}
void NJetsAmplitude::persistentInput(PersistentIStream & is, int) {
is >> colourCorrelatorResults >> spinColourCorrelatorResults;
}
// *** 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<NJetsAmplitude,MatchboxOLPME>
describeHerwigNJetsAmplitude("Herwig::NJetsAmplitude", "HwMatchboxNJet.so");
void NJetsAmplitude::Init() {
static ClassDocumentation<NJetsAmplitude> documentation
("NJetsAmplitude implements an interface to NJets.", "Matrix elements have been calculated using NJet.");
+
+ static Parameter<NJetsAmplitude,string> interfaceNJetsPrefix
+ ("NJetsPrefix",
+ "The prefix for the location of NJets",
+ &NJetAmplitude::NJetsPrefix_, string(NJET_PREFIX),
+ false, false);
}
diff --git a/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.h b/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.h
--- a/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.h
+++ b/MatrixElement/Matchbox/External/NJet/NJetsAmplitude.h
@@ -1,188 +1,193 @@
// -*- C++ -*-
//
// NJetsAmplitude.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_NJetsAmplitude_H
#define Herwig_NJetsAmplitude_H
//
// This is the declaration of the NJetsAmplitude class.
//
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxOLPME.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief NJetsAmplitude implements an interface to NJets
*/
class NJetsAmplitude: public MatchboxOLPME {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
NJetsAmplitude();
/**
* The destructor.
*/
virtual ~NJetsAmplitude();
//@}
public:
/**
* Return true, if this amplitude already includes averaging over
* incoming parton's quantum numbers.
*/
virtual bool hasInitialAverage() const { return false; }
/**
* Return true, if this amplitude already includes symmetry factors
* for identical outgoing particles.
*/
virtual bool hasFinalStateSymmetry() const { return false; }
/**
* Start the one loop provider, if appropriate, giving order and
* contract files
*/
virtual void signOLP(const string&, const string&);
/**
* Start the one loop provider, if appropriate
*/
virtual void startOLP(const string&, int& status);
/**
* 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);
/**
* Call OLP_EvalSubProcess and fill in the results
*/
virtual void evalSubProcess() const;
/**
* Fill in results for the given colour correlator
*/
virtual void evalColourCorrelator(pair<int,int> ij) const;
/**
* Return a positive helicity polarization vector for a gluon of
* momentum p (with reference vector n) to be used when evaluating
* spin correlations.
*/
virtual LorentzVector<Complex> plusPolarization(const Lorentz5Momentum& p,
const Lorentz5Momentum& n,
int id = -1) const;
/**
* Fill in results for the given colour/spin correlator
*/
virtual void evalSpinColourCorrelator(pair<int,int> ij) 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 Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
NJetsAmplitude & operator=(const NJetsAmplitude &);
/**
* Store colour correlator results
*/
mutable vector<double> colourCorrelatorResults;
/**
* Store spin colour correlator results
*/
mutable vector<double> spinColourCorrelatorResults;
+ /**
+ * Location of NJETs
+ */
+ string NJetPrefix_;
+
};
}
#endif /* Herwig_NJetsAmplitude_H */
diff --git a/MatrixElement/Matchbox/External/OpenLoops/Makefile.am b/MatrixElement/Matchbox/External/OpenLoops/Makefile.am
--- a/MatrixElement/Matchbox/External/OpenLoops/Makefile.am
+++ b/MatrixElement/Matchbox/External/OpenLoops/Makefile.am
@@ -1,6 +1,7 @@
pkglib_LTLIBRARIES = HwMatchboxOpenLoops.la
-HwMatchboxOpenLoops_la_LDFLAGS = -module -version-info 11:1:0
-
HwMatchboxOpenLoops_la_SOURCES = \
OpenLoopsAmplitude.h OpenLoopsAmplitude.cc
+HwMatchboxOpenLoops_la_LDFLAGS = -module -version-info 11:1:0 -R$(OPENLOOPSLIBS) -L$(OPENLOOPSLIBS)
+HwMatchboxOpenLoops_la_LIBADD = -lopenloops
+
diff --git a/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc.in b/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc
rename from MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc.in
rename to MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc
--- a/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc.in
+++ b/MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc
@@ -1,430 +1,419 @@
// -*- C++ -*-
//
// OpenLoopsAmplitude.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the OpenLoopsAmplitude class.
//
#include "OpenLoopsAmplitude.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
-#include "ThePEG/Utilities/DynamicLoader.h"
-
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
#include <fstream>
#include <sstream>
#include <string>
#include <cstdlib>
using namespace Herwig;
OpenLoopsAmplitude::OpenLoopsAmplitude():use_cms(true),psp_tolerance(12) {
}
OpenLoopsAmplitude::~OpenLoopsAmplitude() {
}
IBPtr OpenLoopsAmplitude::clone() const {
return new_ptr(*this);
}
IBPtr OpenLoopsAmplitude::fullclone() const {
return new_ptr(*this);
}
extern "C" void OLP_Start(const char*, int* i);
extern "C" void OLP_SetParameter(const char* ,double* ,double*,int*);
extern "C" void ol_setparameter_string(const char*, const char*);
extern "C" void OLP_PrintParameter(const char*);
extern "C" void OLP_EvalSubProcess(int*, double*, double*, double*, double*);
extern "C" void OLP_EvalSubProcess2(int*, double*, double*, double*, double*);
// id ps-point emitter polvec res
extern "C" void ol_evaluate_sc(int, double*, int, double*, double*);
extern "C" void OLP_Polvec(double*,double*,double*);
void OpenLoopsAmplitude::doinitrun() {
MatchboxOLPME::doinitrun();
}
void OpenLoopsAmplitude::startOLP(const string& contract, int& status) {
string tempcontract=contract;
- bool success = DynamicLoader::load("@OPENLOOPSLIBS@/libopenloops.so")
- ||DynamicLoader::load("@OPENLOOPSLIBS@/libopenloops.dylib") ;
-
- if ( !success ) {
- throw Exception() << "Failed to load libopenloops.so/dylib\n"
- << DynamicLoader::lastErrorMessage
- << Exception::abortnow;
- }
-
string stabilityPrefix = factory()->runStorage() + "OpenLoops.StabilityLog";
assert(stabilityPrefix.size() < 256);
ol_setparameter_string("stability_logdir",stabilityPrefix.c_str());
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])->mass()/GeV;
else
mass = it->second/GeV;
double width=getParticleData(part[i])->width()/GeV;
std::stringstream ss;
ss << part[i];
string str = ss.str();
stri="mass("+str+")";
OLP_SetParameter(stri.c_str(),&mass,&null,&a);
stri="width("+str+")";
OLP_SetParameter(stri.c_str(),&width,&null,&a);
}
stri="alphas";
one=SM().alphaS();
OLP_SetParameter( stri.c_str(),&one ,&null,&a);
stri="alpha";
one=SM().alphaEMMZ();
OLP_SetParameter(stri.c_str(),&one ,&null,&a);
OLP_Start(tempcontract.c_str(), &status);
didStartOLP() = true;
}
void OpenLoopsAmplitude::fillOrderFile(const map<pair<Process, int>, int>& procs) {
string orderFileName = factory()->buildStorage() + name() + ".OLPContract.lh";
ofstream orderFile(orderFileName.c_str());
size_t asPower = 100;
size_t minlegs = 100;
size_t maxlegs = 0;
for ( map<pair<Process, int>, int>::const_iterator t = procs.begin() ; t != procs.end() ; ++t ) {
asPower = min(asPower, static_cast<size_t>(t->first.first.orderInAlphaS));
minlegs = min(minlegs, static_cast<size_t>(t->first.first.legs.size()));
maxlegs = max(maxlegs, static_cast<size_t>(t->first.first.legs.size()));
}
orderFile << "# OLP order file created by Herwig++/Matchbox for OpenLoops\n\n";
orderFile << "CorrectionType QCD\n";
orderFile << "IRregularization " << (isDR() ? "DRED" : "CDR") << "\n";
orderFile << "extra answerfile " << (factory()->buildStorage() + name() + ".OLPAnswer.lh") << "\n";
orderFile << "extra psp_tolerance "<<psp_tolerance<<"\n";
orderFile << "extra use_cms "<<(use_cms?"1":"0")<< "\n";
orderFile << "\n";
if (extraOpenLoopsPath!="")
orderFile << "Extra OpenLoopsPath " << extraOpenLoopsPath << "\n";
for ( map<pair<Process, int>, int>::const_iterator p = procs.begin() ; p != procs.end() ; ++p ) {
std::stringstream Processstr;
std::stringstream Typestr;
Processstr << (*p).first.first.legs[0]->id() << " " << (*p).first.first.legs[1]->id() << " -> ";
for ( PDVector::const_iterator o = (*p).first.first.legs.begin() + 2 ; o != (*p).first.first.legs.end() ; ++o )
Processstr << (**o).id() << " ";
if ( (*p).first.second == ProcessType::treeME2 ) {
Typestr << "Tree";
} else if ( (*p).first.second == ProcessType::colourCorrelatedME2 ) {
Typestr << "ccTree";
} else if ( (*p).first.second == ProcessType::spinColourCorrelatedME2 ) {
Typestr << "sctree_polvect";
} else if ( (*p).first.second == ProcessType::oneLoopInterference ) {
Typestr << "Loop";
}
OpenLoopsProcInfo pro = OpenLoopsProcInfo((*p).second, -1, Processstr.str(), Typestr.str());
pro.setOAs(p->first.first.orderInAlphaS);
processmap[(*p).second] = pro;
}
vector < string > types;
types.push_back("Tree");
types.push_back("ccTree");
types.push_back("sctree_polvect");
types.push_back("Loop");
for ( size_t i = asPower ; i != asPower + maxlegs - minlegs + 1 ; i++ ) {
orderFile << "\n\nCouplingPower QCD " << i;
orderFile << "\n\n#AlphasPower " << i;
for ( vector<string>::iterator it = types.begin() ; it != types.end() ; it++ ) {
for ( map<int, OpenLoopsProcInfo>::iterator p = processmap.begin() ; p != processmap.end() ; ++p )
if ( (*p).second.Tstr() == *it && i == (unsigned int) (*p).second.orderAs() ) {
orderFile << "\nAmplitudeType " << *it << "\n";
break;
}
for ( map<int, OpenLoopsProcInfo>::iterator p = processmap.begin() ; p != processmap.end() ; ++p )
if ( (*p).second.Tstr() == *it && i == (unsigned int) (*p).second.orderAs() ) {
orderFile << (*p).second.Pstr() << "\n";
}
}
}
orderFile << flush;
}
bool OpenLoopsAmplitude::checkOLPContract() {
string contractFileName = factory()->buildStorage() + name() + ".OLPAnswer.lh";
ifstream infile(contractFileName.c_str());
string line;
vector < string > contractfile;
while (std::getline(infile, line)) {
contractfile.push_back(line);
}
for ( map<int, OpenLoopsProcInfo>::iterator p = processmap.begin() ; p != processmap.end() ; p++ ) {
bool righttype = false;
for ( vector<string>::iterator linex = contractfile.begin() ; linex != contractfile.end() ; ++linex ) {
if ( (*linex).find("AmplitudeType ")!= std::string::npos ) {
if ( (*linex).find(" " + (*p).second.Tstr() + " ")!= std::string::npos ) {
righttype = true;
} else {
righttype = false;
}
}
if ( righttype ) {
if ( (*linex).find((*p).second.Pstr()) != std::string::npos ){
if( (*p).second.Pstr().length() == (*linex).find("|") ) {
string sub = (*linex).substr((*linex).find("|") + 1, (*linex).find("#") - (*linex).find("|") - 1); // | 1 23 # buggy??
int subint;
int subint2;
istringstream(sub) >> subint >> subint2;
assert(subint==1);
(*p).second.setGID(subint2);
}
}
}
}
}
string ids = factory()->buildStorage() + "OpenLoops.ids.dat";
ofstream IDS(ids.c_str());
for ( map<int, OpenLoopsProcInfo>::iterator p = processmap.begin() ; p != processmap.end() ; p++ ) {
idpair.insert ( std::pair<int,int>((*p).second.HID(),(*p).second.GID()) );
IDS << (*p).second.HID() << " " << (*p).second.GID() << "\n";
if ( (*p).second.GID() == -1 ) return 0;
}
IDS << flush;
return 1;
}
void OpenLoopsAmplitude::getids() const{
string line = factory()->buildStorage() + "OpenLoops.ids.dat";
ifstream infile(line.c_str());
int hid;
int gid;
while (std::getline(infile, line)) {
istringstream(line) >> hid>>gid;
idpair.insert ( std::pair<int,int>(hid,gid) );
}
}
bool OpenLoopsAmplitude::startOLP(const map<pair<Process, int>, int>& procs) {
string contractFileName = factory()->buildStorage() + name() + ".OLPAnswer.lh";
string orderFileName = factory()->buildStorage() + name() + ".OLPContract.lh";
fillOrderFile(procs);
int status = -1;
startOLP(orderFileName, status);
if ( !checkOLPContract() ) {
return false;
}
if ( status != 1 ) return false;
return true;
}
void OpenLoopsAmplitude::evalSubProcess() const {
useMe();
double units = pow(lastSHat() / GeV2, mePartonData().size() - 4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double acc ;
double scale = sqrt(mu2() / GeV2);
if (hasRunningAlphaS()) {
int a=0;double null=0.0;double one=1.0;
string stri="alphas";
one=lastAlphaS();
OLP_SetParameter( stri.c_str(),&one ,&null,&a);
}
double out[7]={};
int id = olpId()[ProcessType::oneLoopInterference] ? olpId()[ProcessType::oneLoopInterference] : olpId()[ProcessType::treeME2];
if ( idpair.size() == 0 ) {
getids();
if ( Debug::level > 1 ) {
string parfile=factory()->runStorage() + name() + ".Parameters.dat";
OLP_PrintParameter(parfile.c_str());
}
}
OLP_EvalSubProcess2(&((*(idpair.find(id))).second), olpMomenta(), &scale, out,&acc );
if ( olpId()[ProcessType::oneLoopInterference] ) {
if(calculateTreeME2())lastTreeME2(out[3] * units);
lastOneLoopInterference((out[2])* units);
lastOneLoopPoles(pair<double, double>(out[0] * units, out[1] * units));
} else if ( olpId()[ProcessType::treeME2] ) {
lastTreeME2(out[0] * units);
}
}
void OpenLoopsAmplitude::evalColourCorrelator(pair<int, int> ) const {
double units = pow(lastSHat() / GeV2, mePartonData().size() - 4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
double acc ;
double scale = sqrt(mu2() / GeV2);
if (hasRunningAlphaS()) {
int a=0;double null=0.0;double one=1.0;
string stri="alphas";
one=lastAlphaS();
OLP_SetParameter( stri.c_str(),&one ,&null,&a);
}
int n = lastXComb().meMomenta().size();
colourCorrelatorResults.resize(n * (n - 1) / 2);
if ( idpair.size() == 0 ) {
getids();
if ( Debug::level > 1 ) {
string parfile=factory()->runStorage() + name() + ".Parameters.dat";
OLP_PrintParameter(parfile.c_str());
}
}
int id = olpId()[ProcessType::colourCorrelatedME2];
OLP_EvalSubProcess2(&((*(idpair.find(id))).second), olpMomenta(), &scale, &colourCorrelatorResults[0],&acc );
for ( int i = 0 ; i < n ; ++i ){
for ( int j = i + 1 ; j < n ; ++j ) {
lastColourCorrelator(make_pair(i, j), colourCorrelatorResults[i+j*(j-1)/2] * units);
}
}
}
void OpenLoopsAmplitude::evalSpinColourCorrelator(pair<int , int > ) const {
assert(false);
}
double OpenLoopsAmplitude::spinColourCorrelatedME2(pair<int,int> ij,
const SpinCorrelationTensor& c) const{
double units = pow(lastSHat() / GeV2, mePartonData().size() - 4.);
fillOLPMomenta(lastXComb().meMomenta(),mePartonData(),reshuffleMasses());
if (hasRunningAlphaS()) {
int a=0;double null=0.0;double one=1.0;
string stri="alphas";
one=lastAlphaS();
OLP_SetParameter( stri.c_str(),&one ,&null,&a);
}
int emitter=ij.first+1;
int n = lastXComb().meMomenta().size();
if ( idpair.size() == 0 ) {
getids();
if ( Debug::level > 1 ) {
string parfile=factory()->runStorage() + name() + ".Parameters.dat";
OLP_PrintParameter(parfile.c_str());
}
}
int id = (*(idpair.find(olpId()[ProcessType::spinColourCorrelatedME2]))).second;
//double * outx =new double[n];
spinColourCorrelatorResults.resize(n);
double polvec[4];
polvec[0]=c.momentum().e()/GeV;
polvec[1]=c.momentum().x()/GeV;
polvec[2]=c.momentum().y()/GeV;
polvec[3]=c.momentum().z()/GeV;
double avg= colourCorrelatedME2(ij)*(-c.diagonal());
ol_evaluate_sc(id, olpMomenta(),emitter,polvec,&spinColourCorrelatorResults[0]);
double corr =-1.*units * spinColourCorrelatorResults[ij.second]/c.scale()*c.momentum().dot(c.momentum());
double Nc = generator()->standardModel()->Nc();
double cfac = 1.;
if ( mePartonData()[ij.first]->iColour() == PDT::Colour8 ) {
cfac = Nc;
} else if ( mePartonData()[ij.first]->iColour() == PDT::Colour3 ||
mePartonData()[ij.first]->iColour() == PDT::Colour3bar ) {
cfac = (sqr(Nc)-1.)/(2.*Nc);
} else assert(false);
return
avg + corr/cfac;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void OpenLoopsAmplitude::persistentOutput(PersistentOStream & os) const {
os << idpair ;
}
void OpenLoopsAmplitude::persistentInput(PersistentIStream & is, int) {
is >> idpair ;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<OpenLoopsAmplitude, MatchboxOLPME> describeHerwigOpenLoopsAmplitude("Herwig::OpenLoopsAmplitude", "HwMatchboxOpenLoops.so");
void OpenLoopsAmplitude::Init() {
static ClassDocumentation<OpenLoopsAmplitude> documentation("OpenLoopsAmplitude implements an interface to OpenLoops.","Matrix elements have been calculated using OpenLoops.");
static Switch<OpenLoopsAmplitude,bool> interfaceUseComplMass
("ComplexMassScheme",
"Switch on or off if Compex Masses.",
&OpenLoopsAmplitude::use_cms, true, false, false);
static SwitchOption interfaceUseComplMassOn
(interfaceUseComplMass,
"True",
"True for Compex Masses.",
true);
static SwitchOption interfaceUseComplMassOff
(interfaceUseComplMass,
"False",
"False for no Compex Masses.",
false);
static Parameter<OpenLoopsAmplitude,int> interfacepsp_tolerance
("PSP_tolerance",
"(Debug)Phase Space Tolerance. Better use e.g.: set OpenLoops:Massless 13",
&OpenLoopsAmplitude::psp_tolerance, 12, 0, 0,
false, false, Interface::lowerlim);
}
diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -1,287 +1,283 @@
dnl Process this file with autoconf to produce a configure script.
AC_PREREQ([2.63])
AC_INIT([Herwig++],[devel],[herwig@projects.hepforge.org],[Herwig++])
AC_CONFIG_SRCDIR([Utilities/HerwigStrategy.cc])
AC_CONFIG_AUX_DIR([Config])
AC_CONFIG_MACRO_DIR([m4])
AC_CONFIG_HEADERS([Config/config.h])
dnl AC_PRESERVE_HELP_ORDER
AC_CANONICAL_HOST
case "${host}" in
*-darwin[[0156]].*)
AC_MSG_ERROR([Herwig++ requires OS X 10.3 or later])
;;
*-darwin7.*)
if test "x$MACOSX_DEPLOYMENT_TARGET" != "x10.3"; then
AC_MSG_ERROR(
[Please add 'MACOSX_DEPLOYMENT_TARGET=10.3' to the configure line.])
fi
;;
esac
dnl === disable debug symbols by default =====
if test "x$CXXFLAGS" = "x"; then
CXXFLAGS=-O3
fi
if test "x$CFLAGS" = "x"; then
CFLAGS=-O3
fi
dnl Looptools manual requires optimization off
dnl if test "x$FCFLAGS" = "x"; then
dnl FCFLAGS=-O0
dnl fi
dnl ==========================================
AC_LANG([C++])
AM_INIT_AUTOMAKE([1.11 subdir-objects gnu dist-bzip2 -Wall])
m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
m4_ifdef([AM_PROG_AR], [AM_PROG_AR])
dnl Checks for C++ compiler. Handle C++11 flags.
AC_PROG_CXX([g++])
dnl this block can go once we decide to have C++11 always on
AC_MSG_CHECKING([whether to include C++11 flag for testing])
AC_ARG_ENABLE(stdcxx11,
AC_HELP_STRING([--enable-stdcxx11],
[turn on C++11 flag (only for testing, do not use in production!)]),
[],
[enable_stdcxx11=no]
)
if test "x$enable_stdcxx11" = "xyes"; then
AC_MSG_RESULT([yes])
dnl remove the wrapper if block and "optional" once we decide to have C++11 always on
AX_CXX_COMPILE_STDCXX_11([noext],[optional])
if test "x$HAVE_CXX11" != "x1"; then
AC_MSG_ERROR([compiler does not recognize requested c++11 option])
fi
else
AC_MSG_RESULT([no])
fi
dnl check for POSIX
AC_CHECK_HEADER([unistd.h],[],
[AC_MSG_ERROR([Herwig++ needs "unistd.h". Non-POSIX systems are not supported.])])
dnl Checks for programs.
AC_PROG_INSTALL
AC_PROG_MAKE_SET
AC_PROG_LN_S
dnl modified search order
AC_PROG_FC([gfortran g95 g77])
dnl xlf95 f95 fort ifort ifc efc pgf95 lf95 ftn xlf90 f90 pgf90 pghpf epcf90 xlf f77 frt pgf77 cf77 fort77 fl32 af77])
AC_LANG_PUSH([Fortran])
AC_MSG_CHECKING([if the Fortran compiler ($FC) works])
AC_COMPILE_IFELSE(
AC_LANG_PROGRAM([],[ print *[,]"Hello"]),
[AC_MSG_RESULT([yes])],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([A Fortran compiler is required to build Herwig++.])
]
)
AC_LANG_POP([Fortran])
LT_PREREQ([2.2.6])
LT_INIT([disable-static dlopen pic-only])
dnl ####################################
dnl ####################################
dnl for Doc/fixinterfaces.pl
AC_PATH_PROG(PERL, perl)
dnl for Models/Feynrules
AM_PATH_PYTHON([2.6],, [:])
AM_CONDITIONAL([HAVE_PYTHON], [test "x$PYTHON" != "x:"])
HERWIG_CHECK_GSL
HERWIG_CHECK_THEPEG
BOOST_REQUIRE([1.41])
BOOST_FIND_HEADER([boost/array.hpp])
BOOST_FIND_HEADER([boost/numeric/ublas/io.hpp])
BOOST_FIND_HEADER([boost/numeric/ublas/matrix.hpp])
BOOST_FIND_HEADER([boost/numeric/ublas/matrix_proxy.hpp])
BOOST_FIND_HEADER([boost/numeric/ublas/matrix_sparse.hpp])
BOOST_FIND_HEADER([boost/numeric/ublas/symmetric.hpp])
BOOST_FIND_HEADER([boost/numeric/ublas/vector.hpp])
BOOST_FIND_HEADER([boost/operators.hpp])
BOOST_FIND_HEADER([boost/progress.hpp])
BOOST_FIND_HEADER([boost/scoped_array.hpp])
BOOST_FIND_HEADER([boost/scoped_ptr.hpp])
BOOST_FIND_HEADER([boost/utility.hpp])
BOOST_FILESYSTEM
HERWIG_CHECK_VBFNLO
HERWIG_CHECK_NLOJET
HERWIG_CHECK_NJET
HERWIG_CHECK_GOSAM
HERWIG_CHECK_GOSAM_CONTRIB
HERWIG_CHECK_OPENLOOPS
HERWIG_CHECK_HEJ
HERWIG_CHECK_MADGRAPH
HERWIG_COMPILERFLAGS
HERWIG_LOOPTOOLS
HERWIG_PDF_PATH
FASTJET_CHECK_FASTJET
HERWIG_VERSIONSTRING
HERWIG_CHECK_ABS_BUG
HERWIG_ENABLE_MODELS
SHARED_FLAG=-shared
AM_CONDITIONAL(NEED_APPLE_FIXES,
[test "xx${host/darwin/foundit}xx" != "xx${host}xx"])
if test "xx${host/darwin/foundit}xx" != "xx${host}xx"; then
APPLE_DSO_FLAGS=-Wl,-undefined,dynamic_lookup
SHARED_FLAG=-bundle
fi
AC_SUBST([APPLE_DSO_FLAGS])
AC_SUBST([SHARED_FLAG])
AC_CONFIG_FILES([UnderlyingEvent/Makefile
Models/Makefile
Models/StandardModel/Makefile
Models/RSModel/Makefile
Models/General/Makefile
Models/Susy/Makefile
Models/Susy/NMSSM/Makefile
Models/Susy/RPV/Makefile
Models/UED/Makefile
Models/LH/Makefile
Models/LHTP/Makefile
Models/Transplanckian/Makefile
Models/Leptoquarks/Makefile
Models/Zprime/Makefile
Models/TTbAsymm/Makefile
Models/Feynrules/Makefile
Models/Feynrules/python/Makefile-FR
Models/ADD/Makefile
Models/Sextet/Makefile
Decay/Makefile
Decay/FormFactors/Makefile
Decay/Tau/Makefile
Decay/Baryon/Makefile
Decay/VectorMeson/Makefile
Decay/Perturbative/Makefile
Decay/ScalarMeson/Makefile
Decay/TensorMeson/Makefile
Decay/WeakCurrents/Makefile
Decay/Partonic/Makefile
Decay/General/Makefile
Decay/Radiation/Makefile
Doc/refman.conf
Doc/refman.h
PDT/Makefile
PDF/Makefile
MatrixElement/Makefile
MatrixElement/General/Makefile
MatrixElement/Lepton/Makefile
MatrixElement/Hadron/Makefile
MatrixElement/DIS/Makefile
MatrixElement/Powheg/Makefile
MatrixElement/Gamma/Makefile
MatrixElement/Matchbox/Makefile
MatrixElement/Matchbox/Base/Makefile
MatrixElement/Matchbox/Utility/Makefile
MatrixElement/Matchbox/Phasespace/Makefile
MatrixElement/Matchbox/Dipoles/Makefile
MatrixElement/Matchbox/InsertionOperators/Makefile
MatrixElement/Matchbox/Matching/Makefile
MatrixElement/Matchbox/Cuts/Makefile
MatrixElement/Matchbox/Scales/Makefile
MatrixElement/Matchbox/Scales/MatchboxScale.cc
MatrixElement/Matchbox/ColorFull/Makefile
MatrixElement/Matchbox/CVolver/Makefile
MatrixElement/Matchbox/Builtin/Makefile
MatrixElement/Matchbox/Builtin/Amplitudes/Makefile
MatrixElement/Matchbox/Tests/Makefile
MatrixElement/Matchbox/External/Makefile
MatrixElement/Matchbox/External/BLHAGeneric/Makefile
MatrixElement/Matchbox/External/NLOJet/Makefile
MatrixElement/Matchbox/External/NLOJet/Amplitudes/Makefile
MatrixElement/Matchbox/External/NLOJet/PP2Jets/Makefile
MatrixElement/Matchbox/External/VBFNLO/Makefile
MatrixElement/Matchbox/External/HEJ/Makefile
MatrixElement/Matchbox/External/NJet/Makefile
- MatrixElement/Matchbox/External/NJet/NJetsAmplitude.cc
MatrixElement/Matchbox/External/GoSam/Makefile
- MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc
MatrixElement/Matchbox/External/OpenLoops/Makefile
- MatrixElement/Matchbox/External/OpenLoops/OpenLoopsAmplitude.cc
MatrixElement/Matchbox/External/MadGraph/Makefile
- MatrixElement/Matchbox/External/MadGraph/MadGraphAmplitude.cc
MatrixElement/Matchbox/External/MadGraph/mg2Matchbox.py
Sampling/Makefile
Sampling/CellGrids/Makefile
Shower/Makefile
Shower/Matching/Makefile
DipoleShower/Makefile
DipoleShower/Base/Makefile
DipoleShower/Kernels/Makefile
DipoleShower/Kinematics/Makefile
DipoleShower/Utility/Makefile
DipoleShower/AlphaS/Makefile
Utilities/Makefile
Utilities/XML/Makefile
Utilities/Statistics/Makefile
Hadronization/Makefile
lib/Makefile
include/Makefile
src/Makefile
src/defaults/Makefile
src/snippets/Makefile
src/Matchbox/Makefile
src/herwig-config
Doc/Makefile
Doc/HerwigDefaults.in
Looptools/Makefile
Analysis/Makefile
src/Makefile-UserModules
src/defaults/Analysis.in
src/defaults/MatchboxDefaults.in
src/defaults/setup.gosam.in
src/Matchbox/LO-DefaultShower.in
src/Matchbox/LO-DipoleShower.in
src/Matchbox/MCatLO-DefaultShower.in
src/Matchbox/MCatLO-DipoleShower.in
src/Matchbox/LO-NoShower.in
src/Matchbox/MCatNLO-DefaultShower.in
src/Matchbox/MCatNLO-DipoleShower.in
src/Matchbox/NLO-NoShower.in
src/Matchbox/Powheg-DefaultShower.in
src/Matchbox/Powheg-DipoleShower.in
Contrib/Makefile
Contrib/make_makefiles.sh
Tests/Makefile
Makefile])
AC_CONFIG_LINKS([Doc/BSMlibs.in:Doc/BSMlibs.in])
AC_CONFIG_FILES([Doc/fixinterfaces.pl],[chmod +x Doc/fixinterfaces.pl])
HERWIG_OVERVIEW
AC_CONFIG_COMMANDS([summary],[cat config.herwig])
AC_OUTPUT

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 1:35 PM (20 h, 5 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022978
Default Alt Text
(147 KB)

Event Timeline