Page MenuHomeHEPForge
Files F7879791

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc b/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc
--- a/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc
+++ b/MatrixElement/Matchbox/External/GoSam/GoSamAmplitude.cc
@@ -1,1065 +1,1063 @@
// -*- C++ -*-
//
// GoSamAmplitude.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the 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 "Herwig/API/Filesystem.h"
-#include <boost/progress.hpp>
-
#include <fstream>
#include <sstream>
#include <string>
#include <cstdlib>
#include <exception>
namespace bfs = Herwig::filesystem;
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 GOSAM_PREFIX
#error Makefile.am needs to define GOSAM_PREFIX
#endif
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
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) {
throw Exception()
<< "--------------------------------------------------------------------------------\n"
<< "The following exception occured:\n\n"
<< " " << e.what() << "\n\n"
<< " -> Please create the parent directory of\n"
<< " " << gosamPath << "\n"
<< " manually!\n"
<< "--------------------------------------------------------------------------------\n"
<< Exception::runerror;
}
}
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 gosam2herwig 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 "+bindir_+"/gosam2herwig ";
cmd+=" --usrinfile="+gosamSetupInFileNameInterface;
cmd+=" --infile="+gosamSetupInFileName+".tbu";
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());
if ( factory()->initVerbose() ) {
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() << "GoSamAmplitude: Ninja reduction needs form optimization!\n" << Exception::runerror;
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;
}
if ( factory()->initVerbose() ) {
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() << "GoSamAmplitude: failed to start GoSam" << Exception::runerror;
}
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() << "GoSamAmplitude: Failed to load GoSam. Please check the log file.\n"
<< Exception::runerror;
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() << "GoSamAmplitude: `Best value' schemes are not supported by GoSam"
<< Exception::runerror;
} else if ( SM().ewScheme() == 4 ) { // EW/Scheme mW (uses mW,GF,sin2thetaW) seems not to be supported by GoSam
throw Exception() << "GoSamAmplitude: `mW' scheme is not supported by GoSam"
<< Exception::runerror;
} else if ( SM().ewScheme() == 1 ) { // EW/Scheme GMuScheme (uses mW,mZ,GF)
double in1=getParticleData(ParticleID::Z0)->hardProcessMass()/GeV;
double in2=getParticleData(ParticleID::Wplus)->hardProcessMass()/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)->hardProcessMass()/GeV;
double in2=getParticleData(ParticleID::Wplus)->hardProcessMass()/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)->hardProcessMass()/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);
} else if ( SM().ewScheme() == 7 ) { // EW/Scheme FeynRulesUFO (uses mZ,GF,alpha(mZ))
double in1=getParticleData(ParticleID::Z0)->hardProcessMass()/GeV;
double in2=SM().alphaEMMZ();
double in3=SM().fermiConstant()*GeV2;
OLP_SetParameter((char *)"mass(23)",&in1,&zero,&pStatus);
OLP_SetParameter((char *)"alpha",&in2,&zero,&pStatus);
OLP_SetParameter((char *)"Gf",&in3,&zero,&pStatus);
}
// hand over mass and width of the Higgs
double wH = getParticleData(25)->hardProcessWidth()/GeV;
double mH = getParticleData(25)->hardProcessMass()/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)->hardProcessMass()/GeV > 0.0) massiveParticles.push_back(i);
// with lepton masses
if (theMassiveLeptons && getParticleData(11)->hardProcessMass()/GeV > 0.0) massiveParticles.push_back(11);
if (theMassiveLeptons && getParticleData(13)->hardProcessMass()/GeV > 0.0) massiveParticles.push_back(13);
if (theMassiveLeptons && getParticleData(15)->hardProcessMass()/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)->hardProcessMass()/GeV;
double width=getParticleData(mInt)->hardProcessWidth()/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)->hardProcessMass()/GeV;
// double width=getParticleData(*mID)->hardProcessWidth()/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)->hardProcessMass()/GeV > 0.0) massiveParticles.push_back(i);
// check if leptons have non-zero masses (iff theMassiveLeptons==true)
if (theMassiveLeptons && getParticleData(11)->hardProcessMass()/GeV > 0.0) massiveParticles.push_back(11);
if (theMassiveLeptons && getParticleData(13)->hardProcessMass()/GeV > 0.0) massiveParticles.push_back(13);
if (theMassiveLeptons && getParticleData(15)->hardProcessMass()/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=" +
gosamSetupInFileName+".tbu" + " --destination=" + gosamSourcePath + " " + order + " > " + cwd + folderMatchboxBuild + "gosam-amplitudes.log 2>&1";
std::system(cmd.c_str());
cmd = "python "+bindir_+"/gosam2herwig ";
cmd += " --makelink ";
// cmd += " --makelinkfrom=contract ";
cmd += " --makelinkfrom="+gosamPath+"/"+name()+".OLPContract.lh";
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() << "GoSamAmplitude: Failed to check contractfile. Please check the logfile.\n"
<< Exception::runerror;
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 ( (olpId()[ProcessType::oneLoopInterference]||olpId()[ProcessType::loopInducedME2]) && acc > accuracyTarget ) {
if ( Debug::level > 1 ) {
accuracyFileStream.open(accuracyFile.c_str(),ios::app);
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
<< bindir_ << pkgdatadir_ << GoSamPrefix_;
}
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
>> bindir_ >> pkgdatadir_ >> GoSamPrefix_;
}
// *** 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 \\cite{Cullen:2011xs}, \\cite{Cullen:2014yla}",
"%\\cite{Cullen:2011xs}\n"
"\\bibitem{Cullen:2011xs}\n"
"G.~Cullen et al.,\n"
"``GoSam: A Program for Automated One-Loop Calculations,''\n"
"arXiv:1111.6534 [hep-ph].\n"
"%%CITATION = ARXIV:1111.6534;%%\n"
"%\\cite{Cullen:2014yla}\n"
"\\bibitem{Cullen:2014yla}\n"
"G.~Cullen et al.,\n"
"``GoSaam-2.0: a tool for automated one-loop calculations within the Standard Model and beyond,''\n"
"arXiv:1404.7096 [hep-ph].\n"
"%%CITATION = ARXIV:1404.7096;%%");
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 interfaceCodeExistsYes
(interfaceCodeExists,
"Yes",
"Switch True if Code already exists.",
true);
static SwitchOption interfaceCodeExistsNo
(interfaceCodeExists,
"No",
"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 interfaceisitDRYes
(interfaceisitDR,
"Yes",
"Switch True.",
true);
static SwitchOption interfaceisitDRNo
(interfaceisitDR,
"No",
"Switch False.",
false);
static Switch<GoSamAmplitude,bool> interfaceFormOpt
("FormOpt",
"Switch On/Off formopt",
&GoSamAmplitude::theFormOpt, true, false, false);
static SwitchOption interfaceFormOptYes
(interfaceFormOpt,
"Yes",
"Yes",
true);
static SwitchOption interfaceFormOptNo
(interfaceFormOpt,
"No",
"No",
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 interfaceNinjaYes
(interfaceNinja,
"Yes",
"Yes",
true);
static SwitchOption interfaceNinjaNo
(interfaceNinja,
"No",
"No",
false);
static Switch<GoSamAmplitude,bool> interfaceHiggsEff
("HiggsEff",
"Switch On/Off for effective higgs model.",
&GoSamAmplitude::theHiggsEff, false, false, false);
static SwitchOption interfaceHiggsEffYes
(interfaceHiggsEff,
"Yes",
"Yes",
true);
static SwitchOption interfaceHiggsEffNo
(interfaceHiggsEff,
"No",
"No",
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/MatchboxFactory.cc b/MatrixElement/Matchbox/MatchboxFactory.cc
--- a/MatrixElement/Matchbox/MatchboxFactory.cc
+++ b/MatrixElement/Matchbox/MatchboxFactory.cc
@@ -1,2114 +1,2106 @@
// -*- C++ -*-
//
// MatchboxFactory.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MatchboxFactory class.
//
#include "MatchboxFactory.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Handlers/SamplerBase.h"
#include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig/MatrixElement/Matchbox/Utility/SU2Helper.h"
#include "Herwig/API/RunDirectories.h"
-#include <boost/progress.hpp>
+#include "Herwig/Utilities/Progress.h"
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
using std::ostream_iterator;
MatchboxFactory::MatchboxFactory()
: SubProcessHandler(), theNLight(0),
theOrderInAlphaS(0), theOrderInAlphaEW(0),
theBornContributions(true), theVirtualContributions(true),
theRealContributions(true), theIndependentVirtuals(false),
theIndependentPKs(false),
theFactorizationScaleFactor(1.0), theRenormalizationScaleFactor(1.0),
theFixedCouplings(false), theFixedQEDCouplings(false), theVetoScales(false),
theDipoleSet(0), theVerbose(false), theInitVerbose(false),
theSubtractionData(""), theSubtractionPlotType(1), theSubtractionScatterPlot(false),
thePoleData(""), theRealEmissionScales(false), theAllProcesses(false),
theMECorrectionsOnly(false), theLoopSimCorrections(false), ranSetup(false),
theFirstPerturbativePDF(true), theSecondPerturbativePDF(true),
inProductionMode(false), theSpinCorrelations(false),theAlphaParameter(1.),
theEnforceChargeConservation(true), theEnforceColourConservation(false),
theEnforceLeptonNumberConservation(false), theEnforceQuarkNumberConservation(false),
theLeptonFlavourDiagonal(false), theQuarkFlavourDiagonal(false) {}
MatchboxFactory::~MatchboxFactory() {}
Ptr<MatchboxFactory>::tptr MatchboxFactory::theCurrentFactory = Ptr<MatchboxFactory>::tptr();
bool& MatchboxFactory::theIsMatchboxRun() {
static bool flag = false;
return flag;
}
IBPtr MatchboxFactory::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxFactory::fullclone() const {
return new_ptr(*this);
}
void MatchboxFactory::prepareME(Ptr<MatchboxMEBase>::ptr me) {
Ptr<MatchboxAmplitude>::ptr amp =
dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>((*me).amplitude());
me->matchboxAmplitude(amp);
if ( phasespace() && !me->phasespace() )
me->phasespace(phasespace());
if ( scaleChoice() && !me->scaleChoice() )
me->scaleChoice(scaleChoice());
if ( !reweighters().empty() ) {
for ( vector<ReweightPtr>::const_iterator rw = reweighters().begin();
rw != reweighters().end(); ++rw )
me->addReweighter(*rw);
}
if ( !preweighters().empty() ) {
for ( vector<ReweightPtr>::const_iterator rw = preweighters().begin();
rw != preweighters().end(); ++rw )
me->addPreweighter(*rw);
}
}
string pid(const PDVector& key) {
ostringstream res;
res << "[" << key[0]->PDGName() << ","
<< key[1]->PDGName() << "->";
for ( PDVector::const_iterator k =
key.begin() + 2; k != key.end(); ++k )
res << (**k).PDGName() << (k != --key.end() ? "," : "");
res << "]";
return res.str();
}
vector<Ptr<MatchboxMEBase>::ptr> MatchboxFactory::
makeMEs(const vector<string>& proc, unsigned int orderas, bool virt) {
generator()->log() << "determining subprocesses for ";
copy(proc.begin(),proc.end(),ostream_iterator<string>(generator()->log()," "));
generator()->log() << "\n" << flush;
map<Ptr<MatchboxAmplitude>::ptr,set<Process> > ampProcs;
map<Process,set<Ptr<MatchboxAmplitude>::ptr> > procAmps;
set<PDVector> processes = makeSubProcesses(proc);
set<PDVector> colouredProcesses;
for ( set<PDVector>::const_iterator pr = processes.begin();
pr != processes.end(); ++pr ) {
for ( PDVector::const_iterator pp = pr->begin();
pp != pr->end(); ++pp ) {
if ( (**pp).coloured() ) {
colouredProcesses.insert(*pr);
break;
}
}
}
if ( colouredProcesses.size() != processes.size() &&
(virtualContributions() || realContributions()) ) {
// NLO not working for non coloured legs
throw Exception()
<< "Found processes without coloured legs.\n"
<< "We currently do not support NLO corrections for those processes.\n"
<< "Please switch to a setup for LO production."
<< Exception::runerror;
}
// detect external particles with non-zero width for the hard process
bool trouble = false;
string troubleMaker;
for ( set<PDVector>::const_iterator pr = processes.begin();
pr != processes.end(); ++pr ) {
for ( PDVector::const_iterator pp = pr->begin();
pp != pr->end(); ++pp ) {
if ( (**pp).hardProcessWidth() != ZERO ) {
trouble = true;
troubleMaker = (**pp).PDGName();
break;
}
}
}
if ( trouble ) {
throw Exception()
<< "MatchboxFactory::makeMEs(): Particle '"
<< troubleMaker << "' appears as external\nprocess leg with non-zero "
<< "width to be used in the hard process calculation.\n"
<< "Please check your setup and consider setting HardProcessWidth to zero."
<< Exception::runerror;
}
vector<Ptr<MatchboxAmplitude>::ptr> matchAmplitudes;
unsigned int lowestAsOrder =
allProcesses() ? 0 : orderas;
unsigned int highestAsOrder = orderas;
unsigned int lowestAeOrder =
allProcesses() ? 0 : orderInAlphaEW();
unsigned int highestAeOrder = orderInAlphaEW();
for ( unsigned int oas = lowestAsOrder; oas <= highestAsOrder; ++oas ) {
for ( unsigned int oae = lowestAeOrder; oae <= highestAeOrder; ++oae ) {
for ( vector<Ptr<MatchboxAmplitude>::ptr>::const_iterator amp
= amplitudes().begin(); amp != amplitudes().end(); ++amp ) {
if ( !theSelectedAmplitudes.empty() ) {
if ( find(theSelectedAmplitudes.begin(),theSelectedAmplitudes.end(),*amp)
== theSelectedAmplitudes.end() )
continue;
}
if ( !theDeselectedAmplitudes.empty() ) {
if ( find(theDeselectedAmplitudes.begin(),theDeselectedAmplitudes.end(),*amp)
!= theDeselectedAmplitudes.end() )
continue;
}
(**amp).orderInGs(oas);
(**amp).orderInGem(oae);
if ( (**amp).orderInGs() != oas ||
(**amp).orderInGem() != oae ) {
continue;
}
matchAmplitudes.push_back(*amp);
}
}
}
size_t combinations = processes.size()*matchAmplitudes.size();
size_t procCount = 0;
generator()->log() << "building matrix elements." << flush;
- boost::progress_display * progressBar =
- new boost::progress_display(combinations,generator()->log());
+ progress_display progressBar{ combinations, generator()->log() };
for ( unsigned int oas = lowestAsOrder; oas <= highestAsOrder; ++oas ) {
for ( unsigned int oae = lowestAeOrder; oae <= highestAeOrder; ++oae ) {
for ( vector<Ptr<MatchboxAmplitude>::ptr>::const_iterator amp
= matchAmplitudes.begin(); amp != matchAmplitudes.end(); ++amp ) {
(**amp).orderInGs(oas);
(**amp).orderInGem(oae);
for ( set<PDVector>::const_iterator p = processes.begin();
p != processes.end(); ++p ) {
- ++(*progressBar);
+ ++progressBar;
if ( !(**amp).canHandle(*p,this,virt) )
continue;
if ( (**amp).isExternal() )
externalAmplitudes().insert(*amp);
++procCount;
Process proc(*p,oas,oae);
ampProcs[*amp].insert(proc);
procAmps[proc].insert(*amp);
}
}
}
}
- delete progressBar;
generator()->log() << flush;
bool clash = false;
for ( map<Process,set<Ptr<MatchboxAmplitude>::ptr> >::const_iterator check =
procAmps.begin(); check != procAmps.end(); ++check ) {
if ( check->second.size() > 1 ) {
clash = true;
generator()->log() << "Several different amplitudes have been found for: "
<< check->first.legs[0]->PDGName() << " "
<< check->first.legs[1]->PDGName() << " -> ";
for ( PDVector::const_iterator p = check->first.legs.begin() + 2;
p != check->first.legs.end(); ++p )
generator()->log() << (**p).PDGName() << " ";
generator()->log() << "at alpha_s^" << check->first.orderInAlphaS
<< " and alpha_ew^" << check->first.orderInAlphaEW
<< "\n";
generator()->log() << "The following amplitudes claim responsibility:\n";
for ( set<Ptr<MatchboxAmplitude>::ptr>::const_iterator a = check->second.begin();
a != check->second.end(); ++a ) {
generator()->log() << (**a).name() << " ";
}
generator()->log() << "\n";
}
}
if ( clash ) {
throw Exception() << "MatchboxFactory: Ambiguous amplitude setup - please check your input files.\n"
<< "To avoid this problem use the SelectAmplitudes or DeselectAmplitudes interfaces.\n"
<< Exception::runerror;
}
bool canDoSpinCorrelations = true;
vector<Ptr<MatchboxMEBase>::ptr> res;
for ( map<Ptr<MatchboxAmplitude>::ptr,set<Process> >::const_iterator
ap = ampProcs.begin(); ap != ampProcs.end(); ++ap ) {
canDoSpinCorrelations &= ap->first->canFillRhoMatrix();
for ( set<Process>::const_iterator m = ap->second.begin();
m != ap->second.end(); ++m ) {
Ptr<MatchboxMEBase>::ptr me = ap->first->makeME(m->legs);
me->subProcess() = *m;
me->amplitude(ap->first);
me->matchboxAmplitude(ap->first);
prepareME(me);
string pname = "ME" + pid(m->legs);
if ( ! (generator()->preinitRegister(me,pname) ) )
throw Exception() << "MatchboxFactory: Matrix element " << pname << " already existing."
<< Exception::runerror;
if ( me->diagrams().empty() )continue;
res.push_back(me);
if ( theFirstPerturbativePDF )
theIncoming.insert(m->legs[0]->id());
if ( theSecondPerturbativePDF )
theIncoming.insert(m->legs[1]->id());
}
}
if ( spinCorrelations() && !canDoSpinCorrelations ) {
generator()->log() << "Warning: Spin correlations have been requested, but no amplitude is "
<< "capable of performing these.\n";
theSpinCorrelations = false;
}
generator()->log() << "created "
<< procCount << " subprocesses.\n";
generator()->log() << "---------------------------------------------------\n"
<< flush;
return res;
}
int MatchboxFactory::orderOLPProcess(const Process& proc,
Ptr<MatchboxAmplitude>::tptr amp,
int type) {
map<pair<Process,int>,int>& procs =
olpProcesses()[amp];
map<pair<Process,int>,int>::const_iterator it =
procs.find(make_pair(proc,type));
if ( it != procs.end() )
return it->second;
int id = procs.size();
procs[make_pair(proc,type)] = id + 1;
return id + 1;
}
void MatchboxFactory::productionMode() {
if ( inProductionMode )
return;
if ( !bornContributions() && !virtualContributions() && !realContributions() )
throw Exception() << "MatchboxFactory: At least one cross section contribution needs to be enabled.\n"
<< "Please check your setup.\n"
<< Exception::runerror;
bool needTrueVirtuals =
virtualContributions() && !meCorrectionsOnly() && !loopSimCorrections();
for ( vector<Ptr<MatchboxAmplitude>::ptr>::iterator amp
= amplitudes().begin(); amp != amplitudes().end(); ++amp ) {
if ( !needTrueVirtuals && (**amp).oneLoopAmplitude() ) {
Repository::clog() << "One-loop contributions from '"
<< (**amp).name()
<< "' are not required and will be disabled.\n"
<< flush;
(**amp).disableOneLoop();
}
}
if ( subtractionData() != "" && !subProcessGroups() ) {
throw Exception() << "MatchboxFactory: Plain NLO settings are required for subtraction checks.\n"
<< "Please check your setup.\n"
<< Exception::runerror;
}
if ( showerApproximation() && !virtualContributions() && !realContributions() ) {
Repository::clog() << "Warning: Matching requested for LO run. Matching disabled.\n" << flush;
showerApproximation(Ptr<ShowerApproximation>::tptr());
}
if ( showerApproximation() && (subtractionData() != "" || subProcessGroups()) ) {
Repository::clog() << "Warning: Matching requested for plain NLO run. Matching disabled.\n" << flush;
showerApproximation(Ptr<ShowerApproximation>::tptr());
}
if ( showerApproximation() ) {
if ( spinCorrelations() && !showerApproximation()->hasSpinCorrelations() ) {
Repository::clog() << "Warning: Spin correlations have been requested but the matching "
<< "object is not capable of these. Spin correlations will be turned of.\n"
<< flush;
theSpinCorrelations = false;
}
}
inProductionMode = true;
}
void MatchboxFactory::setup() {
useMe();
if ( !ranSetup ) {
if ( !inProductionMode )
throw Exception() << "MatchboxFactory: The MatchboxFactory object '"
<< name() << "' has not been switched to production mode.\n"
<< "Did you use 'do "
<< name() << ":ProductionMode' before isolating the event generator?\n"
<< Exception::runerror;
olpProcesses().clear();
externalAmplitudes().clear();
theHighestVirtualsize = 0;
theIncoming.clear();
bool needTrueVirtuals =
virtualContributions() && !meCorrectionsOnly() && !loopSimCorrections();
if ( bornMEs().empty() ) {
if ( particleGroups().find("j") == particleGroups().end() )
throw Exception() << "MatchboxFactory: Could not find a jet particle group named 'j'"
<< Exception::runerror;
// rebind the particle data objects
for ( map<string,PDVector>::iterator g = particleGroups().begin();
g != particleGroups().end(); ++g )
for ( PDVector::iterator p = g->second.begin();
p != g->second.end(); ++p ) {
#ifndef NDEBUG
long checkid = (**p).id();
#endif
*p = getParticleData((**p).id());
assert((**p).id() == checkid);
}
const PDVector& partons = particleGroups()["j"];
unsigned int nl = 0;
for ( PDVector::const_iterator p = partons.begin();
p != partons.end(); ++p ) {
if ( abs((**p).id()) < 7 && (**p).hardProcessMass() == ZERO )
++nl;
if ( (**p).id() > 0 && (**p).id() < 7 && (**p).hardProcessMass() == ZERO )
nLightJetVec( (**p).id() );
if ( (**p).id() > 0 && (**p).id() < 7 && (**p).hardProcessMass() != ZERO )
nHeavyJetVec( (**p).id() );
}
nLight(nl/2);
if ( particleGroups().find("p") == particleGroups().end() )
throw Exception() << "MatchboxFactory: Could not find a hadron particle group named 'p'"
<< Exception::runerror;
const PDVector& partonsInP = particleGroups()["p"];
for ( PDVector::const_iterator pip = partonsInP.begin();
pip != partonsInP.end(); ++pip ) {
if ( (**pip).id() > 0 && (**pip).id() < 7 && (**pip).hardProcessMass() == ZERO )
nLightProtonVec( (**pip).id() );
}
vector<Ptr<MatchboxMEBase>::ptr> mes;
for ( vector<vector<string> >::const_iterator p = processes.begin();
p != processes.end(); ++p ) {
if( needTrueVirtuals ) {
theHighestVirtualsize = max(theHighestVirtualsize,(int((*p).size())));
}
mes = makeMEs(*p,orderInAlphaS(),needTrueVirtuals);
copy(mes.begin(),mes.end(),back_inserter(bornMEs()));
if ( realContributions() || meCorrectionsOnly() ||
(showerApproximation() && virtualContributions()) ||
(showerApproximation() && loopSimCorrections()) ) {
if ( realEmissionProcesses.empty() ) {
vector<string> rproc = *p;
rproc.push_back("j");
mes = makeMEs(rproc,orderInAlphaS()+1,false);
copy(mes.begin(),mes.end(),back_inserter(realEmissionMEs()));
}
}
}
if ( realContributions() || meCorrectionsOnly() ||
(showerApproximation() && virtualContributions()) ||
(showerApproximation() && loopSimCorrections()) ) {
if ( !realEmissionProcesses.empty() ) {
for ( vector<vector<string> >::const_iterator q =
realEmissionProcesses.begin(); q != realEmissionProcesses.end(); ++q ) {
mes = makeMEs(*q,orderInAlphaS()+1,false);
copy(mes.begin(),mes.end(),back_inserter(realEmissionMEs()));
}
}
}
}
if ( loopInducedMEs().empty() ) {
for ( vector<vector<string> >::const_iterator p = loopInducedProcesses.begin();
p != loopInducedProcesses.end(); ++p ) {
vector<Ptr<MatchboxMEBase>::ptr> mes = makeMEs(*p,orderInAlphaS(),false);
copy(mes.begin(),mes.end(),back_inserter(loopInducedMEs()));
}
}
if( bornMEs().empty() && realEmissionMEs().empty() && loopInducedMEs().empty() )
throw Exception() << "MatchboxFactory: No matrix elements have been found.\n\
Please check if your order of Alpha_s and Alpha_ew have the right value.\n"
<< Exception::runerror;
// check if we have virtual contributions
bool haveVirtuals = true;
// check DR conventions of virtual contributions
bool virtualsAreDR = false;
bool virtualsAreCDR = false;
// check finite term conventions of virtual contributions
bool virtualsAreCS = false;
bool virtualsAreBDK = false;
bool virtualsAreExpanded = false;
// renormalization scheme
bool virtualsAreDRbar = false;
// check and prepare the Born and virtual matrix elements
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
prepareME(*born);
haveVirtuals &= (**born).haveOneLoop();
if ( needTrueVirtuals ) {
if ( (**born).haveOneLoop() ) {
virtualsAreDRbar |= (**born).isDRbar();
virtualsAreDR |= (**born).isDR();
virtualsAreCDR |= !(**born).isDR();
virtualsAreCS |= (**born).isCS();
virtualsAreBDK |= (**born).isBDK();
virtualsAreExpanded |= (**born).isExpanded();
}
}
}
// prepare the loop induced matrix elements
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator looped
= loopInducedMEs().begin(); looped != loopInducedMEs().end(); ++looped ) {
prepareME(*looped);
}
if ( needTrueVirtuals ) {
// check the additional insertion operators
if ( !virtuals().empty() )
haveVirtuals = true;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= virtuals().begin(); virt != virtuals().end(); ++virt ) {
virtualsAreDRbar |= (**virt).isDRbar();
virtualsAreDR |= (**virt).isDR();
virtualsAreCDR |= !(**virt).isDR();
virtualsAreCS |= (**virt).isCS();
virtualsAreBDK |= (**virt).isBDK();
virtualsAreExpanded |= (**virt).isExpanded();
}
// check for consistent conventions on virtuals, if we are to include them
if ( virtualContributions() ) {
if ( !haveVirtuals ) {
throw Exception() << "MatchboxFactory: Could not find amplitudes for all virtual contributions needed.\n"
<< Exception::runerror;
}
if ( virtualsAreDR && virtualsAreCDR ) {
throw Exception() << "MatchboxFactory: Virtual corrections use inconsistent regularization schemes.\n"
<< Exception::runerror;
}
if ( (virtualsAreCS && virtualsAreBDK) ||
(virtualsAreCS && virtualsAreExpanded) ||
(virtualsAreBDK && virtualsAreExpanded) ||
(!virtualsAreCS && !virtualsAreBDK && !virtualsAreExpanded) ) {
throw Exception() << "MatchboxFactory: Virtual corrections use inconsistent conventions on finite terms.\n"
<< Exception::runerror;
}
}
}
// prepare the real emission matrix elements
if ( realContributions() || meCorrectionsOnly() ||
(showerApproximation() && virtualContributions()) ||
(showerApproximation() && loopSimCorrections()) ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator real
= realEmissionMEs().begin(); real != realEmissionMEs().end(); ++real ) {
prepareME(*real);
}
}
// start creating matrix elements
MEs().clear();
// setup born and virtual contributions
if ( bornContributions() || virtualContributions() ) {
generator()->log() << "preparing Born"
<< (virtualContributions() ? " and virtual" : "")
<< " matrix elements.\n" << flush;
}
if ( (bornContributions() && !virtualContributions()) ||
(bornContributions() && meCorrectionsOnly()) ||
(bornContributions() && virtualContributions() && independentVirtuals()) ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
if ( (**born).onlyOneLoop() )
continue;
Ptr<MatchboxMEBase>::ptr bornme = (**born).cloneMe();
string pname = fullName() + "/" + (**born).name();
if ( virtualContributions() && independentVirtuals() )
pname += ".Born";
if ( ! (generator()->preinitRegister(bornme,pname) ) )
throw Exception() << "MatchboxFactory: Matrix element " << pname << " already existing."
<< Exception::runerror;
if ( bornme->isOLPTree() ) {
int id = orderOLPProcess(bornme->subProcess(),
(**born).matchboxAmplitude(),
ProcessType::treeME2);
bornme->olpProcess(ProcessType::treeME2,id);
}
bornme->needsNoCorrelations();
bornme->cloneDependencies();
MEs().push_back(bornme);
}
}
if ( bornContributions() && !loopInducedMEs().empty() ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator looped
= loopInducedMEs().begin(); looped != loopInducedMEs().end(); ++looped ) {
Ptr<MatchboxMEBase>::ptr loopme = (**looped).cloneMe();
string pname = fullName() + "/" + (**looped).name() + ".LoopInduced";
if ( ! (generator()->preinitRegister(loopme,pname) ) )
throw Exception() << "MatchboxFactory: Matrix element " << pname << " already existing."
<< Exception::runerror;
if ( loopme->isOLPTree() ) {
int id = orderOLPProcess(loopme->subProcess(),
(**looped).matchboxAmplitude(),
ProcessType::loopInducedME2);
loopme->olpProcess(ProcessType::loopInducedME2,id);
}
loopme->needsNoCorrelations();
loopme->cloneDependencies();
MEs().push_back(loopme);
}
}
if ( needTrueVirtuals ) {
bornVirtualMEs().clear();
- boost::progress_display * progressBar =
- new boost::progress_display(bornMEs().size(),generator()->log());
+ progress_display progressBar{ bornMEs().size(), generator()->log() };
if ( thePoleData != "" )
if ( thePoleData[thePoleData.size()-1] != '/' )
thePoleData += "/";
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
Ptr<MatchboxMEBase>::ptr nlo = (**born).cloneMe();
string pname = fullName() + "/" + (**born).name();
if ( !independentVirtuals() && !(!bornContributions() && virtualContributions()) )
pname += ".BornVirtual";
else if ( independentPKs() && !nlo->onlyOneLoop() )
pname += ".VirtualVI";
else
pname += ".Virtual";
if ( ! (generator()->preinitRegister(nlo,pname) ) )
throw Exception() << "MatchboxFactory: NLO ME " << pname << " already existing."
<< Exception::runerror;
nlo->virtuals().clear();
if ( !nlo->onlyOneLoop() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= virtuals().begin(); virt != virtuals().end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlo->virtuals().push_back(*virt);
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= DipoleRepository::insertionIOperators(dipoleSet()).begin();
virt != DipoleRepository::insertionIOperators(dipoleSet()).end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlo->virtuals().push_back(*virt);
}
if ( !independentVirtuals() || ( independentVirtuals() && !independentPKs() ) ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= DipoleRepository::insertionPKOperators(dipoleSet()).begin();
virt != DipoleRepository::insertionPKOperators(dipoleSet()).end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlo->virtuals().push_back(*virt);
}
}
if ( nlo->virtuals().empty() )
throw Exception() << "MatchboxFactory: No insertion operators have been found for "
<< (**born).name() << "\n"
<< Exception::runerror;
if ( checkPoles() ) {
if ( !virtualsAreExpanded ) {
throw Exception()
<< "MatchboxFactory: Cannot check epsilon poles if virtuals are not in `expanded' convention.\n"
<< Exception::runerror;
}
}
}
if ( !bornContributions() || independentVirtuals() ) {
nlo->doOneLoopNoBorn();
} else {
nlo->doOneLoop();
}
if ( nlo->isOLPLoop() ) {
int id = orderOLPProcess(nlo->subProcess(),
(**born).matchboxAmplitude(),
ProcessType::oneLoopInterference);
nlo->olpProcess(ProcessType::oneLoopInterference,id);
if ( !nlo->onlyOneLoop() && nlo->needsOLPCorrelators() ) {
id = orderOLPProcess(nlo->subProcess(),
(**born).matchboxAmplitude(),
ProcessType::colourCorrelatedME2);
nlo->olpProcess(ProcessType::colourCorrelatedME2,id);
}
}
nlo->needsCorrelations();
nlo->cloneDependencies();
bornVirtualMEs().push_back(nlo);
MEs().push_back(nlo);
if ( independentVirtuals() && independentPKs() && !nlo->onlyOneLoop() ) {
Ptr<MatchboxMEBase>::ptr nlopk = (**born).cloneMe();
string pnamepk = fullName() + "/" + (**born).name();
pnamepk += ".VirtualPK";
if ( ! (generator()->preinitRegister(nlopk,pnamepk) ) )
throw Exception() << "MatchboxFactory: NLO ME " << pnamepk << " already existing."
<< Exception::runerror;
nlopk->virtuals().clear();
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= DipoleRepository::insertionPKOperators(dipoleSet()).begin();
virt != DipoleRepository::insertionPKOperators(dipoleSet()).end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlopk->virtuals().push_back(*virt);
}
if ( !nlopk->virtuals().empty() ) {
nlopk->doOneLoopNoBorn();
nlopk->doOneLoopNoLoops();
if ( nlopk->isOLPLoop() ) {
int id = orderOLPProcess(nlopk->subProcess(),
(**born).matchboxAmplitude(),
ProcessType::treeME2);
nlopk->olpProcess(ProcessType::treeME2,id);
if ( nlopk->needsOLPCorrelators() ) {
id = orderOLPProcess(nlopk->subProcess(),
(**born).matchboxAmplitude(),
ProcessType::colourCorrelatedME2);
nlopk->olpProcess(ProcessType::colourCorrelatedME2,id);
}
}
nlopk->needsCorrelations();
nlopk->cloneDependencies();
bornVirtualMEs().push_back(nlopk);
MEs().push_back(nlopk);
}
}
- ++(*progressBar);
+ ++progressBar;
}
- delete progressBar;
-
generator()->log() << "---------------------------------------------------\n"
<< flush;
}
theSplittingDipoles.clear();
set<cPDVector> bornProcs;
if ( showerApproximation() ) {
if ( showerApproximation()->needsSplittingGenerator() ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born )
for ( MEBase::DiagramVector::const_iterator d = (**born).diagrams().begin();
d != (**born).diagrams().end(); ++d )
bornProcs.insert((**d).partons());
}
}
if ( realContributions() || meCorrectionsOnly() ||
(showerApproximation() && virtualContributions()) ||
(showerApproximation() && loopSimCorrections()) ) {
generator()->log() << "preparing subtracted matrix elements.\n" << flush;
if ( theSubtractionData != "" )
if ( theSubtractionData[theSubtractionData.size()-1] != '/' )
theSubtractionData += "/";
subtractedMEs().clear();
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
if ( (**born).onlyOneLoop() )
continue;
(**born).needsCorrelations();
if ( (**born).isOLPTree() ) {
int id = orderOLPProcess((**born).subProcess(),
(**born).matchboxAmplitude(),
ProcessType::colourCorrelatedME2);
(**born).olpProcess(ProcessType::colourCorrelatedME2,id);
bool haveGluon = false;
for ( PDVector::const_iterator p = (**born).subProcess().legs.begin();
p != (**born).subProcess().legs.end(); ++p )
if ( (**p).id() == 21 ) {
haveGluon = true;
break;
}
if ( haveGluon ) {
id = orderOLPProcess((**born).subProcess(),
(**born).matchboxAmplitude(),
ProcessType::spinColourCorrelatedME2);
(**born).olpProcess(ProcessType::spinColourCorrelatedME2,id);
}
if ( showerApproximation() ) {
id = orderOLPProcess((**born).subProcess(),
(**born).matchboxAmplitude(),
ProcessType::treeME2);
(**born).olpProcess(ProcessType::treeME2,id);
}
}
}
- boost::progress_display * progressBar =
- new boost::progress_display(realEmissionMEs().size(),generator()->log());
+ progress_display progressBar{ realEmissionMEs().size(), generator()->log() };
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator real
= realEmissionMEs().begin(); real != realEmissionMEs().end(); ++real ) {
Ptr<SubtractedME>::ptr sub = new_ptr(SubtractedME());
string pname = fullName() + "/" + (**real).name() + ".SubtractedReal";
if ( ! (generator()->preinitRegister(sub,pname) ) )
throw Exception() << "MatchboxFactory: Subtracted ME " << pname << " already existing."
<< Exception::runerror;
(**real).needsNoCorrelations();
if ( (**real).isOLPTree() ) {
int id = orderOLPProcess((**real).subProcess(),
(**real).matchboxAmplitude(),
ProcessType::treeME2);
(**real).olpProcess(ProcessType::treeME2,id);
}
sub->head(*real);
sub->dependent().clear();
sub->getDipoles();
if ( sub->dependent().empty() ) {
// finite real contribution
if ( realContributions() ) {
Ptr<MatchboxMEBase>::ptr fme =
dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(sub->head())->cloneMe();
string qname = fullName() + "/" + (**real).name() + ".FiniteReal";
if ( ! (generator()->preinitRegister(fme,qname) ) )
throw Exception() << "MatchboxFactory: ME " << qname << " already existing."
<< Exception::runerror;
MEs().push_back(fme);
finiteRealMEs().push_back(fme);
}
sub->head(tMEPtr());
- ++(*progressBar);
+ ++progressBar;
continue;
}
if ( realEmissionScales() )
sub->doRealEmissionScales();
subtractedMEs().push_back(sub);
if ( realContributions() )
if ( !showerApproximation() || (showerApproximation() && showerApproximation()->hasHEvents()) )
MEs().push_back(sub);
if ( showerApproximation() ) {
if ( virtualContributions() && !meCorrectionsOnly() && !loopSimCorrections() ) {
Ptr<SubtractedME>::ptr subv = new_ptr(*sub);
string vname = sub->fullName() + ".SubtractionIntegral";
if ( ! (generator()->preinitRegister(subv,vname) ) )
throw Exception() << "MatchboxFactory: Subtracted ME " << vname << " already existing."
<< Exception::runerror;
subv->cloneDependencies(vname);
subv->doVirtualShowerSubtraction();
subtractedMEs().push_back(subv);
MEs().push_back(subv);
}
if ( loopSimCorrections() ) {
Ptr<SubtractedME>::ptr subv = new_ptr(*sub);
string vname = sub->fullName() + ".SubtractionIntegral";
if ( ! (generator()->preinitRegister(subv,vname) ) )
throw Exception() << "MatchboxFactory: Subtracted ME " << vname << " already existing."
<< Exception::runerror;
subv->cloneDependencies(vname);
subv->doLoopSimSubtraction();
subtractedMEs().push_back(subv);
MEs().push_back(subv);
}
sub->doRealShowerSubtraction();
if ( showerApproximation()->needsSplittingGenerator() )
for ( set<cPDVector>::const_iterator p = bornProcs.begin();
p != bornProcs.end(); ++p ) {
vector<Ptr<SubtractionDipole>::ptr> sdip = sub->splitDipoles(*p);
set<Ptr<SubtractionDipole>::ptr>& dips = theSplittingDipoles[*p];
copy(sdip.begin(),sdip.end(),inserter(dips,dips.begin()));
}
}
- ++(*progressBar);
+ ++progressBar;
}
- delete progressBar;
-
generator()->log() << "---------------------------------------------------\n"
<< flush;
}
if ( !theSplittingDipoles.empty() ) {
map<Ptr<SubtractionDipole>::ptr,Ptr<SubtractionDipole>::ptr> cloneMap;
for ( map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >::const_iterator sd = theSplittingDipoles.begin();
sd != theSplittingDipoles.end(); ++sd ) {
for ( set<Ptr<SubtractionDipole>::ptr>::const_iterator d = sd->second.begin();
d != sd->second.end(); ++d ) {
cloneMap[*d] = Ptr<SubtractionDipole>::ptr();
}
}
for ( map<Ptr<SubtractionDipole>::ptr,Ptr<SubtractionDipole>::ptr>::iterator cd =
cloneMap.begin(); cd != cloneMap.end(); ++cd ) {
Ptr<SubtractionDipole>::ptr cloned = cd->first->cloneMe();
string dname = cd->first->fullName() + ".splitting";
if ( ! (generator()->preinitRegister(cloned,dname)) )
throw Exception() << "MatchboxFactory: Dipole '" << dname << "' already existing."
<< Exception::runerror;
cloned->cloneDependencies();
cloned->showerApproximation(Ptr<ShowerApproximation>::tptr());
cloned->doSplitting();
cd->second = cloned;
}
for ( map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >::iterator sd = theSplittingDipoles.begin();
sd != theSplittingDipoles.end(); ++sd ) {
set<Ptr<SubtractionDipole>::ptr> cloned;
for ( set<Ptr<SubtractionDipole>::ptr>::iterator d = sd->second.begin();
d != sd->second.end(); ++d ) {
cloned.insert(cloneMap[*d]);
}
sd->second = cloned;
}
}
if ( !externalAmplitudes().empty() ) {
generator()->log() << "Initializing external amplitudes.\n" << flush;
for ( set<Ptr<MatchboxAmplitude>::tptr>::const_iterator ext =
externalAmplitudes().begin(); ext != externalAmplitudes().end(); ++ext ) {
if ( !(**ext).initializeExternal() ) {
throw Exception() << "Failed to initialize amplitude '" << (**ext).name() << "'\n"
<< Exception::runerror;
}
}
generator()->log() << "---------------------------------------------------\n"
<< flush;
}
if ( !olpProcesses().empty() ) {
generator()->log() << "Initializing one-loop provider(s).\n" << flush;
map<Ptr<MatchboxAmplitude>::tptr,map<pair<Process,int>,int> > olps;
for ( map<Ptr<MatchboxAmplitude>::tptr,map<pair<Process,int>,int> >::const_iterator
oit = olpProcesses().begin(); oit != olpProcesses().end(); ++oit ) {
olps[oit->first] = oit->second;
}
for ( map<Ptr<MatchboxAmplitude>::tptr,map<pair<Process,int>,int> >::const_iterator
olpit = olps.begin(); olpit != olps.end(); ++olpit ) {
if ( !olpit->first->startOLP(olpit->second) ) {
throw Exception() << "MatchboxFactory: Failed to start OLP for amplitude '" << olpit->first->name() << "'\n"
<< Exception::runerror;
}
}
generator()->log() << "---------------------------------------------------\n"
<< flush;
}
generator()->log() << "Process setup finished.\n" << flush;
ranSetup = true;
}
}
void MatchboxFactory::SplittingChannel::print(ostream& os) const {
os << "--- SplittingChannel setup -----------------------------------------------------\n";
os << " Born process ";
const StandardXComb& bxc = *bornXComb;
os << bxc.mePartonData()[0]->PDGName() << " "
<< bxc.mePartonData()[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = bxc.mePartonData().begin() + 2;
p != bxc.mePartonData().end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "\n";
os << " to real emission process ";
const StandardXComb& rxc = *realXComb;
os << rxc.mePartonData()[0]->PDGName() << " "
<< rxc.mePartonData()[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = rxc.mePartonData().begin() + 2;
p != rxc.mePartonData().end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "\n";
os << " with dipole:\n";
dipole->print(os);
os << "---------------------------------------------------\n";
os << flush;
}
list<MatchboxFactory::SplittingChannel>
MatchboxFactory::getSplittingChannels(tStdXCombPtr xcptr) const {
if ( xcptr->lastProjector() )
xcptr = xcptr->lastProjector();
const StandardXComb& xc = *xcptr;
cPDVector proc = xc.mePartonData();
map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >::const_iterator splitEntries
= splittingDipoles().find(proc);
list<SplittingChannel> res;
if ( splitEntries == splittingDipoles().end() )
return res;
const set<Ptr<SubtractionDipole>::ptr>& splitDipoles = splitEntries->second;
SplittingChannel channel;
if ( !splitDipoles.empty() ) {
Ptr<MatchboxMEBase>::tptr bornME =
const_ptr_cast<Ptr<MatchboxMEBase>::tptr>((**splitDipoles.begin()).underlyingBornME());
channel.bornXComb =
bornME->makeXComb(xc.maxEnergy(),xc.particles(),xc.eventHandlerPtr(),
const_ptr_cast<tSubHdlPtr>(xc.subProcessHandler()),
xc.pExtractor(),xc.CKKWHandler(),
xc.partonBins(),xc.cuts(),xc.diagrams(),xc.mirror(),
PartonPairVec());
}
for ( set<Ptr<SubtractionDipole>::ptr>::const_iterator sd =
splitDipoles.begin(); sd != splitDipoles.end(); ++sd ) {
channel.dipole = *sd;
vector<StdXCombPtr> realXCombs = (**sd).makeRealXCombs(channel.bornXComb);
for ( vector<StdXCombPtr>::const_iterator rxc = realXCombs.begin();
rxc != realXCombs.end(); ++rxc ) {
channel.realXComb = *rxc;
if ( showerApproximation()->needsTildeXCombs() ) {
channel.tildeXCombs.clear();
assert(!channel.dipole->partnerDipoles().empty());
for ( vector<Ptr<SubtractionDipole>::tptr>::const_iterator p =
channel.dipole->partnerDipoles().begin();
p != channel.dipole->partnerDipoles().end(); ++p ) {
StdXCombPtr txc = channel.dipole->makeBornXComb(channel.realXComb);
if ( txc )
channel.tildeXCombs.push_back(txc);
}
}
res.push_back(channel);
}
}
if ( initVerbose() ) {
generator()->log()
<< "--- MatchboxFactory splitting channels ----------------------------------------------\n";
const StandardXComb& bxc = *xcptr;
generator()->log() << " hard process handled is: ";
generator()->log() << bxc.mePartonData()[0]->PDGName() << " "
<< bxc.mePartonData()[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = bxc.mePartonData().begin() + 2;
p != bxc.mePartonData().end(); ++p ) {
generator()->log() << (**p).PDGName() << " ";
}
generator()->log() << "\n";
for ( list<MatchboxFactory::SplittingChannel>::const_iterator sp =
res.begin(); sp != res.end(); ++sp ) {
sp->print(generator()->log());
}
generator()->log()
<< "--------------------------------------------------------\n"
<< flush;
}
return res;
}
void MatchboxFactory::print(ostream& os) const {
os << "--- MatchboxFactory setup -----------------------------------------------------------\n";
if ( !amplitudes().empty() ) {
os << " generated Born matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator m = bornMEs().begin();
m != bornMEs().end(); ++m ) {
(**m).print(os);
}
os << flush;
os << " generated real emission matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator m = realEmissionMEs().begin();
m != realEmissionMEs().end(); ++m ) {
(**m).print(os);
}
os << flush;
}
os << " generated Born+virtual matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator bv
= bornVirtualMEs().begin(); bv != bornVirtualMEs().end(); ++bv ) {
(**bv).print(os);
}
os << " generated subtracted matrix elements:\n";
for ( vector<Ptr<SubtractedME>::ptr>::const_iterator sub
= subtractedMEs().begin(); sub != subtractedMEs().end(); ++sub ) {
os << " '" << (**sub).name() << "'\n";
}
os << "---------------------------------------------------\n";
os << flush;
}
void MatchboxFactory::summary(ostream& os) const {
os << "\n\n================================================================================\n"
<< " Matchbox hard process summary\n"
<< "================================================================================\n\n";
os << " Electro-weak parameter summary:\n"
<< "---------------------------------------------------\n\n";
os << " Electro-weak scheme : ";
switch ( SM().ewScheme() ) {
case 0: os << "Default"; break;
case 1: os << "GMuScheme"; break;
case 2: os << "alphaMZScheme"; break;
case 3: os << "NoMass"; break;
case 4: os << "mW"; break;
case 5: os << "mZ"; break;
case 6: os << "Independent"; break;
case 7: os << "FeynRulesUFO"; break;
default: assert(false);
}
os << "\n";
os << " alphaEM is "
<< (SM().ewScheme() == 0 && !theFixedQEDCouplings ? "running" : "fixed at alphaEM(m(Z))") << "\n";
if ( SM().ewScheme() == 0 && !theFixedQEDCouplings )
os << " alphaEM is running at " << SM().alphaEMPtr()->nloops()
<< " loops\n\n";
else
os << "\n";
os << (SM().ewScheme() != 0 ? " Tree level relations " : " Best values ")
<< "yield:\n\n"
<< " m(Z)/GeV = "
<< getParticleData(ParticleID::Z0)->hardProcessMass()/GeV
<< "\n"
<< " g(Z)/GeV = "
<< getParticleData(ParticleID::Z0)->hardProcessWidth()/GeV
<< "\n"
<< " m(W)/GeV = "
<< getParticleData(ParticleID::Wplus)->hardProcessMass()/GeV
<< "\n"
<< " g(W)/GeV = "
<< getParticleData(ParticleID::Wplus)->hardProcessWidth()/GeV
<< "\n"
<< " m(H)/GeV = "
<< getParticleData(ParticleID::h0)->hardProcessMass()/GeV
<< "\n"
<< " g(H)/GeV = "
<< getParticleData(ParticleID::h0)->hardProcessWidth()/GeV
<< "\n"
<< " alphaEM(m(Z)) = "
<< SM().alphaEMME(sqr(getParticleData(ParticleID::Z0)->hardProcessMass())) << "\n"
<< " sin^2(theta) = " << SM().sin2ThetaW()
<< "\n"
<< " GeV^2 GF = " << GeV2*SM().fermiConstant()
<< "\n\n";
os << " Quark masses and widths are:\n"
<< "---------------------------------------------------\n\n"
<< " m(u)/GeV = " << getParticleData(ParticleID::u)->hardProcessMass()/GeV << "\n"
<< " m(d)/GeV = " << getParticleData(ParticleID::d)->hardProcessMass()/GeV << "\n"
<< " m(c)/GeV = " << getParticleData(ParticleID::c)->hardProcessMass()/GeV << "\n"
<< " m(s)/GeV = " << getParticleData(ParticleID::s)->hardProcessMass()/GeV << "\n"
<< " m(t)/GeV = " << getParticleData(ParticleID::t)->hardProcessMass()/GeV << "\n"
<< " g(t)/GeV = " << getParticleData(ParticleID::t)->hardProcessWidth()/GeV << "\n"
<< " m(b)/GeV = " << getParticleData(ParticleID::b)->hardProcessMass()/GeV << "\n\n";
os << " Lepton masses and widths are:\n"
<< "---------------------------------------------------\n\n"
<< " m(n_e)/GeV = " << getParticleData(ParticleID::nu_e)->hardProcessMass()/GeV << "\n"
<< " m(e)/GeV = " << getParticleData(ParticleID::eminus)->hardProcessMass()/GeV << "\n"
<< " m(n_mu)/GeV = " << getParticleData(ParticleID::nu_mu)->hardProcessMass()/GeV << "\n"
<< " m(mu)/GeV = " << getParticleData(ParticleID::muminus)->hardProcessMass()/GeV << "\n"
<< " m(nu_tau)/GeV = " << getParticleData(ParticleID::nu_tau)->hardProcessMass()/GeV << "\n"
<< " m(tau)/GeV = " << getParticleData(ParticleID::tauminus)->hardProcessMass()/GeV << "\n\n";
os << " Strong coupling summary:\n"
<< "---------------------------------------------------\n\n";
os << " alphaS is";
if ( !theFixedCouplings ) {
os << " running at " << SM().alphaSPtr()->nloops()
<< " loops with\n"
<< " alphaS(m(Z)) = " << SM().alphaSPtr()->value(sqr(getParticleData(ParticleID::Z0)->mass()))
<< "\n\n";
} else {
os << " fixed at "
<< SM().alphaS()
<< "\n\n";
}
if ( !theFixedCouplings ) {
os << " flavour thresholds are matched at\n";
for ( long id = 1; id <= 6; ++id ) {
os << " m(" << id << ")/GeV = "
<< (SM().alphaSPtr()->quarkMasses().empty() ?
getParticleData(id)->mass()/GeV :
SM().alphaSPtr()->quarkMasses()[id-1]/GeV)
<< "\n";
}
}
os << "\n\n" << flush;
}
void MatchboxFactory::doinit() {
theIsMatchboxRun() = true;
theCurrentFactory = Ptr<MatchboxFactory>::tptr(this);
if ( RunDirectories::empty() )
RunDirectories::pushRunId(generator()->runName());
setup();
if ( theShowerApproximation )
theShowerApproximation->init();
if ( initVerbose() && !ranSetup )
print(Repository::clog());
Ptr<StandardEventHandler>::tptr eh =
dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(generator()->eventHandler());
assert(eh);
if ( initVerbose() && !ranSetup ) {
assert(standardModel());
standardModel()->init();
summary(Repository::clog());
}
SubProcessHandler::doinit();
}
void MatchboxFactory::doinitrun() {
theIsMatchboxRun() = true;
theCurrentFactory = Ptr<MatchboxFactory>::tptr(this);
if ( theShowerApproximation )
theShowerApproximation->initrun();
Ptr<StandardEventHandler>::tptr eh =
dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(generator()->eventHandler());
assert(eh);
SubProcessHandler::doinitrun();
}
const string& MatchboxFactory::buildStorage() {
return RunDirectories::buildStorage();
}
const string& MatchboxFactory::runStorage() {
return RunDirectories::runStorage();
}
void MatchboxFactory::persistentOutput(PersistentOStream & os) const {
os << theDiagramGenerator << theProcessData
<< theNLight
<< theNLightJetVec << theNHeavyJetVec << theNLightProtonVec
<< theOrderInAlphaS << theOrderInAlphaEW
<< theBornContributions << theVirtualContributions
<< theRealContributions << theIndependentVirtuals << theIndependentPKs
<< thePhasespace << theScaleChoice
<< theFactorizationScaleFactor << theRenormalizationScaleFactor
<< theFixedCouplings << theFixedQEDCouplings << theVetoScales
<< theAmplitudes
<< theBornMEs << theVirtuals << theRealEmissionMEs << theLoopInducedMEs
<< theBornVirtualMEs << theSubtractedMEs << theFiniteRealMEs
<< theVerbose << theInitVerbose << theSubtractionData << theSubtractionPlotType
<< theSubtractionScatterPlot << thePoleData
<< theParticleGroups << processes << loopInducedProcesses << realEmissionProcesses
<< theShowerApproximation << theSplittingDipoles
<< theRealEmissionScales << theAllProcesses
<< theOLPProcesses << theExternalAmplitudes
<< theSelectedAmplitudes << theDeselectedAmplitudes
<< theDipoleSet << theReweighters << thePreweighters
<< theMECorrectionsOnly<< theLoopSimCorrections<<theHighestVirtualsize << ranSetup
<< theIncoming << theFirstPerturbativePDF << theSecondPerturbativePDF
<< inProductionMode << theSpinCorrelations << theAlphaParameter
<< theEnforceChargeConservation << theEnforceColourConservation
<< theEnforceLeptonNumberConservation << theEnforceQuarkNumberConservation
<< theLeptonFlavourDiagonal << theQuarkFlavourDiagonal;
}
void MatchboxFactory::persistentInput(PersistentIStream & is, int) {
is >> theDiagramGenerator >> theProcessData
>> theNLight
>> theNLightJetVec >> theNHeavyJetVec >> theNLightProtonVec
>> theOrderInAlphaS >> theOrderInAlphaEW
>> theBornContributions >> theVirtualContributions
>> theRealContributions >> theIndependentVirtuals >> theIndependentPKs
>> thePhasespace >> theScaleChoice
>> theFactorizationScaleFactor >> theRenormalizationScaleFactor
>> theFixedCouplings >> theFixedQEDCouplings >> theVetoScales
>> theAmplitudes
>> theBornMEs >> theVirtuals >> theRealEmissionMEs >> theLoopInducedMEs
>> theBornVirtualMEs >> theSubtractedMEs >> theFiniteRealMEs
>> theVerbose >> theInitVerbose >> theSubtractionData >> theSubtractionPlotType
>> theSubtractionScatterPlot >> thePoleData
>> theParticleGroups >> processes >> loopInducedProcesses >> realEmissionProcesses
>> theShowerApproximation >> theSplittingDipoles
>> theRealEmissionScales >> theAllProcesses
>> theOLPProcesses >> theExternalAmplitudes
>> theSelectedAmplitudes >> theDeselectedAmplitudes
>> theDipoleSet >> theReweighters >> thePreweighters
>> theMECorrectionsOnly>> theLoopSimCorrections>>theHighestVirtualsize >> ranSetup
>> theIncoming >> theFirstPerturbativePDF >> theSecondPerturbativePDF
>> inProductionMode >> theSpinCorrelations >> theAlphaParameter
>> theEnforceChargeConservation >> theEnforceColourConservation
>> theEnforceLeptonNumberConservation >> theEnforceQuarkNumberConservation
>> theLeptonFlavourDiagonal >> theQuarkFlavourDiagonal;
}
string MatchboxFactory::startParticleGroup(string name) {
particleGroupName = StringUtils::stripws(name);
particleGroup.clear();
return "";
}
string MatchboxFactory::endParticleGroup(string) {
if ( particleGroup.empty() )
throw Exception() << "MatchboxFactory: Empty particle group."
<< Exception::runerror;
particleGroups()[particleGroupName] = particleGroup;
particleGroup.clear();
return "";
}
vector<string> MatchboxFactory::parseProcess(string in) {
vector<string> process = StringUtils::split(in);
if ( process.size() < 3 )
throw Exception() << "MatchboxFactory: Invalid process."
<< Exception::runerror;
for ( vector<string>::iterator p = process.begin();
p != process.end(); ++p ) {
*p = StringUtils::stripws(*p);
}
vector<string> pprocess;
for ( vector<string>::const_iterator p = process.begin();
p != process.end(); ++p ) {
if ( *p == "->" )
continue;
pprocess.push_back(*p);
}
return pprocess;
}
string MatchboxFactory::doProcess(string in) {
processes.push_back(parseProcess(in));
return "";
}
string MatchboxFactory::doLoopInducedProcess(string in) {
loopInducedProcesses.push_back(parseProcess(in));
return "";
}
string MatchboxFactory::doSingleRealProcess(string in) {
realEmissionProcesses.push_back(parseProcess(in));
return "";
}
struct SortPID {
inline bool operator()(PDPtr a, PDPtr b) const {
return a->id() < b->id();
}
};
//
// @TODO
//
// SP: After improving this for standard model process building this should
// actually got into a separate process builder class or something along these
// lines to have it better factored for use with BSM models.
//
//
set<PDVector> MatchboxFactory::
makeSubProcesses(const vector<string>& proc) const {
if ( proc.empty() )
throw Exception() << "MatchboxFactory: No process specified."
<< Exception::runerror;
vector<PDVector> groups;
typedef map<string,PDVector>::const_iterator GroupIterator;
for ( vector<string>::const_iterator gr = proc.begin();
gr != proc.end(); ++gr ) {
GroupIterator git = particleGroups().find(*gr);
if ( git == particleGroups().end() ) {
throw Exception() << "MatchboxFactory: Particle group '"
<< *gr << "' not defined." << Exception::runerror;
}
groups.push_back(git->second);
}
vector<size_t> counts(groups.size(),0);
PDVector proto(groups.size());
set<PDVector> allProcs;
while ( true ) {
for ( size_t k = 0; k < groups.size(); ++k )
proto[k] = groups[k][counts[k]];
int charge = 0;
int colour = 0;
int nleptons = 0;
int nquarks = 0;
int ncolour = 0;
int nleptonsGen[4];
int nquarksGen[4];
for ( size_t i = 0; i < 4; ++i ) {
nleptonsGen[i] = 0;
nquarksGen[i] = 0;
}
for ( size_t k = 0; k < proto.size(); ++k ) {
int sign = k > 1 ? 1 : -1;
charge += sign * proto[k]->iCharge();
colour += sign * proto[k]->iColour();
if ( abs(proto[k]->id()) <= 8 ) {
int generation = (abs(proto[k]->id()) - 1)/2;
nquarks += sign * ( proto[k]->id() < 0 ? -1 : 1);
nquarksGen[generation] += sign * ( proto[k]->id() < 0 ? -1 : 1);
}
if ( abs(proto[k]->id()) > 10 &&
abs(proto[k]->id()) <= 18 ) {
int generation = (abs(proto[k]->id()) - 11)/2;
nleptons += sign * ( proto[k]->id() < 0 ? -1 : 1);
nleptonsGen[generation] += sign * ( proto[k]->id() < 0 ? -1 : 1);
}
if ( proto[k]->coloured() )
++ncolour;
}
bool pass = true;
if ( theEnforceChargeConservation )
pass &= (charge == 0);
if ( theEnforceColourConservation )
pass &= (colour % 8 == 0);
if ( theEnforceLeptonNumberConservation ) {
pass &= (nleptons == 0);
if ( theLeptonFlavourDiagonal ) {
for ( size_t i = 0; i < 4; ++i )
pass &= (nleptonsGen[i] == 0);
}
}
if ( theEnforceQuarkNumberConservation ) {
pass &= (nquarks == 0);
if ( theQuarkFlavourDiagonal ) {
for ( size_t i = 0; i < 4; ++i )
pass &= (nquarksGen[i] == 0);
}
}
if ( pass ) {
for ( int i = 0; i < 2; ++i ) {
if ( proto[i]->coloured() &&
proto[i]->hardProcessMass() != ZERO )
throw Exception()
<< "Inconsistent flavour scheme detected with massive incoming "
<< proto[i]->PDGName() << ". Check your setup."
<< Exception::runerror;
}
sort(proto.begin()+2,proto.end(),SortPID());
allProcs.insert(proto);
}
vector<size_t>::reverse_iterator c = counts.rbegin();
vector<PDVector>::const_reverse_iterator g = groups.rbegin();
while ( c != counts.rend() ) {
if ( ++(*c) == g->size() ) {
*c = 0;
++c; ++g;
} else {
break;
}
}
if ( c == counts.rend() )
break;
}
return allProcs;
}
void MatchboxFactory::Init() {
static ClassDocumentation<MatchboxFactory> documentation
("MatchboxFactory",
"NLO QCD corrections have been calculated "
"using Matchbox \\cite{Platzer:2011bc}, \\cite{Matchbox:2015}",
"%\\cite{Platzer:2011bc}\n"
"\\bibitem{Platzer:2011bc}\n"
"S.~Platzer and S.~Gieseke,\n"
"``Dipole Showers and Automated NLO Matching in Herwig,''\n"
"arXiv:1109.6256 [hep-ph].\n"
"%%CITATION = ARXIV:1109.6256;%%\n"
"%\\cite{Matchbox:2015}\n"
"\\bibitem{Matchbox:2015}\n"
"Herwig collaboration,\n"
"``Precision LHC Event Generation with Herwig,''\n"
"in preparation.");
static Reference<MatchboxFactory,Tree2toNGenerator> interfaceDiagramGenerator
("DiagramGenerator",
"Set the diagram generator.",
&MatchboxFactory::theDiagramGenerator, false, false, true, true, false);
interfaceDiagramGenerator.rank(-1);
static Reference<MatchboxFactory,ProcessData> interfaceProcessData
("ProcessData",
"Set the process data object to be used.",
&MatchboxFactory::theProcessData, false, false, true, true, false);
interfaceProcessData.rank(-1);
static Parameter<MatchboxFactory,unsigned int> interfaceOrderInAlphaS
("OrderInAlphaS",
"The order in alpha_s to consider.",
&MatchboxFactory::theOrderInAlphaS, 0, 0, 0,
false, false, Interface::lowerlim);
static Parameter<MatchboxFactory,unsigned int> interfaceOrderInAlphaEW
("OrderInAlphaEW",
"The order in alpha_EW to consider.",
&MatchboxFactory::theOrderInAlphaEW, 2, 0, 0,
false, false, Interface::lowerlim);
static Switch<MatchboxFactory,bool> interfaceBornContributions
("BornContributions",
"Switch on or off the Born contributions.",
&MatchboxFactory::theBornContributions, true, false, false);
static SwitchOption interfaceBornContributionsYes
(interfaceBornContributions,
"Yes",
"Switch on Born contributions.",
true);
static SwitchOption interfaceBornContributionsNo
(interfaceBornContributions,
"No",
"Switch off Born contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceVirtualContributions
("VirtualContributions",
"Switch on or off the virtual contributions.",
&MatchboxFactory::theVirtualContributions, true, false, false);
static SwitchOption interfaceVirtualContributionsYes
(interfaceVirtualContributions,
"Yes",
"Switch on virtual contributions.",
true);
static SwitchOption interfaceVirtualContributionsNo
(interfaceVirtualContributions,
"No",
"Switch off virtual contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceRealContributions
("RealContributions",
"Switch on or off the real contributions.",
&MatchboxFactory::theRealContributions, true, false, false);
static SwitchOption interfaceRealContributionsYes
(interfaceRealContributions,
"Yes",
"Switch on real contributions.",
true);
static SwitchOption interfaceRealContributionsNo
(interfaceRealContributions,
"No",
"Switch off real contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceIndependentVirtuals
("IndependentVirtuals",
"Switch on or off virtual contributions as separate subprocesses.",
&MatchboxFactory::theIndependentVirtuals, true, false, false);
static SwitchOption interfaceIndependentVirtualsYes
(interfaceIndependentVirtuals,
"Yes",
"Switch on virtual contributions as separate subprocesses.",
true);
static SwitchOption interfaceIndependentVirtualsNo
(interfaceIndependentVirtuals,
"No",
"Switch off virtual contributions as separate subprocesses.",
false);
static Switch<MatchboxFactory,bool> interfaceIndependentPKs
("IndependentPKOperators",
"Switch on or off PK oeprators as separate subprocesses.",
&MatchboxFactory::theIndependentPKs, true, false, false);
static SwitchOption interfaceIndependentPKsYes
(interfaceIndependentPKs,
"Yes",
"Switch on PK operators as separate subprocesses.",
true);
static SwitchOption interfaceIndependentPKsNo
(interfaceIndependentPKs,
"No",
"Switch off PK operators as separate subprocesses.",
false);
static Reference<MatchboxFactory,MatchboxPhasespace> interfacePhasespace
("Phasespace",
"Set the phasespace generator.",
&MatchboxFactory::thePhasespace, false, false, true, true, false);
static Reference<MatchboxFactory,MatchboxScaleChoice> interfaceScaleChoice
("ScaleChoice",
"Set the scale choice object.",
&MatchboxFactory::theScaleChoice, false, false, true, true, false);
static Parameter<MatchboxFactory,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&MatchboxFactory::theFactorizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<MatchboxFactory,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&MatchboxFactory::theRenormalizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Switch<MatchboxFactory,bool> interfaceFixedCouplings
("FixedCouplings",
"Switch on or off fixed couplings.",
&MatchboxFactory::theFixedCouplings, true, false, false);
static SwitchOption interfaceFixedCouplingsYes
(interfaceFixedCouplings,
"Yes",
"Yes",
true);
static SwitchOption interfaceFixedCouplingsNo
(interfaceFixedCouplings,
"No",
"No",
false);
interfaceFixedCouplings.rank(-1);
static Switch<MatchboxFactory,bool> interfaceFixedQEDCouplings
("FixedQEDCouplings",
"Switch on or off fixed QED couplings.",
&MatchboxFactory::theFixedQEDCouplings, true, false, false);
static SwitchOption interfaceFixedQEDCouplingsYes
(interfaceFixedQEDCouplings,
"Yes",
"Yes",
true);
static SwitchOption interfaceFixedQEDCouplingsNo
(interfaceFixedQEDCouplings,
"No",
"No",
false);
interfaceFixedQEDCouplings.rank(-1);
// @TDOO SP to remove this in the code as well
/*
static Switch<MatchboxFactory,bool> interfaceVetoScales
("VetoScales",
"Switch on or setting veto scales.",
&MatchboxFactory::theVetoScales, false, false, false);
static SwitchOption interfaceVetoScalesYes
(interfaceVetoScales,
"Yes",
"Yes",
true);
static SwitchOption interfaceVetoScalesNo
(interfaceVetoScales,
"No",
"No",
false);
*/
static RefVector<MatchboxFactory,MatchboxAmplitude> interfaceAmplitudes
("Amplitudes",
"The amplitude objects.",
&MatchboxFactory::theAmplitudes, -1, false, false, true, true, false);
static Switch<MatchboxFactory,bool> interfaceVerbose
("Verbose",
"Print full infomation on each evaluated phase space point.",
&MatchboxFactory::theVerbose, false, false, false);
static SwitchOption interfaceVerboseYes
(interfaceVerbose,
"Yes",
"Yes",
true);
static SwitchOption interfaceVerboseNo
(interfaceVerbose,
"No",
"No",
false);
interfaceVerbose.rank(-1);
static Switch<MatchboxFactory,bool> interfaceInitVerbose
("InitVerbose",
"Print setup information.",
&MatchboxFactory::theInitVerbose, false, false, false);
static SwitchOption interfaceInitVerboseYes
(interfaceInitVerbose,
"Yes",
"Yes",
true);
static SwitchOption interfaceInitVerboseNo
(interfaceInitVerbose,
"No",
"No",
false);
interfaceInitVerbose.rank(-1);
static Parameter<MatchboxFactory,string> interfaceSubtractionData
("SubtractionData",
"Prefix for subtraction check data.",
&MatchboxFactory::theSubtractionData, "",
false, false);
static Switch<MatchboxFactory,int> interfaceSubtractionPlotType
("SubtractionPlotType",
"Switch for controlling what kind of plot is generated for checking the subtraction",
&MatchboxFactory::theSubtractionPlotType, 1, false, false);
static SwitchOption interfaceSubtractionPlotTypeLinearRatio
(interfaceSubtractionPlotType,
"LinRatio",
"Switch on the linear plot of the ratio",
1);
static SwitchOption interfaceSubtractionPlotTypeLogRelDiff
(interfaceSubtractionPlotType,
"LogRelDiff",
"Switch on the logarithmic plot of the relative difference",
2);
static Switch<MatchboxFactory,bool> interfaceSubtractionScatterPlot
("SubtractionScatterPlot",
"Switch for controlling whether subtraction data should be plotted for each phase space point individually",
&MatchboxFactory::theSubtractionScatterPlot, false, false, false);
static SwitchOption interfaceSubtractionScatterPlotNo
(interfaceSubtractionScatterPlot,
"No", "Switch off the scatter plot", false);
static SwitchOption interfaceSubtractionScatterPlotYes
(interfaceSubtractionScatterPlot,
"Yes", "Switch on the scatter plot", true);
static Parameter<MatchboxFactory,string> interfacePoleData
("PoleData",
"Prefix for subtraction check data.",
&MatchboxFactory::thePoleData, "",
false, false);
static RefVector<MatchboxFactory,ParticleData> interfaceParticleGroup
("ParticleGroup",
"The particle group just started.",
&MatchboxFactory::particleGroup, -1, false, false, true, false, false);
static Command<MatchboxFactory> interfaceStartParticleGroup
("StartParticleGroup",
"Start a particle group.",
&MatchboxFactory::startParticleGroup, false);
static Command<MatchboxFactory> interfaceEndParticleGroup
("EndParticleGroup",
"End a particle group.",
&MatchboxFactory::endParticleGroup, false);
static Command<MatchboxFactory> interfaceProcess
("Process",
"Set the process(es) to consider.",
&MatchboxFactory::doProcess, false);
static Command<MatchboxFactory> interfaceLoopInducedProcess
("LoopInducedProcess",
"Set the loop induced process(es) to consider.",
&MatchboxFactory::doLoopInducedProcess, false);
static Command<MatchboxFactory> interfaceSingleRealProcess
("SingleRealProcess",
"Set the real emission process(es) to consider.",
&MatchboxFactory::doSingleRealProcess, false);
static Reference<MatchboxFactory,ShowerApproximation> interfaceShowerApproximation
("ShowerApproximation",
"Set the shower approximation to be considered.",
&MatchboxFactory::theShowerApproximation, false, false, true, true, false);
static Switch<MatchboxFactory,bool> interfaceRealEmissionScales
("RealEmissionScales",
"Switch on or off calculation of subtraction scales from real emission kinematics.",
&MatchboxFactory::theRealEmissionScales, false, false, false);
static SwitchOption interfaceRealEmissionScalesYes
(interfaceRealEmissionScales,
"Yes",
"Yes",
true);
static SwitchOption interfaceRealEmissionScalesNo
(interfaceRealEmissionScales,
"No",
"No",
false);
interfaceRealEmissionScales.rank(-1);
static Switch<MatchboxFactory,bool> interfaceAllProcesses
("AllProcesses",
"Consider all processes up to a maximum coupling order specified by the coupling order interfaces.",
&MatchboxFactory::theAllProcesses, false, false, false);
static SwitchOption interfaceAllProcessesYes
(interfaceAllProcesses,
"Yes",
"Include all processes.",
true);
static SwitchOption interfaceAllProcessesNo
(interfaceAllProcesses,
"No",
"Only consider processes matching the exact order in the couplings.",
false);
interfaceAllProcesses.rank(-1);
static RefVector<MatchboxFactory,MatchboxAmplitude> interfaceSelectAmplitudes
("SelectAmplitudes",
"The amplitude objects to be favoured in clashing responsibilities.",
&MatchboxFactory::theSelectedAmplitudes, -1, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxAmplitude> interfaceDeselectAmplitudes
("DeselectAmplitudes",
"The amplitude objects to be disfavoured in clashing responsibilities.",
&MatchboxFactory::theDeselectedAmplitudes, -1, false, false, true, true, false);
static Switch<MatchboxFactory,int> interfaceDipoleSet
("DipoleSet",
"The set of subtraction terms to be considered.",
&MatchboxFactory::theDipoleSet, 0, false, false);
static SwitchOption interfaceDipoleSetCataniSeymour
(interfaceDipoleSet,
"CataniSeymour",
"Use default Catani-Seymour dipoles.",
0);
interfaceDipoleSet.rank(-1);
static RefVector<MatchboxFactory,ReweightBase> interfaceReweighters
("Reweighters",
"Reweight objects for matrix elements.",
&MatchboxFactory::theReweighters, -1, false, false, true, false, false);
static RefVector<MatchboxFactory,ReweightBase> interfacePreweighters
("Preweighters",
"Preweight objects for matrix elements.",
&MatchboxFactory::thePreweighters, -1, false, false, true, false, false);
static Switch<MatchboxFactory,bool> interfaceMECorrectionsOnly
("MECorrectionsOnly",
"Prepare only ME corrections, but no NLO calculation.",
&MatchboxFactory::theMECorrectionsOnly, false, false, false);
static SwitchOption interfaceMECorrectionsOnlyYes
(interfaceMECorrectionsOnly,
"Yes",
"Produce only ME corrections.",
true);
static SwitchOption interfaceMECorrectionsOnlyNo
(interfaceMECorrectionsOnly,
"No",
"Produce full NLO.",
false);
static Switch<MatchboxFactory,bool> interfaceLoopSimCorrections
("LoopSimCorrections",
"Prepare LoopSim corrections.",
&MatchboxFactory::theLoopSimCorrections, false, false, false);
static SwitchOption interfaceLoopSimCorrectionsYes
(interfaceLoopSimCorrections,
"Yes",
"Produce loopsim corrections.",
true);
static SwitchOption interfaceLoopSimCorrectionsNo
(interfaceLoopSimCorrections,
"No",
"Produce full NLO.",
false);
static Switch<MatchboxFactory,bool> interfaceFirstPerturbativePDF
("FirstPerturbativePDF",
"",
&MatchboxFactory::theFirstPerturbativePDF, true, false, false);
static SwitchOption interfaceFirstPerturbativePDFYes
(interfaceFirstPerturbativePDF,
"Yes",
"",
true);
static SwitchOption interfaceFirstPerturbativePDFNo
(interfaceFirstPerturbativePDF,
"No",
"",
false);
interfaceFirstPerturbativePDF.rank(-1);
static Switch<MatchboxFactory,bool> interfaceSecondPerturbativePDF
("SecondPerturbativePDF",
"",
&MatchboxFactory::theSecondPerturbativePDF, true, false, false);
static SwitchOption interfaceSecondPerturbativePDFYes
(interfaceSecondPerturbativePDF,
"Yes",
"",
true);
static SwitchOption interfaceSecondPerturbativePDFNo
(interfaceSecondPerturbativePDF,
"No",
"",
false);
interfaceSecondPerturbativePDF.rank(-1);
static Command<MatchboxFactory> interfaceProductionMode
("ProductionMode",
"Switch this factory to production mode.",
&MatchboxFactory::doProductionMode, false);
static Switch<MatchboxFactory,bool> interfaceSpinCorrelations
("SpinCorrelations",
"Fill information for the spin correlations, if possible.",
&MatchboxFactory::theSpinCorrelations, false, false, false);
static SwitchOption interfaceSpinCorrelationsYes
(interfaceSpinCorrelations,
"Yes",
"",
true);
static SwitchOption interfaceSpinCorrelationsNo
(interfaceSpinCorrelations,
"No",
"",
false);
static Parameter<MatchboxFactory,double> interfaceAlphaParameter
("AlphaParameter",
"Nagy-AlphaParameter.",
&MatchboxFactory::theAlphaParameter, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Switch<MatchboxFactory,bool> interfaceEnforceChargeConservation
("EnforceChargeConservation",
"Enforce charge conservation while generating the hard process.",
&MatchboxFactory::theEnforceChargeConservation, true, false, false);
static SwitchOption interfaceEnforceChargeConservationYes
(interfaceEnforceChargeConservation,
"Yes",
"Enforce charge conservation.",
true);
static SwitchOption interfaceEnforceChargeConservationNo
(interfaceEnforceChargeConservation,
"No",
"Do not enforce charge conservation.",
false);
interfaceEnforceChargeConservation.rank(-1);
static Switch<MatchboxFactory,bool> interfaceEnforceColourConservation
("EnforceColourConservation",
"Enforce colour conservation while generating the hard process.",
&MatchboxFactory::theEnforceColourConservation, false, false, false);
static SwitchOption interfaceEnforceColourConservationYes
(interfaceEnforceColourConservation,
"Yes",
"Enforce colour conservation.",
true);
static SwitchOption interfaceEnforceColourConservationNo
(interfaceEnforceColourConservation,
"No",
"Do not enforce colour conservation.",
false);
interfaceEnforceColourConservation.rank(-1);
static Switch<MatchboxFactory,bool> interfaceEnforceLeptonNumberConservation
("EnforceLeptonNumberConservation",
"Enforce lepton number conservation while generating the hard process.",
&MatchboxFactory::theEnforceLeptonNumberConservation, false, false, false);
static SwitchOption interfaceEnforceLeptonNumberConservationYes
(interfaceEnforceLeptonNumberConservation,
"Yes",
"Enforce lepton number conservation.",
true);
static SwitchOption interfaceEnforceLeptonNumberConservationNo
(interfaceEnforceLeptonNumberConservation,
"No",
"Do not enforce lepton number conservation.",
false);
interfaceEnforceLeptonNumberConservation.rank(-1);
static Switch<MatchboxFactory,bool> interfaceEnforceQuarkNumberConservation
("EnforceQuarkNumberConservation",
"Enforce quark number conservation while generating the hard process.",
&MatchboxFactory::theEnforceQuarkNumberConservation, false, false, false);
static SwitchOption interfaceEnforceQuarkNumberConservationYes
(interfaceEnforceQuarkNumberConservation,
"Yes",
"Enforce quark number conservation.",
true);
static SwitchOption interfaceEnforceQuarkNumberConservationNo
(interfaceEnforceQuarkNumberConservation,
"No",
"Do not enforce quark number conservation.",
false);
interfaceEnforceQuarkNumberConservation.rank(-1);
static Switch<MatchboxFactory,bool> interfaceLeptonFlavourDiagonal
("LeptonFlavourDiagonal",
"Assume that lepton interactions are flavour diagonal while generating the hard process.",
&MatchboxFactory::theLeptonFlavourDiagonal, false, false, false);
static SwitchOption interfaceLeptonFlavourDiagonalYes
(interfaceLeptonFlavourDiagonal,
"Yes",
"Assume that lepton interactions are flavour diagonal.",
true);
static SwitchOption interfaceLeptonFlavourDiagonalNo
(interfaceLeptonFlavourDiagonal,
"No",
"Do not assume that lepton interactions are flavour diagonal.",
false);
interfaceLeptonFlavourDiagonal.rank(-1);
static Switch<MatchboxFactory,bool> interfaceQuarkFlavourDiagonal
("QuarkFlavourDiagonal",
"Assume that quark interactions are flavour diagonal while generating the hard process.",
&MatchboxFactory::theQuarkFlavourDiagonal, false, false, false);
static SwitchOption interfaceQuarkFlavourDiagonalYes
(interfaceQuarkFlavourDiagonal,
"Yes",
"Assume that quark interactions are flavour diagonal.",
true);
static SwitchOption interfaceQuarkFlavourDiagonalNo
(interfaceQuarkFlavourDiagonal,
"No",
"Do not assume that quark interactions are flavour diagonal.",
false);
interfaceQuarkFlavourDiagonal.rank(-1);
}
// *** 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<MatchboxFactory,SubProcessHandler>
describeHerwigMatchboxFactory("Herwig::MatchboxFactory", "Herwig.so");
diff --git a/Sampling/BinSampler.cc b/Sampling/BinSampler.cc
--- a/Sampling/BinSampler.cc
+++ b/Sampling/BinSampler.cc
@@ -1,730 +1,730 @@
// -*- C++ -*-
//
// BinSampler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the BinSampler class.
//
#include "BinSampler.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/Repository/Repository.h"
#include "ThePEG/Utilities/ColourOutput.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Handlers/StandardEventHandler.h"
#include "ThePEG/Handlers/StandardXComb.h"
-#include <boost/progress.hpp>
+#include "Herwig/Utilities/Progress.h"
#include "GeneralSampler.h"
using namespace Herwig;
BinSampler::BinSampler()
: MultiIterationStatistics(),
theBias(1.),
theWeighted(false),
theInitialPoints(1000000),
theNIterations(1),
theEnhancementFactor(1.0),
theNonZeroInPresampling(false),
theHalfPoints(false),
theMaxNewMax(30),
theReferenceWeight(1.0),
theBin(-1),
theInitialized(false),
theRemapperPoints(0),
theRemapChannelDimension(false),
theLuminosityMapperBins(0),
theGeneralMapperBins(0),
theRemapperMinSelection(0.00001),
theIntegrated(false),
theRemappersFilled(false),
theHasGrids(false),
theKappa(1.){}
BinSampler::~BinSampler() {}
IBPtr BinSampler::clone() const {
return new_ptr(*this);
}
IBPtr BinSampler::fullclone() const {
return new_ptr(*this);
}
void BinSampler::sampler(Ptr<GeneralSampler>::tptr s) {
theSampler = s;
}
Ptr<GeneralSampler>::tptr BinSampler::sampler() const {
return theSampler;
}
string BinSampler::process() const {
ostringstream os("");
const StandardEventHandler& eh = *theEventHandler;
const StandardXComb& xc = *eh.xCombs()[theBin];
os << xc.matrixElement()->name() << " : ";
os << xc.mePartonData()[0]->PDGName() << " "
<< xc.mePartonData()[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator pid =
xc.mePartonData().begin() + 2;
pid != xc.mePartonData().end(); ++pid )
os << (**pid).PDGName() << " ";
return os.str();
}
string BinSampler::shortprocess() const {
ostringstream os("");
const StandardEventHandler& eh = *theEventHandler;
const StandardXComb& xc = *eh.xCombs()[theBin];
os << xc.mePartonData()[0]->id() << " "
<< xc.mePartonData()[1]->id() << " : ";
for ( cPDVector::const_iterator pid =
xc.mePartonData().begin() + 2;
pid != xc.mePartonData().end(); ++pid )
os << (**pid).id() << " ";
return os.str();
}
string BinSampler::id() const {
ostringstream os("");
const StandardEventHandler& eh = *theEventHandler;
const StandardXComb& xc = *eh.xCombs()[theBin];
string name = xc.matrixElement()->name();
string::size_type i = name.find_first_of("[");
string nameFirst = name.substr(0,i);
i = name.find_first_of("]");
string nameSecond = name.substr(i+1);
os << nameFirst << nameSecond << ":";
for ( cPDVector::const_iterator pid =
xc.mePartonData().begin();
pid != xc.mePartonData().end(); ++pid )
os << (**pid).id() << (pid != (--xc.mePartonData().end()) ? "," : "");
return os.str();
}
double BinSampler::evaluate(vector<double> p,
bool remap) {
double w = 1.0;
if ( remap && !remappers.empty() ) {
for ( size_t k = 0; k < p.size(); ++k ) {
map<size_t,Remapper>::const_iterator r =
remappers.find(k);
if ( r != remappers.end() ) {
pair<double,double> f = r->second.generate(p[k]);
p[k] = f.first;
w /= f.second;
}
}
}
try {
w *= eventHandler()->dSigDR(p) / nanobarn;
} catch (Veto&) {
w = 0.0;
} catch (...) {
throw;
}
if (randomNumberString()!="")
for ( size_t k = 0; k < p.size(); ++k ) {
RandomNumberHistograms[RandomNumberIndex(id(),k)].first.book(p[k],w);
RandomNumberHistograms[RandomNumberIndex(id(),k)].second+=w;
}
return w;
}
double BinSampler::generate() {
double w = 1.;
for ( size_t k = 0; k < lastPoint().size(); ++k ) {
lastPoint()[k] = UseRandom::rnd();
}
try {
w = evaluate(lastPoint());
} catch (Veto&) {
w = 0.0;
} catch (...) {
throw;
}
if ( !weighted() && initialized() ) {
double p = min(abs(w),kappa()*referenceWeight())/(kappa()*referenceWeight());
double sign = w >= 0. ? 1. : -1.;
if ( p < 1 && UseRandom::rnd() > p )
w = 0.;
else
w = sign*max(abs(w),referenceWeight()*kappa());
}
select(w);
if ( w != 0.0 )
accept();
assert(kappa()==1.||sampler()->almostUnweighted());
return w;
}
void BinSampler::fillRemappers(bool progress) {
if ( remappers.empty() )
return;
unsigned long nanPoints = 0;
- boost::progress_display* progressBar = 0;
+ progress_display* progressBar = nullptr;
if ( progress ) {
Repository::clog() << "warming up " << ANSI::red << process() << ANSI::reset;
- progressBar = new boost::progress_display(theRemapperPoints,Repository::clog());
+ progressBar = new progress_display{ theRemapperPoints, Repository::clog() };
}
unsigned long countzero =0;
for ( unsigned long k = 0; k < theRemapperPoints; ++k,++countzero ) {
if (countzero>=theRemapperPoints)break;
double w = 1.;
for ( size_t j = 0; j < lastPoint().size(); ++j ) {
lastPoint()[j] = UseRandom::rnd();
}
try {
w = evaluate(lastPoint(),false);
} catch (Veto&) {
w = 0.0;
} catch (...) {
throw;
}
if ( ! isfinite(w) )
++nanPoints;
if ( theNonZeroInPresampling && w==0. ){
k--;
continue;
}
if ( w != 0.0 ) {
countzero=0;
for ( map<size_t,Remapper>::iterator r = remappers.begin();
r != remappers.end(); ++r )
r->second.fill(lastPoint()[r->first],w);
}
if ( progressBar )
++(*progressBar);
}
if ( progressBar ) {
delete progressBar;
}
if ( nanPoints ) {
Repository::clog() << "Warning: " << nanPoints
<< " out of " << theRemapperPoints << " points with nan or inf "
<< "weight encountered while filling remappers.\n" << flush;
}
}
void BinSampler::saveIntegrationData() const {
XML::Element stats = MultiIterationStatistics::toXML();
stats.appendAttribute("process",id());
sampler()->grids().append(stats);
}
void BinSampler::readIntegrationData() {
if ( theIntegrated )
return;
bool haveStats = false;
list<XML::Element>::iterator sit = sampler()->grids().children().begin();
for ( ; sit != sampler()->grids().children().end(); ++sit ) {
if ( sit->type() != XML::ElementTypes::Element )
continue;
if ( sit->name() != "MultiIterationStatistics" )
continue;
string proc;
sit->getFromAttribute("process",proc);
if ( proc == id() ) {
haveStats = true;
break;
}
}
if ( haveStats ) {
MultiIterationStatistics::fromXML(*sit);
sampler()->grids().erase(sit);
theIntegrated = true;
} else {
throw Exception()
<< "\n---------------------------------------------------\n\n"
<< "Expected integration data.\n\n"
<< "* When using the build setup make sure the integrate command has been run.\n\n"
<< "* Check the [EventGenerator].log file for further information.\n\n"
<< "* Make sure that the Herwig folder can be found and that it contains a HerwigGrids.xml file.\n\n"
<< "* If you have split the integration jobs, make sure that each integration job was finished.\n"
<< " Afterwards delete the global HerwigGrids.xml file in the Herwig subfolder\n"
<< " to automatically create an updated version of the global HerwigGrids.xml file.\n\n"
<< "---------------------------------------------------\n"
<< Exception::abortnow;
}
}
void BinSampler::saveRemappers() const {
if ( remappers.empty() )
return;
XML::Element maps(XML::ElementTypes::Element,"Remappers");
maps.appendAttribute("process",id());
for ( map<size_t,Remapper>::const_iterator r = remappers.begin();
r != remappers.end(); ++r ) {
XML::Element rmap = r->second.toXML();
rmap.appendAttribute("dimension",r->first);
maps.append(rmap);
}
sampler()->grids().append(maps);
}
void BinSampler::setupRemappers(bool progress) {
if ( !theRemapperPoints )
return;
if ( theRemappersFilled )
return;
lastPoint().resize(dimension());
bool haveGrid = false;
list<XML::Element>::iterator git = sampler()->grids().children().begin();
for ( ; git != sampler()->grids().children().end(); ++git ) {
if ( git->type() != XML::ElementTypes::Element )
continue;
if ( git->name() != "Remappers" )
continue;
string proc;
git->getFromAttribute("process",proc);
if ( proc == id() ) {
haveGrid = true;
break;
}
}
if ( haveGrid ) {
for ( list<XML::Element>::iterator cit = git->children().begin();
cit != git->children().end(); ++cit ) {
if ( cit->type() != XML::ElementTypes::Element )
continue;
if ( cit->name() != "Remapper" )
continue;
size_t dimension = 0;
cit->getFromAttribute("dimension",dimension);
remappers[dimension].fromXML(*cit);
}
sampler()->grids().erase(git);
}
if ( !haveGrid ) {
const StandardEventHandler& eh = *eventHandler();
const StandardXComb& xc = *eh.xCombs()[bin()];
const pair<int,int>& pdims = xc.partonDimensions();
set<int> remapped;
if ( theRemapChannelDimension && xc.diagrams().size() > 1 &&
dimension() > pdims.first + pdims.second ) {
remappers[pdims.first] = Remapper(xc.diagrams().size(),theRemapperMinSelection,false);
remapped.insert(pdims.first);
}
if ( theLuminosityMapperBins > 1 && dimension() >= pdims.first + pdims.second ) {
for ( int n = 0; n < pdims.first; ++n ) {
remappers[n] = Remapper(theLuminosityMapperBins,theRemapperMinSelection,true);
remapped.insert(n);
}
for ( int n = dimension() - pdims.second; n < dimension(); ++n ) {
remappers[n] = Remapper(theLuminosityMapperBins,theRemapperMinSelection,true);
remapped.insert(n);
}
}
if ( theGeneralMapperBins > 1 ) {
for ( int n = 0; n < dimension(); n++ ) {
if ( remapped.find(n) == remapped.end() ) {
remappers[n] = Remapper(theGeneralMapperBins,theRemapperMinSelection,true);
remapped.insert(n);
}
}
}
fillRemappers(progress);
for ( map<size_t,Remapper>::iterator r = remappers.begin();
r != remappers.end(); ++r ) {
r->second.finalize();
}
}
theRemappersFilled = true;
}
void BinSampler::runIteration(unsigned long points, bool progress) {
- boost::progress_display* progressBar = 0;
+ progress_display* progressBar = 0;
if ( progress ) {
Repository::clog() << "integrating " << ANSI::red << process()<< ANSI::reset << ", iteration "
<< (iterations().size() + 1);
- progressBar = new boost::progress_display(points,Repository::clog());
+ progressBar = new progress_display(points,Repository::clog());
}
double w=0.;
double maxweight=0;
int numlastmax=0;
unsigned long countzero =0;
int newmax=0;
for ( unsigned long k = 0; k < points; ++k,++countzero ) {
if (countzero>=points)break;
w=abs(generate());
if(theNonZeroInPresampling && w==0.0){
k--;
continue;
}
if (w!=0.0)
countzero =0;
numlastmax++;
if (theHalfPoints&&maxweight<w&&
numlastmax<(int)(points/2.)){
if(++newmax>theMaxNewMax){
throw Exception()
<< "\n---------------------------------------------------\n\n"
<< "To many new Maxima.\n\n"
<< "* With the option:\n\n"
<< "* set Sampler:BinSampler:HalfPoints Yes\n\n"
<< "* for every new maximum weight found until the half of the persampling points\n"
<< "* the counter is set to zero. We count the number of new maxima.\n"
<< "* You have reached: "<<newmax<<"\n"
<< "* Did you apply reasonable cuts to the process?\n"
<< "* You can set the maximum allowed new maxima by:"
<< "* set Sampler:BinSampler:MaxNewMax N\n\n"
<< "---------------------------------------------------\n"
<< Exception::abortnow;
}
maxweight=w;
k=0;
numlastmax=0;
}
if ( progress ) {
++(*progressBar);
}
}
if ( progress ) {
Repository::clog() << "integrated ( " << ANSI::yellow
<< averageWeight() << " +/- " << sqrt(averageWeightVariance())
<< ANSI::reset << " ) nb\nepsilon = "
<< (abs(maxWeight()) != 0. ? averageAbsWeight()/abs(maxWeight()) : 0.);
if ( !iterations().empty() )
Repository::clog() << " chi2 = " << chi2();
Repository::clog() << "\n";
Repository::clog() << "---------------------------------------------------\n";
}
if ( progressBar )
delete progressBar;
}
void BinSampler::initialize(bool progress) {
lastPoint().resize(dimension());
if (randomNumberString()!="")
for(size_t i=0;i<lastPoint().size();i++){
RandomNumberHistograms[RandomNumberIndex(id(),i)] = make_pair( RandomNumberHistogram(),0.);
}
if ( initialized() )
return;
if ( !sampler()->grids().children().empty() ) {
nIterations(1);
}
if ( !integrated() ) {
unsigned long points = initialPoints();
for ( unsigned long k = 0; k < nIterations(); ++k ) {
runIteration(points,progress);
if ( k < nIterations() - 1 ) {
points = (unsigned long)(points*enhancementFactor());
adapt();
nextIteration();
}
}
}
isInitialized();
}
void BinSampler::finalize(bool){
if (theRandomNumbers!="")
for ( map<RandomNumberIndex,pair<RandomNumberHistogram,double> >::
const_iterator b = RandomNumberHistograms.begin();
b != RandomNumberHistograms.end(); ++b ) {
b->second.first.dump(randomNumberString(), b->first.first,shortprocess(),b->first.second);
}
}
BinSampler::RandomNumberHistogram::
RandomNumberHistogram(double low,
double up,
unsigned int nbins)
: lower(low) {
nbins = nbins + 1;
double c = up / (nbins-1.);
for ( unsigned int k = 1; k < nbins; ++k ) {
bins[low+c*k] = 0.;
binsw1[low+c*k] = 0.;
}
}
void BinSampler::RandomNumberHistogram::
dump(const std::string& folder,const std::string& prefix, const std::string& process,
const int NR) const {
ostringstream fname("");
std::string prefix2;
std::string prefix3=prefix;
std::remove_copy(prefix.begin(), prefix.end(), std::back_inserter(prefix2), '.');
prefix3=prefix2;prefix2.clear();
std::remove_copy(prefix3.begin(), prefix3.end(), std::back_inserter(prefix2), ':');
prefix3=prefix2;prefix2.clear();
std::remove_copy(prefix3.begin(), prefix3.end(), std::back_inserter(prefix2), ',');
fname << "RN-"<< NR ;
ofstream out((folder+"/"+prefix2+fname.str()+".dat").c_str());
double sumofweights=0.;
for ( map<double,double >::const_iterator b = bins.begin();b != bins.end(); ++b )
sumofweights+=b->second;
double sumofweights2=0.;
for ( map<double,double >::const_iterator b = binsw1.begin();b != binsw1.end(); ++b )
sumofweights2+=b->second;
map<double,double >::const_iterator b2 = binsw1.begin();
if ( sumofweights == 0 ) {
cerr << "Not enough statistic accumulated for "
<< process << " skipping random number diagnostic.\n"
<< flush;
return;
}
for ( map<double,double >::const_iterator b = bins.begin();
b != bins.end(); ++b, ++b2) {
out << " " << b->first
<< " " << b->second/sumofweights*100.
<< " " << b2->second/sumofweights2*100.
<< "\n" << flush;
}
double xmin = -0.01;
double xmax = 1.01;
ofstream gpout((folder+"/"+prefix2+fname.str()+".gp").c_str());
gpout << "set terminal epslatex color solid\n"
<< "set output '" << prefix2+fname.str() << "-plot.tex'\n"
<< "set xrange [" << xmin << ":" << xmax << "]\n";
gpout << "set xlabel 'rn "<<NR <<"' \n";
gpout << "set size 0.5,0.6\n";
gpout << "plot '" << prefix2+fname.str()
<< ".dat' u ($1):($2) w boxes lc rgbcolor \"blue\" t '{\\tiny "<<process <<" }',";
gpout << " '" << prefix2+fname.str();
gpout << ".dat' u ($1):($3) w boxes lc rgbcolor \"red\" t '';";
gpout << "reset\n";
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void BinSampler::persistentOutput(PersistentOStream & os) const {
MultiIterationStatistics::put(os);
os << theBias << theWeighted << theInitialPoints << theNIterations
<< theEnhancementFactor << theNonZeroInPresampling << theHalfPoints
<< theMaxNewMax << theReferenceWeight
<< theBin << theInitialized << theLastPoint
<< theEventHandler << theSampler << theRandomNumbers
<< theRemapperPoints << theRemapChannelDimension
<< theLuminosityMapperBins << theGeneralMapperBins << theKappa;
}
void BinSampler::persistentInput(PersistentIStream & is, int) {
MultiIterationStatistics::get(is);
is >> theBias >> theWeighted >> theInitialPoints >> theNIterations
>> theEnhancementFactor >> theNonZeroInPresampling >> theHalfPoints
>> theMaxNewMax >> theReferenceWeight
>> theBin >> theInitialized >> theLastPoint
>> theEventHandler >> theSampler >> theRandomNumbers
>> theRemapperPoints >> theRemapChannelDimension
>> theLuminosityMapperBins >> theGeneralMapperBins >> theKappa;
}
// *** 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<BinSampler,MultiIterationStatistics>
describeHerwigBinSampler("Herwig::BinSampler", "HwSampling.so");
void BinSampler::Init() {
static ClassDocumentation<BinSampler> documentation
("BinSampler samples XCombs bins. This default implementation performs flat MC integration.");
static Parameter<BinSampler,unsigned long> interfaceInitialPoints
("InitialPoints",
"The number of points to use for initial integration.",
&BinSampler::theInitialPoints, 1000000, 1, 0,
false, false, Interface::lowerlim);
static Parameter<BinSampler,size_t> interfaceNIterations
("NIterations",
"The number of iterations to perform initially.",
&BinSampler::theNIterations, 1, 1, 0,
false, false, Interface::lowerlim);
static Parameter<BinSampler,double> interfaceEnhancementFactor
("EnhancementFactor",
"The enhancement factor for the number of points in the next iteration.",
&BinSampler::theEnhancementFactor, 2.0, 1.0, 0,
false, false, Interface::lowerlim);
static Switch<BinSampler,bool> interfaceNonZeroInPresampling
("NonZeroInPresampling",
"Switch on to count only non zero weights in presampling.",
&BinSampler::theNonZeroInPresampling, true, false, false);
static SwitchOption interfaceNonZeroInPresamplingYes
(interfaceNonZeroInPresampling,
"Yes",
"",
true);
static SwitchOption interfaceNonZeroInPresamplingNo
(interfaceNonZeroInPresampling,
"No",
"",
false);
static Switch<BinSampler,bool> interfaceHalfPoints
("HalfPoints",
"Switch on to reset the counter of points if new maximumis was found in the first 1/2 points.",
&BinSampler::theHalfPoints, true, false, false);
static SwitchOption interfaceHalfPointsYes
(interfaceHalfPoints,
"Yes",
"",
true);
static SwitchOption interfaceHalfPointsNo
(interfaceHalfPoints,
"No",
"",
false);
static Parameter<BinSampler,int> interfaceMaxNewMax
("MaxNewMax",
"The maximum number of allowed new maxima in combination with the HalfPoints option.",
&BinSampler::theMaxNewMax, 30, 1, 0,
false, false, Interface::lowerlim);
static Parameter<BinSampler,string> interfaceRandomNumbers
("RandomNumbers",
"Prefix for distributions of the random numbers.",
&BinSampler::theRandomNumbers, "",
false, false);
static Parameter<BinSampler,unsigned long> interfaceRemapperPoints
("RemapperPoints",
"The number of points to be used for filling remappers.",
&BinSampler::theRemapperPoints, 10000, 0, 0,
false, false, Interface::lowerlim);
static Switch<BinSampler,bool> interfaceRemapChannelDimension
("RemapChannelDimension",
"Switch on remapping of the channel dimension.",
&BinSampler::theRemapChannelDimension, true, false, false);
static SwitchOption interfaceRemapChannelDimensionYes
(interfaceRemapChannelDimension,
"Yes",
"",
true);
static SwitchOption interfaceRemapChannelDimensionNo
(interfaceRemapChannelDimension,
"No",
"",
false);
static Parameter<BinSampler,unsigned long> interfaceLuminosityMapperBins
("LuminosityMapperBins",
"The number of bins to be used for remapping parton luminosities.",
&BinSampler::theLuminosityMapperBins, 0, 0, 0,
false, false, Interface::lowerlim);
static Parameter<BinSampler,unsigned long> interfaceGeneralMapperBins
("GeneralMapperBins",
"The number of bins to be used for remapping other phase space dimensions.",
&BinSampler::theGeneralMapperBins, 0, 0, 0,
false, false, Interface::lowerlim);
static Parameter<BinSampler,double> interfaceRemapperMinSelection
("RemapperMinSelection",
"The minimum bin selection probability for remappers.",
&BinSampler::theRemapperMinSelection, 0.00001, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<BinSampler,double> interfaceKappa
("Kappa",
"In AllmostUnweighted mode unweight to Kappa ReferenceWeight.",
&BinSampler::theKappa, 1., 0.000001, 1.0,
false, false, Interface::limited);
}
diff --git a/Sampling/CellGrids/CellGridSampler.cc b/Sampling/CellGrids/CellGridSampler.cc
--- a/Sampling/CellGrids/CellGridSampler.cc
+++ b/Sampling/CellGrids/CellGridSampler.cc
@@ -1,360 +1,360 @@
// -*- C++ -*-
//
// CellGridSampler.cpp is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the CellGridSampler class.
//
#include "CellGridSampler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Handlers/StandardEventHandler.h"
#include "ThePEG/Handlers/StandardXComb.h"
-#include <boost/progress.hpp>
+#include "Herwig/Utilities/Progress.h"
#include "CellGridSampler.h"
#include "Herwig/Sampling/GeneralSampler.h"
using namespace Herwig;
using namespace ExSample;
CellGridSampler::CellGridSampler()
: BinSampler(), SimpleCellGrid(),
theExplorationPoints(1000), theExplorationSteps(8),
theGain(0.3), theEpsilon(0.01),
theMinimumSelection(0.0001), theLuminositySplits(0),
theChannelSplits(0), theAllChannelSplits(false),
theUnweightCells(true) {}
CellGridSampler::~CellGridSampler() {}
IBPtr CellGridSampler::clone() const {
return new_ptr(*this);
}
IBPtr CellGridSampler::fullclone() const {
return new_ptr(*this);
}
double CellGridSampler::generate() {
UseRandom rnd;
double w = SimpleCellGrid::sample(rnd,*this,lastPoint(),
!weighted() && initialized() && theUnweightCells,
!initialized());
if ( !weighted() && initialized() ) {
double p = min(abs(w),kappa()*referenceWeight())/(kappa()*referenceWeight());
double sign = w >= 0. ? 1. : -1.;
if ( p < 1 && UseRandom::rnd() > p )
w = 0.;
else
w = sign*max(abs(w),referenceWeight()*kappa());
}
select(w);
if ( w != 0.0 )
accept();
assert(kappa()==1.||sampler()->almostUnweighted());
return w;
}
void CellGridSampler::adapt() {
UseRandom rnd;
set<SimpleCellGrid*> newCells;
SimpleCellGrid::adapt(theGain,theEpsilon,newCells);
SimpleCellGrid::explore(theExplorationPoints,rnd,*this,newCells,Repository::clog());
SimpleCellGrid::setWeights();
SimpleCellGrid::updateIntegral();
SimpleCellGrid::minimumSelection(theMinimumSelection);
}
void CellGridSampler::saveGrid() const {
XML::Element grid = SimpleCellGrid::toXML();
grid.appendAttribute("process",id());
sampler()->grids().append(grid);
}
bool CellGridSampler::existsGrid() const {
list<XML::Element>::iterator git = sampler()->grids().children().begin();
for ( ; git != sampler()->grids().children().end(); ++git ) {
if ( git->type() != XML::ElementTypes::Element )
continue;
if ( git->name() != "CellGrid" )
continue;
string proc;
git->getFromAttribute("process",proc);
if ( proc == id() )
return true;
}
return false;
}
void CellGridSampler::initialize(bool progress) {
bool haveGrid = false;
list<XML::Element>::iterator git = sampler()->grids().children().begin();
for ( ; git != sampler()->grids().children().end(); ++git ) {
if ( git->type() != XML::ElementTypes::Element )
continue;
if ( git->name() != "CellGrid" )
continue;
string proc;
git->getFromAttribute("process",proc);
if ( proc == id() ) {
haveGrid = true;
break;
}
}
if ( haveGrid ) {
SimpleCellGrid::fromXML(*git);
sampler()->grids().erase(git);
didReadGrids();
}
lastPoint().resize(dimension());
if (randomNumberString()!="")
for(size_t i=0;i<lastPoint().size();i++){
RandomNumberHistograms[RandomNumberIndex(id(),i)] = make_pair( RandomNumberHistogram(),0.);
}
if ( initialized() ) {
if ( !hasGrids() )
throw Exception() << "CellGridSampler: Require existing grid when starting to run.\n"
<< "Did you miss setting --setupfile?"
<< Exception::abortnow;
return;
}
if ( haveGrid ) {
if ( !integrated() )
runIteration(initialPoints(),progress);
isInitialized();
return;
}
SimpleCellGrid::boundaries(vector<double>(dimension(),0.0),vector<double>(dimension(),1.0));
SimpleCellGrid::weightInformation().resize(dimension());
UseRandom rnd;
- boost::progress_display* progressBar = 0;
+ progress_display* progressBar = nullptr;
if ( progress ) {
Repository::clog() << "exploring " << process();
- progressBar = new boost::progress_display(theExplorationSteps,Repository::clog());
+ progressBar = new progress_display{ theExplorationSteps, Repository::clog() };
}
std::set<SimpleCellGrid*> newCells;
if ( pre_adaption_splits().empty() &&
(theLuminositySplits || theChannelSplits || theAllChannelSplits) ) {
const StandardEventHandler& eh = *eventHandler();
const StandardXComb& xc = *eh.xCombs()[bin()];
the_pre_adaption_splits.resize(dimension(),0);
const pair<int,int>& pdims = xc.partonDimensions();
if ( theLuminositySplits && dimension() >= pdims.first + pdims.second ) {
for ( int n = 0; n < pdims.first; ++n )
the_pre_adaption_splits[n] = theLuminositySplits;
for ( int n = dimension() - pdims.second; n < dimension(); ++n )
the_pre_adaption_splits[n] = theLuminositySplits;
}
if ( theChannelSplits && xc.diagrams().size() &&
dimension() > pdims.first + pdims.second ) {
the_pre_adaption_splits[pdims.first] = theChannelSplits;
}
if ( theAllChannelSplits && xc.diagrams().size() > 1 &&
dimension() > pdims.first + pdims.second ) {
the_pre_adaption_splits[pdims.first] = xc.diagrams().size() - 1;
}
}
for(int splitdim=0; splitdim<min(dimension(),(int)pre_adaption_splits().size());splitdim++)
SimpleCellGrid::splitter(splitdim,pre_adaption_splits()[splitdim]);
SimpleCellGrid::explore(theExplorationPoints,rnd,*this,newCells,Repository::clog());
bool notAll = false;
for ( std::size_t step = 1; step < theExplorationSteps; ++step ) {
newCells.clear();
SimpleCellGrid::adapt(theGain,theEpsilon,newCells);
if ( progressBar )
++(*progressBar);
if ( newCells.empty() ) {
notAll = true;
break;
}
SimpleCellGrid::explore(theExplorationPoints,rnd,*this,newCells,Repository::clog());
}
if ( progressBar )
++(*progressBar);
SimpleCellGrid::setWeights();
SimpleCellGrid::updateIntegral();
SimpleCellGrid::minimumSelection(theMinimumSelection);
if ( progressBar ) {
if ( notAll )
cout << "\n" << flush;
delete progressBar;
}
unsigned long points = initialPoints();
for ( unsigned long k = 0; k < nIterations(); ++k ) {
runIteration(points,progress);
if ( k < nIterations() - 1 ) {
points = (unsigned long)(points*enhancementFactor());
adapt();
nextIteration();
}
}
didReadGrids();
isInitialized();
}
void CellGridSampler::finalize(bool) {
XML::Element grid = SimpleCellGrid::toXML();
grid.appendAttribute("process",id());
sampler()->grids().append(grid);
if (randomNumberString()!="")
for ( map<RandomNumberIndex,pair<RandomNumberHistogram,double> >::
const_iterator b = RandomNumberHistograms.begin();
b != RandomNumberHistograms.end(); ++b ) {
b->second.first.dump(randomNumberString(), b->first.first,shortprocess(),b->first.second);
}
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void CellGridSampler::persistentOutput(PersistentOStream & os) const {
os << theExplorationPoints << theExplorationSteps
<< theGain << theEpsilon << theMinimumSelection
<< the_pre_adaption_splits
<< theLuminositySplits << theChannelSplits
<< theAllChannelSplits << theUnweightCells;
}
void CellGridSampler::persistentInput(PersistentIStream & is, int) {
is >> theExplorationPoints >> theExplorationSteps
>> theGain >> theEpsilon >> theMinimumSelection
>> the_pre_adaption_splits
>> theLuminositySplits >> theChannelSplits
>> theAllChannelSplits >> theUnweightCells;
}
// *** 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<CellGridSampler,BinSampler>
describeHerwigCellGridSampler("Herwig::CellGridSampler", "HwSampling.so");
void CellGridSampler::Init() {
static ClassDocumentation<CellGridSampler> documentation
("CellGridSampler samples XCombs bins using CellGrids.");
static Parameter<CellGridSampler,size_t> interfaceExplorationPoints
("ExplorationPoints",
"The number of points to use for cell exploration.",
&CellGridSampler::theExplorationPoints, 1000, 1, 0,
false, false, Interface::lowerlim);
static Parameter<CellGridSampler,size_t> interfaceExplorationSteps
("ExplorationSteps",
"The number of exploration steps to perform.",
&CellGridSampler::theExplorationSteps, 8, 1, 0,
false, false, Interface::lowerlim);
static Parameter<CellGridSampler,double> interfaceGain
("Gain",
"The gain factor used for adaption.",
&CellGridSampler::theGain, 0.3, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<CellGridSampler,double> interfaceEpsilon
("Epsilon",
"The efficieny threshold used for adaption.",
&CellGridSampler::theEpsilon, 0.01, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<CellGridSampler,double> interfaceMinimumSelection
("MinimumSelection",
"The minimum cell selection probability.",
&CellGridSampler::theMinimumSelection, 0.0001, 0.0, 1.0,
false, false, Interface::limited);
static ParVector<CellGridSampler,int> interfacethe_pre_adaption_splits
("preadaptionsplit",
"The splittings for each dimension befor adaption.",
&CellGridSampler::the_pre_adaption_splits, 1., -1, 0.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<CellGridSampler,int> interfaceLuminositySplits
("LuminositySplits",
"",
&CellGridSampler::theLuminositySplits, 0, 0, 0,
false, false, Interface::lowerlim);
static Parameter<CellGridSampler,int> interfaceChannelSplits
("ChannelSplits",
"",
&CellGridSampler::theChannelSplits, 0, 0, 0,
false, false, Interface::lowerlim);
static Switch<CellGridSampler,bool> interfaceAllChannelSplits
("AllChannelSplits",
"",
&CellGridSampler::theAllChannelSplits, false, false, false);
static SwitchOption interfaceAllChannelSplitsYes
(interfaceAllChannelSplits,
"Yes",
"",
true);
static SwitchOption interfaceAllChannelSplitsNo
(interfaceAllChannelSplits,
"No",
"",
false);
static Switch<CellGridSampler,bool> interfaceUnweightCells
("UnweightCells",
"",
&CellGridSampler::theUnweightCells, true, false, false);
static SwitchOption interfaceUnweightCellsYes
(interfaceUnweightCells,
"Yes",
"",
true);
static SwitchOption interfaceUnweightCellsNo
(interfaceUnweightCells,
"No",
"",
false);
}
diff --git a/Sampling/GeneralSampler.cc b/Sampling/GeneralSampler.cc
--- a/Sampling/GeneralSampler.cc
+++ b/Sampling/GeneralSampler.cc
@@ -1,1074 +1,1072 @@
// -*- C++ -*-
//
// GeneralSampler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the GeneralSampler class.
//
#include "GeneralSampler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Utilities/LoopGuard.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Utilities/ColourOutput.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Handlers/StandardEventHandler.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig/API/RunDirectories.h"
#include "Herwig/API/Filesystem.h"
#include "Herwig/Utilities/XML/ElementIO.h"
+#include "Herwig/Utilities/Progress.h"
-#include <boost/progress.hpp>
#include <cstdlib>
#include <sstream>
using namespace Herwig;
GeneralSampler::GeneralSampler()
: theVerbose(false),
theIntegratedXSec(ZERO), theIntegratedXSecErr(ZERO),
theUpdateAfter(1), crossSectionCalls(0), gotCrossSections(false),
theSumWeights(0.), theSumWeights2(0.),
theAttempts(0), theAccepts(0),
theMaxWeight(0.0), theAddUpSamplers(false),
theGlobalMaximumWeight(true), theFlatSubprocesses(false),
isSampling(false), theMinSelection(0.01), runCombinationData(false),
theAlmostUnweighted(false), maximumExceeds(0),
maximumExceededBy(0.), correctWeights(0.),theMaxEnhancement(1.05), didReadGrids(false),
theParallelIntegration(false),
theIntegratePerJob(0), theIntegrationJobs(0), theIntegrationJobsCreated(0),
justAfterIntegrate(false), theWriteGridsOnFinish(false) {}
GeneralSampler::~GeneralSampler() {}
IBPtr GeneralSampler::clone() const {
return new_ptr(*this);
}
IBPtr GeneralSampler::fullclone() const {
return new_ptr(*this);
}
double sign(double x) {
return x >= 0. ? 1. : -1.;
}
void GeneralSampler::initialize() {
if ( theParallelIntegration &&
runLevel() == ReadMode )
throw Exception()
<< "\n----------------------------------------------------\n\n"
<< "Parallel integration is only supported\n in the build/integrate/run mode\n\n"
<< "----------------------------------------------------\n"
<< Exception::abortnow;
if ( runLevel() == ReadMode ||
runLevel() == IntegrationMode ) {
assert(theSamplers.empty());
}
if ( theParallelIntegration ) {
if ( !theIntegratePerJob && !theIntegrationJobs )
throw Exception()
<< "Please specify the number of subprocesses per integration job or the "
<< "number of integration jobs to be created."
<< Exception::abortnow;
if ( theIntegrationJobs ) {
unsigned int nintegrate = eventHandler()->nBins()/theIntegrationJobs;
if ( eventHandler()->nBins() % theIntegrationJobs != 0 )
++nintegrate;
theIntegratePerJob = max(theIntegratePerJob,nintegrate);
}
unsigned int jobCount = 0;
ofstream* jobList = 0;
generator()->log()
<< "----------------------------------------------------\n"
<< "preparing integration jobs ...\n" << flush;
vector<int> randomized;
vector<int> pickfrom;
for ( int b = 0; b < eventHandler()->nBins(); ++b )
pickfrom.push_back(b);
//set<int> check;
while ( !pickfrom.empty() ) {
size_t idx = UseRandom::irnd(pickfrom.size());
randomized.push_back(pickfrom[idx]);
pickfrom.erase(pickfrom.begin() + idx);
}
int b = 0;
for ( vector<int>::const_iterator bx = randomized.begin();
bx != randomized.end(); ++bx, ++b ) {
if ( b == 0 || b % theIntegratePerJob == 0 ) {
if ( jobList ) {
jobList->close();
delete jobList;
jobList = 0;
}
ostringstream name;
string prefix = RunDirectories::runStorage();
if ( prefix.empty() )
prefix = "./";
else if ( *prefix.rbegin() != '/' )
prefix += "/";
name << prefix << "integrationJob" << jobCount<<"List";
++jobCount;
string fname = name.str();
jobList = new ofstream(fname.c_str());
if ( !*jobList ) {
delete jobList;
throw Exception() << "Failed to write integration job list"
<< Exception::abortnow;
}
}
*jobList << *bx << " ";
}
theIntegrationJobsCreated = jobCount;
generator()->log()
<< "---------------------------------------------------\n\n"
<< "Wrote " << jobCount << " integration jobs\n"
<< "Please submit integration jobs with the\nintegrate --jobid=x\ncommand for job ids "
<< "from 0 to " << (jobCount-1) << "\n\n"
<< "e.g.:\n\n" << ANSI::yellow
<< " for i in $(seq 0 "<< (jobCount-1) <<");do Herwig integrate --jobid=$i "<<generator()->runName()<<".run & done \n\n" << ANSI::reset
<< "---------------------------------------------------\n"
<< flush;
if ( jobList ) {
jobList->close();
delete jobList;
jobList = 0;
}
theParallelIntegration = false;
return;
}
if ( runLevel() == BuildMode )
return;
if ( !samplers().empty() )
return;
if ( binSampler()->adaptsOnTheFly() ) {
if ( !theAddUpSamplers ) {
Repository::clog() << "Warning: On-the-fly adapting samplers require cross section calculation from "
<< "adding up individual samplers. The AddUpSamplers flag will be switched on.";
}
theAddUpSamplers = true;
}
if ( !weighted() && !binSampler()->canUnweight() )
throw Exception() << "Unweighted events requested from weighted bin sampler object.";
if ( theFlatSubprocesses && !theGlobalMaximumWeight ) {
Repository::clog() << "Warning: Can only use a global maximum weight when selecting subprocesses "
<< "uniformly. The GlobalMaximumWeight flag will be switched on.";
theGlobalMaximumWeight = true;
}
set<int> binsToIntegrate;
if ( integrationList() != "" ) {
string prefix = RunDirectories::runStorage();
assert ( !prefix.empty() );
string fname = prefix.substr(0, prefix.size()-1) + "List";
ifstream jobList(fname.c_str());
if ( jobList ) {
int b = 0;
while ( jobList >> b )
binsToIntegrate.insert(b);
} else {
Repository::clog()
<< "Job list '"
<< fname << "' not found.\n"
<< "Assuming empty integration job\n" << flush;
return;
}
}
if ( binsToIntegrate.empty() ) {
for ( int b = 0; b < eventHandler()->nBins(); ++b )
binsToIntegrate.insert(b);
}
- boost::progress_display* progressBar = 0;
+ progress_display* progressBar = nullptr;
if ( !theVerbose && !justAfterIntegrate ) {
Repository::clog() << "integrating subprocesses";
- progressBar = new boost::progress_display(binsToIntegrate.size(),Repository::clog());
+ progressBar = new progress_display{ binsToIntegrate.size(), Repository::clog() };
}
int count=0;
bool reuseGrid = false;
bool missingGrid = false;
for ( set<int>::const_iterator bit = binsToIntegrate.begin(); bit != binsToIntegrate.end(); ++bit ) {
count++;
if(theVerbose&&
(runLevel() == ReadMode ||
runLevel() == IntegrationMode))
cout<<"\nIntegrate "<< count <<" of "<<binsToIntegrate.size() <<":\n"<<flush;
Ptr<BinSampler>::ptr s = theBinSampler->cloneMe();
s->eventHandler(eventHandler());
s->sampler(this);
s->bin(*bit);
lastSampler(s);
s->doWeighted(eventHandler()->weighted());
s->setupRemappers(theVerbose);
if ( justAfterIntegrate )
s->readIntegrationData();
reuseGrid = reuseGrid || s->existsGrid();
missingGrid = missingGrid || ( ! s->existsGrid() );
s->initialize(theVerbose);
samplers()[*bit] = s;
if ( !theVerbose && !justAfterIntegrate )
++(*progressBar);
if ( s->nanPoints() && theVerbose ) {
Repository::clog() << "warning: "
<< s->nanPoints() << " of "
<< s->allPoints() << " points with nan or inf weight.\n"
<< flush;
}
}
- if ( progressBar ) {
- delete progressBar;
- progressBar = 0;
- }
-
+ delete progressBar;
+ progressBar = nullptr;
+
if ( missingGrid && runLevel() == RunMode )
generator()->log()
<< "\n----------------------------------------------------\n\n"
<< "Warning:No grid file could be found at the start of\n"
<< "this run.\n\n"
<< "* For a read/run setup intented to be used with \n"
<< " --setupfile \n"
<< " please consider using the build/integrate/run setup.\n"
<< "* For a build/integrate/run setup to be used with\n"
<< " --setupfile\n"
<< " please ensure that the same setupfile is provided\n"
<< " to both the integrate and run steps.\n\n"
<< "---------------------------------------------------\n" << flush;
if ( runLevel() == ReadMode ||
runLevel() == IntegrationMode ) {
if ( reuseGrid )
Repository::clog()
<< "----------------------------------------------------\n\n"
<< "Re-using an existing grid as starting point for grid\n"
<< "optimization. Please consider removing the grid files \n"
<< "and re-running the grid adaption when there have\n"
<< "been significant changes to parameters, cuts, etc.\n\n"
<< "----------------------------------------------------\n"
<< flush;
}
if ( runLevel() == IntegrationMode ) {
theGrids = XML::Element(XML::ElementTypes::Element,"Grids");
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s ) {
s->second->saveGrid();
s->second->saveRemappers();
s->second->saveIntegrationData();
}
writeGrids();
return;
}
if ( theVerbose ) {
bool oldAdd = theAddUpSamplers;
theAddUpSamplers = true;
try {
Repository::clog() << "estimated total cross section is ( "
<< integratedXSec()/nanobarn << " +/- "
<< integratedXSecErr()/nanobarn << " ) nb\n" << flush;
} catch (...) {
theAddUpSamplers = oldAdd;
throw;
}
theAddUpSamplers = oldAdd;
}
updateSamplers();
if ( samplers().empty() ) {
throw Exception() << "No processes with non-zero cross section present."
<< Exception::abortnow;
}
if ( !justAfterIntegrate ) {
theGrids = XML::Element(XML::ElementTypes::Element,"Grids");
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s ) {
s->second->saveGrid();
s->second->saveRemappers();
}
writeGrids();
}
}
double GeneralSampler::generate() {
long excptTries = 0;
gotCrossSections = false;
lastSampler(samplers().upper_bound(UseRandom::rnd())->second);
double weight = 0.;
while ( true ) {
try {
weight = 1.0;
double p = lastSampler()->referenceWeight()/lastSampler()->bias()/theMaxWeight;
if ( weighted() )
weight *= p;
else if ( p < UseRandom::rnd() ){
weight = 0.0;
// The lastSampler was picked according to the bias of the process.
--excptTries;
}
if ( weight != 0.0 )
weight *= lastSampler()->generate()/lastSampler()->referenceWeight();
} catch(BinSampler::NextIteration) {
updateSamplers();
lastSampler(samplers().upper_bound(UseRandom::rnd())->second);
if ( ++excptTries == eventHandler()->maxLoop() )
break;
continue;
} catch (...) {
throw;
}
if ( ! isfinite(lastSampler()->lastWeight()) ) {
lastSampler() = samplers().upper_bound(UseRandom::rnd())->second;
if ( ++excptTries == eventHandler()->maxLoop() )
break;
continue;
}
theAttempts += 1;
if ( abs(weight) == 0.0 ) {
lastSampler(samplers().upper_bound(UseRandom::rnd())->second);
if ( ++excptTries == eventHandler()->maxLoop() )
break;
continue;
}
if ( !eventHandler()->weighted() && !theAlmostUnweighted ) {
if ( abs(weight) > 1. ) {
++maximumExceeds;
maximumExceededBy += abs(weight)-1.;
}
correctWeights+=weight;
if ( weight > 0.0 )
weight = 1.;
else
weight = -1.;
}
break;
}
theAccepts += 1;
if ( excptTries == eventHandler()->maxLoop() )
throw Exception()
<< "GeneralSampler::generate() : Maximum number of tries to re-run event "
<< "selection reached. Aborting now." << Exception::runerror;
lastPoint() = lastSampler()->lastPoint();
lastSampler()->accept();
theSumWeights += weight;
theSumWeights2 += sqr(weight);
return weight;
}
void GeneralSampler::rejectLast() {
if ( !lastSampler() )
return;
double w = 0.0;
if ( weighted() )
w = lastSampler()->lastWeight()/lastSampler()->bias()/theMaxWeight;
else
w = lastSampler()->lastWeight()/lastSampler()->referenceWeight();
lastSampler()->reject();
theSumWeights -= w;
theSumWeights2 -= sqr(w);
theAttempts -= 1;
theAccepts -= 1;
}
void GeneralSampler::updateSamplers() {
map<double,Ptr<BinSampler>::ptr> checkedSamplers;
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s ) {
if ( s->second->averageAbsWeight() == 0.0 ) {
generator()->log() << "Warning: no phase space points with non-zero cross section\n"
<< "could be obtained for the process: "
<< s->second->process() << "\n"
<< "This process will not be considered. Try increasing InitialPoints.\n"
<< flush;
if ( s->second->nanPoints() ) {
generator()->log() << "Warning: "
<< s->second->nanPoints() << " of "
<< s->second->allPoints() << " points with nan or inf weight\n"
<< "in " << s->second->process() << "\n" << flush;
}
continue;
}
checkedSamplers.insert(*s);
}
theSamplers = checkedSamplers;
if ( samplers().empty() )
return;
double allMax = 0.0;
double sumbias = 0.;
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s ) {
double bias = 1.;
if ( !theFlatSubprocesses )
bias *= s->second->averageAbsWeight();
s->second->bias(bias);
sumbias += bias;
allMax = max(allMax,s->second->maxWeight()*theMaxEnhancement);
}
double nsumbias = 0.0;
bool needAdjust = false;
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s ) {
needAdjust |= s->second->bias()/sumbias < theMinSelection;
s->second->bias(max(s->second->bias()/sumbias,theMinSelection));
nsumbias += s->second->bias();
}
if ( nsumbias == 0.0 ) {
samplers().clear();
return;
}
if ( needAdjust ) {
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s ) {
s->second->bias(s->second->bias()/nsumbias);
}
}
theMaxWeight = 0.0;
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s ) {
double wref = theGlobalMaximumWeight ? allMax :
s->second->maxWeight()*theMaxEnhancement;
s->second->referenceWeight(wref);
theMaxWeight = max(theMaxWeight,wref/s->second->bias());
if ( (isSampling && s->second == lastSampler()) ||
!isSampling )
s->second->nextIteration();
}
map<double,Ptr<BinSampler>::ptr> newSamplers;
double current = 0.;
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s ) {
if ( s->second->bias() == 0.0 )
continue;
current += s->second->bias();
newSamplers[current] = s->second;
}
samplers() = newSamplers;
}
void GeneralSampler::currentCrossSections() const {
if ( !theAddUpSamplers ) {
double n = attempts();
if ( n > 1 ) {
theIntegratedXSec = sumWeights()*maxXSec()/attempts();
double sw = sumWeights(); double sw2 = sumWeights2();
theIntegratedXSecErr = maxXSec()*sqrt(abs(sw2/n-sqr(sw/n))/(n-1));
} else {
theIntegratedXSec = ZERO;
theIntegratedXSecErr = ZERO;
}
return;
}
if ( gotCrossSections )
return;
if ( crossSectionCalls > 0 ) {
if ( ++crossSectionCalls == theUpdateAfter ) {
crossSectionCalls = 0;
} else return;
}
++crossSectionCalls;
gotCrossSections = true;
theIntegratedXSec = ZERO;
double var = 0.0;
for ( map<double,Ptr<BinSampler>::ptr>::const_iterator s = samplers().begin();
s != samplers().end(); ++s ) {
theIntegratedXSec += s->second->integratedXSec();
var += sqr(s->second->integratedXSecErr()/nanobarn);
}
theIntegratedXSecErr = sqrt(var)*nanobarn;
}
void GeneralSampler::prepare() {
readGrids();
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void GeneralSampler::doinit() {
if ( RunDirectories::empty() )
RunDirectories::pushRunId(generator()->runName());
if ( integratePerJob() || integrationJobs() ) {
theParallelIntegration = true;
theIntegratePerJob = integratePerJob();
theIntegrationJobs = integrationJobs();
}
readGrids();
bool missingGrid = false;
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s )
missingGrid = missingGrid || ( ! s->second->existsGrid() );
if ( missingGrid && runLevel() == RunMode )
generator()->log()
<< "\n---------------------------------------------------\n\n"
<< "Warning: No grid file could be found at the start of this run.\n\n"
<< "* For a read/run setup intented to be used with --setupfile please consider\n"
<< " using the build/integrate/run setup.\n"
<< "* For a build/integrate/run setup to be used with --setupfile please ensure\n"
<< " that the same setupfile is provided to both, the integrate and run steps.\n\n"
<< "---------------------------------------------------\n" << flush;
if ( samplers().empty() && runLevel() == RunMode )
justAfterIntegrate = true;
SamplerBase::doinit();
}
void GeneralSampler::dofinish() {
set<string> compensating;
for ( map<double,Ptr<BinSampler>::ptr>::const_iterator s =
samplers().begin(); s != samplers().end(); ++s ) {
if ( s->second->compensating() ) {
compensating.insert(s->second->process());
}
if ( s->second->nanPoints() ) {
generator()->log() << "warning: "
<< s->second->nanPoints() << " of "
<< s->second->allPoints() << " points with nan or inf weight\n"
<< "in " << s->second->process() << "\n" << flush;
}
s->second->finalize(theVerbose);
}
if ( theVerbose ) {
if ( !compensating.empty() ) {
generator()->log() << "warning: sampling for the following processes is still compensating:\n";
for ( set<string>::const_iterator c = compensating.begin();
c != compensating.end(); ++c )
generator()->log() << *c << "\n";
}
generator()->log() << "final integrated cross section is ( "
<< integratedXSec()/nanobarn << " +/- "
<< integratedXSecErr()/nanobarn << " ) nb\n" << flush;
}
if ( !compensating.empty() ) {
generator()->log() << "Warning: Some samplers are still in compensating mode.\n" << flush;
}
if ( maximumExceeds != 0 ) {
//generator()->log() << maximumExceeds << " of " << theAttempts
// << " attempted points exceeded the guessed maximum weight\n"
// << "with an average relative deviation of "
// << maximumExceededBy/maximumExceeds << "\n\n" << flush;
generator()->log() <<"\n\n\nNote: In this run "<<maximumExceeds<<" of the "<<theAccepts<<" accepted events\n"
<<"were found with a weight W larger than the expected Wmax.\n";
generator()->log() <<"This corresponds to a cross section difference between:\n"
<<" UnitWeights: "<< theMaxWeight*theSumWeights/theAttempts<<"nb\n"
<<" AlmostUnweighted: "<< theMaxWeight*correctWeights/theAttempts<< "nb\n"<<
" use 'set Sampler:AlmostUnweighted Yes' to switch to non-unit weights.\n\n";
generator()->log() <<"The maximum weight determined in the read/integrate step has been enhanced by \n"<<
" set /Herwig/Samplers/Sampler:MaxEnhancement "<< theMaxEnhancement<<
".\nIf the rate of excessions ("<<(double)maximumExceeds*100/(double)theAccepts<<
"%) or the change of the cross section is large,\nyou can try to:\n\n"<<
"Enhance the number of points used in the read/integrate step\n"<<
" set /Herwig/Samplers/Sampler:BinSampler:InitialPoints ...\n\n"<<
"and/or enhance the reference weight found in the read/integrate step\n"<<
" set /Herwig/Samplers/Sampler:MaxEnhancement 1.x\n\n"<<
"If this does not help (and your process is well defined by cuts)\n"<<
"don't hesitate to contact herwig@projects.hepforge.org.\n\n";
}
if ( runCombinationData ) {
string dataName = RunDirectories::runStorage();
if ( dataName.empty() )
dataName = "./";
else if ( *dataName.rbegin() != '/' )
dataName += "/";
dataName += "HerwigSampling.dat";
ofstream data(dataName.c_str());
double runXSec =
theMaxWeight*theSumWeights/theAttempts;
double runXSecErr =
sqr(theMaxWeight)*(1./theAttempts)*(1./(theAttempts-1.))*
abs(theSumWeights2 - sqr(theSumWeights)/theAttempts);
data << setprecision(17);
data << "CrossSectionCombined "
<< (integratedXSec()/nanobarn) << " +/- "
<< (integratedXSecErr()/nanobarn) << "\n"
<< "CrossSectionRun "
<< runXSec << " +/- " << sqrt(runXSecErr) << "\n"
<< "PointsAttempted " << theAttempts << "\n"
<< "PointsAccepted " << theAccepts << "\n"
<< "SumWeights " << theSumWeights*theMaxWeight << "\n"
<< "SumWeights2 " << theSumWeights2*sqr(theMaxWeight) << "\n"
<< flush;
}
theGrids = XML::Element(XML::ElementTypes::Element,"Grids");
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s ) {
s->second->saveGrid();
s->second->saveRemappers();
if ( justAfterIntegrate )
s->second->saveIntegrationData();
}
if ( theWriteGridsOnFinish )
writeGrids();
SamplerBase::dofinish();
}
void GeneralSampler::doinitrun() {
readGrids();
if ( samplers().empty() ) {
justAfterIntegrate = true;
if ( !hasSetupFile() )
initialize();
} else {
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s ) {
s->second->setupRemappers(theVerbose);
if ( justAfterIntegrate )
s->second->readIntegrationData();
s->second->initialize(theVerbose);
}
}
bool missingGrid = false;
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers().begin();
s != samplers().end(); ++s )
missingGrid = missingGrid || ( ! s->second->existsGrid() );
if ( missingGrid && !didReadGrids )
generator()->log()
<< "\n----------------------------------------------------\n\n"
<< "Warning:No grid file could be found at the start of\n"
<< "this run.\n\n"
<< "* For a read/run setup intented to be used with \n"
<< " --setupfile \n"
<< " please consider using the build/integrate/run setup.\n"
<< "* For a build/integrate/run setup to be used with\n"
<< " --setupfile\n"
<< " please ensure that the same setupfile is provided\n"
<< " to both the integrate and run steps.\n\n"
<< "---------------------------------------------------\n" << flush;
isSampling = true;
SamplerBase::doinitrun();
}
void GeneralSampler::rebind(const TranslationMap & trans) {
for ( map<double,Ptr<BinSampler>::ptr>::iterator s =
samplers().begin(); s != samplers().end(); ++s )
s->second = trans.translate(s->second);
SamplerBase::rebind(trans);
}
IVector GeneralSampler::getReferences() {
IVector ret = SamplerBase::getReferences();
for ( map<double,Ptr<BinSampler>::ptr>::iterator s =
samplers().begin(); s != samplers().end(); ++s )
ret.push_back(s->second);
return ret;
}
void GeneralSampler::writeGrids() const {
if ( theGrids.children().empty() )
return;
string dataName = RunDirectories::runStorage();
if ( dataName.empty() )
dataName = "./";
else if ( *dataName.rbegin() != '/' )
dataName += "/";
dataName += "HerwigGrids.xml";
ofstream out(dataName.c_str());
XML::ElementIO::put(theGrids,out);
}
void GeneralSampler::readGrids() {
// return if grids were already read
if ( didReadGrids )
return;
// check for global HerwigGrids.xml file or combine integration jobs to a global HerwigGrids.xml file
// Show messages of integration job combination only in the first run (if no global HerwigGrids.xml file is found in one of the directories)
// or in case of an error
// Check if a globalHerwigGridsFileFound was found and keep messages in a stringstream buffer beforehand
bool globalHerwigGridsFileFound = false;
bool integrationJobCombinationSuccessful = true;
std::stringstream messageBuffer;
RunDirectories directories;
while ( directories && !didReadGrids ) {
string dataName = directories.nextRunStorage();
if ( dataName.empty() )
dataName = "./";
else if ( *dataName.rbegin() != '/' )
dataName += "/";
string directoryName = dataName;
dataName += "HerwigGrids.xml";
ifstream in(dataName.c_str());
if ( in ) {
theGrids = XML::ElementIO::get(in);
didReadGrids = true;
// Set to true if in any of the directories a global HerwigGrid.xml file was found
globalHerwigGridsFileFound = true;
}
else {
// Check if integrationJob was split and try to merge single integrationJobs together
// integrationJobsCreated() == 0 indicates that parallel integration has not been
// requested, while the parallel integration parameters may well yield a single job
if(integrationJobsCreated() >= 1 && runLevel() == RunMode) {
messageBuffer << "\n\n* Global HerwigGrids.xml file does not exist yet"
<< "\n and integration jobs were split into " << integrationJobsCreated() << " integration jobs."
<< "\n Trying to combine single integration jobs to a global HerwigGrids.xml file"
<< "\n using the following directory " << directoryName << ".";
theGrids = XML::Element(XML::ElementTypes::Element,"Grids");
integrationJobCombinationSuccessful = true;
for(unsigned int currentProcessedIntegrationJobNum = 0; currentProcessedIntegrationJobNum < integrationJobsCreated(); ++currentProcessedIntegrationJobNum) {
ostringstream currentProcessedIntegrationJob;
currentProcessedIntegrationJob << directoryName << "integrationJob" << currentProcessedIntegrationJobNum << "/HerwigGrids.xml";
if(filesystem::exists(currentProcessedIntegrationJob.str())) {
ifstream localGridFileIN(currentProcessedIntegrationJob.str().c_str());
if(localGridFileIN) {
theGrids = theGrids + XML::ElementIO::get(localGridFileIN);
messageBuffer << "\n* Added integration job " << currentProcessedIntegrationJobNum << " to global HerwigGrids.xml file.";
}
else {
integrationJobCombinationSuccessful = false;
messageBuffer << "\n* Could not open/add integration job " << currentProcessedIntegrationJobNum << " to global HerwigGrids.xml file.";
}
}
else {
integrationJobCombinationSuccessful = false;
messageBuffer << "\n* Could not find integration job " << currentProcessedIntegrationJob.str();
}
}
if(integrationJobCombinationSuccessful) {
string globalGridFile = directoryName + "HerwigGrids.xml";
ofstream globalGridFileOF(globalGridFile.c_str());
XML::ElementIO::put(theGrids,globalGridFileOF);
messageBuffer << "\n* Global HerwigGrids.xml file was created, the integration jobs 0 to " << integrationJobsCreated()-1
<< " were combined."
<< "\n* If previous warnings in regards to the HerwigGrids.xml file occured, these can be safely ignored."
<< "\n* Note: This message will occur only in the first run and will be suppressed in further runs.\n"
<< flush;
didReadGrids = true;
}
else {
messageBuffer << "\n* Global HerwigGrids.xml file could not be created due to failed combination of integration jobs."
<< "\n Please check the above-mentioned missing/failed integration jobs which are needed for the combination."
<< "\n* Note: It can be that the HerwigGrids.xml file is searched and can be found in further directories."
<< "\n In this case you can ignore this warning message.\n" << flush;
}
}
}
}
// Show messages if global HerwigGrids.xml file was not found or first combination run
if (!globalHerwigGridsFileFound && (theVerbose || !integrationJobCombinationSuccessful))
BaseRepository::cout() << messageBuffer.str() << "\n" << flush;
if ( !didReadGrids )
theGrids = XML::Element(XML::ElementTypes::Element,"Grids");
}
void GeneralSampler::persistentOutput(PersistentOStream & os) const {
os << theVerbose << theBinSampler << theSamplers << theLastSampler
<< theUpdateAfter << crossSectionCalls << gotCrossSections
<< ounit(theIntegratedXSec,nanobarn)
<< ounit(theIntegratedXSecErr,nanobarn)
<< theSumWeights << theSumWeights2
<< theAttempts << theAccepts << theMaxWeight
<< theAddUpSamplers << theGlobalMaximumWeight
<< theFlatSubprocesses << isSampling << theMinSelection
<< runCombinationData << theAlmostUnweighted << maximumExceeds
<< maximumExceededBy << correctWeights << theMaxEnhancement
<< theParallelIntegration
<< theIntegratePerJob << theIntegrationJobs
<< theIntegrationJobsCreated << theWriteGridsOnFinish;
}
void GeneralSampler::persistentInput(PersistentIStream & is, int) {
is >> theVerbose >> theBinSampler >> theSamplers >> theLastSampler
>> theUpdateAfter >> crossSectionCalls >> gotCrossSections
>> iunit(theIntegratedXSec,nanobarn)
>> iunit(theIntegratedXSecErr,nanobarn)
>> theSumWeights >> theSumWeights2
>> theAttempts >> theAccepts >> theMaxWeight
>> theAddUpSamplers >> theGlobalMaximumWeight
>> theFlatSubprocesses >> isSampling >> theMinSelection
>> runCombinationData >> theAlmostUnweighted >> maximumExceeds
>> maximumExceededBy >> correctWeights >> theMaxEnhancement
>> theParallelIntegration
>> theIntegratePerJob >> theIntegrationJobs
>> theIntegrationJobsCreated >> theWriteGridsOnFinish;
}
// *** 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<GeneralSampler,SamplerBase>
describeHerwigGeneralSampler("Herwig::GeneralSampler", "HwSampling.so");
void GeneralSampler::Init() {
static ClassDocumentation<GeneralSampler> documentation
("A GeneralSampler class");
static Reference<GeneralSampler,BinSampler> interfaceBinSampler
("BinSampler",
"The bin sampler to be used.",
&GeneralSampler::theBinSampler, false, false, true, false, false);
static Parameter<GeneralSampler,size_t> interfaceUpdateAfter
("UpdateAfter",
"Update cross sections every number of events.",
&GeneralSampler::theUpdateAfter, 1, 1, 0,
false, false, Interface::lowerlim);
static Switch<GeneralSampler,bool> interfaceVerbose
("Verbose",
"",
&GeneralSampler::theVerbose, false, false, false);
static SwitchOption interfaceVerboseYes
(interfaceVerbose,
"Yes",
"",
true);
static SwitchOption interfaceVerboseNo
(interfaceVerbose,
"No",
"",
false);
static Switch<GeneralSampler,bool> interfaceAddUpSamplers
("AddUpSamplers",
"Calculate cross sections from adding up individual samplers.",
&GeneralSampler::theAddUpSamplers, false, false, false);
static SwitchOption interfaceAddUpSamplersYes
(interfaceAddUpSamplers,
"Yes",
"",
true);
static SwitchOption interfaceAddUpSamplersNo
(interfaceAddUpSamplers,
"No",
"",
false);
static Switch<GeneralSampler,bool> interfaceGlobalMaximumWeight
("GlobalMaximumWeight",
"Use a global maximum weight instead of partial unweighting.",
&GeneralSampler::theGlobalMaximumWeight, true, false, false);
static SwitchOption interfaceGlobalMaximumWeightYes
(interfaceGlobalMaximumWeight,
"Yes",
"",
true);
static SwitchOption interfaceGlobalMaximumWeightNo
(interfaceGlobalMaximumWeight,
"No",
"",
false);
static Parameter<GeneralSampler,double> interfaceMaxEnhancement
("MaxEnhancement",
"Enhance the maximum reference weight found in the read step.",
&GeneralSampler::theMaxEnhancement, 1.1, 1.0, 1.5,
false, false, Interface::limited);
static Switch<GeneralSampler,bool> interfaceFlatSubprocesses
("FlatSubprocesses",
"[debug] Perform a flat subprocess selection.",
&GeneralSampler::theFlatSubprocesses, false, false, false);
static SwitchOption interfaceFlatSubprocessesYes
(interfaceFlatSubprocesses,
"Yes",
"",
true);
static SwitchOption interfaceFlatSubprocessesNo
(interfaceFlatSubprocesses,
"No",
"",
false);
static Parameter<GeneralSampler,double> interfaceMinSelection
("MinSelection",
"A minimum subprocess selection probability.",
&GeneralSampler::theMinSelection, 0.01, 0.0, 1.0,
false, false, Interface::limited);
static Switch<GeneralSampler,bool> interfaceRunCombinationData
("RunCombinationData",
"",
&GeneralSampler::runCombinationData, false, false, false);
static SwitchOption interfaceRunCombinationDataYes
(interfaceRunCombinationData,
"Yes",
"",
true);
static SwitchOption interfaceRunCombinationDataNo
(interfaceRunCombinationData,
"No",
"",
false);
static Switch<GeneralSampler,bool> interfaceAlmostUnweighted
("AlmostUnweighted",
"",
&GeneralSampler::theAlmostUnweighted, false, false, false);
static SwitchOption interfaceAlmostUnweightedYes
(interfaceAlmostUnweighted,
"Yes",
"",
true);
static SwitchOption interfaceAlmostUnweightedNo
(interfaceAlmostUnweighted,
"No",
"",
false);
static Switch<GeneralSampler,bool> interfaceParallelIntegration
("ParallelIntegration",
"Prepare parallel jobs for integration.",
&GeneralSampler::theParallelIntegration, false, false, false);
static SwitchOption interfaceParallelIntegrationYes
(interfaceParallelIntegration,
"Yes",
"",
true);
static SwitchOption interfaceParallelIntegrationNo
(interfaceParallelIntegration,
"No",
"",
false);
static Parameter<GeneralSampler,unsigned int> interfaceIntegratePerJob
("IntegratePerJob",
"The number of subprocesses to integrate per job.",
&GeneralSampler::theIntegratePerJob, 0, 0, 0,
false, false, Interface::lowerlim);
static Parameter<GeneralSampler,unsigned int> interfaceIntegrationJobs
("IntegrationJobs",
"The maximum number of integration jobs to create.",
&GeneralSampler::theIntegrationJobs, 0, 0, 0,
false, false, Interface::lowerlim);
static Parameter<GeneralSampler,unsigned int> interfaceIntegrationJobsCreated
("IntegrationJobsCreated",
"The number of integration jobs which were actually created.",
&GeneralSampler::theIntegrationJobsCreated, 1, 1, 0,
false, false, Interface::lowerlim);
static Switch<GeneralSampler,bool> interfaceWriteGridsOnFinish
("WriteGridsOnFinish",
"Write grids on finishing a run.",
&GeneralSampler::theWriteGridsOnFinish, false, false, false);
static SwitchOption interfaceWriteGridsOnFinishYes
(interfaceWriteGridsOnFinish,
"Yes",
"",
true);
static SwitchOption interfaceWriteGridsOnFinishNo
(interfaceWriteGridsOnFinish,
"No",
"",
false);
}
diff --git a/Sampling/MonacoSampler.cc b/Sampling/MonacoSampler.cc
--- a/Sampling/MonacoSampler.cc
+++ b/Sampling/MonacoSampler.cc
@@ -1,414 +1,412 @@
// -*- C++ -*-
//
// MonacoSampler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MonacoSampler class.
//
#include "MonacoSampler.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/Repository/Repository.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Handlers/StandardEventHandler.h"
#include "ThePEG/Handlers/StandardXComb.h"
-#include <boost/progress.hpp>
-
#include "MonacoSampler.h"
#include "Herwig/Sampling/GeneralSampler.h"
using namespace Herwig;
MonacoSampler::MonacoSampler()
: BinSampler(),
theAlpha(0.875),
theGridDivisions(48),
theIterationPoints(0) {}
MonacoSampler::~MonacoSampler() {}
IBPtr MonacoSampler::clone() const {
return new_ptr(*this);
}
IBPtr MonacoSampler::fullclone() const {
return new_ptr(*this);
}
double MonacoSampler::generate() {
double w = 1.;
// cout<<"\npoint: ";
std::valarray<int> upperb(dimension());
for ( int k = 0; k < dimension(); ++k ) {
double div = (1 - UseRandom::rnd()) * theGridDivisions;
upperb[k] = static_cast<int>(div);
double gupper, glower;
if ( upperb[k] <= 0 ) {
upperb[k] = 0;
glower = 0.;
gupper = theGrid(k,0);
} else if (upperb[k] >= static_cast<int>(theGridDivisions)) {
upperb[k] = theGridDivisions-1;
glower = theGrid(k,theGridDivisions-2);
gupper = theGrid(k,theGridDivisions-1);
} else {
glower = theGrid(k,upperb[k]-1);
gupper = theGrid(k,upperb[k]);
}
double gdiff = gupper - glower;
lastPoint()[k] = glower + (div-upperb[k])*gdiff;
w *= gdiff * theGridDivisions;
}
// cout<<lastPoint()[k]<<" ";
try {
w *= eventHandler()->dSigDR(lastPoint()) / nanobarn;
} catch (Veto&) {
w = 0.0;
} catch (...) {
throw;
}
// only store numbers
double wgt = w;
if ( ! isfinite(wgt) ) wgt = 0;
// save results for later grid optimization
theIterationPoints++;
for ( int k = 0; k < dimension(); ++k ) {
theGridData(k,upperb[k]) += wgt*wgt;
}
if (randomNumberString()!="")
for ( size_t k = 0; k < lastPoint().size(); ++k ) {
RandomNumberHistograms[RandomNumberIndex(id(),k)].first.book(lastPoint()[k],wgt);
RandomNumberHistograms[RandomNumberIndex(id(),k)].second+=wgt;
}
if ( !weighted() && initialized() ) {
double p = min(abs(w),kappa()*referenceWeight())/(kappa()*referenceWeight());
double sign = w >= 0. ? 1. : -1.;
if ( p < 1 && UseRandom::rnd() > p )
w = 0.;
else
w = sign*max(abs(w),kappa()*referenceWeight());
}
select(w);
assert(kappa()==1.||sampler()->almostUnweighted());
if ( w != 0.0 )
accept();
return w;
}
void MonacoSampler::saveGrid() const {
XML::Element grid = toXML();
grid.appendAttribute("process",id());
sampler()->grids().append(grid);
}
bool MonacoSampler::existsGrid() const {
list<XML::Element>::iterator git = sampler()->grids().children().begin();
for ( ; git != sampler()->grids().children().end(); ++git ) {
if ( git->type() != XML::ElementTypes::Element )
continue;
if ( git->name() != "Monaco" )
continue;
string proc;
git->getFromAttribute("process",proc);
if ( proc == id() )
return true;
}
return false;
}
void MonacoSampler::initialize(bool progress) {
//read in grid
bool haveGrid = false;
list<XML::Element>::iterator git = sampler()->grids().children().begin();
for ( ; git != sampler()->grids().children().end(); ++git ) {
if ( git->type() != XML::ElementTypes::Element )
continue;
if ( git->name() != "Monaco" )
continue;
string proc;
git->getFromAttribute("process",proc);
if ( proc == id() ) {
haveGrid = true;
break;
}
}
if ( haveGrid ) {
fromXML(*git);
sampler()->grids().erase(git);
didReadGrids();
} else {
// flat grid
theGrid.resize(dimension(),theGridDivisions);
for (int k = 0; k < dimension(); k++)
for (size_t l = 0; l < theGridDivisions; l++)
theGrid(k,l) = (l+1)/static_cast<double>(theGridDivisions);
theGridData = boost::numeric::ublas::zero_matrix<double>(dimension(),theGridDivisions);
theIterationPoints = 0;
}
lastPoint().resize(dimension());
if (randomNumberString()!="")
for(size_t i=0;i<lastPoint().size();i++){
RandomNumberHistograms[RandomNumberIndex(id(),i)] = make_pair( RandomNumberHistogram(),0.);
}
if ( initialized() ) {
if ( !hasGrids() )
throw Exception() << "MonacoSampler: Require existing grid when starting to run.\n"
<< "Did you miss setting --setupfile?"
<< Exception::abortnow;
return;
}
if ( haveGrid ) {
if ( !integrated() ) {
runIteration(initialPoints(),progress);
adapt();
}
isInitialized();
return;
}
// if ( !sampler()->grids().children().empty() ) {
// nIterations(1);
// }
unsigned long points = initialPoints();
for ( unsigned long k = 0; k < nIterations(); ++k ) {
runIteration(points,progress);
if ( k < nIterations() - 1 ) {
points = (unsigned long)(points*enhancementFactor());
adapt();
nextIteration();
}
}
adapt();
didReadGrids();
isInitialized();
}
void MonacoSampler::adapt() {
int dim = dimension();
// refine grid
std::valarray<double> gridcumul(dim);
for (int k=0; k<dim; ++k) {
double gridold = theGridData(k,0);
double gridnew = theGridData(k,1);
theGridData(k,0) = (gridold + gridnew) / 2.0;
gridcumul[k] = theGridData(k,0);
for (size_t l=1; l<theGridDivisions-1; ++l) {
theGridData(k,l) = gridold + gridnew;
gridold = gridnew;
gridnew = theGridData(k,l+1);
theGridData(k,l) = (theGridData(k,l) + gridnew) / 3.0;
gridcumul[k] += theGridData(k,l);
}
theGridData(k,theGridDivisions-1) = (gridnew + gridold) / 2.0;
gridcumul[k] += theGridData(k,theGridDivisions-1);
}
for (int k=0; k<dim; ++k) {
double rc = 0.;
std::valarray<double> ri(theGridDivisions);
for (size_t l=0; l<theGridDivisions; ++l) {
ri[l] = 0.;
if ((theGridData(k,l) >= 0) && (gridcumul[k] != 0)) {
theGridData(k,l) = max( 1.0e-30, theGridData(k,l) );
double gpart = gridcumul[k] / theGridData(k,l);
ri[l] = pow( (gpart - 1.0) / (gpart * log( gpart )), theAlpha);
} else {
ri[l] = pow( 1. / log( 1e30 ), theAlpha);
}
rc += ri[l];
}
rc /= theGridDivisions;
double gridold = 0, gridnew = 0.;
double deltar = 0.;
unsigned int m = 0;
std::valarray<double> theGridRowNew(theGridDivisions);
for (size_t l = 0; l < theGridDivisions; ++l) {
deltar += ri[l];
gridold = gridnew;
gridnew = theGrid(k,l);
for (; deltar > rc; m++) {
deltar -= rc;
theGridRowNew[m] = gridnew - (gridnew - gridold) * deltar / ri[l];
}
}
for (size_t l = 0; l < theGridDivisions-1; ++l) {
theGrid(k,l) = theGridRowNew[l];
}
theGrid(k,theGridDivisions-1) = 1.0;
}
theGridData = boost::numeric::ublas::zero_matrix<double>(dimension(),theGridDivisions);
theIterationPoints = 0;
}
void MonacoSampler::finalize(bool) {
// save grid
adapt();
XML::Element grid = MonacoSampler::toXML();
grid.appendAttribute("process",id());
sampler()->grids().append(grid);
if (randomNumberString()!="")
for ( map<RandomNumberIndex,pair<RandomNumberHistogram,double> >::
const_iterator b = RandomNumberHistograms.begin();
b != RandomNumberHistograms.end(); ++b ) {
b->second.first.dump(randomNumberString(), b->first.first,shortprocess(),b->first.second);
}
}
void MonacoSampler::fromXML(const XML::Element& grid) {
int dim = 0;
grid.getFromAttribute("Dimension",dim);
if ( dim != dimension() ) {
throw std::runtime_error("[MonacoSampler] Number of dimensions in grid file does not match expectation.");
}
size_t griddivisions = 0;
grid.getFromAttribute("GridDivisions",griddivisions);
boost::numeric::ublas::matrix<double> tmpgrid(dim,griddivisions);
pair<multimap<pair<int,string>,list<XML::Element>::iterator>::const_iterator,multimap<pair<int,string>,list<XML::Element>::iterator>::const_iterator> cit;
cit = grid.findAll(XML::ElementTypes::Element,"GridVector");
if ( cit.first->second == grid.children().end() )
throw std::runtime_error("[MonacoSampler] Expected a GridVector element.");
for (multimap<pair<int,string>,list<XML::Element>::iterator>::const_iterator iit=cit.first; iit!=cit.second; ++iit) {
const XML::Element& gridvector = *iit->second;
int k = 0;
gridvector.getFromAttribute("Index",k);
if ( k >= dim ) {
throw std::runtime_error("[MonacoSampler] Index of grid dimension larger than grid size.");
} else {
list<XML::Element>::const_iterator git;
git = gridvector.findFirst(XML::ElementTypes::ParsedCharacterData,"");
if ( git == gridvector.children().end() )
throw std::runtime_error("[MonacoSampler] Expected grid data.");
istringstream bdata(git->content());
for ( size_t l = 0; l < griddivisions; ++l ) {
bdata >> tmpgrid(k,l);
}
}
}
// store back into main variable
// if griddivisions do not match, rebin preserving bin density
theGrid.resize(dim,theGridDivisions);
theIterationPoints = 0;
double divratio = griddivisions / static_cast<double>(theGridDivisions);
for (int k = 0; k < dim; k++) {
double xold = 0, xnew = 0, deltar = 0;
size_t l = 0;
for (size_t m = 0; m < griddivisions; m++) {
deltar += 1;
xold = xnew;
xnew = tmpgrid(k,m);
for (; deltar > divratio; l++) {
deltar -= divratio;
theGrid(k,l) = xnew - (xnew - xold) * deltar;
}
}
theGrid(k,theGridDivisions-1) = 1.0;
}
theGridData = boost::numeric::ublas::zero_matrix<double>(dimension(),theGridDivisions);
}
XML::Element MonacoSampler::toXML() const {
XML::Element grid(XML::ElementTypes::Element,"Monaco");
grid.appendAttribute("Dimension",dimension());
grid.appendAttribute("GridDivisions",theGridDivisions);
for ( int k = 0; k < dimension(); ++k ) {
XML::Element gridvector(XML::ElementTypes::Element,"GridVector");
gridvector.appendAttribute("Index",k);
ostringstream bdata;
bdata << setprecision(17);
for ( size_t l = 0; l < theGridDivisions; ++l )
bdata << theGrid(k,l) << " ";
XML::Element belem(XML::ElementTypes::ParsedCharacterData,bdata.str());
gridvector.append(belem);
grid.append(gridvector);
}
return grid;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void MonacoSampler::persistentOutput(PersistentOStream & os) const {
BinSampler::put(os);
os << theAlpha << theGridDivisions;
}
void MonacoSampler::persistentInput(PersistentIStream & is, int) {
BinSampler::get(is);
is >> theAlpha >> theGridDivisions;
}
// *** 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<MonacoSampler,BinSampler>
describeHerwigMonacoSampler("Herwig::MonacoSampler", "HwSampling.so");
void MonacoSampler::Init() {
static ClassDocumentation<MonacoSampler> documentation
("MonacoSampler samples XCombs bins. This implementation performs weighted MC integration using Monaco, an adapted Vegas algorithm.");
static Parameter<MonacoSampler,double> interfaceAlpha
("Alpha",
"Rate of grid modification (0 for no modification).",
&MonacoSampler::theAlpha, 0.875, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<MonacoSampler,size_t> interfaceGridDivisions
("GridDivisions",
"The number of divisions per grid dimension.",
&MonacoSampler::theGridDivisions, 48, 1, 0,
false, false, Interface::lowerlim);
}
diff --git a/Shower/Dipole/Merging/MergingFactory.cc b/Shower/Dipole/Merging/MergingFactory.cc
--- a/Shower/Dipole/Merging/MergingFactory.cc
+++ b/Shower/Dipole/Merging/MergingFactory.cc
@@ -1,650 +1,646 @@
// -*- C++ -*-
//
// MergeboxFactory.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MergeboxFactory class.
//
#include "MergingFactory.h"
#include "Node.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Utilities/ColourOutput.h"
using namespace Herwig;
using std::ostream_iterator;
IBPtr MergingFactory::clone() const {
return new_ptr(*this);
}
IBPtr MergingFactory::fullclone() const {
return new_ptr(*this);
}
void MergingFactory::doinit(){
MatchboxFactory::doinit();
if (subProcessGroups()) {
throw InitException() << "There are no subprocess groups in merging!";
}
}
void MergingFactory::productionMode() {
if(M()<0)
for ( vector<Ptr<MatchboxAmplitude>::ptr>::iterator amp
= amplitudes().begin(); amp != amplitudes().end(); ++amp ) {
Repository::clog() << "One-loop contributions from '"
<< (**amp).name()
<< "' are not required and will be disabled.\n"
<< flush;
(**amp).disableOneLoop();
}
MatchboxFactory::productionMode();
}
void MergingFactory::fillMEsMap() {
olpProcesses().clear();
assert( getProcesses().size() == 1 );
processMap[0] = getProcesses()[0];
if ( MH()->M() >= 0 )
setHighestVirt(processMap[0].size()+MH()->M());
MH()->N0(processMap[0].size());
for ( int i = 1 ; i <= MH()->N() ; ++i ) {
processMap[i] = processMap[i - 1];
processMap[i].push_back("j");
}
for ( int i = 0 ; i <= MH()->N() ; ++i ) {
const bool below_maxNLO = i < MH()->M() + 1;
vector<MatchboxMEBasePtr> ames
= makeMEs(processMap[i], orderInAlphaS() + i, below_maxNLO );
copy(ames.begin(), ames.end(), back_inserter(pureMEsMap()[i]));
}
}
#include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h"
void MergingFactory::prepare_BV(int i) {
// check if we have virtual contributions
bool haveVirtuals = true;
for ( auto born : pureMEsMap()[i]) {
prepareME(born);
if ( born->isOLPTree() ) {
int id = orderOLPProcess(born->subProcess(),
born->matchboxAmplitude(),
ProcessType::treeME2);
born->olpProcess(ProcessType::treeME2,id);
id = orderOLPProcess(born->subProcess(),
born->matchboxAmplitude(),
ProcessType::colourCorrelatedME2);
born->olpProcess(ProcessType::colourCorrelatedME2,id);
bool haveGluon = false;
for ( const auto & p : born->subProcess().legs )
if ( p->id() == 21 ) {
haveGluon = true;
break;
}
if ( haveGluon ) {
id = orderOLPProcess(born->subProcess(),
born->matchboxAmplitude(),
ProcessType::spinColourCorrelatedME2);
born->olpProcess(ProcessType::spinColourCorrelatedME2,id);
}
}
if ( born->isOLPLoop() && i <= MH()->M() ) {
int id = orderOLPProcess(born->subProcess(),
born->matchboxAmplitude(),
ProcessType::oneLoopInterference);
born->olpProcess(ProcessType::oneLoopInterference,id);
if ( !born->onlyOneLoop() && born->needsOLPCorrelators() ) {
id = orderOLPProcess(born->subProcess(),
born->matchboxAmplitude(),
ProcessType::colourCorrelatedME2);
born->olpProcess(ProcessType::colourCorrelatedME2,id);
}
}
haveVirtuals &= born->haveOneLoop();
}
// check for consistent conventions on virtuals, if we are to include MH()->M()
if (!(i > MH()->M()||haveVirtuals))
throw InitException()
<< MH()->M()
<< " NLO corrections requested,\n"
<< "but no virtual contributions are found.";
}
void MergingFactory::prepare_R(int i) {
for ( auto real : pureMEsMap()[i])
prepareME(real);
}
#include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h"
void MergingFactory::getVirtuals(MatchboxMEBasePtr nlo, bool clone){
const auto & partons = nlo->diagrams().front()->partons();
for ( auto I : DipoleRepository::insertionIOperators(dipoleSet()) )
if ( I->apply(partons) ){
auto myI = I;
if ( clone ) myI = I->cloneMe();
nlo->virtuals().push_back(myI);
}
for ( auto PK : DipoleRepository::insertionPKOperators(dipoleSet()) )
if ( PK->apply(partons) ){
auto myPK = PK;
if ( clone ) myPK = PK->cloneMe();
nlo->virtuals().push_back(myPK);
}
}
void MergingFactory::pushB(MatchboxMEBasePtr born, int i) {
MatchboxMEBasePtr bornme = born->cloneMe();
bornme->maxMultCKKW(1);
bornme->minMultCKKW(0);
string pname = fullName() + "/" + bornme->name() + ".Born";
if ( !(generator()->preinitRegister(bornme, pname)) )
throw InitException()
<< "Born ME "<< pname << " already existing.";
if (MH()->gamma()!=1.)
getVirtuals(bornme,false);
NodePtr clusternode = new_ptr(Node(bornme, 0, MH()));
clusternode->deepHead(clusternode);
MH()->firstNodeMap(bornme,clusternode);
bornme->merger(MH());
vector<NodePtr> current = {{clusternode}};
vector<NodePtr> children;
unsigned int k = 1;
while ( ! thePureMEsMap[i - k].empty() ) {
for ( auto tmp : current ){//j
tmp->birth(thePureMEsMap[i - k]);
for ( auto tmpchild : tmp->children() ) {//m
children.push_back(tmpchild);
}
}
current = children;
children.clear();
++k;
}
if ( MH()->N() > i )
bornme->needsCorrelations();
else
bornme->needsNoCorrelations();
bornme->cloneDependencies();
MEs().push_back(bornme);
}
void MergingFactory::pushV(MatchboxMEBasePtr born, int i) {
MatchboxMEBasePtr nlo = born->cloneMe();
string pname = fullName() + "/" + nlo->name() + ".Virtual";
if ( !(generator()->preinitRegister(nlo, pname)) )
throw InitException()
<< "Virtual ME "<< pname << " already existing.";
////////////////////////////////////NLO///////////////////////////
nlo->virtuals().clear();
getVirtuals(nlo , false);
if ( nlo->virtuals().empty() )
throw InitException()
<< "No insertion operators have been found for "
<< born->name() << ".\n";
nlo->doOneLoopNoBorn();
////////////////////////////////////NLO///////////////////////////
NodePtr clusternode = new_ptr(Node(nlo, 0,MH()));
clusternode->deepHead(clusternode);
clusternode->virtualContribution(true);
MH()->firstNodeMap(nlo,clusternode);
nlo->merger(MH());
vector<NodePtr> current = {{clusternode}};
vector<NodePtr> children;
unsigned int k = 1;
while ( ! thePureMEsMap[i - k].empty() ) {
for ( auto tmp : current ){
tmp->birth(thePureMEsMap[i - k]);
for ( auto tmpchild : tmp->children())
children.push_back(tmpchild);
}
current = children;
children.clear();
++k;
}
if ( nlo->isOLPLoop() ) {
int id = orderOLPProcess(nlo->subProcess(),
born->matchboxAmplitude(),
ProcessType::oneLoopInterference);
nlo->olpProcess(ProcessType::oneLoopInterference,id);
if ( !nlo->onlyOneLoop() && nlo->needsOLPCorrelators() ) {
id = orderOLPProcess(nlo->subProcess(),
born->matchboxAmplitude(),
ProcessType::colourCorrelatedME2);
nlo->olpProcess(ProcessType::colourCorrelatedME2,id);
}
}
nlo->needsCorrelations();
nlo->cloneDependencies();
MEs().push_back(nlo);
}
void MergingFactory::pushR(MatchboxMEBasePtr born, int i) {
MatchboxMEBasePtr bornme = born->cloneMe();
string pname = fullName() + "/" + bornme->name() + ".Real";
if ( !(generator()->preinitRegister(bornme, pname)) )
throw InitException()
<< "Subtracted ME " << pname << " already existing.";
NodePtr clusternode = new_ptr(Node(bornme, 1, MH()));
clusternode->deepHead(clusternode);
clusternode->subtractedReal(true);
MH()->firstNodeMap(bornme,clusternode);
bornme->merger(MH());
vector<NodePtr> current = {{clusternode}};
vector<NodePtr> children;
unsigned int k = 1;
while ( ! thePureMEsMap[i - k].empty() ) {
for ( auto tmp : current ){
tmp->birth(thePureMEsMap[i - k]);
for ( auto tmpchild : tmp->children())
children.push_back(tmpchild);
}
current = children;
children.clear();
++k;
}
if(clusternode->children().empty()){
// This is a finite real contribution.
// This process is included in the LO merging.
return;
}
if ( MH()->N() > i ) bornme->needsCorrelations();
else bornme->needsNoCorrelations();
bornme->cloneDependencies(pname);
MEs().push_back(bornme);
}
// MergingFactory should never order OLPs here,
// they're done elsewhere.
void MergingFactory::orderOLPs() {}
#include "ThePEG/Utilities/StringUtils.h"
vector<string> MergingFactory::parseProcess(string in) {
vector<string> process = StringUtils::split(in);
if ( process.size() < 3 )
throw Exception()
<< "MatchboxFactory: Invalid process."<< Exception::runerror;
for ( string & p : process) {
p = StringUtils::stripws(p);
}
theN = 0;
bool prodprocess = true;
vector<string> result;
for ( const string & p : process ) {
if ( p == "->" )
continue;
if (p=="[") {
prodprocess = false;
} else if (p=="]") {
prodprocess = false;
// TODO what if there's stuff after the bracket?
assert( p == process.back() );
break;
} else if (p=="[j") {
prodprocess = false;
++theN;
} else if (p=="j" && !prodprocess) {
++theN;
prodprocess = false;
} else if (p=="j]") {
++theN;
prodprocess = false;
// TODO what if there's stuff after the bracket?
assert( p == process.back() );
break;
} else if ( prodprocess ) {
result.push_back(p);
} else {
throw InitException()
<< "Unknown particle class \"" << p << '"'
<< " in the process definition merging bracket.\n"
<< "Only jets (\"j\") are supported at the moment.";
}
}
return result;
}
-#include <boost/progress.hpp>
+#include "Herwig/Utilities/Progress.h"
void MergingFactory::setup() {
useMe();
DipoleShowerHandlerPtr dsh=dynamic_ptr_cast<DipoleShowerHandlerPtr>(this->CKKWHandler());
if(! dsh )throw InitException() << "The showerhandlerfor the MergingFactory must be the DipoleShower. ";
dsh->setMerger(MH());
MH()->setFactory(this);
MH()->setDipoleShower(dsh);
if(!ransetup){
generator()->log() <<"\nStarting merging setup.\n\n";
olpProcesses().clear();
externalAmplitudes().clear();
// We set the couplings in the ME to be fixed
// and reweight in the history weight for this.
setFixedCouplings(true);
setFixedQEDCouplings(true);
// rebind the particle data objects, can't use rebind() function
for ( auto & g : particleGroups()) {
for ( auto & p : g.second) {
p = getParticleData(p->id());
}
}
const PDVector& partons = particleGroups()["j"];
unsigned int nl = 0;
for ( const auto p : partons ) {
const Energy mass = p->hardProcessMass();
const long pid = p->id();
if ( abs(pid) < 7 && mass == ZERO )
++nl;
if ( pid > 0 && pid < 7 && mass == ZERO )
nLightJetVec( pid );
if ( pid > 0 && pid < 7 && mass != ZERO )
nHeavyJetVec( pid );
}
nLight(nl/2);
const PDVector& partonsInP = particleGroups()["p"];
for ( const auto pip : partonsInP )
if ( pip->id() > 0 && pip->id() < 7 && pip->hardProcessMass() == ZERO )
nLightProtonVec( pip->id() );
// fill the amplitudes
if ( !amplitudes().empty() ) fillMEsMap();
// Use the colour basis of the first element of amplitudes
// to set the large N colour basis for the MergingHelper
assert(!amplitudes().empty() );
if ( !amplitudes()[0]->colourBasis() )
throw Exception() << "MergingFactory::setup(): Expecting a colour basis object."
<< Exception::runerror;
auto largeNBasis =
amplitudes()[0]->colourBasis()->cloneMe();
largeNBasis->clear();
largeNBasis->doLargeN();
MH()->largeNBasis(largeNBasis);
// prepare the Born and virtual matrix elements
for ( int i = 0 ; i <= max(0, MH()->N()) ; ++i ) prepare_BV(i);
// prepare the real emission matrix elements
for ( int i = 0 ; i <= MH()->N() ; ++i ) prepare_R(i);
if (MH()->N()<=MH()->M()) {
throw InitException() << "Merging: The number of NLOs need to be"
<< "\nsmaller than the number of LO processes.\n";
}
// Order the external Amplitudes.
orderOLPs();
// start creating matrix elements
MEs().clear();
// count the subprocesses
size_t numb = 0;
size_t numv = 0;
size_t numr = 0;
for (int i = 0; i <= max(0, MH()->N()) ; ++i ) {
for ( auto born : thePureMEsMap[i] ) {
if (bornContributions()
) {
numb++;
}
}
}
for (int i = 0 ; i <=max(0, MH()->N()); ++i )
for ( auto virt : thePureMEsMap[i] )
if ( virtualContributions() && i <= MH()->M()) {
numv++;
}
for (int i = 1; i <= max(0, MH()->N()) ; ++i )
for ( auto real : thePureMEsMap[i] )
if (realContributions() && i <= MH()->M() + 1 ){
numr++;
}
generator()->log() << ANSI::red << "Preparing Merging: ";
generator()->log() << ANSI::green << numb << " x Born " << ANSI::red;
if (MH()->M()>-1) {
generator()->log() << ANSI::yellow << numv << " x Virtual ";
generator()->log() << ANSI::blue << numr << " x Real " << ANSI::red << flush;
}
- boost::progress_display * progressBar = new boost::progress_display(numb+numv+numr,generator()->log());
+ progress_display progressBar{ numb+numv+numr, generator()->log() };
for (int i = 0; i <= max(0, MH()->N()) ; ++i ){
for ( auto born : thePureMEsMap[i]){
if (bornContributions() ){
pushB(born, i);
generator()->log() << ANSI::green;
- ++(*progressBar);
+ ++progressBar;
generator()->log() << ANSI::reset;
}
}
}
for (int i = 0 ; i <=max(0, MH()->N()); ++i )
for ( auto virt : thePureMEsMap[i])
if ( virtualContributions() && i <= MH()->M()){
pushV(virt, i);
generator()->log() << ANSI::yellow;
- ++(*progressBar);
+ ++progressBar;
}
for (int i = 1; i <= max(0, MH()->N()) ; ++i )
for ( auto real : thePureMEsMap[i] )
if (realContributions()&& i <= MH()->M() + 1 ){
pushR(real, i);
generator()->log() << ANSI::blue;
- ++(*progressBar);
+ ++progressBar;
}
generator()->log() << ANSI::reset;
- delete progressBar;
if ( !externalAmplitudes().empty() ) {
generator()->log() << "Initializing external amplitudes." << endl;
- boost::progress_display * progressBar =
- new boost::progress_display(externalAmplitudes().size(),generator()->log());
- for ( const auto ext : externalAmplitudes()) {
+ progress_display progressBar{ externalAmplitudes().size(), generator()->log() };
+ for ( const auto ext : externalAmplitudes() ) {
if ( ! ext->initializeExternal() ) {
throw InitException()
<< "error: failed to initialize amplitude '" << ext->name() << "'\n";
}
- ++(*progressBar);
+ ++progressBar;
}
- delete progressBar;
generator()->log()
<< "---------------------------------------------------" << endl;
}
if ( !olpProcesses().empty() ) {
generator()->log() << "Initializing one-loop provider(s)." << endl;
map<Ptr<MatchboxAmplitude>::tptr, map<pair<Process, int>, int> > olps;
for (const auto oit : olpProcesses()) {
olps[oit.first] = oit.second;
}
- boost::progress_display * progressBar = new boost::progress_display(olps.size(), generator()->log());
+ progress_display progressBar{ olps.size(), generator()->log() };
for ( const auto olpit : olps ) {
if ( !olpit.first->startOLP(olpit.second) ) {
throw InitException() << "error: failed to start OLP for amplitude '" << olpit.first->name() << "'\n";
}
- ++(*progressBar);
+ ++progressBar;
}
- delete progressBar;
generator()->log()
<< "---------------------------------------------------\n" << flush;
}
generator()->log() <<"\nGenerated "<<MEs().size()<<" Subprocesses.\n"<<flush;
generator()->log()
<< "---------------------------------------------------\n" << flush;
generator()->log() <<"\n\n" << ANSI::red
<<"Note: Due to the unitarization of the higher "
<<"\nmultiplicities, the individual cross sections "
<<"\ngiven in the integration and run step are not"
<<"\nmeaningful without merging." << ANSI::reset << endl;
ransetup=true;
}
}
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
void MergingFactory::persistentOutput(PersistentOStream & os) const {
os
<< theonlymulti
<< ransetup
<< processMap << theMergingHelper <<theM<<theN<<theNonQCDCuts;
}
void MergingFactory::persistentInput(PersistentIStream & is, int) {
is
>> theonlymulti
>> ransetup
>> processMap >> theMergingHelper >>theM>>theN>>theNonQCDCuts;
}
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Command.h"
void MergingFactory::Init() {
static Parameter<MergingFactory, int> interfaceonlymulti("onlymulti",
"Calculate only the ME with k additional partons.",
&MergingFactory::theonlymulti, -1, -1, 0,
false, false, Interface::lowerlim);
static Switch<MergingFactory, bool> interface_Unitarized("Unitarized",
"Unitarize the cross section (default is unitarised. NLO merging must be unitarised).",
&MergingFactory::unitarized, true, false, false);
static SwitchOption interface_UnitarizedYes(interface_Unitarized, "Yes",
"Switch on the unitarized cross section.", true);
static SwitchOption interface_UnitarizedNo(interface_Unitarized, "No",
"Switch off the unitarized cross section.", false);
static Reference<MergingFactory,Merger> interfaceMergingHelper("MergingHelper",
"Pointer to the Merging Helper.",
&MergingFactory::theMergingHelper, false, false, true, true, false);
static Parameter<MergingFactory, int> interfaceaddNLOLegs("NLOProcesses",
"Set the number of virtual corrections to consider. 0 is default for no virtual correction.",
&MergingFactory::theM, 0, 0, 0, false, false, Interface::lowerlim);
static Reference<MergingFactory, Cuts > interfaceNonQcdCuts("NonQCDCuts",
"Cut on non-QCD modified observables. Be carefull!",
&MergingFactory::theNonQCDCuts, false, false, true, true, false);
}
// *** 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<MergingFactory, Herwig::MatchboxFactory>
describeHerwigMergingFactory("Herwig::MergingFactory", "HwDipoleShower.so");
diff --git a/Utilities/Makefile.am b/Utilities/Makefile.am
--- a/Utilities/Makefile.am
+++ b/Utilities/Makefile.am
@@ -1,44 +1,45 @@
SUBDIRS = XML Statistics
noinst_LTLIBRARIES = libHwUtils.la
libHwUtils_la_SOURCES = \
EnumParticles.h \
Interpolator.tcc Interpolator.h \
Kinematics.cc Kinematics.h \
+Progress.h Progress.cc \
Maths.h Maths.cc \
StandardSelectors.cc StandardSelectors.h\
Histogram.cc Histogram.fh Histogram.h \
GaussianIntegrator.cc GaussianIntegrator.h \
GaussianIntegrator.tcc \
Statistic.h HerwigStrategy.cc HerwigStrategy.h \
GSLIntegrator.h GSLIntegrator.tcc \
GSLBisection.h GSLBisection.tcc GSLHelper.h \
expm-1.h
nodist_libHwUtils_la_SOURCES = hgstamp.inc
BUILT_SOURCES = hgstamp.inc
CLEANFILES = hgstamp.inc
HGVERSION := $(shell hg -R $(top_srcdir) parents --template '"Herwig {node|short} ({branch})"' 2> /dev/null || echo \"$(PACKAGE_STRING)\" || true )
.PHONY: update_hgstamp
hgstamp.inc: update_hgstamp
@[ -f $@ ] || touch $@
@echo '$(HGVERSION)' | cmp -s $@ - || echo '$(HGVERSION)' > $@
libHwUtils_la_LIBADD = \
XML/libHwXML.la \
Statistics/libHwStatistics.la
check_PROGRAMS = utilities_test
utilities_test_SOURCES = \
tests/utilitiesTestsMain.cc \
tests/utilitiesTestsGlobalFixture.h \
tests/utilitiesTestsKinematics.h \
tests/utilitiesTestMaths.h \
tests/utilitiesTestsStatistic.h
utilities_test_LDADD = $(BOOST_UNIT_TEST_FRAMEWORK_LIBS) $(THEPEGLIB) -ldl libHwUtils.la
utilities_test_LDFLAGS = $(AM_LDFLAGS) -export-dynamic $(BOOST_UNIT_TEST_FRAMEWORK_LDFLAGS) $(THEPEGLDFLAGS)
utilities_test_CPPFLAGS = $(AM_CPPFLAGS) $(BOOST_CPPFLAGS)
TESTS = utilities_test
diff --git a/Utilities/Progress.cc b/Utilities/Progress.cc
new file mode 100644
--- /dev/null
+++ b/Utilities/Progress.cc
@@ -0,0 +1,90 @@
+// -*- C++ -*-
+//
+// Progress.h is a part of Herwig - A multi-purpose Monte Carlo event generator
+//
+
+// boost progress.hpp header file ------------------------------------------//
+
+// Copyright Beman Dawes 1994-99. Distributed under the Boost
+// Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// See http://www.boost.org/libs/timer for documentation.
+
+// Revision History
+
+// 2017-10-14 Modified for Herwig (D. Grellscheid)
+
+// 1 Dec 01 Add leading progress display strings (suggested by Toon Knapen)
+// 20 May 01 Introduce several static_casts<> to eliminate warning messages
+// (Fixed by Beman, reported by Herve Bronnimann)
+// 12 Jan 01 Change to inline implementation to allow use without library
+// builds. See docs for more rationale. (Beman Dawes)
+// 22 Jul 99 Name changed to .hpp
+// 16 Jul 99 Second beta
+// 6 Jul 99 Initial boost version
+
+#include "Progress.h"
+#include <iostream> // for ostream, cout, etc
+
+namespace Herwig {
+
+// progress_display --------------------------------------------------------//
+
+// progress_display displays an appropriate indication of
+// progress at an appropriate place in an appropriate form.
+
+// NOTE: (Jan 12, 2001) Tried to change unsigned long to boost::uintmax_t, but
+// found some compilers couldn't handle the required conversion to double.
+// Reverted to unsigned long until the compilers catch up.
+
+progress_display::progress_display( unsigned long expected_count,
+ std::ostream & os,
+ const std::string & s1, //leading strings
+ const std::string & s2,
+ const std::string & s3 )
+ // os is hint; implementation may ignore, particularly in embedded systems
+ : m_os(os), m_s1(s1), m_s2(s2), m_s3(s3) { restart(expected_count); }
+
+
+ void progress_display::restart( unsigned long expected_count )
+ // Effects: display appropriate scale
+ // Postconditions: count()==0, expected_count()==expected_count
+ {
+ _count = _next_tic_count = _tic = 0;
+ _expected_count = expected_count;
+
+ m_os << m_s1 << "0% 10 20 30 40 50 60 70 80 90 100%\n"
+ << m_s2 << "|----|----|----|----|----|----|----|----|----|----|"
+ << std::endl // endl implies flush, which ensures display
+ << m_s3;
+ if ( !_expected_count ) _expected_count = 1; // prevent divide by zero
+ } // restart
+
+ unsigned long progress_display::operator+=( unsigned long increment )
+ // Effects: Display appropriate progress tic if needed.
+ // Postconditions: count()== original count() + increment
+ // Returns: count().
+ {
+ if ( (_count += increment) >= _next_tic_count ) { display_tic(); }
+ return _count;
+ }
+
+ void progress_display::display_tic()
+ {
+ // use of floating point ensures that both large and small counts
+ // work correctly. static_cast<>() is also used several places
+ // to suppress spurious compiler warnings.
+ unsigned int tics_needed =
+ static_cast<unsigned int>(
+ (static_cast<double>(_count)/_expected_count)*50.0 );
+ do { m_os << '*' << std::flush; } while ( ++_tic < tics_needed );
+ _next_tic_count =
+ static_cast<unsigned long>((_tic/50.0)*_expected_count);
+ if ( _count == _expected_count ) {
+ if ( _tic < 51 ) m_os << '*';
+ m_os << std::endl;
+ }
+ } // display_tic
+
+} // namespace Herwig
diff --git a/Utilities/Progress.h b/Utilities/Progress.h
new file mode 100644
--- /dev/null
+++ b/Utilities/Progress.h
@@ -0,0 +1,75 @@
+// -*- C++ -*-
+//
+// Progress.h is a part of Herwig - A multi-purpose Monte Carlo event generator
+//
+
+// boost progress.hpp header file ------------------------------------------//
+
+// Copyright Beman Dawes 1994-99. Distributed under the Boost
+// Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// See http://www.boost.org/libs/timer for documentation.
+
+// Revision History
+
+// 2017-10-14 Modified for Herwig (D. Grellscheid)
+
+// 1 Dec 01 Add leading progress display strings (suggested by Toon Knapen)
+// 20 May 01 Introduce several static_casts<> to eliminate warning messages
+// (Fixed by Beman, reported by Herve Bronnimann)
+// 12 Jan 01 Change to inline implementation to allow use without library
+// builds. See docs for more rationale. (Beman Dawes)
+// 22 Jul 99 Name changed to .hpp
+// 16 Jul 99 Second beta
+// 6 Jul 99 Initial boost version
+
+#ifndef HERWIG_PROGRESS_H
+#define HERWIG_PROGRESS_H
+
+#include <iosfwd> // for ostream, cout, etc
+#include <string> // for string
+
+namespace Herwig {
+
+// progress_display --------------------------------------------------------//
+
+// progress_display displays an appropriate indication of
+// progress at an appropriate place in an appropriate form.
+
+// NOTE: (Jan 12, 2001) Tried to change unsigned long to boost::uintmax_t, but
+// found some compilers couldn't handle the required conversion to double.
+// Reverted to unsigned long until the compilers catch up.
+
+class progress_display
+{
+ public:
+ explicit progress_display( unsigned long expected_count,
+ std::ostream & os,
+ const std::string & s1 = "\n", //leading strings
+ const std::string & s2 = "",
+ const std::string & s3 = "" );
+ progress_display( const progress_display& ) = delete;
+ progress_display& operator=( const progress_display& ) = delete;
+ void restart( unsigned long expected_count );
+
+ unsigned long operator+=( unsigned long increment );
+
+ unsigned long operator++() { return operator+=( 1 ); }
+ unsigned long count() const { return _count; }
+ unsigned long expected_count() const { return _expected_count; }
+
+ private:
+ std::ostream & m_os; // may not be present in all imps
+ const std::string m_s1; // string is more general, safer than
+ const std::string m_s2; // const char *, and efficiency or size are
+ const std::string m_s3; // not issues
+
+ unsigned long _count, _expected_count, _next_tic_count;
+ unsigned int _tic;
+ void display_tic();
+};
+
+} // namespace Herwig
+
+#endif // HERWIG_PROGRESS_H
diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -1,243 +1,242 @@
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
dnl === disable debug symbols by default =====
if test "x$CXXFLAGS" = "x"; then
CXXFLAGS=-O2
fi
if test "x$CFLAGS" = "x"; then
CFLAGS=-O2
fi
AC_LANG([C++])
AM_INIT_AUTOMAKE([1.11 subdir-objects gnu dist-bzip2 no-dist-gzip -Wall -Wno-portability])
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
AX_CXX_COMPILE_STDCXX([11],[noext],[mandatory])
dnl check for POSIX
AC_CHECK_HEADER([unistd.h],[],
[AC_MSG_ERROR([Herwig needs "unistd.h". Non-POSIX systems are not supported.])])
AC_CHECK_HEADER([sys/stat.h],[],
[AC_MSG_ERROR([Herwig needs "sys/stat.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/numeric/ublas/io.hpp])
dnl Boost 1.64 is missing a required header to make these work
dnl we just assume they're there if io.hpp has been found OK above
dnl BOOST_FIND_HEADER([boost/numeric/ublas/matrix.hpp])
dnl BOOST_FIND_HEADER([boost/numeric/ublas/matrix_proxy.hpp])
dnl BOOST_FIND_HEADER([boost/numeric/ublas/matrix_sparse.hpp])
dnl BOOST_FIND_HEADER([boost/numeric/ublas/symmetric.hpp])
dnl BOOST_FIND_HEADER([boost/numeric/ublas/vector.hpp])
BOOST_FIND_HEADER([boost/operators.hpp])
-BOOST_FIND_HEADER([boost/progress.hpp])
BOOST_TEST()
HERWIG_CHECK_VBFNLO
HERWIG_CHECK_NJET
HERWIG_CHECK_GOSAM
HERWIG_CHECK_GOSAM_CONTRIB
HERWIG_CHECK_OPENLOOPS
HERWIG_CHECK_MADGRAPH
HERWIG_CHECK_EVTGEN
HERWIG_CHECK_PYTHIA
HERWIG_COMPILERFLAGS
HERWIG_LOOPTOOLS
FASTJET_CHECK_FASTJET
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
Decay/EvtGen/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/Reweighters/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/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/VBFNLO/Makefile
MatrixElement/Matchbox/External/NJet/Makefile
MatrixElement/Matchbox/External/GoSam/Makefile
MatrixElement/Matchbox/External/OpenLoops/Makefile
MatrixElement/Matchbox/External/MadGraph/Makefile
MatrixElement/Matchbox/External/MadGraph/mg2herwig.py
Sampling/Makefile
Sampling/CellGrids/Makefile
Shower/Makefile
Shower/QTilde/Makefile
Shower/QTilde/Matching/Makefile
Shower/Dipole/Makefile
Shower/Dipole/Base/Makefile
Shower/Dipole/Kernels/Makefile
Shower/Dipole/Kinematics/Makefile
Shower/Dipole/Utility/Makefile
Shower/Dipole/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
API/Makefile
src/Makefile-UserModules
src/defaults/Analysis.in
src/defaults/MatchboxDefaults.in
src/defaults/Decays.in
src/defaults/decayers.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
src/Merging/Makefile
Shower/Dipole/Merging/Makefile
src/defaults/MatchboxMergingDefaults.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
Tue, Nov 19, 9:02 PM (1 d, 2 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806195
Default Alt Text
(239 KB)

Event Timeline