file ColliderBit/getPy8Collider.hpp
[No description available] More…
Namespaces
| Name |
|---|
| Gambit TODO: see if we can use this one: |
| Gambit::ColliderBit |
Defines
| Name | |
|---|---|
| DEBUG_PREFIX | |
| GET_SPECTRUM_AND_DECAYS_FOR_PYTHIA_SUSY(NAME, SPECTRUM) | |
| GET_SPECTRUM_AND_DECAYS_FOR_PYTHIA_NONSUSY(NAME, SPECTRUM) | |
| HEPMC_TYPE(PYTHIA_NS) Work out last template arg of Py8Collider depending on whether we are using HepMC. | |
| GET_SPECIFIC_PYTHIA(NAME, PYTHIA_NS, MODEL_EXTENSION) Retrieve a specific Pythia hard-scattering Monte Carlo simulation. | |
| GET_SPECIFIC_PYTHIA_SLHA(NAME, PYTHIA_NS, MODEL_EXTENSION) | |
| GET_PYTHIA_AS_BASE_COLLIDER(NAME) Get a specific Pythia hard-scattering sim as a generator-independent pointer-to-BaseCollider. |
Detailed Description
Author:
- Abram Krislock (a.m.b.krislock@fys.uio.no)
- Aldo Saavedra
- Andy Buckley
- Chris Rogan (crogan@cern.ch)
- Pat Scott (p.scott@imperial.ac.uk)
- Anders Kvellestad (anders.kvellestad@fys.uio.no)
Date:
- 2014 Aug
- 2015 May
- 2015 Jul
- 2018 Jan
- 2019 Jan, Feb, May
- 2017 March
- 2018 Jan
- 2018 May
ColliderBit event loop functions returning collider Monte Carlo event simulators.
Authors (add name and date if you modify):
Macros Documentation
define DEBUG_PREFIX
#define DEBUG_PREFIX "DEBUG: OMP thread " << omp_get_thread_num() << ": "
define GET_SPECTRUM_AND_DECAYS_FOR_PYTHIA_SUSY
#define GET_SPECTRUM_AND_DECAYS_FOR_PYTHIA_SUSY(
NAME,
SPECTRUM
)
define GET_SPECTRUM_AND_DECAYS_FOR_PYTHIA_NONSUSY
#define GET_SPECTRUM_AND_DECAYS_FOR_PYTHIA_NONSUSY(
NAME,
SPECTRUM
)
void NAME(SLHAstruct& result) \
{ \
using namespace Pipes::NAME; \
result.clear(); \
/* Get decays */ \
result = Dep::decay_rates->getSLHAea(2); \
/* Get spectrum */ \
SLHAstruct slha_spectrum = Dep::SPECTRUM->getSLHAea(2); \
result.insert(result.begin(), slha_spectrum.begin(), slha_spectrum.end()); \
}
define HEPMC_TYPE
#define HEPMC_TYPE(
PYTHIA_NS
)
PYTHIA_NS::Pythia8::GAMBIT_hepmc_writer
Work out last template arg of Py8Collider depending on whether we are using HepMC.
define GET_SPECIFIC_PYTHIA
#define GET_SPECIFIC_PYTHIA(
NAME,
PYTHIA_NS,
MODEL_EXTENSION
)
void NAME(Py8Collider<PYTHIA_NS::Pythia8::Pythia, \
PYTHIA_NS::Pythia8::Event, \
HEPMC_TYPE(PYTHIA_NS)> &result) \
{ \
using namespace Pipes::NAME; \
static SLHAstruct slha; \
\
if (*Loop::iteration == BASE_INIT) \
{ \
/* SLHAea object constructed from dependencies on the spectrum and decays. */ \
slha.clear(); \
slha = *Dep::SpectrumAndDecaysForPythia; \
} \
\
getPy8Collider(result, *Dep::RunMC, slha, #MODEL_EXTENSION, \
*Loop::iteration, Loop::wrapup, *runOptions); \
}
Retrieve a specific Pythia hard-scattering Monte Carlo simulation.
define GET_SPECIFIC_PYTHIA_SLHA
#define GET_SPECIFIC_PYTHIA_SLHA(
NAME,
PYTHIA_NS,
MODEL_EXTENSION
)
void NAME(Py8Collider<PYTHIA_NS::Pythia8::Pythia, \
PYTHIA_NS::Pythia8::Event, \
HEPMC_TYPE(PYTHIA_NS)> &result) \
{ \
using namespace Pipes::NAME; \
static SLHAstruct slha; \
\
if (*Loop::iteration == COLLIDER_INIT) \
{ \
const pair_str_SLHAstruct& filename_content_pair = *Dep::SLHAFileNameAndContent; \
if (filename_content_pair.first.empty()) \
{ \
piped_invalid_point.request("Got empty SLHA filename. Will invalidate point."); \
} \
} \
\
getPy8Collider(result, *Dep::RunMC, Dep::SLHAFileNameAndContent->second, \
#MODEL_EXTENSION, *Loop::iteration, Loop::wrapup, *runOptions); \
}
Retrieve a specific Pythia hard-scattering Monte Carlo simulation from reading a SLHA file rather than getting a Spectrum + DecayTable
define GET_PYTHIA_AS_BASE_COLLIDER
#define GET_PYTHIA_AS_BASE_COLLIDER(
NAME
)
void NAME(const BaseCollider* &result) \
{ \
result = &(*Pipes::NAME::Dep::HardScatteringSim); \
} \
Get a specific Pythia hard-scattering sim as a generator-independent pointer-to-BaseCollider.
Source code
// GAMBIT: Global and Modular BSM Inference Tool
// *********************************************
/// \file
///
/// ColliderBit event loop functions returning
/// collider Monte Carlo event simulators.
///
/// *********************************************
///
/// Authors (add name and date if you modify):
///
/// \author Abram Krislock
/// (a.m.b.krislock@fys.uio.no)
///
/// \author Aldo Saavedra
///
/// \author Andy Buckley
///
/// \author Chris Rogan
/// (crogan@cern.ch)
/// \date 2014 Aug
/// \date 2015 May
///
/// \author Pat Scott
/// (p.scott@imperial.ac.uk)
/// \date 2015 Jul
/// \date 2018 Jan
/// \date 2019 Jan, Feb, May
///
/// \author Anders Kvellestad
/// (anders.kvellestad@fys.uio.no)
/// \date 2017 March
/// \date 2018 Jan
/// \date 2018 May
///
/// *********************************************
#include "gambit/ColliderBit/ColliderBit_eventloop.hpp"
// #define COLLIDERBIT_DEBUG
#define DEBUG_PREFIX "DEBUG: OMP thread " << omp_get_thread_num() << ": "
namespace Gambit
{
namespace ColliderBit
{
/// Retrieve a Pythia hard-scattering Monte Carlo simulation
template<typename PythiaT, typename EventT, typename hepmc_writerT>
void getPy8Collider(Py8Collider<PythiaT, EventT, hepmc_writerT>& result,
const MCLoopInfo& RunMC,
const SLHAstruct& slha,
const str model_suffix,
const int iteration,
void(*wrapup)(),
const Options& runOptions)
{
static bool first = true;
static str pythia_doc_path;
static double xsec_veto_fb;
if (iteration == BASE_INIT)
{
// Setup the Pythia documentation path and print the banner once
if (first)
{
const str be = "Pythia" + model_suffix;
const str ver = Backends::backendInfo().default_version(be);
pythia_doc_path = Backends::backendInfo().path_dir(be, ver) + "/../share/Pythia8/xmldoc/";
result.banner(pythia_doc_path);
first = false;
}
}
// To make sure that the Pythia instance on each OMP thread gets all the
// options it should, all the options parsing and initialisation happens in
// COLLIDER_INIT_OMP (OMP parallel) rather than COLLIDER_INIT (only thread 0).
// We may want to split this up, so that all the yaml options are parsed in
// COLLIDER_INIT (by thread 0), and used to initialize the 'result' instance
// of each thread within COLLIDER_INIT_OMP.
//
// else if (iteration == COLLIDER_INIT)
// {
// // Do the option parsing here?
// }
else if (iteration == COLLIDER_INIT_OMP)
{
std::vector<str> pythiaOptions;
// By default we tell Pythia to be quiet. (Can be overridden from yaml settings)
pythiaOptions.push_back("Print:quiet = on");
pythiaOptions.push_back("SLHA:verbose = 0");
// Get options from yaml file.
const double xsec_veto_default = 0.0;
const bool partonOnly_default = false;
if (runOptions.hasKey(RunMC.current_collider()))
{
YAML::Node colNode = runOptions.getValue<YAML::Node>(RunMC.current_collider());
Options colOptions(colNode);
xsec_veto_fb = colOptions.getValueOrDef<double>(xsec_veto_default, "xsec_veto");
result.partonOnly = colOptions.getValueOrDef<bool>(partonOnly_default, "partonOnly");
// Fill the jet collection settings.
result.all_jet_collection_settings.clear();
result.jetcollection_taus = "";
if (colOptions.hasKey("jet_collections"))
{
YAML::Node all_jetcollections_node = colOptions.getValue<YAML::Node>("jet_collections");
Options all_jetcollection_options(all_jetcollections_node);
std::vector<str> jetcollection_names = all_jetcollection_options.getNames();
for (str key : jetcollection_names)
{
YAML::Node current_jc_node = all_jetcollection_options.getValue<YAML::Node>(key);
Options current_jc_options(current_jc_node);
str algorithm = current_jc_options.getValue<str>("algorithm");
double R = current_jc_options.getValue<double>("R");
str recombination_scheme = current_jc_options.getValue<str>("recombination_scheme");
str strategy = current_jc_options.getValue<str>("strategy");
(result.all_jet_collection_settings).push_back({key, algorithm, R, recombination_scheme, strategy});
}
result.jetcollection_taus = colOptions.getValue<str>("jet_collection_taus");
// Throw an error if the "jet_collection_taus" setting does not match an entry in "jet_collections".
if (std::find(jetcollection_names.begin(), jetcollection_names.end(), result.jetcollection_taus) == jetcollection_names.end())
{
ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections.");
}
}
else
{
ColliderBit_error().raise(LOCAL_INFO,"Could not find jet_collections option for collider " + RunMC.current_collider() + ". Please provide this in the YAML file.");
}
if (colOptions.hasKey("pythia_settings"))
{
std::vector<str> addPythiaOptions = colNode["pythia_settings"].as<std::vector<str> >();
pythiaOptions.insert(pythiaOptions.end(), addPythiaOptions.begin(), addPythiaOptions.end());
}
}
else
{
ColliderBit_error().raise(LOCAL_INFO,"Could not find runOptions for collider " + RunMC.current_collider() + ".");
}
// We need showProcesses for the xsec veto.
pythiaOptions.push_back("Init:showProcesses = on");
// We need "SLHA:file = slhaea" for the SLHAea interface.
pythiaOptions.push_back("SLHA:file = slhaea");
// If the collider energy is not given in the list of Pythia options, we set it to 13 TeV by default.
// We search for the substring "Beams:e", meaning that if any of the Pythia options "Beams:eCM", "Beams:eA"
// or "Beams:eB" are present we don't apply the default.
bool has_beam_energy_option = std::any_of(pythiaOptions.begin(), pythiaOptions.end(), [](const str& s){ return s.find("Beams:e") != str::npos; });
if (!has_beam_energy_option)
{
pythiaOptions.push_back("Beams:eCM = 13000");
logger() << LogTags::debug << "Could not find a beam energy in the list of Pythia settings. Will add the setting 'Beams:eCM = 13000'." << EOM;
}
// Variables needed for the xsec veto
std::stringstream processLevelOutput;
str _junk, readline;
int code, nxsec;
double xsec, totalxsec;
// Each thread needs an independent Pythia instance at the start
// of each event generation loop.
// Thus, the actual Pythia initialization is
// *after* COLLIDER_INIT, within omp parallel.
result.clear();
// Add the thread-specific seed to the Pythia options
str seed = std::to_string(int(Random::draw() * 899990000.));
pythiaOptions.push_back("Random:seed = " + seed);
#ifdef COLLIDERBIT_DEBUG
cout << DEBUG_PREFIX << "getPythia"+model_suffix+": My Pythia seed is: " << seed << endl;
#endif
try
{
result.init(pythia_doc_path, pythiaOptions, &slha, processLevelOutput);
}
catch (typename Py8Collider<PythiaT,EventT,hepmc_writerT>::InitializationError& e)
{
// Append new seed to override the previous one
int newSeedBase = int(Random::draw() * 899990000.);
pythiaOptions.push_back("Random:seed = " + std::to_string(newSeedBase));
try
{
result.init(pythia_doc_path, pythiaOptions, &slha, processLevelOutput);
}
catch (typename Py8Collider<PythiaT,EventT,hepmc_writerT>::InitializationError& e)
{
#ifdef COLLIDERBIT_DEBUG
cout << DEBUG_PREFIX << "Py8Collider::InitializationError caught in getPy8Collider. Will discard this point." << endl;
#endif
piped_invalid_point.request("Bad point: Pythia can't initialize");
wrapup();
return;
}
}
// Should we apply the xsec veto and skip event generation?
// - Get the upper limt xsec as estimated by Pythia
code = -1;
nxsec = 0;
totalxsec = 0.;
while(true)
{
std::getline(processLevelOutput, readline);
std::istringstream issPtr(readline);
issPtr.seekg(47, issPtr.beg);
issPtr >> code;
if (!issPtr.good() && nxsec > 0) break;
issPtr >> _junk >> xsec;
if (issPtr.good())
{
totalxsec += xsec;
nxsec++;
}
}
#ifdef COLLIDERBIT_DEBUG
cout << DEBUG_PREFIX << "totalxsec [fb] = " << totalxsec * 1e12 << ", veto limit [fb] = " << xsec_veto_fb << endl;
#endif
// - Check for NaN xsec
if (Utils::isnan(totalxsec))
{
#ifdef COLLIDERBIT_DEBUG
cout << DEBUG_PREFIX << "Got NaN cross-section estimate from Pythia." << endl;
#endif
piped_invalid_point.request("Got NaN cross-section estimate from Pythia.");
wrapup();
return;
}
// - Wrap up loop if veto applies
if (totalxsec * 1e12 < xsec_veto_fb)
{
#ifdef COLLIDERBIT_DEBUG
cout << DEBUG_PREFIX << "Cross-section veto applies. Will now call Loop::wrapup() to skip event generation for this collider." << endl;
#endif
wrapup();
} else {
// Create a dummy event to make Pythia fill its internal list of process codes
EventT dummy_pythia_event;
try
{
result.nextEvent(dummy_pythia_event);
}
catch (typename Py8Collider<PythiaT,EventT,hepmc_writerT>::EventGenerationError& e)
{
// Try again...
try
{
result.nextEvent(dummy_pythia_event);
}
catch (typename Py8Collider<PythiaT,EventT,hepmc_writerT>::EventGenerationError& e)
{
piped_invalid_point.request("Failed to generate dummy test event. Will invalidate point.");
#ifdef COLLIDERBIT_DEBUG
cout << DEBUG_PREFIX << "Failed to generate dummy test event during COLLIDER_INIT_OMP in getPy8Collider. Check the logs for event details." << endl;
#endif
#pragma omp critical (pythia_event_failure)
{
std::stringstream ss;
dummy_pythia_event.list(ss, 1);
logger() << LogTags::debug << "Failed to generate dummy test event during COLLIDER_INIT_OMP iteration in getPy8Collider. Pythia record for the event that failed:\n" << ss.str() << EOM;
}
}
}
}
}
}
// Get SLHAea object with spectrum and decays for Pythia -- SUSY version
#define GET_SPECTRUM_AND_DECAYS_FOR_PYTHIA_SUSY(NAME, SPECTRUM) \
void NAME(SLHAstruct& result) \
{ \
using namespace Pipes::NAME; \
static const int slha_version = runOptions->getValueOrDef<int>(2, "slha_version"); \
static const bool write_summary_to_log = \
runOptions->getValueOrDef<bool>(false, "write_summary_to_log"); \
\
if ((slha_version != 1) && (slha_version != 2)) \
{ \
ColliderBit_error().raise(LOCAL_INFO, \
"The option 'slha_version' must be set to 1 or 2 (default)."); \
} \
result.clear(); \
/* Get decays */ \
result = Dep::decay_rates->getSLHAea(slha_version, false, *Dep::SLHA_pseudonyms); \
/* Get spectrum */ \
SLHAstruct slha_spectrum = Dep::SPECTRUM->getSLHAea(slha_version); \
result.insert(result.begin(), slha_spectrum.begin(), slha_spectrum.end()); \
/* Add MODSEL block if not found */ \
if(result.find("MODSEL") == result.end()) \
{ \
SLHAea::Block block("MODSEL"); \
block.push_back("BLOCK MODSEL # Model selection"); \
SLHAea::Line line; \
line << 1 << 0 << "# Tell Pythia that this is a SUSY model."; \
block.push_back(line); \
result.push_front(block); \
} \
\
if (write_summary_to_log) \
{ \
std::stringstream SLHA_log_output; \
SLHA_log_output << "SLHA" << slha_version << " input to Pythia:\n" << result.str() \
<< "\n"; \
logger() << SLHA_log_output.str() << EOM; \
} \
}
// Get SLHAea object with spectrum and decays for Pythia -- non-SUSY version
#define GET_SPECTRUM_AND_DECAYS_FOR_PYTHIA_NONSUSY(NAME, SPECTRUM) \
void NAME(SLHAstruct& result) \
{ \
using namespace Pipes::NAME; \
result.clear(); \
/* Get decays */ \
result = Dep::decay_rates->getSLHAea(2); \
/* Get spectrum */ \
SLHAstruct slha_spectrum = Dep::SPECTRUM->getSLHAea(2); \
result.insert(result.begin(), slha_spectrum.begin(), slha_spectrum.end()); \
}
/// Work out last template arg of Py8Collider depending on whether we are using HepMC
#ifdef EXCLUDE_HEPMC
#define HEPMC_TYPE(PYTHIA_NS) void
#else
#define HEPMC_TYPE(PYTHIA_NS) PYTHIA_NS::Pythia8::GAMBIT_hepmc_writer
#endif
/// Retrieve a specific Pythia hard-scattering Monte Carlo simulation
#define GET_SPECIFIC_PYTHIA(NAME, PYTHIA_NS, MODEL_EXTENSION) \
void NAME(Py8Collider<PYTHIA_NS::Pythia8::Pythia, \
PYTHIA_NS::Pythia8::Event, \
HEPMC_TYPE(PYTHIA_NS)> &result) \
{ \
using namespace Pipes::NAME; \
static SLHAstruct slha; \
\
if (*Loop::iteration == BASE_INIT) \
{ \
/* SLHAea object constructed from dependencies on the spectrum and decays. */ \
slha.clear(); \
slha = *Dep::SpectrumAndDecaysForPythia; \
} \
\
getPy8Collider(result, *Dep::RunMC, slha, #MODEL_EXTENSION, \
*Loop::iteration, Loop::wrapup, *runOptions); \
}
/// Retrieve a specific Pythia hard-scattering Monte Carlo simulation
/// from reading a SLHA file rather than getting a Spectrum + DecayTable
#define GET_SPECIFIC_PYTHIA_SLHA(NAME, PYTHIA_NS, MODEL_EXTENSION) \
void NAME(Py8Collider<PYTHIA_NS::Pythia8::Pythia, \
PYTHIA_NS::Pythia8::Event, \
HEPMC_TYPE(PYTHIA_NS)> &result) \
{ \
using namespace Pipes::NAME; \
static SLHAstruct slha; \
\
if (*Loop::iteration == COLLIDER_INIT) \
{ \
const pair_str_SLHAstruct& filename_content_pair = *Dep::SLHAFileNameAndContent; \
if (filename_content_pair.first.empty()) \
{ \
piped_invalid_point.request("Got empty SLHA filename. Will invalidate point."); \
} \
} \
\
getPy8Collider(result, *Dep::RunMC, Dep::SLHAFileNameAndContent->second, \
#MODEL_EXTENSION, *Loop::iteration, Loop::wrapup, *runOptions); \
}
/// Get a specific Pythia hard-scattering sim as a generator-independent pointer-to-BaseCollider
#define GET_PYTHIA_AS_BASE_COLLIDER(NAME) \
void NAME(const BaseCollider* &result) \
{ \
result = &(*Pipes::NAME::Dep::HardScatteringSim); \
} \
}
}
Updated on 2024-07-18 at 13:53:34 +0000