gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
getColliderPythia.hpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
37 
38 
40 
41  #define COLLIDERBIT_DEBUG
42 
43 namespace Gambit
44 {
45 
46  namespace ColliderBit
47  {
48 
50  template<typename PythiaT, typename EventT>
52  const MCLoopInfo& RunMC,
53  const Spectrum& spectrum,
54  const DecayTable& decay_rates,
55  const str model_suffix,
56  const int iteration,
57  void(*wrapup)(),
58  const Options& runOptions,
59  bool is_SUSY)
60  {
61  static bool first = true;
62  static std::vector<str> filenames;
63  static str pythia_doc_path;
64  static std::vector<str> pythiaCommonOptions;
65  static SLHAstruct slha;
66  static SLHAstruct slha_spectrum;
67  static double xsec_veto_fb;
68  static unsigned int fileCounter = 0;
69 
70  if (iteration == BASE_INIT)
71  {
72  // Setup the Pythia documentation path and print the banner once
73  if (first)
74  {
75  const str be = "Pythia" + model_suffix;
76  const str ver = Backends::backendInfo().default_version(be);
77  pythia_doc_path = Backends::backendInfo().path_dir(be, ver) + "/../share/Pythia8/xmldoc/";
78  result.banner(pythia_doc_path);
79  if (runOptions.hasKey("SLHA_filenames"))
80  {
81  filenames = runOptions.getValue<std::vector<str> >("SLHA_filenames");
82  }
83  first = false;
84  }
85 
86  if (filenames.empty())
87  {
88  // SLHAea object constructed from dependencies on the spectrum and decays.
89  slha.clear();
90  slha_spectrum.clear();
91  slha = decay_rates.getSLHAea(2);
92  // SLHAea in SLHA2 format, please.
93  slha_spectrum = spectrum.getSLHAea(2);
94  slha.insert(slha.begin(), slha_spectrum.begin(), slha_spectrum.end());
95  if (is_SUSY)
96  {
97  SLHAea::Block block("MODSEL");
98  block.push_back("BLOCK MODSEL # Model selection");
99  SLHAea::Line line;
100  line << 1 << 0 << "# Tell Pythia that this is a SUSY model.";
101  block.push_back(line);
102  slha.push_front(block);
103  }
104  }
105  else
106  {
107  if (filenames.size() <= fileCounter) invalid_point().raise("No more SLHA files. My work is done.");
108  }
109 
110  }
111 
112  else if (iteration == COLLIDER_INIT)
113  {
114  // Collect Pythia options that are common across all OMP threads
115  pythiaCommonOptions.clear();
116 
117  // By default we tell Pythia to be quiet. (Can be overridden from yaml settings)
118  pythiaCommonOptions.push_back("Print:quiet = on");
119  pythiaCommonOptions.push_back("SLHA:verbose = 0");
120 
121  // Get options from yaml file.
122  double xsec_veto_default = 0.0;
123  if (runOptions.hasKey(RunMC.current_collider()))
124  {
125  YAML::Node colNode = runOptions.getValue<YAML::Node>(RunMC.current_collider());
126  Options colOptions(colNode);
127  xsec_veto_fb = colOptions.getValueOrDef<double>(xsec_veto_default, "xsec_veto");
128 
129  if (colOptions.hasKey("pythia_settings"))
130  {
131  std::vector<str> addPythiaOptions = colNode["pythia_settings"].as<std::vector<str> >();
132  pythiaCommonOptions.insert(pythiaCommonOptions.end(), addPythiaOptions.begin(), addPythiaOptions.end());
133  }
134  }
135  else xsec_veto_fb = xsec_veto_default;
136 
137  // We need showProcesses for the xsec veto.
138  pythiaCommonOptions.push_back("Init:showProcesses = on");
139 
140  // We need "SLHA:file = slhaea" for the SLHAea interface, and the filename for the SLHA interface.
141  str slha_string = (filenames.empty() ? "slhaea" : filenames.at(fileCounter));
142  pythiaCommonOptions.push_back("SLHA:file = " + slha_string);
143  }
144 
145  else if (iteration == START_SUBPROCESS)
146  {
147  // Variables needed for the xsec veto
148  std::stringstream processLevelOutput;
149  str _junk, readline;
150  int code, nxsec;
151  double xsec, totalxsec;
152 
153  // Each thread needs an independent Pythia instance at the start
154  // of each event generation loop.
155  // Thus, the actual Pythia initialization is
156  // *after* COLLIDER_INIT, within omp parallel.
157 
158  result.clear();
159 
160  if (not filenames.empty() and omp_get_thread_num() == 0)
161  logger() << "Reading SLHA file: " << filenames.at(fileCounter) << EOM;
162 
163  // Get the Pythia options that are common across all OMP threads ('pythiaCommonOptions')
164  // and then add the thread-specific seed
165  std::vector<str> pythiaOptions = pythiaCommonOptions;
166  str seed = std::to_string(int(Random::draw() * 899990000.));
167  pythiaOptions.push_back("Random:seed = " + seed);
168 
169  #ifdef COLLIDERBIT_DEBUG
170  cout << debug_prefix() << "getPythia"+model_suffix+": My Pythia seed is: " << seed << endl;
171  #endif
172 
173  try
174  {
175  if (filenames.empty())
176  {
177  result.init(pythia_doc_path, pythiaOptions, &slha, processLevelOutput);
178  }
179  else
180  {
181  result.init(pythia_doc_path, pythiaOptions, processLevelOutput);
182  }
183  }
185  {
186  // Append new seed to override the previous one
187  int newSeedBase = int(Random::draw() * 899990000.);
188  pythiaOptions.push_back("Random:seed = " + std::to_string(newSeedBase));
189  try
190  {
191  if (filenames.empty())
192  {
193  result.init(pythia_doc_path, pythiaOptions, &slha, processLevelOutput);
194  }
195  else
196  {
197  result.init(pythia_doc_path, pythiaOptions, processLevelOutput);
198  }
199  }
201  {
202  #ifdef COLLIDERBIT_DEBUG
203  cout << debug_prefix() << "ColliderPythia::InitializationError caught in getColliderPythia. Will discard this point." << endl;
204  #endif
205  piped_invalid_point.request("Bad point: Pythia can't initialize");
206  wrapup();
207  return;
208  }
209  }
210 
211  // Should we apply the xsec veto and skip event generation?
212 
213  // - Get the upper limt xsec as estimated by Pythia
214  code = -1;
215  nxsec = 0;
216  totalxsec = 0.;
217  while(true)
218  {
219  std::getline(processLevelOutput, readline);
220  std::istringstream issPtr(readline);
221  issPtr.seekg(47, issPtr.beg);
222  issPtr >> code;
223  if (!issPtr.good() && nxsec > 0) break;
224  issPtr >> _junk >> xsec;
225  if (issPtr.good())
226  {
227  totalxsec += xsec;
228  nxsec++;
229  }
230  }
231 
232  #ifdef COLLIDERBIT_DEBUG
233  cout << debug_prefix() << "totalxsec [fb] = " << totalxsec * 1e12 << ", veto limit [fb] = " << xsec_veto_fb << endl;
234  #endif
235 
236  // - Check for NaN xsec
237  if (Utils::isnan(totalxsec))
238  {
239  #ifdef COLLIDERBIT_DEBUG
240  cout << debug_prefix() << "Got NaN cross-section estimate from Pythia." << endl;
241  #endif
242  piped_invalid_point.request("Got NaN cross-section estimate from Pythia.");
243  wrapup();
244  return;
245  }
246 
247  // - Wrap up loop if veto applies
248  if (totalxsec * 1e12 < xsec_veto_fb)
249  {
250  #ifdef COLLIDERBIT_DEBUG
251  cout << debug_prefix() << "Cross-section veto applies. Will now call Loop::wrapup() to skip event generation for this collider." << endl;
252  #endif
253  wrapup();
254  }
255 
256  }
257 
258  else if (iteration == BASE_FINALIZE and not filenames.empty()) fileCounter++;
259 
260  }
261 
263  #define IS_SUSY true
264  #define NOT_SUSY false
265  #define GET_SPECIFIC_PYTHIA(NAME, PYTHIA_NS, SPECTRUM, MODEL_EXTENSION, SUSY_FLAG)\
266  void NAME(ColliderPythia<PYTHIA_NS::Pythia8::Pythia, \
267  PYTHIA_NS::Pythia8::Event> &result) \
268  { \
269  using namespace Pipes::NAME; \
270  getColliderPythia(result, *Dep::RunMC, *Dep::SPECTRUM, \
271  *Dep::decay_rates, #MODEL_EXTENSION, *Loop::iteration, \
272  Loop::wrapup, *runOptions, SUSY_FLAG); \
273  }
274 
276  #define GET_PYTHIA_AS_BASE_COLLIDER(NAME) \
277  void NAME(const BaseCollider* &result) \
278  { \
279  result = &(*Pipes::NAME::Dep::HardScatteringSim); \
280  } \
281 
282  }
283 
284 }
void clear()
Reset this instance for reuse, avoiding the need for "new" or "delete".
A specializable, recyclable class interfacing ColliderBit and Pythia.
SLHAstruct getSLHAea(int) const
SLHAea object getter First constructs an SLHAea object from the SMINPUTS object, then adds the info f...
Definition: spectrum.cpp:409
void init()
Initialize with no settings (error): override version.
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:524
Declarations commont to all ColliderBit event loop functions.
bool hasKey(const args &... keys) const
Getters for key/value pairs (which is all the options node should contain)
TYPE getValue(const args &... keys) const
SLHAea::Coll SLHAstruct
Less confusing name for SLHAea container class.
static double draw()
Draw a single uniform random deviate from the interval (0,1) using the chosen RNG engine...
const str & current_collider() const
Get the current collider.
Definition: MCLoopInfo.cpp:95
virtual void raise(const std::string &)
Raise the exception, i.e. throw it.
Definition: exceptions.cpp:422
void getColliderPythia(ColliderPythia< PythiaT, EventT > &result, const MCLoopInfo &RunMC, const Spectrum &spectrum, const DecayTable &decay_rates, const str model_suffix, const int iteration, void(*wrapup)(), const Options &runOptions, bool is_SUSY)
Retrieve a Pythia hard-scattering Monte Carlo simulation.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
DS5_MSPCTM DS_INTDOF int
GAMBIT native decay table class.
Definition: decay_table.hpp:35
void request(std::string message)
Request an exception.
Definition: exceptions.cpp:463
void banner(const std::string pythiaDocPath)
Create a useless Pythia instance just to print the banner.
invalid_point_exception & invalid_point()
Invalid point exceptions.
str debug_prefix()
Debug prefix string giving thread number.
A class for holding cross-section info within ColliderBit.
Definition: xsec.hpp:28
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry decay_rates
SLHAstruct getSLHAea(int SLHA_version, bool include_zero_bfs=false, const mass_es_pseudonyms &psn=mass_es_pseudonyms()) const
Output entire decay table as an SLHAea file full of DECAY blocks.
An exception for when Pythia fails to initialize.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
A small wrapper object for &#39;options&#39; nodes.
Container for event loop status data and settings.
Definition: MCLoopInfo.hpp:31
"Standard Model" (low-energy) plus high-energy model container class
Definition: spectrum.hpp:110