hepmc - Rev 35

Subversion Repositories:
Rev:
//////////////////////////////////////////////////////////////////////////
// testHepMC.cc.in
//
// garren@fnal.gov, March 2006
// Use Matt's example_EventSelection to check HepMC.
// 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"
//////////////////////////////////////////////////////////////////////////
//

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

class IsGoodEvent {
    // event selection predicate. returns true if the event contains
    // a photon with pT > 50 GeV
public:
    bool operator()( const HepMC::GenEvent* evt ) { 
        for ( HepMC::GenEvent::particle_const_iterator p 
                  = evt->particles_begin(); p != evt->particles_end(); ++p ){
            if ( (*p)->pdg_id() == 22 && (*p)->momentum().perp() > 25. ) {
                //std::cout << "Event " << evt->event_number()
                //     << " is a good event." << std::endl;
                //(*p)->print();
                return 1;
            }
        }
        return 0;
    }
};

int main() { 
    // declare an input strategy to read the data produced with the 
    // example_MyPythia
    HepMC::IO_Ascii ascii_in("@srcdir@/testHepMC.input",std::ios::in);
    // declare another IO_Ascii for writing out the good events
    HepMC::IO_Ascii 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;
        }
        delete evt;
        ascii_in >> evt;
    }
    //........................................PRINT RESULT
    std::cout << num_good_events << " out of " << icount 
              << " processed events passed the cuts. Finished." << std::endl;
}