main62.cc

// main62.cc is a part of the PYTHIA event generator.
// Copyright (C) 2012 Mikhail Kirsanov, Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

// This is a simple test program.
// It illustrates how a file with HepMC events can be generated by Pythia8.
// It is similar to main42, except that it is compiled with LHAPDF
// and allows for several subruns, e.g. from related LHEF.
// Input and output files are specified on the command line, e.g. like
// ./main62.exe main62.cmnd hepmcout62.dat > out
// The main program contains no analysis; this is intended to happen later.
// It therefore "never" has to be recompiled to handle different tasks.

// WARNING: typically one needs 25 MB/100 events at the LHC.
// Therefore large event samples may be impractical.

#include "Pythia.h"
#include "HepMCInterface.h"

#include "HepMC/GenEvent.h"   
#include "HepMC/IO_GenEvent.h"

// Following line is a deprecated alternative, removed in recent versions.
//#include "HepMC/IO_Ascii.h"
//#include "HepMC/IO_AsciiParticles.h"

// Following line to be used with HepMC 2.04 onwards.
#ifdef HEPMC_HAS_UNITS
#include "HepMC/Units.h"
#endif

using namespace Pythia8; 

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

  // Check that correct number of command-line arguments
  if (argc != 3) {
    cerr << " Unexpected number of command-line arguments. \n You are"
         << " expected to provide one input and one output file name. \n"
         << " Program stopped! " << endl;
    return 1;
  }

  // Check that the provided input name corresponds to an existing file.
  ifstream is(argv[1]);  
  if (!is) {
    cerr << " Command-line file " << argv[1] << " was not found. \n"
         << " Program stopped! " << endl;
    return 1;
  }

  // Confirm that external files will be used for input and output.
  cout << "\n >>> PYTHIA settings will be read from file " << argv[1] 
       << " <<< \n >>> HepMC events will be written to file " 
       << argv[2] << " <<< \n" << endl;

  // Interface for conversion from Pythia8::Event to HepMC one. 
  HepMC::I_Pythia8 ToHepMC;
  //  ToHepMC.set_crash_on_problem();

  // Specify file where HepMC events will be stored.
  HepMC::IO_GenEvent ascii_io(argv[2], std::ios::out);
  // Following line is a deprecated alternative, removed in recent versions
  // HepMC::IO_Ascii ascii_io("hepmcout32.dat", std::ios::out);
  // Line below is an eye-readable one-way output, uncomment the include above
  // HepMC::IO_AsciiParticles ascii_io("hepmcout32.dat", std::ios::out);
 
  // Generator. 
  Pythia pythia;

  // Read in subrun-independent commands from external file.
  pythia.readFile( argv[1]);    
                          
  // Extract data to be used in main program. Set counters.
  int nSubrun = pythia.mode("Main:numberOfSubruns"); 
  int nAbort  = pythia.mode("Main:timesAllowErrors");
  int iAbort  = 0;

  // Begin loop over subruns.
  for (int iSubrun = 1; iSubrun <= nSubrun; ++iSubrun) {
 
    // Read in subrun-specific data from external file.
    pythia.readFile( argv[1], iSubrun); 

    // Initialization.
    pythia.init();

    // Print name of Les Houches Event File.
    string lheFile = pythia.settings.word("Beams:LHEF");
    cout << "\n >>> Now begin subrun " << iSubrun 
         << " with events from file " << lheFile << " <<< \n"
         << endl; 

    // Begin infinite event loop - to be exited at end of file.
    for (int iEvent = 0; ; ++iEvent) {

      // Generate event. 
      if (!pythia.next()) {

        // Leave event loop if at end of file.
        if (pythia.info.atEndOfFile()) break;

        // First few failures write off as "acceptable" errors, then quit.
        if (++iAbort < nAbort) continue;
        cout << " Event generation aborted prematurely, owing to error!\n"; 
        break;
      }

      // Construct new empty HepMC event. 
#ifdef HEPMC_HAS_UNITS
      // This form with arguments is only meaningful for HepMC 2.04 onwards, 
      // and even then unnecessary if HepMC was built with GeV and mm as units.
      HepMC::GenEvent* hepmcevt = new HepMC::GenEvent( 
        HepMC::Units::GEV, HepMC::Units::MM);
#else
      // This form is needed for backwards compatibility. 
      // In HepMCInterface.cc a conversion from GeV to MeV will be done.
      HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
#endif

      // Fill HepMC event, including PDF info.
      ToHepMC.fill_next_event( pythia, hepmcevt );
      // This alternative older method fills event, without PDF info.
      // ToHepMC.fill_next_event( pythia.event, hepmcevt );

      // Write the HepMC event to file. Done with it.
      ascii_io << hepmcevt;
      delete hepmcevt;

    // End of event loop.
    }
  
  // End of subrun loop.
  }

  // Statistics. 
  pythia.stat();

  // Done.
  return 0;
}