include/Pythia8Plugins/LHAHelaconia.h

0001 // LHAHelaconia.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005
0006 // Author: Philip Ilten, December 2017.
0007
0008 #ifndef Pythia8_LHAHelaconia_H
0009 #define Pythia8_LHAHelaconia_H
0010
0011 #include "Pythia8/Pythia.h"
0012 #include <unistd.h>
0013 #include <sys/stat.h>
0014 #include <limits>
0015
0016 namespace Pythia8 {
0017
0018 //==========================================================================
0019
0020 // A derived class from LHAup which generates events with HelacOnia.
0021
0022 // This class automatically generates hard processes with HelacOnia
0023 // and reads in the LHEF file output.
0024
0025 // The user can send commands to HelacOnia via the readString
0026 // method. The launch command, random seed, and shower choice are
0027 // automatically handled. For example, the following will produce
0028 // J/psi events from 13 TeV proton proton collisions:
0029
0030 //    readString("generate u u~ > cc~(3S11) g")
0031
0032 // The number of events generated per HelacOnia run is controlled by
0033 // setEvents, while the random seed is controlled by setSeed.
0034
0035 class LHAupHelaconia : public LHAup {
0036
0037 public:
0038
0039   // Constructor.
0040   LHAupHelaconia(Pythia *pythiaIn, string dirIn = "helaconiarun",
0041                  string exeIn = "ho_cluster");
0042
0043   // Destructor: Print error statistics before exiting. Printing code
0044   // basically copied from Info class.
0045   ~LHAupHelaconia();
0046
0047   // Read a HelacOnia command string.
0048   bool readString(string line);
0049
0050   // Set the number of events to generate per run.
0051   void setEvents(int eventsIn);
0052
0053   // Set the random seed and maximum runs.
0054   bool setSeed(int seedIn, int runsIn = 30081);
0055
0056   // Set the initialization information.
0057   bool setInit();
0058
0059   // Set the event information.
0060   bool setEvent(int = 0);
0061
0062   // Note: The functions below have been made public to ease the generation
0063   // of Python bindings.
0064   //protected:
0065
0066   // Execute a system command.
0067   bool execute(string line);
0068
0069   // Run HelacOnia.
0070   bool run(int eventsIn, int seedIn = -1);
0071
0072   // Create the LHEF reader.
0073   bool reader(bool init);
0074
0075   // Convert the color octet HelacOnia ID to a Pythia 8 ID.
0076   int convert(int idIn);
0077
0078   // Print a message the first few times. Insert in database.
0079   void errorMsg(string messageIn) {
0080     // Recover number of times message occured. Also inserts new string.
0081     int times = messages[messageIn];
0082     ++messages[messageIn];
0083     // Print message the first few times.
0084     if (times < TIMESTOPRINT) cout << " PYTHIA " << messageIn << endl;
0085   }
0086
0087 private:
0088
0089   // The PYTHIA object and LHEF file reader and matching hook.
0090   Pythia *pythia;
0091   LHAupLHEF *lhef;
0092
0093   // Stored members.
0094   int events, seed, runs, nRuns, nId, nQ, nR, nL, nJ;
0095   string dir, exe, lhegz;
0096   double mQ;
0097
0098   // The HelacOnia commands.
0099   vector<string> lines;
0100
0101   // Map for all error messages.
0102   map<string, int> messages;
0103   // Number of times the same error message is repeated, unless overridden.
0104   static const int TIMESTOPRINT = 1;
0105
0106 };
0107
0108 //--------------------------------------------------------------------------
0109
0110 // Constructor.
0111
0112 LHAupHelaconia::LHAupHelaconia(Pythia *pythiaIn, string dirIn, string exeIn) :
0113   pythia(pythiaIn), lhef(0), events(10000), seed(-1), runs(30081),
0114   nRuns(0), nId(443), nQ(4), nR(0), nL(0), nJ(3),
0115   dir(dirIn), exe(exeIn), lhegz(dirIn + "/events.lhe"), mQ(-2) {
0116   mkdir(dir.c_str(), 0777);
0117   if (pythia) pythia->readString("Beams:frameType = 5");
0118   pythia->settings.addMode("Onia:state", -1, false, false, 0, 0);
0119 }
0120
0121 //--------------------------------------------------------------------------
0122
0123 // Destructor.
0124
0125 LHAupHelaconia::~LHAupHelaconia() {
0126
0127   if (lhef) delete lhef;
0128
0129   // Header.
0130   cout << "\n *-------  LHAupHelaconia Error and Warning Messages Statistics"
0131        << "  --------------------------------------------------* \n"
0132        << " |                                                       "
0133        << "                                                          | \n"
0134        << " |  times   message                                      "
0135        << "                                                          | \n"
0136        << " |                                                       "
0137        << "                                                          | \n";
0138
0139   // Loop over all messages
0140   map<string, int>::iterator messageEntry = messages.begin();
0141   if (messageEntry == messages.end())
0142     cout << " |      0   no errors or warnings to report              "
0143          << "                                                          | \n";
0144   while (messageEntry != messages.end()) {
0145     // Message printout.
0146     string temp = messageEntry->first;
0147     int len = temp.length();
0148     temp.insert( len, max(0, 102 - len), ' ');
0149     cout << " | " << setw(6) << messageEntry->second << "   "
0150          << temp << " | \n";
0151     ++messageEntry;
0152   }
0153
0154   // Done.
0155   cout << " |                                                       "
0156        << "                                                          | \n"
0157        << " *-------  End LHAupHelaconia Error and Warning Messages "
0158        << "Statistics  ----------------------------------------------* "
0159        << endl;
0160 }
0161
0162 //--------------------------------------------------------------------------
0163
0164 // Read a HelacOnia command string.
0165
0166 // The special command "set state = <pid>" is not passed to HelacOnia,
0167 // but rather is used to set the color singlet state being produced
0168 // from the color octet state (if a color octet state is being
0169 // produced). If color octet production is enabled then the
0170 // appropriate quark mass is modified to half the mass of the color
0171 // octet state plus the color octet mass splitting.
0172
0173 bool LHAupHelaconia::readString(string line) {
0174
0175   size_t n = line.find("state");
0176   if (line.find("8)") != string::npos) mQ = -1;
0177   if (n != string::npos && pythia) {
0178     pythia->settings.readString("Onia:" + line.substr(n));
0179     nId = abs(pythia->settings.mode("Onia:state"));
0180     nQ  = int(nId/1e2) % 10;
0181     nR  = int(nId/1e5) % 10;
0182     nL  = int(nId/1e4) % 10;
0183     nJ  = int(nId/1e0) % 10;
0184   } else lines.push_back(line);
0185   return true;
0186
0187 }
0188
0189 //--------------------------------------------------------------------------
0190
0191 // Set the random seed and maximum allowed runs.
0192
0193 // If the random seed is negative (default of -1), then the HelacOnia
0194 // seed is taken as the Pythia parameter "Random:seed", which must be
0195 // greater than 0. If the maximum number of allowed runs is exceeded
0196 // (default of 30081) an error is thrown. The seed for a HelacOnia run
0197 // is set as:
0198
0199 //    (random seed - 1) * (maximum runs) + (number of runs) + 1
0200
0201 // HelacOnia can only handle random seeds up to 30081 * 30081. So, with
0202 // this strategy, one can generate Pythia jobs with seeds from 0 to
0203 // 30081, with each job running HelacOnia less than 30081 times, and
0204 // ensure a fully statistically independent sample. If more than 30081
0205 // jobs are needed, then the maximum allowed runs can be lowered
0206 // accordingly, and if need be, setEvents can be used to increase the
0207 // number of events generated per run.
0208
0209 bool LHAupHelaconia::setSeed(int seedIn, int runsIn) {
0210
0211   if (!pythia) return false;
0212   seed = seedIn;
0213   if (seed < 0) {
0214     seed = pythia->settings.mode("Random:seed");
0215     if (seed < 1) {
0216       errorMsg("Error from LHAupHelaconia::setSeed: the given "
0217                "Pythia seed is less than 1."); return false;}
0218   }
0219   runs = runsIn;
0220   if (seed * runs > 30081 * 30081) {
0221     errorMsg("Error from LHAupHelaconia::setSeed: the given seed "
0222              "exceeds the HelacOnia limit."); return false;}
0223   nRuns = 0;
0224   return true;
0225
0226 }
0227
0228 //--------------------------------------------------------------------------
0229
0230 // Set the number of events to generate per HelacOnia run; default is 10000.
0231
0232 void LHAupHelaconia::setEvents(int eventsIn) {events = eventsIn;}
0233
0234 //--------------------------------------------------------------------------
0235
0236 // Execute a system command.
0237
0238 bool LHAupHelaconia::execute(string line) {return system(line.c_str()) != -1;}
0239
0240 //--------------------------------------------------------------------------
0241
0242 // Run HelacOnia.
0243
0244 bool LHAupHelaconia::run(int eventsIn, int seedIn) {
0245
0246   // Set up run and seed.
0247   if (!pythia) return false;
0248   if (nRuns >= runs) {
0249     errorMsg("Error from LHAupHelaconia::run: maximum number "
0250              "of allowed runs exceeded."); return false;}
0251   if (seed < 0 && !setSeed(seed, runs)) return false;
0252   if (seedIn < 0) seedIn = (seed - 1) * runs + nRuns + 1;
0253
0254   // Determine the heavy quark mass.
0255   if (mQ == -1)
0256     mQ = (pythia->particleData.m0(nId)
0257           + pythia->settings.parm("Onia:massSplit"))/2.0;
0258
0259   // Write the generation file.
0260   if (!pythia) return false;
0261   fstream config((dir + "/generate.py").c_str(), ios::out);
0262   for (int iLine = 0; iLine < (int)lines.size(); ++iLine)
0263     config << lines[iLine] << "\n";
0264   config << "set seed = " << seedIn << "\n"
0265          << "set unwgt = T\n"
0266          << "set unwevt = " << eventsIn << "\n"
0267          << "set preunw = " << 3.0/2.0*eventsIn << "\n";
0268   if (mQ > 0) config << "set " << (nQ == 4 ? "c" : "b") << "mass = " << mQ
0269                      << "\n";
0270   config << "launch\n";
0271   config.close();
0272
0273   // Create the event shuffler.
0274   fstream shuffle((dir + "/shuffle.py").c_str(), ios::out);
0275   shuffle <<
0276     "import random, os\n"
0277     "random.seed(" << seedIn << "); tag, pre, post, events = '', [], [], []\n"
0278     "for line in open('events.lhe').readlines():\n"
0279     "    if line.strip().startswith('<'):\n"
0280     "        tag = line.strip()\n"
0281     "        if tag == '<event>':  events += ['<event>\\n']; continue\n"
0282     "        if tag == '</event>': events[-1] += '</event>\\n'; continue\n"
0283     "    if tag == '<event>': events[-1] += line\n"
0284     "    elif len(events) == 0: pre += [line]\n"
0285     "    else: post += [line]\n"
0286     "random.shuffle(events); os.unlink('events.lhe')\n"
0287     "open('events.lhe', 'w').writelines(pre + events + post)\n";
0288   shuffle.close();
0289
0290   // Execute the run commands.
0291   if (!execute("rm -rf " + dir + "/PROC* " + lhegz)) return false;
0292   if (!execute("cd " + dir + "; cat generate.py | " + exe)) return false;
0293   if (!execute("cd " + dir + "; ln -s PROC_HO_0/P0_calc_0/output/*.lhe "
0294                "events.lhe;# python shuffle.py")) return false;
0295   if (access(lhegz.c_str(), F_OK) == -1) return false;
0296   ++nRuns;
0297   return true;
0298
0299 }
0300
0301 //--------------------------------------------------------------------------
0302
0303 // Create the LHEF reader.
0304
0305 bool LHAupHelaconia::reader(bool init) {
0306
0307   // Check valid LHE file.
0308   if (!pythia) return false;
0309   if (lhef) delete lhef;
0310   bool setScales(pythia->settings.flag("Beams:setProductionScalesFromLHEF"));
0311   lhef = new LHAupLHEF(infoPtr, lhegz.c_str(), nullptr, false, setScales);
0312   if (!lhef->setInit()) {
0313     errorMsg("Error from LHAupHelaconia::reader: failed to "
0314              "initialize the LHEF reader"); return false;}
0315   if (lhef->sizeProc() != 1) {
0316     errorMsg("Error from LHAupHelaconia::reader: number of "
0317              "processes is not 1"); return false;}
0318
0319   if (init) {
0320
0321     // Determine the cross-section (if needed).
0322     double sig(lhef->xSec(0)), err(lhef->xErr(0));
0323
0324     // Set the info.
0325     setBeamA(lhef->idBeamA(), lhef->eBeamA(), lhef->pdfGroupBeamA(),
0326              lhef->pdfSetBeamA());
0327     setBeamB(lhef->idBeamB(), lhef->eBeamB(), lhef->pdfGroupBeamB(),
0328              lhef->pdfSetBeamB());
0329     setStrategy(lhef->strategy());
0330     addProcess(lhef->idProcess(0), sig, err, lhef->xMax(0));
0331     xSecSumSave = sig; xErrSumSave = err;
0332   }
0333   return true;
0334
0335 }
0336
0337 //--------------------------------------------------------------------------
0338
0339 // Convert the color octet HelacOnia ID to a Pythia 8 ID.
0340
0341 int LHAupHelaconia::convert(int idIn) {
0342
0343     if (abs(idIn) < 9900000) return idIn;
0344     else idIn = abs(idIn) - 9900000;
0345     int nS = 2;
0346     if (idIn == nQ*110 + 3) nS = 0;
0347     else if (idIn == nQ*110 + 1) nS = 1;
0348     return 9900000 + 10000*nQ + 1000*nS + 100*nR + 10*nL + nJ;
0349
0350 }
0351
0352 //--------------------------------------------------------------------------
0353
0354 // Set the initialization information.
0355
0356 bool LHAupHelaconia::setInit() {
0357
0358   // Create the LHEF LHAup object and run setInit.
0359   if (!pythia) return false;
0360   if (!run(events)) return false;
0361   if (!reader(true)) return false;
0362   listInit();
0363   return true;
0364
0365 }
0366
0367 //--------------------------------------------------------------------------
0368
0369 // Set the event information.
0370
0371 bool LHAupHelaconia::setEvent(int) {
0372
0373   // Run setEvent from the LHEF object and launch HelacOnia if failed.
0374   if (!pythia) return false;
0375   if (!lhef) {
0376     errorMsg("Error from LHAupHelaconia::setEvent: LHAupLHEF "
0377              "object not correctly initialized"); return false;}
0378   if (!lhef->fileFound()) {
0379     errorMsg("Error from LHAupHelaconia::setEvent: LHEF "
0380              "event file was not found"); return false;}
0381   if (!lhef->setEvent()) {
0382     if (!run(events)) return false;
0383     if (!reader(false)) return false;
0384     lhef->setEvent();
0385   }
0386
0387   // Read the event from the LHEF object.
0388   particlesSave.clear();
0389   int mom1, mom2;
0390   for (int ip = 1; ip < lhef->sizePart(); ++ip) {
0391     mom1 = lhef->mother1(ip);
0392     mom2 = lhef->mother2(ip);
0393     particlesSave.push_back
0394       (LHAParticle(convert(lhef->id(ip)),
0395                    lhef->status(ip), mom1, mom2, lhef->col1(ip),
0396                    lhef->col2(ip), lhef->px(ip), lhef->py(ip), lhef->pz(ip),
0397                    lhef->e(ip), lhef->m(ip), lhef->tau(ip), lhef->spin(ip),
0398                    lhef->scale(ip)));
0399     if (mom1 > 0 && mom1 < (int)particlesSave.size() && mom2 == 0)
0400       particlesSave[mom1 - 1].statusPart = 2;
0401   }
0402
0403   // Write the event.
0404   setProcess(lhef->idProcess(), lhef->weight(), lhef->scale(),
0405     lhef->alphaQED(), lhef->alphaQCD());
0406   for (int ip = 0; ip < (int)particlesSave.size(); ++ip)
0407     addParticle(particlesSave[ip]);
0408   setIdX(lhef->id1(), lhef->id2(), lhef->x1(), lhef->x2());
0409   setPdf(lhef->id1pdf(), lhef->id2pdf(), lhef->x1pdf(), lhef->x2pdf(),
0410          lhef->scalePDF(), lhef->pdf1(), lhef->pdf2(), lhef->pdfIsSet());
0411   return true;
0412
0413 }
0414
0415 //==========================================================================
0416
0417 } // end namespace Pythia8
0418
0419 #endif // Pythia8_LHAHelaconia_H