main144
Back to index.
// main144.cc is a part of the PYTHIA event generator.
// Copyright (C) 2025 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.
// Authors:
// Christian Bierlich
// Keywords:
// Analysis
// Hepmc
// Command file
// Command line option
// Root
// Rivet
// Tuning
// Streamlined event generation with possibility to output ROOT files,
// output HepMC files and run RIVET analyses, all by specifying output modes
// in a cmnd file, where also the event generator settings are specified.
// The example is run with command line options, run ./main144 -h to see a
// full list. See ROOT Usage for information about ROOT output, RIVET Usage
// for information about RIVET and HepMC Interface for information about HepMC.
#include "Pythia8/Pythia.h"
#include "Pythia8/HeavyIons.h"
#include "Pythia8Plugins/InputParser.h"
#include <chrono>
#ifdef RIVET
#include "Pythia8Plugins/Pythia8Rivet.h"
#endif
#ifdef PY8ROOT
#include "TSystem.h"
#include "TTree.h"
#include "TFile.h"
#endif
// Use the Pythia namespace.
using namespace Pythia8;
//==========================================================================
// Define filling ROOT particle and events if using ROOT.
#ifdef PY8ROOT
#include "main144Dct.h"
// Fill a ROOT particle.
RootParticle::RootParticle(Pythia8::Particle &prt) {
phi = prt.phi();
eta = prt.eta();
y = prt.y();
pT = prt.pT();
pid = prt.id();
}
// Fill a ROOT event to a TTree.
void RootEvent::fill(const Pythia8::Info &infoIn, vector<RootParticle> &prtsIn,
TTree *treeIn) {
weight = infoIn.weight();
particles = prtsIn;
treeIn->Fill();
}
#endif
//==========================================================================
// Implement a user supplied UserHooks derived class inside this
// wrapper, which will allow you to give settings that can be supplied
// in command files.
class UserHooksWrapper : public UserHooks {
public:
// Add the settings you want available in the run card in this method.
void additionalSettings(Settings* settingsIn) {
settings = settingsIn;
settings->addFlag("UserHooks:doMPICut", false);
settings->addMode("UserHooks:nMPICut", 0, true, false, 0, 0);
}
// Check if parton level can be vetoed.
bool canVetoPartonLevel() final {
return settings->flag("UserHooks:doMPICut");}
// Check if parton level should be vetoed.
bool doVetoPartonLevel(const Event&) final {
return infoPtr->nMPI() < settings->mode("UserHooks:nMPICut");}
private:
Settings* settings{};
};
//==========================================================================
// Main example execution.
int main(int argc, char* argv[]) {
// Parser object for command line input.
InputParser ip("Run Pythia with cmnd file input, and get Rivet, HepMC or"
" standard Pythia output.",
{"./main144 [options]", "./main144 -c main144.cmnd -n 1000 -o myoutput"},
"Additional options in cmnd file:\n"
"\tMain:writeLog = on\n\t\tRedirect output to <-o prefix>.log.\n"
"\tMain:writeHepMC = on \n\t\tWrite HepMC output, requires HepMC linked.\n"
"\tMain:writeRoot = on \n\t\tWrite a ROOT tree declared in "
"RootEvent.h, requires ROOT linked.\n"
"\tMain:runRivet = on \n\t\tRun Rivet analyses, requires Rivet linked.\n"
"\tMain:rivetAnalyses = {ANALYSIS1,ANALYSIS2,...}\n "
"\t\tComma separated list of Rivet analyses to run.\n"
"\t\tAnalysis names can be post-fixed with analysis parameters.\n"
"\t\tANALYSIS:parm=value:parm2=value2:...\n"
"\tMain:rivetRunName = STRING \n\t\tAdd an optional run name to "
"the Rivet analysis.\n"
"\tMain:rivetIgnoreBeams = on\n\t\tIgnore beams in Rivet. \n"
"\tMain:rivetDumpPeriod = NUMBER\n\t\tDump Rivet histograms "
"to file evert NUMBER of events.\n"
"\tMain:rivetDumpFile = STRING\n\t\t Specify alternative "
"name for Rivet dump file. Default = OUT.\n");
// Set up command line options.
ip.require("c", "User-written command file, can use multiple times.",
{"-cmnd"});
ip.add("s", "-1", "Specify seed for the random number generator.",
{"-seed"});
ip.add("o", "main144", "Output prefix for log file, Rivet, HepMC, and ROOT.",
{"-out"});
ip.add("n", "-1", "Number of events. Overrides the command files.",
{"-nevents"});
ip.add("l", "false", "Silence the splash screen.");
ip.add("t", "false", "Time event generation.", {"-time"});
ip.add("v", "false", "Print Pythia version number and exit.", {"-version"});
// Initialize the parser and exit if necessary.
InputParser::Status status = ip.init(argc, argv);
if (status != InputParser::Valid) return status;
// Print version number and exit.
if (ip.get<bool>("v")) {
cout << "PYTHIA version: " << PYTHIA_VERSION << endl;
return 0;
}
// Get the command files.
vector<string> cmnds = ip.getVector<string>("c");
if (cmnds.size() == 0) {
cout << "Please provide one or more command files with the -c option."
<< endl;
return 1;
}
// Random number seed.
string seed = ip.get<string>("s");
// Output filename.
string out = ip.get<string>("o");
// Time event generation.
bool writeTime = ip.get<bool>("t");
// Command line number of event, overrides the one set in input .cmnd file.
int nev = ip.get<int>("n");
// Catch the splash screen in a buffer.
stringstream splashBuf;
std::streambuf* sBuf = cout.rdbuf();
cout.rdbuf(splashBuf.rdbuf());
// The Pythia object.
Pythia pythia;
// Direct cout back.
cout.rdbuf(sBuf);
// UserHooks wrapper.
shared_ptr<UserHooksWrapper> userHooksWrapper =
make_shared<UserHooksWrapper>();
userHooksWrapper->additionalSettings(&pythia.settings);
pythia.setUserHooksPtr(userHooksWrapper);
// Some extra parameters.
pythia.settings.addFlag("Main:writeLog", false);
pythia.settings.addFlag("Main:writeHepMC", false);
pythia.settings.addFlag("Main:writeRoot", false);
pythia.settings.addFlag("Main:runRivet", false);
pythia.settings.addFlag("Main:rivetIgnoreBeams", false);
pythia.settings.addMode("Main:rivetDumpPeriod", -1, true, false, -1, 0);
pythia.settings.addWord("Main:rivetDumpFile", "");
pythia.settings.addWord("Main:rivetRunName", "");
pythia.settings.addWVec("Main:rivetAnalyses", {});
pythia.settings.addWVec("Main:rivetPreload", {});
// Read the command files.
for (int iCmnd = 0; iCmnd < (int)cmnds.size(); ++iCmnd)
if (!cmnds[iCmnd].empty()) pythia.readFile(cmnds[iCmnd]);
// Set seed after reading input.
if(seed != "-1") {
pythia.readString("Random:setSeed = on");
pythia.readString("Random:seed = "+seed);
}
// Read the extra parameters.
if (nev > -1) pythia.settings.mode("Main:numberOfEvents", nev);
int nEvent = pythia.mode("Main:numberOfEvents");;
int nError = pythia.mode("Main:timesAllowErrors");
bool writeLog = pythia.flag("Main:writeLog");
bool writeHepmc = pythia.flag("Main:writeHepMC");
bool writeRoot = pythia.flag("Main:writeRoot");
bool runRivet = pythia.flag("Main:runRivet");
bool countErrors = nError > 0;
// Check if Rivet, HepMC, and ROOT are requested and available.
bool valid = true;
#ifndef RIVET
valid = valid && !runRivet && !writeHepmc;
if (runRivet)
cout << "Option Main::runRivet = on requires the Rivet library.\n";
if (writeHepmc)
cout << "Option Main::writeHepMC = on requires the HepMC library.\n";
#endif
#ifndef PY8ROOT
valid = valid && !writeRoot;
if (writeRoot)
cout << "Option Main::writeRoot = on requires the ROOT library.\n";
#endif
if (!valid) return 1;
// Rivet and HepMC initialization.
#ifdef RIVET
// Initialize HepMC.
Pythia8ToHepMC hepmc;
if (writeHepmc) hepmc.setNewFile(out + ".hepmc");
// Initialize Rivet.
Pythia8Rivet rivet(pythia, out + ".yoda");
rivet.ignoreBeams(pythia.flag("Main:rivetIgnoreBeams"));
rivet.dump(pythia.settings.mode("Main:rivetDumpPeriod"),
pythia.settings.word("Main:rivetDumpFile"));
// Load the analyses.
vector<string> rivetAnalyses = pythia.settings.wvec("Main:rivetAnalyses");
for (int iAna = 0; iAna < (int)rivetAnalyses.size(); ++iAna)
rivet.addAnalysis(rivetAnalyses[iAna]);
// Pre-load the YODA histograms.
vector<string> rivetPreload = pythia.settings.wvec("Main:rivetPreload");
for (int iYoda = 0; iYoda < (int)rivetPreload.size(); ++iYoda)
rivet.addPreload(rivetPreload[iYoda]);
// Add the run name.
rivet.addRunName(pythia.settings.word("Main:rivetRunName"));
#endif
// ROOT initialization.
#ifdef PY8ROOT
// Create the ROOT TFile and TTree.
TFile *file;
TTree *tree;
RootEvent *evt;
if (writeRoot) {
// Open the ROOT file.
file = TFile::Open((out + ".root").c_str(), "recreate" );
tree = new TTree("t", "Pythia8 event tree");
evt = new RootEvent();
// Set the TTree branch to the ROOT event.
tree->Branch("events", &evt);
}
#endif
// Logfile initialization.
ofstream logBuf;
streambuf *oldCout;
if (writeLog) {
oldCout = cout.rdbuf(logBuf.rdbuf());
logBuf.open(out + ".log");
}
// Remove splash screen, if requested.
ostream cnull(NULL);
if (ip.get<bool>("l")) cnull << splashBuf.str();
else cout << splashBuf.str();
// If Pythia fails to initialize, exit with error.
if (!pythia.init()) return 1;
// Make a sanity check of initialized Rivet analyses.
#ifdef RIVET
if (!runRivet && rivetAnalyses.size() > 0 )
cout << "Rivet analyses are set with Main:rivetAnalyses, "
<< "but Main:runRivet = off.\n";
#endif
// Loop over events.
auto startAllEvents = std::chrono::high_resolution_clock::now();
for ( int iEvent = 0; iEvent < nEvent; ++iEvent ) {
auto startThisEvent = std::chrono::high_resolution_clock::now();
// Exit if too many failures.
if (!pythia.next()) {
if (countErrors && --nError < 0) {
pythia.stat();
cout << " \n *------- PYTHIA STOPPED! -----------------------*\n"
<< " | Event generation failed due to too many errors. |\n"
<< " *-------------------------------------------------*\n";
return 1;
}
continue;
}
// Calculate the event time.
auto stopThisEvent = std::chrono::high_resolution_clock::now();
auto eventTime = std::chrono::duration_cast<std::chrono::milliseconds>
(stopThisEvent - startThisEvent);
double tt = eventTime.count();
// Run the Rivet analyses.
#ifdef RIVET
if (runRivet) {
if (writeTime) rivet.addAttribute("EventTime", tt);
rivet();
}
if (writeHepmc) hepmc.writeNextEvent(pythia);
#endif
// Write to ROOT file output.
#ifdef PY8ROOT
if (writeRoot) {
vector<RootParticle> prts;
for (int iPrt = 0; iPrt < pythia.event.size(); ++iPrt) {
Particle& prt = pythia.event[iPrt];
// Any particle cuts can be placed here. Here, only final
// state particles are kept.
if (!prt.isFinal()) continue;
// Push back the ROOT particle.
prts.push_back(RootParticle(prt));
}
// Fill the ROOT event and tree.
evt->fill(pythia.info, prts, tree);
}
#endif
}
// Finalize.
pythia.stat();
#ifdef PY8ROOT
if (writeRoot) {
tree->Print();
tree->Write();
delete file, tree, evt;
}
#endif
// Print timing.
auto stopAllEvents = std::chrono::high_resolution_clock::now();
auto durationAll = std::chrono::duration_cast<std::chrono::milliseconds>
(stopAllEvents - startAllEvents);
if (writeTime) {
cout << " \n *------- Generation time -----------------------*\n"
<< " | Event generation, analysis and writing to files |\n"
<< " | took: " << double(durationAll.count()) << " ms or "
<< double(durationAll.count())/double(nEvent)
<< " ms per event |\n"
<< " *-------------------------------------------------*\n";
}
// Put cout back in its place.
if (writeLog) cout.rdbuf(oldCout);
return 0;
}