source: HiSusy/trunk/hepmc/x86_64-slc5-gcc41-opt/share/HepMC/examples/fio/example_PythiaStreamIO.cc @ 1

//////////////////////////////////////////////////////////////////////////
// example_PythiaStreamIO.cc
//
// garren@fnal.gov, May 2009
// 
//////////////////////////////////////////////////////////////////////////
/// example of generating events with Pythia using HepMC/PythiaWrapper.h 
/// Events are read into the HepMC event record from the FORTRAN HEPEVT 
/// common block using the IO_HEPEVT strategy 
///
/// To Compile: go to the HepMC example directory and type:
/// make example_PythiaStreamIO.exe
///
/// This example uses streaming I/O 
/// writePythiaStreamIO() sets the cross section in GenRun
/// readPythiaStreamIO() reads the file written by writePythiaStreamIO()
///
//////////////////////////////////////////////////////////////////////////


#include <fstream>
#include <iostream>
#include "HepMC/PythiaWrapper.h"
#include "HepMC/IO_HEPEVT.h"
#include "HepMC/GenEvent.h"
#include "PythiaHelper.h"

void writePythiaStreamIO();
void readPythiaStreamIO();

int main() { 

    writePythiaStreamIO();
    readPythiaStreamIO();

    return 0;
}
   

void writePythiaStreamIO() {
    // example to generate events and write output
    std::cout << std::endl;
    std::cout << "Begin pythia_out()" << std::endl;
   //........................................HEPEVT
    // Pythia 6.1 uses HEPEVT with 4000 entries and 8-byte floating point
    //  numbers. We need to explicitly pass this information to the 
    //  HEPEVT_Wrapper.
    //
    HepMC::HEPEVT_Wrapper::set_max_number_entries(4000);
    HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
    //
    //........................................PYTHIA INITIALIZATIONS
    initPythia();

    //........................................HepMC INITIALIZATIONS
    //
    // Instantiate an IO strategy for reading from HEPEVT.
    HepMC::IO_HEPEVT hepevtio;
    //
    { // begin scope of ascii_io
        // declare an output stream
        const char outfile[] = "example_PythiaStreamIO_write.dat";
        std::ofstream ascii_io( outfile );
        if( !ascii_io ) {
          std::cerr << "cannot open " << outfile << std::endl;
          exit(-1);
        }
        // use the default IO_GenEvent precision
        ascii_io.precision(16);
        // write the line that defines the beginning of a GenEvent block
        HepMC::write_HepMC_IO_block_begin( ascii_io );
        //
        //........................................EVENT LOOP
        for ( int i = 1; i <= 100; i++ ) {
            if ( i%50==1 ) std::cout << "Processing Event Number " 
                                     << i << std::endl;
            call_pyevnt();      // generate one event with Pythia
            // pythia pyhepc routine converts common PYJETS in common HEPEVT
            call_pyhepc( 1 );
            HepMC::GenEvent* evt = hepevtio.read_next_event();
            // define the units (Pythia uses GeV and mm)
            evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
            // add some information to the event
            evt->set_event_number(i);
            evt->set_signal_process_id(20);
            // set number of multi parton interactions
            evt->set_mpi( pypars.msti[31-1] );
            // set cross section information
            evt->set_cross_section( HepMC::getPythiaCrossSection() );
            // write the event out to the ascii files
            ascii_io << (*evt);;
            // we also need to delete the created event from memory
            delete evt;
        }
        // write the line that defines the end of a GenEvent block
        HepMC::write_HepMC_IO_block_end( ascii_io );
        //........................................TERMINATION
        // write out some information from Pythia to the screen
        call_pystat( 1 );    
    } // end scope of ascii_io
}

void readPythiaStreamIO() {
    // example to read events written by writePythiaStreamIO
    // and write them back out
    std::cout << std::endl;
    // input units are GeV and mm
    const char infile[] = "example_PythiaStreamIO_write.dat";
    std::ifstream is( infile );
    if( !is ) {
      std::cerr << "cannot open " << infile << std::endl;
      exit(-1);
    }
    //
    { // begin scope of ascii_io
        // declare an output stream
        const char outfile[] = "example_PythiaStreamIO_read.dat";
        std::ofstream ascii_io( outfile );
        if( !ascii_io ) {
          std::cerr << "cannot open " << outfile << std::endl;
          exit(-1);
        }
        ascii_io.precision(16);
        HepMC::write_HepMC_IO_block_begin( ascii_io );
        //
        //........................................EVENT LOOP
        HepMC::GenEvent evt;
        int i = 0;
        while ( is ) {
            evt.read( is );
            // make sure we have a valid event
            if( evt.is_valid() ) {
                ++i;
                if ( i%50==1 ) std::cout << "Processing Event Number " 
                                         << i << std::endl;
                if ( i%25==2 ) {
                    // write the cross section if it exists
                    if( evt.cross_section() ) {
                        std::cout << "cross section at event " << i << " is " 
                                  << evt.cross_section()->cross_section()
                                  << std::endl;
                    }
                }
                // write the event out to the ascii files
                evt.write( ascii_io );
            }
        }
        //........................................TERMINATION
        HepMC::write_HepMC_IO_block_end( ascii_io );
    } // end scope of ascii_io
}