HepMC Reference Documentation HepMC

testHepMC.cc.in

based on example_EventSelection Apply an event selection to the events in testHepMC.input Events containing a photon of pT > 25 GeV pass the selection and are written to "testHepMC.out" Add arbitrary PDF information to the good events Also write events using IO_AsciiParticles

// testHepMC.cc.in
//
// garren@fnal.gov, March 2006
// based on example_EventSelection
// Apply an event selection to the events in testHepMC.input
// Events containing a photon of pT > 25 GeV pass the selection and are
// written to "testHepMC.out"
// Add arbitrary PDF information to the good events
// Also write events using IO_AsciiParticles 
//

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

// define methods and classes used by this test
#include "IsGoodEvent.h"

void read_testIOGenEvent();
void read_testAscii();
void read_variousFormats();

double findPiZero( HepMC::GenEvent* );

int main() { 
    read_testIOGenEvent();
    read_testAscii();
    read_variousFormats();
    return 0;
}

void read_testIOGenEvent()
{
    std::cout << std::endl;
    std::cout << "basic IO_GenEvent input and output" << std::endl;
    // declare an input strategy to read the data produced with the 
    // example_MyPythia - units are GeV and mm
    HepMC::IO_GenEvent ascii_in("@srcdir@/testIOGenEvent.input",std::ios::in);
    ascii_in.use_input_units( HepMC::Units::GEV, HepMC::Units::MM );
    // declare another IO_GenEvent for writing out the good events
    HepMC::IO_GenEvent ascii_out("testHepMC.out",std::ios::out);
    // declare an IO_AsciiParticle for output
    HepMC::IO_AsciiParticles particle_out("testHepMCParticle.out",std::ios::out);
    // declare an instance of the event selection predicate
    IsGoodEvent is_good_event;
    //........................................EVENT LOOP
    int icount=0;
    int num_good_events=0;
    HepMC::GenEvent* evt = ascii_in.read_next_event();
    while ( evt ) {
        icount++;
        if ( icount%50==1 ) std::cout << "Processing Event Number " << icount
                                      << " its # " << evt->event_number() 
                                      << std::endl;
        if ( is_good_event(evt) ) {
            ascii_out << evt;
            particle_out << evt;
            ++num_good_events;
        }
        
        // clean up and get next event
        delete evt;
        ascii_in >> evt;
    }
    //........................................PRINT RESULT
    std::cout << num_good_events << " out of " << icount 
              << " processed events passed the cuts. Finished." << std::endl;
}

void read_testAscii()
{
    std::cout << std::endl;
    std::cout << "read file written with IO_Ascii" << std::endl;
    // declare an input strategy to read the data produced with the 
    // example_MyPythia - units are GeV and mm
    HepMC::IO_GenEvent ascii_in("@srcdir@/testAscii.input",std::ios::in);
    ascii_in.use_input_units( HepMC::Units::GEV, HepMC::Units::MM );
    if ( ascii_in.rdstate() == std::ios::failbit ) {
        std::cerr << "ERROR input file @srcdir@/testAscii.input is needed "
                  << "and does not exist.  Exit." << std::endl;
        return;
    }
    // use IO_Ascii because we are checking against the original input file 
    HepMC::IO_Ascii ascii_out("testIOAscii.dat",std::ios::out);
    if ( ascii_out.rdstate() == std::ios::failbit ) {
        std::cerr << "ERROR opening output file testAscii.dat.  Exit."
                  << std::endl;
        return;
    }
    //........................................EVENT LOOP
    int icount=0;
    HepMC::GenEvent* evt = ascii_in.read_next_event();
    while ( evt ) {
        icount++;
        if ( icount%50==1 ) std::cout << "Processing Event Number " << icount
                                      << " its # " << evt->event_number() 
                                      << std::endl;
        ascii_out << evt;

        // clean up and get next event
        delete evt;
        ascii_in >> evt;
    }
    //........................................PRINT RESULT
    std::cout << icount << " events processed. Finished." << std::endl;
}

void read_variousFormats()
{
    std::cout << std::endl;
    std::cout << "process varied input" << std::endl;
    // declare an input strategy 
    HepMC::IO_GenEvent ascii_in("@srcdir@/testHepMCVarious.input",std::ios::in);
    ascii_in.use_input_units( HepMC::Units::GEV, HepMC::Units::MM );
    // declare another IO_GenEvent for writing out the good events
    HepMC::IO_GenEvent ascii_out("testHepMCVarious.out",std::ios::out);
    //........................................EVENT LOOP
    int icount=0;
    HepMC::GenEvent* evt = ascii_in.read_next_event();
    while ( evt ) {
        icount++;
        double pim;
        std::cout << "Processing Event Number " << icount
                  << " its # " << evt->event_number() 
                  << std::endl;
        ascii_out << evt;
        // units should be unknown
        evt->write_units();
        pim = findPiZero(evt);
        std::cout << " pizero mass: " << pim << std::endl;
        // set units to GeV and mm
        evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
        evt->write_units();
        pim = findPiZero(evt);
        std::cout << " pizero mass: " << pim 
                  << " " << HepMC::Units::name( evt->momentum_unit() ) << std::endl;
        // convert units to MeV
        evt->use_units(HepMC::Units::MEV, HepMC::Units::MM);
        evt->write_units();
        pim = findPiZero(evt);
        std::cout << " pizero mass: " << pim 
                  << " " << HepMC::Units::name( evt->momentum_unit() ) << std::endl;
        // clean up and get next event
        delete evt;
        ascii_in >> evt;
    }
    //........................................PRINT RESULT
    std::cout << icount << " events processed. Finished." << std::endl;
}

double findPiZero( HepMC::GenEvent* evt )
{
    for ( HepMC::GenEvent::particle_const_iterator p 
              = evt->particles_begin(); p != evt->particles_end(); ++p ){
        if ( (*p)->pdg_id() == 111 ) {
            return (*p)->generated_mass();
        }
    }
    return 0.;
}

---