main134

// main134.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:
//            Stefan Prestel

// Contact: Christian Preuss <preuss@uni-wuppertal.de>

// Keywords:
//            LHE file
//            Hepmc

// This program illustrates how a file with HepMC2 or HepMC3 events
// can be generated by Pythia8 using LHEF input. Input and output files are
// specified on the command line, e.g. like
//     ./main134 -c main134.cmnd -o main134.hepmc > main134.log

#include "Pythia8/Pythia.h"
#include "Pythia8Plugins/InputParser.h"
#ifndef HEPMC2
#include "Pythia8Plugins/HepMC3.h"
#else
#include "Pythia8Plugins/HepMC2.h"
#endif
#include <unistd.h>

using namespace Pythia8;

//==========================================================================

// Example main programm to illustrate merging.

int main( int argc, char* argv[] ) {

  // Set up command line options.
  InputParser ip("This program illustrates how HepMC files can be written by"
    " Pythia8.", {"./main133 -c main133.cmnd -o main133.hepmc"});
  ip.require("c", "Use this user-written command file.", {"-cmnd"});
  ip.require("o", "Specify HepMC output filename.", {"-out"});

  // Initialize the parser and exit if necessary.
  InputParser::Status status = ip.init(argc, argv);
  if (status != InputParser::Valid) return status;

  // Confirm that external files will be used for input and output.
  string cmnd(ip.get<string>("c")), out(ip.get<string>("o"));
  cout << "\n >>> PYTHIA settings will be read from file '" << cmnd
       << "' <<< \n >>> HepMC events will be written to file '"
       << out << "' <<< \n";

  // Input parameters.
  Pythia pythia;
  pythia.readFile(cmnd, 0);

  // Interface for conversion from Pythia8::Event to HepMC one.
  // Specify file where HepMC events will be stored.
  Pythia8ToHepMC toHepMC(out);

  // Allow abort of run if many errors.
  int  nAbort  = pythia.mode("Main:timesAllowErrors");
  int  iAbort  = 0;
  bool doAbort = false;

  // Read in loop parameters.
  cout << endl << endl << endl;
  cout << "Start generating events" << endl;
  long nEvent = pythia.settings.mode("Main:numberOfEvents");
  int  nRuns  = pythia.mode("Main:numberOfSubruns");
  double sigmaTotal(0.), errorTotal(0.);

  // Loop over subruns with varying number of jets.
  for (int iRuns = 0; iRuns < nRuns; ++iRuns) {
    double sigmaSample = 0., errorSample = 0.;

    // Read in name of LHE file for current subrun and initialize.
    pythia.readFile(cmnd, iRuns);

    // If Pythia fails to initialize, exit with error.
    if (!pythia.init()) return 1;

    // Get the inclusive x-section by summing over all process x-sections.
    double xs = 0.;
    for (int i=0; i < pythia.info.nProcessesLHEF(); ++i)
      xs += pythia.info.sigmaLHEF(i);

    // Start generation loop.
    while( pythia.info.nSelected() < nEvent ) {

      // Generate next event.
      if( !pythia.next() ) {
        if ( pythia.info.atEndOfFile() ) break;
        else if (++iAbort > nAbort) {doAbort = true; break;}
        else continue;
      }

      // Get event weight(s).
      double evtweight = pythia.info.weight();

      // Do not print zero-weight events.
      if ( evtweight == 0. ) continue;

      // Inform HepMC3 about the naming of the weights.
#ifndef HEPMC2
      toHepMC.setWeightNames(pythia.info.weightNameVector());
#endif

      // Work with weighted (LHA strategy=-4) events.
      double normhepmc = 1.;
      if (abs(pythia.info.lhaStrategy()) == 4)
        normhepmc = 1. / double(1e9*nEvent);
      // Work with unweighted events.
      else
        normhepmc = xs / double(1e9*nEvent);

      // Set event weights (optional).
#ifndef HEPMC2
      // hepmcevt.weights().push_back(evtweight*normhepmc);
#endif

      // Fill a new HepMC event.
      toHepMC.fillNextEvent( pythia );

      // Add the weight of the current event to the cross section.
      sigmaTotal  += evtweight*normhepmc;
      sigmaSample += evtweight*normhepmc;
      errorTotal  += pow2(evtweight*normhepmc);
      errorSample += pow2(evtweight*normhepmc);

      // Report cross section to HepMC.
      toHepMC.setXSec( sigmaTotal*1e9, errorTotal*1e9 );

      // Write the HepMC event to file. Done with it.
      toHepMC.writeEvent();

    } // End loop over events to generate.
    if (doAbort) break;

    // Print cross section and errors.
    pythia.stat();
    cout << endl << " Contribution of sample " << iRuns
         << " to the inclusive cross section : "
         << scientific << setprecision(8)
         << sigmaSample << "  +-  " << sqrt(errorSample)  << endl;
  }

  // Abort warning.
  cout << endl << endl << endl;
  if (doAbort)
    cout << " Run was not completed owing to too many aborted events" << endl;
  cout << endl << endl << endl;

  // Done
  return 0;

}