PYTHIA  8.312
LHAHelaconia.h
1 // LHAHelaconia.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Author: Philip Ilten, December 2017.
7 
8 #ifndef Pythia8_LHAHelaconia_H
9 #define Pythia8_LHAHelaconia_H
10 
11 #include "Pythia8/Pythia.h"
12 #include <unistd.h>
13 #include <sys/stat.h>
14 #include <limits>
15 
16 namespace Pythia8 {
17 
18 //==========================================================================
19 
20 // A derived class from LHAup which generates events with HelacOnia.
21 
22 // This class automatically generates hard processes with HelacOnia
23 // and reads in the LHEF file output.
24 
25 // The user can send commands to HelacOnia via the readString
26 // method. The launch command, random seed, and shower choice are
27 // automatically handled. For example, the following will produce
28 // J/psi events from 13 TeV proton proton collisions:
29 
30 // readString("generate u u~ > cc~(3S11) g")
31 
32 // The number of events generated per HelacOnia run is controlled by
33 // setEvents, while the random seed is controlled by setSeed.
34 
35 class LHAupHelaconia : public LHAup {
36 
37 public:
38 
39  // Constructor.
40  LHAupHelaconia(Pythia *pythiaIn, string dirIn = "helaconiarun",
41  string exeIn = "ho_cluster");
42 
43  // Destructor: Print error statistics before exiting. Printing code
44  // basically copied from Info class.
46 
47  // Read a HelacOnia command string.
48  bool readString(string line);
49 
50  // Set the number of events to generate per run.
51  void setEvents(int eventsIn);
52 
53  // Set the random seed and maximum runs.
54  bool setSeed(int seedIn, int runsIn = 30081);
55 
56  // Set the initialization information.
57  bool setInit();
58 
59  // Set the event information.
60  bool setEvent(int = 0);
61 
62  // Note: The functions below have been made public to ease the generation
63  // of Python bindings.
64  //protected:
65 
66  // Execute a system command.
67  bool execute(string line);
68 
69  // Run HelacOnia.
70  bool run(int eventsIn, int seedIn = -1);
71 
72  // Create the LHEF reader.
73  bool reader(bool init);
74 
75  // Convert the color octet HelacOnia ID to a Pythia 8 ID.
76  int convert(int idIn);
77 
78  // Print a message the first few times. Insert in database.
79  void errorMsg(string messageIn) {
80  // Recover number of times message occured. Also inserts new string.
81  int times = messages[messageIn];
82  ++messages[messageIn];
83  // Print message the first few times.
84  if (times < TIMESTOPRINT) cout << " PYTHIA " << messageIn << endl;
85  }
86 
87 private:
88 
89  // The PYTHIA object and LHEF file reader and matching hook.
90  Pythia *pythia;
91  LHAupLHEF *lhef;
92 
93  // Stored members.
94  int events, seed, runs, nRuns, nId, nQ, nR, nL, nJ;
95  string dir, exe, lhegz;
96  double mQ;
97 
98  // The HelacOnia commands.
99  vector<string> lines;
100 
101  // Map for all error messages.
102  map<string, int> messages;
103  // Number of times the same error message is repeated, unless overridden.
104  static const int TIMESTOPRINT = 1;
105 
106 };
107 
108 //--------------------------------------------------------------------------
109 
110 // Constructor.
111 
112 LHAupHelaconia::LHAupHelaconia(Pythia *pythiaIn, string dirIn, string exeIn) :
113  pythia(pythiaIn), lhef(0), events(10000), seed(-1), runs(30081),
114  nRuns(0), nId(443), nQ(4), nR(0), nL(0), nJ(3),
115  dir(dirIn), exe(exeIn), lhegz(dirIn + "/events.lhe"), mQ(-2) {
116  mkdir(dir.c_str(), 0777);
117  if (pythia) pythia->readString("Beams:frameType = 5");
118  pythia->settings.addMode("Onia:state", -1, false, false, 0, 0);
119 }
120 
121 //--------------------------------------------------------------------------
122 
123 // Destructor.
124 
126 
127  if (lhef) delete lhef;
128 
129  // Header.
130  cout << "\n *------- LHAupHelaconia Error and Warning Messages Statistics"
131  << " --------------------------------------------------* \n"
132  << " | "
133  << " | \n"
134  << " | times message "
135  << " | \n"
136  << " | "
137  << " | \n";
138 
139  // Loop over all messages
140  map<string, int>::iterator messageEntry = messages.begin();
141  if (messageEntry == messages.end())
142  cout << " | 0 no errors or warnings to report "
143  << " | \n";
144  while (messageEntry != messages.end()) {
145  // Message printout.
146  string temp = messageEntry->first;
147  int len = temp.length();
148  temp.insert( len, max(0, 102 - len), ' ');
149  cout << " | " << setw(6) << messageEntry->second << " "
150  << temp << " | \n";
151  ++messageEntry;
152  }
153 
154  // Done.
155  cout << " | "
156  << " | \n"
157  << " *------- End LHAupHelaconia Error and Warning Messages "
158  << "Statistics ----------------------------------------------* "
159  << endl;
160 }
161 
162 //--------------------------------------------------------------------------
163 
164 // Read a HelacOnia command string.
165 
166 // The special command "set state = <pid>" is not passed to HelacOnia,
167 // but rather is used to set the color singlet state being produced
168 // from the color octet state (if a color octet state is being
169 // produced). If color octet production is enabled then the
170 // appropriate quark mass is modified to half the mass of the color
171 // octet state plus the color octet mass splitting.
172 
173 bool LHAupHelaconia::readString(string line) {
174 
175  size_t n = line.find("state");
176  if (line.find("8)") != string::npos) mQ = -1;
177  if (n != string::npos && pythia) {
178  pythia->settings.readString("Onia:" + line.substr(n));
179  nId = abs(pythia->settings.mode("Onia:state"));
180  nQ = int(nId/1e2) % 10;
181  nR = int(nId/1e5) % 10;
182  nL = int(nId/1e4) % 10;
183  nJ = int(nId/1e0) % 10;
184  } else lines.push_back(line);
185  return true;
186 
187 }
188 
189 //--------------------------------------------------------------------------
190 
191 // Set the random seed and maximum allowed runs.
192 
193 // If the random seed is negative (default of -1), then the HelacOnia
194 // seed is taken as the Pythia parameter "Random:seed", which must be
195 // greater than 0. If the maximum number of allowed runs is exceeded
196 // (default of 30081) an error is thrown. The seed for a HelacOnia run
197 // is set as:
198 
199 // (random seed - 1) * (maximum runs) + (number of runs) + 1
200 
201 // HelacOnia can only handle random seeds up to 30081 * 30081. So, with
202 // this strategy, one can generate Pythia jobs with seeds from 0 to
203 // 30081, with each job running HelacOnia less than 30081 times, and
204 // ensure a fully statistically independent sample. If more than 30081
205 // jobs are needed, then the maximum allowed runs can be lowered
206 // accordingly, and if need be, setEvents can be used to increase the
207 // number of events generated per run.
208 
209 bool LHAupHelaconia::setSeed(int seedIn, int runsIn) {
210 
211  if (!pythia) return false;
212  seed = seedIn;
213  if (seed < 0) {
214  seed = pythia->settings.mode("Random:seed");
215  if (seed < 1) {
216  errorMsg("Error from LHAupHelaconia::setSeed: the given "
217  "Pythia seed is less than 1."); return false;}
218  }
219  runs = runsIn;
220  if (seed * runs > 30081 * 30081) {
221  errorMsg("Error from LHAupHelaconia::setSeed: the given seed "
222  "exceeds the HelacOnia limit."); return false;}
223  nRuns = 0;
224  return true;
225 
226 }
227 
228 //--------------------------------------------------------------------------
229 
230 // Set the number of events to generate per HelacOnia run; default is 10000.
231 
232 void LHAupHelaconia::setEvents(int eventsIn) {events = eventsIn;}
233 
234 //--------------------------------------------------------------------------
235 
236 // Execute a system command.
237 
238 bool LHAupHelaconia::execute(string line) {return system(line.c_str()) != -1;}
239 
240 //--------------------------------------------------------------------------
241 
242 // Run HelacOnia.
243 
244 bool LHAupHelaconia::run(int eventsIn, int seedIn) {
245 
246  // Set up run and seed.
247  if (!pythia) return false;
248  if (nRuns >= runs) {
249  errorMsg("Error from LHAupHelaconia::run: maximum number "
250  "of allowed runs exceeded."); return false;}
251  if (seed < 0 && !setSeed(seed, runs)) return false;
252  if (seedIn < 0) seedIn = (seed - 1) * runs + nRuns + 1;
253 
254  // Determine the heavy quark mass.
255  if (mQ == -1)
256  mQ = (pythia->particleData.m0(nId)
257  + pythia->settings.parm("Onia:massSplit"))/2.0;
258 
259  // Write the generation file.
260  if (!pythia) return false;
261  fstream config((dir + "/generate.py").c_str(), ios::out);
262  for (int iLine = 0; iLine < (int)lines.size(); ++iLine)
263  config << lines[iLine] << "\n";
264  config << "set seed = " << seedIn << "\n"
265  << "set unwgt = T\n"
266  << "set unwevt = " << eventsIn << "\n"
267  << "set preunw = " << 3.0/2.0*eventsIn << "\n";
268  if (mQ > 0) config << "set " << (nQ == 4 ? "c" : "b") << "mass = " << mQ
269  << "\n";
270  config << "launch\n";
271  config.close();
272 
273  // Create the event shuffler.
274  fstream shuffle((dir + "/shuffle.py").c_str(), ios::out);
275  shuffle <<
276  "import random, os\n"
277  "random.seed(" << seedIn << "); tag, pre, post, events = '', [], [], []\n"
278  "for line in open('events.lhe').readlines():\n"
279  " if line.strip().startswith('<'):\n"
280  " tag = line.strip()\n"
281  " if tag == '<event>': events += ['<event>\\n']; continue\n"
282  " if tag == '</event>': events[-1] += '</event>\\n'; continue\n"
283  " if tag == '<event>': events[-1] += line\n"
284  " elif len(events) == 0: pre += [line]\n"
285  " else: post += [line]\n"
286  "random.shuffle(events); os.unlink('events.lhe')\n"
287  "open('events.lhe', 'w').writelines(pre + events + post)\n";
288  shuffle.close();
289 
290  // Execute the run commands.
291  if (!execute("rm -rf " + dir + "/PROC* " + lhegz)) return false;
292  if (!execute("cd " + dir + "; cat generate.py | " + exe)) return false;
293  if (!execute("cd " + dir + "; ln -s PROC_HO_0/P0_calc_0/output/*.lhe "
294  "events.lhe;# python shuffle.py")) return false;
295  if (access(lhegz.c_str(), F_OK) == -1) return false;
296  ++nRuns;
297  return true;
298 
299 }
300 
301 //--------------------------------------------------------------------------
302 
303 // Create the LHEF reader.
304 
305 bool LHAupHelaconia::reader(bool init) {
306 
307  // Check valid LHE file.
308  if (!pythia) return false;
309  if (lhef) delete lhef;
310  bool setScales(pythia->settings.flag("Beams:setProductionScalesFromLHEF"));
311  lhef = new LHAupLHEF(infoPtr, lhegz.c_str(), nullptr, false, setScales);
312  if (!lhef->setInit()) {
313  errorMsg("Error from LHAupHelaconia::reader: failed to "
314  "initialize the LHEF reader"); return false;}
315  if (lhef->sizeProc() != 1) {
316  errorMsg("Error from LHAupHelaconia::reader: number of "
317  "processes is not 1"); return false;}
318 
319  if (init) {
320 
321  // Determine the cross-section (if needed).
322  double sig(lhef->xSec(0)), err(lhef->xErr(0));
323 
324  // Set the info.
325  setBeamA(lhef->idBeamA(), lhef->eBeamA(), lhef->pdfGroupBeamA(),
326  lhef->pdfSetBeamA());
327  setBeamB(lhef->idBeamB(), lhef->eBeamB(), lhef->pdfGroupBeamB(),
328  lhef->pdfSetBeamB());
329  setStrategy(lhef->strategy());
330  addProcess(lhef->idProcess(0), sig, err, lhef->xMax(0));
331  xSecSumSave = sig; xErrSumSave = err;
332  }
333  return true;
334 
335 }
336 
337 //--------------------------------------------------------------------------
338 
339 // Convert the color octet HelacOnia ID to a Pythia 8 ID.
340 
341 int LHAupHelaconia::convert(int idIn) {
342 
343  if (abs(idIn) < 9900000) return idIn;
344  else idIn = abs(idIn) - 9900000;
345  int nS = 2;
346  if (idIn == nQ*110 + 3) nS = 0;
347  else if (idIn == nQ*110 + 1) nS = 1;
348  return 9900000 + 10000*nQ + 1000*nS + 100*nR + 10*nL + nJ;
349 
350 }
351 
352 //--------------------------------------------------------------------------
353 
354 // Set the initialization information.
355 
357 
358  // Create the LHEF LHAup object and run setInit.
359  if (!pythia) return false;
360  if (!run(events)) return false;
361  if (!reader(true)) return false;
362  listInit();
363  return true;
364 
365 }
366 
367 //--------------------------------------------------------------------------
368 
369 // Set the event information.
370 
372 
373  // Run setEvent from the LHEF object and launch HelacOnia if failed.
374  if (!pythia) return false;
375  if (!lhef) {
376  errorMsg("Error from LHAupHelaconia::setEvent: LHAupLHEF "
377  "object not correctly initialized"); return false;}
378  if (!lhef->fileFound()) {
379  errorMsg("Error from LHAupHelaconia::setEvent: LHEF "
380  "event file was not found"); return false;}
381  if (!lhef->setEvent()) {
382  if (!run(events)) return false;
383  if (!reader(false)) return false;
384  lhef->setEvent();
385  }
386 
387  // Read the event from the LHEF object.
388  particlesSave.clear();
389  int mom1, mom2;
390  for (int ip = 1; ip < lhef->sizePart(); ++ip) {
391  mom1 = lhef->mother1(ip);
392  mom2 = lhef->mother2(ip);
393  particlesSave.push_back
394  (LHAParticle(convert(lhef->id(ip)),
395  lhef->status(ip), mom1, mom2, lhef->col1(ip),
396  lhef->col2(ip), lhef->px(ip), lhef->py(ip), lhef->pz(ip),
397  lhef->e(ip), lhef->m(ip), lhef->tau(ip), lhef->spin(ip),
398  lhef->scale(ip)));
399  if (mom1 > 0 && mom1 < (int)particlesSave.size() && mom2 == 0)
400  particlesSave[mom1 - 1].statusPart = 2;
401  }
402 
403  // Write the event.
404  setProcess(lhef->idProcess(), lhef->weight(), lhef->scale(),
405  lhef->alphaQED(), lhef->alphaQCD());
406  for (int ip = 0; ip < (int)particlesSave.size(); ++ip)
407  addParticle(particlesSave[ip]);
408  setIdX(lhef->id1(), lhef->id2(), lhef->x1(), lhef->x2());
409  setPdf(lhef->id1pdf(), lhef->id2pdf(), lhef->x1pdf(), lhef->x2pdf(),
410  lhef->scalePDF(), lhef->pdf1(), lhef->pdf2(), lhef->pdfIsSet());
411  return true;
412 
413 }
414 
415 //==========================================================================
416 
417 } // end namespace Pythia8
418 
419 #endif // Pythia8_LHAHelaconia_H
void setEvents(int eventsIn)
Set the number of events to generate per run.
Definition: LHAHelaconia.h:232
void setBeamA(int idIn, double eIn, int pdfGroupIn=0, int pdfSetIn=0)
Input beam info.
Definition: LesHouches.h:222
void addProcess(int idProcIn, double xSecIn=1., double xErrIn=0., double xMaxIn=1.)
Input process info.
Definition: LesHouches.h:233
bool flag(string keyIn)
Give back current value, with check that key exists.
Definition: Settings.cc:1478
bool fileFound()
Confirm that file was found and opened as expected.
Definition: LesHouches.h:408
void setStrategy(int strategyIn)
Input process weight strategy.
Definition: LesHouches.h:230
ParticleData particleData
ParticleData: the particle data table/database.
Definition: Pythia.h:335
bool readString(string line, bool warn=true)
Read in one update from a single line.
Definition: Settings.cc:374
void setIdX(int id1In, int id2In, double x1In, double x2In)
Input info on flavour and x values of hard-process initiators.
Definition: LesHouches.h:263
bool setSeed(int seedIn, int runsIn=30081)
Set the random seed and maximum runs.
Definition: LHAHelaconia.h:209
A class for the particles stored in LHAup.
Definition: LesHouches.h:48
Settings settings
Settings: databases of flags/modes/parms/words to control run.
Definition: Pythia.h:332
int id1() const
Give back info on flavour and x values of hard-process initiators.
Definition: LesHouches.h:158
Definition: LesHouches.h:78
void addParticle(LHAParticle particleIn)
Input particle info, one particle at the time.
Definition: LesHouches.h:252
void setProcess(int idProcIn=0, double weightIn=1., double scaleIn=0., double alphaQEDIn=0.0073, double alphaQCDIn=0.12)
Input info on the selected process.
Definition: LesHouches.h:243
A derived class with information read from a Les Houches Event File.
Definition: LesHouches.h:346
void setPdf(int id1pdfIn, int id2pdfIn, double x1pdfIn, double x2pdfIn, double scalePDFIn, double pdf1In, double pdf2In, bool pdfIsSetIn)
Optionally input info on parton density values of event.
Definition: LesHouches.h:267
int convert(int idIn)
Convert the color octet HelacOnia ID to a Pythia 8 ID.
Definition: LHAHelaconia.h:341
A derived class from LHAup which generates events with HelacOnia.
Definition: LHAHelaconia.h:35
int sizePart() const
Give back info on separate particle.
Definition: LesHouches.h:141
void errorMsg(string messageIn)
Print a message the first few times. Insert in database.
Definition: LHAHelaconia.h:79
bool readString(string line)
Read a HelacOnia command string.
Definition: LHAHelaconia.h:173
LHAupHelaconia(Pythia *pythiaIn, string dirIn="helaconiarun", string exeIn="ho_cluster")
Constructor.
Definition: LHAHelaconia.h:112
bool setEvent(int=0)
Set the event information.
Definition: LHAHelaconia.h:371
bool execute(string line)
Execute a system command.
Definition: LHAHelaconia.h:238
~LHAupHelaconia()
Destructor.
Definition: LHAHelaconia.h:125
bool pdfIsSet() const
Optional: give back info on parton density values of event.
Definition: LesHouches.h:164
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
The Pythia class contains the top-level routines to generate an event.
Definition: Pythia.h:71
bool setEvent(int=0)
Routine for doing the job of reading and setting info on next event.
Definition: LesHouches.h:420
Info * infoPtr
Pointer to various information on the generation.
Definition: LesHouches.h:218
bool reader(bool init)
Create the LHEF reader.
Definition: LHAHelaconia.h:305
bool readString(string, bool warn=true, int subrun=SUBRUNDEFAULT)
Read in one update for a setting or particle data from a single line.
Definition: Pythia.cc:365
bool run(int eventsIn, int seedIn=-1)
Run HelacOnia.
Definition: LHAHelaconia.h:244
void listInit()
Print the initialization info; useful to check that setting it worked.
Definition: LesHouches.cc:33
bool setInit()
Set the initialization information.
Definition: LHAHelaconia.h:356