hepmc - Rev 87

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"
// Also write events using IO_AsciiParticles and IO_ExtendedAscii
//////////////////////////////////////////////////////////////////////////
//

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

bool findW( HepMC::GenEvent* evt, std::ofstream& os);

/// \class  IsGoodEvent
/// event selection predicate. returns true if the event contains
/// a photon with pT > 50 GeV
class IsGoodEvent {
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;
    }
};

/// \class  IsW_Boson
/// this predicate returns true if the input particle is a W+/W-
class IsW_Boson {
public:
    /// returns true if the GenParticle is a W
    bool operator()( const HepMC::GenParticle* p ) { 
        if ( fabs(p->pdg_id()) == 24 ) 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 IO_ExtendedAscii for output
    HepMC::IO_ExtendedAscii xout("testHepMCExtended.out",std::ios::out);
    // declare an instance of the event selection predicate
    IsGoodEvent is_good_event;
    // define an output stream
    std::ofstream os( "testIterator.out" );
    //........................................EVENT LOOP
    int icount=0;
    int num_good_events=0;
    double x=0., y=0., z=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 ( icount==100 ) std::cout << "Processing Event Number " << icount
                                      << " its # " << evt->event_number() 
                                      << std::endl;
        if ( is_good_event(evt) ) {
            // add some arbitrary PDF information
            x = 0.1 * icount;
            y = 0.13 * icount;
            z = 0.012 * icount;
            HepMC::PdfInfo* pdf = new HepMC::PdfInfo( 11, 12, x, y, z, 0.11, 0.34);
            evt->set_pdf_info(pdf);
            ascii_out << evt;
            particle_out << evt;
            xout << evt;
            ++num_good_events;
            // test iterators
            findW( evt, os );
            delete pdf;
        }
        
        // 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;
}



bool findW( HepMC::GenEvent* evt, std::ofstream& os )
{
    int num_W=0;
        // use GenEvent::particle_iterator to find all W's in the event,
        // then 
        // (1) for each W user the GenVertex::particle_iterator with a range of
        //     parents to return and print the immediate mothers of these W's.
        // (2) for each W user the GenVertex::particle_iterator with a range of
        //     descendants to return and print all descendants of these W's.
        IsW_Boson isw;
        for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin();
              p != evt->particles_end(); ++p ) {
            if ( isw(*p) ) {
                ++num_W;
                os << "A W boson has been found in event: " << evt->event_number() << std::endl;
                //std::cout << "A W boson has been found in event: " << evt->event_number() << std::endl;
                (*p)->print( os );
                // return all parents
                // we do this by pointing to the production vertex of the W 
                // particle and asking for all particle parents of that vertex
                os << "\t Its parents are: " << std::endl;
                if ( (*p)->production_vertex() ) {
                    for ( HepMC::GenVertex::particle_iterator mother 
                              = (*p)->production_vertex()->
                              particles_begin(HepMC::parents);
                          mother != (*p)->production_vertex()->
                              particles_end(HepMC::parents); 
                          ++mother ) {
                        os << "\t";
                        (*mother)->print( os );
                    }
                }
                // return all descendants
                // we do this by pointing to the end vertex of the W 
                // particle and asking for all particle descendants of that vertex
                os << "\t\t Its descendants are: " << std::endl;
                if ( (*p)->end_vertex() ) {
                    for ( HepMC::GenVertex::particle_iterator des 
                              =(*p)->end_vertex()->
                              particles_begin(HepMC::descendants);
                          des != (*p)->end_vertex()->
                              particles_end(HepMC::descendants);
                          ++des ) {
                        os << "\t\t";
                        (*des)->print( os );
                    }
                }
            }
        }
        return true;
}