HepMC Reference Documentation HepMC

GenEvent.h

Go to the documentation of this file.

//--------------------------------------------------------------------------
#ifndef HEPMC_GEN_EVENT_H
#define HEPMC_GEN_EVENT_H

// Matt.Dobbs@Cern.CH, September 1999, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// Event record for MC generators (for use at any stage of generation)
//
// This class is intended as both a "container class" ( to store a MC
//  event for interface between MC generators and detector simulation )
//  and also as a "work in progress class" ( that could be used inside
//  a generator and modified as the event is built ).
//
// Iterators are provided which allow the user to easily obtain a
//  list of particles or vertices in an event --- this list can be filled
//  subject to some sort of selection criteria. Examples are given below
//  ( see HepMC::copy_if and std::copy )

namespace HepMC {

    // To create a list from an iterator, use: (i.e. for a list of particles);
    // #include <algorithm>
    //     list<GenParticle*> thelist;
    //     copy( evt->particles_begin(), evt->particles_end(), 
    //           back_inserter(thelist) );
    // to create a list subject to a condition (predicate) use:
    //     list<GenParticle*> thelist;
    //     HepMC::copy_if( evt->particles_begin(), evt->particles_end(), 
    //                     back_inserter(thelist), is_photon() );
    // where is_photon() is a predicate like:
    //     class is_photon {
    //       public:
    //         bool operator() ( GenParticle const * p ) {
    //             if ( p && p->pdg_id() == 22 ) return true;
    //             return false;
    //         }
    //     };
    // which the user defines herself.

    template <class InputIterator, class OutputIterator, class Predicate>
    void copy_if( InputIterator first, InputIterator last, OutputIterator out,
                  Predicate pred ) {
        for ( ; first != last; ++first ) { if ( pred(*first) ) out = *first; }
    }
} // HepMC

// Since a container of all vertices in the event is maintained, the time
//  required to loop over all vertices (or particles) is very fast -- and 
//  the user does not gain much by first making his own list.
//  (this is not true for the GenVertex:: versions of these iterators, which
//   allow you to specify the vertex starting point and range)

// Data Members:
// signal_process_id()   The integer ID that uniquely specifies this signal
//                       process, i.e. MSUB in Pythia. It is necessary to
//                       package this with each event rather than with the run
//                       because many processes may be generated within one
//                       run.
// event_number()        Strictly speaking we cannot think of any reason that
//                       an event would need to know its own event number, it
//                       is more likely something that would be assigned by
//                       a database. It is included anyway (tradition?) since
//                       we expect it may be useful for debugging. It can
//                       be reset later by a database.
// mpi()                 The number of multi parton interactions in the event.
//                       This is NOT beam pileup.  Set to -1 by default.
// beam_particles()      A pair of pointers to the incoming beam particles.
// signal_process_vertex() pointer to the vertex containing the signal process
// weights()             Vector of doubles which specify th weight of the evnt,
//                       the first entry will be the "event weight" used for
//                       hit and miss etc., but a general vector is used to
//                       allow for reweighting etc. We envision a list of
//                       WeightTags to be included with a run class which
//                       would specify the meaning of the Weights .
// random_states()       Vector of integers which specify the random number 
//                       generator's state for this event. It is left to the
//                       generator to make use of this. We envision a vector of
//                       RndmStatesTags to be included with a run class which
//                       would specify the meaning of the random_states.
//
// Memory allocation //
// -When a vertex (particle) is added to a event (vertex), it is "adopted" 
//  and becomes the responsibility of the event (vertex) to delete that 
//  particle. 
// -objects responsible for deleting memory:
//    -events delete included vertices
//    -each vertex deletes its outgoing particles which do not have decay
//     vertices
//    -each vertex deletes its incoming particles which do not
//     have creation vertices 
//
// About the Barcodes //
// - each vertex or particle has a barcode, which is just an integer which
//   uniquely identifies it inside the event (i.e. there is a one to one
//   mapping between particle memory addresses and particle barcodes... and 
//   the same applied for vertices)
// - The value of a barcode has NO MEANING and NO ORDER!
//   For the user's convenience, when an event is read in via an IO_method
//   from an indexed list (like the HEPEVT common block), then the index will
//   become the barcode for that particle.
// - particle barcodes are always positive integers
//   vertex barcodes are always negative integers
//   The barcodes are chosen and set automatically when a vertex or particle
//   comes under the ownership of an event (i.e. it is contained in an event).
// - You can tell when a particle or vertex is owned, because its 
//   parent_event() return value will return a pointer to the event which owns
//   it (or null if its an orphan).
// 

#include "HepMC/GenVertex.h"
#include "HepMC/GenParticle.h"
#include "HepMC/WeightContainer.h"
#include "HepMC/GenCrossSection.h"
#include "HepMC/HeavyIon.h"
#include "HepMC/PdfInfo.h"
#include "HepMC/Units.h"
#include "HepMC/HepMCDefs.h"
#include <map>
#include <string>
#include <vector>
#include <algorithm>
#include <iostream>

namespace HepMC {


    class GenEvent {
        friend class GenParticle;
        friend class GenVertex;  
    public:
        GenEvent( int signal_process_id = 0, int event_number = 0,
                  GenVertex* signal_vertex = 0,
                  const WeightContainer& weights = std::vector<double>(),
                  const std::vector<long>& randomstates = std::vector<long>(),
                  Units::MomentumUnit = Units::default_momentum_unit(), 
                  Units::LengthUnit = Units::default_length_unit() );
        GenEvent( int signal_process_id, int event_number,
                  GenVertex* signal_vertex, const WeightContainer& weights,
                  const std::vector<long>& randomstates,
                  const HeavyIon& ion, const PdfInfo& pdf,
                  Units::MomentumUnit = Units::default_momentum_unit(), 
                  Units::LengthUnit = Units::default_length_unit() );
        GenEvent( Units::MomentumUnit, Units::LengthUnit,
                  int signal_process_id = 0, int event_number = 0,
                  GenVertex* signal_vertex = 0,
                  const WeightContainer& weights = std::vector<double>(),
                  const std::vector<long>& randomstates = std::vector<long>() );
        GenEvent( Units::MomentumUnit, Units::LengthUnit,
                  int signal_process_id, int event_number,
                  GenVertex* signal_vertex, const WeightContainer& weights,
                  const std::vector<long>& randomstates,
                  const HeavyIon& ion, const PdfInfo& pdf );
        GenEvent( const GenEvent& inevent );          
        GenEvent& operator=( const GenEvent& inevent ); 
        virtual ~GenEvent(); 

        void swap( GenEvent & other );  
    
        void print( std::ostream& ostr = std::cout ) const; 
        void print_version( std::ostream& ostr = std::cout ) const; 

        GenParticle* barcode_to_particle( int barCode ) const;
        GenVertex*   barcode_to_vertex(   int barCode ) const;

        // access methods //

        int signal_process_id() const; 
        int event_number() const; 
        int mpi() const;          
        double event_scale() const; 
        double alphaQCD() const; 
        double alphaQED() const; 

        GenVertex* signal_process_vertex() const;
        bool valid_beam_particles() const;
        std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beam_particles() const;
        bool is_valid() const;

        WeightContainer&        weights(); 
        const WeightContainer&  weights() const; 

        GenCrossSection const *     cross_section() const;
        GenCrossSection*            cross_section();
        HeavyIon const *         heavy_ion() const;
        HeavyIon*                heavy_ion();
        PdfInfo const *          pdf_info() const;
        PdfInfo*                 pdf_info();

        const std::vector<long>& random_states() const;

        int     particles_size() const;
        bool    particles_empty() const;
        int     vertices_size() const;
        bool    vertices_empty() const;

        void write_units( std::ostream & os = std::cout ) const; 

        Units::MomentumUnit momentum_unit() const;
        Units::LengthUnit   length_unit()   const;
        
        std::ostream& write(std::ostream&);
        std::istream& read(std::istream&);

        // mutator methods //

        bool    add_vertex( GenVertex* vtx );    
        bool    remove_vertex( GenVertex* vtx ); 
        void    clear();                         
        
        void set_signal_process_id( int id ); 
        void set_event_number( int eventno ); 
        void set_mpi( int  ); 
        void set_event_scale( double scale ); 
        void set_alphaQCD( double a ); 
        void set_alphaQED( double a ); 

        void set_signal_process_vertex( GenVertex* );
        bool set_beam_particles(GenParticle*, GenParticle*);
        bool set_beam_particles(std::pair<HepMC::GenParticle*,HepMC::GenParticle*> const &);
        void set_random_states( const std::vector<long>& randomstates );

        void set_cross_section( const GenCrossSection& );
        void set_heavy_ion( const HeavyIon& ion );
        void set_pdf_info( const PdfInfo& p );
        
        void use_units( Units::MomentumUnit, Units::LengthUnit );
        void use_units( std::string&, std::string& );

    public:
        // vertex_iterators          //
        // Note:  the XXX_iterator is "resolvable" as XXX_const_iterator, but 
        //  not the reverse, which is consistent with STL, 
        //  see Musser, Derge, Saini 2ndEd. p. 69,70.


        class vertex_const_iterator :
          public std::iterator<std::forward_iterator_tag,HepMC::GenVertex*,ptrdiff_t>{
            // Iterates over all vertices in this event
        public:
            vertex_const_iterator(
                const 
                std::map<int,HepMC::GenVertex*,std::greater<int> >::const_iterator& i)
                : m_map_iterator(i) {}
            vertex_const_iterator() {}
            vertex_const_iterator( const vertex_const_iterator& i )
                { *this = i; }
            virtual ~vertex_const_iterator() {}
            vertex_const_iterator&  operator=( const vertex_const_iterator& i )
                { m_map_iterator = i.m_map_iterator; return *this; }
            GenVertex* operator*(void) const { return m_map_iterator->second; }
            vertex_const_iterator&  operator++(void)  //Pre-fix increment 
                { ++m_map_iterator; return *this; }
            vertex_const_iterator   operator++(int)   //Post-fix increment
                { vertex_const_iterator out(*this); ++(*this); return out; }
            bool  operator==( const vertex_const_iterator& a ) const
                { return m_map_iterator == a.m_map_iterator; }
            bool  operator!=( const vertex_const_iterator& a ) const
                { return !(m_map_iterator == a.m_map_iterator); }
        protected:
            std::map<int,HepMC::GenVertex*,std::greater<int> >::const_iterator 
                                                                m_map_iterator;
        private:
            vertex_const_iterator&  operator--(void);
            vertex_const_iterator   operator--(int);
        };
        friend class vertex_const_iterator;
        vertex_const_iterator      vertices_begin() const
            { return GenEvent::vertex_const_iterator( 
                m_vertex_barcodes.begin() ); }
        vertex_const_iterator      vertices_end() const
            { return GenEvent::vertex_const_iterator(
                m_vertex_barcodes.end() ); }



        class vertex_iterator :
          public std::iterator<std::forward_iterator_tag,HepMC::GenVertex*,ptrdiff_t>{
            // Iterates over all vertices in this event
        public:
            vertex_iterator( 
                const 
                std::map<int,HepMC::GenVertex*,std::greater<int> >::iterator& i )
                : m_map_iterator( i ) {}
            vertex_iterator() {}
            vertex_iterator( const vertex_iterator& i ) { *this = i; }
            virtual ~vertex_iterator() {}
            vertex_iterator&  operator=( const vertex_iterator& i ) {
                m_map_iterator = i.m_map_iterator;
                return *this;
            }
            operator vertex_const_iterator() const
                { return vertex_const_iterator(m_map_iterator); }
            GenVertex*        operator*(void) const
                { return m_map_iterator->second; }
            vertex_iterator&  operator++(void)  //Pre-fix increment 
                { ++m_map_iterator;     return *this; }
            vertex_iterator   operator++(int)   //Post-fix increment
                { vertex_iterator out(*this); ++(*this); return out; }
            bool              operator==( const vertex_iterator& a ) const
                { return m_map_iterator == a.m_map_iterator; }
            bool              operator!=( const vertex_iterator& a ) const
                { return !(m_map_iterator == a.m_map_iterator); }
        protected:
            std::map<int,HepMC::GenVertex*,std::greater<int> >::iterator 
                                                               m_map_iterator;
        private:
            vertex_iterator&  operator--(void);
            vertex_iterator   operator--(int);

        };
        friend class vertex_iterator;
        vertex_iterator            vertices_begin() 
            { return GenEvent::vertex_iterator( 
                m_vertex_barcodes.begin() ); }
        vertex_iterator            vertices_end()
            { return GenEvent::vertex_iterator(
                m_vertex_barcodes.end() ); }

    public:
        // particle_iterator         //
        // Example of iterating over all particles in the event:
        //      for ( GenEvent::particle_const_iterator p = particles_begin();
        //            p != particles_end(); ++p ) {
        //         (*p)->print();
        //      }
        //


        class particle_const_iterator :
          public std::iterator<std::forward_iterator_tag,HepMC::GenParticle*,ptrdiff_t>{
            // Iterates over all vertices in this event
        public:
            particle_const_iterator(
                const std::map<int,HepMC::GenParticle*>::const_iterator& i )
                : m_map_iterator(i) {}
            particle_const_iterator() {}
            particle_const_iterator( const particle_const_iterator& i )
                { *this = i; }
            virtual ~particle_const_iterator() {}
            particle_const_iterator& operator=(
                const particle_const_iterator& i )
                { m_map_iterator = i.m_map_iterator; return *this; }
            GenParticle*        operator*(void) const
                { return m_map_iterator->second; }
            particle_const_iterator&  operator++(void)  //Pre-fix increment 
                { ++m_map_iterator; return *this; }
            particle_const_iterator   operator++(int)   //Post-fix increment
                { particle_const_iterator out(*this); ++(*this); return out; }
            bool  operator==( const particle_const_iterator& a ) const
                { return m_map_iterator == a.m_map_iterator; }
            bool  operator!=( const particle_const_iterator& a ) const
                { return !(m_map_iterator == a.m_map_iterator); }
        protected:
            std::map<int,HepMC::GenParticle*>::const_iterator m_map_iterator;
        private:
            particle_const_iterator&  operator--(void);
            particle_const_iterator   operator--(int);
        };      
        friend class particle_const_iterator;
        particle_const_iterator      particles_begin() const
            { return GenEvent::particle_const_iterator( 
                m_particle_barcodes.begin() ); }
        particle_const_iterator      particles_end() const
            { return GenEvent::particle_const_iterator(
                m_particle_barcodes.end() ); }


        class particle_iterator :
          public std::iterator<std::forward_iterator_tag,HepMC::GenParticle*,ptrdiff_t>{
            // Iterates over all vertices in this event
        public:
            particle_iterator( const std::map<int,HepMC::GenParticle*>::iterator& i )
                : m_map_iterator( i ) {}
            particle_iterator() {}
            particle_iterator( const particle_iterator& i ) { *this = i; }
            virtual ~particle_iterator() {}
            particle_iterator&  operator=( const particle_iterator& i ) {
                m_map_iterator = i.m_map_iterator;
                return *this;
            }
            operator particle_const_iterator() const
                { return particle_const_iterator(m_map_iterator); }
            GenParticle*        operator*(void) const
                { return m_map_iterator->second; }
            particle_iterator&  operator++(void) 
                { ++m_map_iterator;     return *this; }
            particle_iterator   operator++(int)   
                { particle_iterator out(*this); ++(*this); return out; }
            bool              operator==( const particle_iterator& a ) const
                { return m_map_iterator == a.m_map_iterator; }
            bool              operator!=( const particle_iterator& a ) const
                { return !(m_map_iterator == a.m_map_iterator); }
        protected:
            std::map<int,HepMC::GenParticle*>::iterator m_map_iterator;
        private:
            particle_iterator&  operator--(void);
            particle_iterator   operator--(int);
        };
        friend class particle_iterator;
        particle_iterator particles_begin() 
            { return GenEvent::particle_iterator(
                m_particle_barcodes.begin() ); }
        particle_iterator particles_end()
            { return GenEvent::particle_iterator(
                m_particle_barcodes.end() ); }

    protected:
        //
        // Following methods intended for use by GenParticle/Vertex classes:
        // In general there is no reason they should be used elsewhere.
        bool         set_barcode( GenParticle* p, int suggested_barcode =false );
        bool         set_barcode( GenVertex*   v, int suggested_barcode =false );
        void         remove_barcode( GenParticle* p );
        void         remove_barcode( GenVertex*   v );

        void delete_all_vertices(); 

     private: // methods
        bool use_momentum_unit( Units::MomentumUnit );
        bool use_momentum_unit( std::string& );
        bool use_length_unit( Units::LengthUnit );
        bool use_length_unit( std::string& );
        
        // the following internal methods are used by read() and write()

        std::ostream & write_beam_particles( std::ostream &, 
                             std::pair<HepMC::GenParticle *,HepMC::GenParticle *> );
        std::ostream & write_vertex( std::ostream &, GenVertex const * );
        std::ostream & write_particle( std::ostream&, GenParticle const * );
        std::istream & find_file_type( std::istream & );
        std::istream & find_end_key( std::istream &, int & );
        std::istream & read_units( std::istream & );

    private: // data members
        int                   m_signal_process_id;
        int                   m_event_number;  
        int                   m_mpi;        // number of multi paricle interactions
        double                m_event_scale;// energy scale, see hep-ph/0109068
        double                m_alphaQCD;   // QCD coupling, see hep-ph/0109068
        double                m_alphaQED;   // QED coupling, see hep-ph/0109068
        GenVertex*            m_signal_process_vertex;
        GenParticle*          m_beam_particle_1;
        GenParticle*          m_beam_particle_2;
        WeightContainer       m_weights; // weights for this event first weight
                                         // is used by default for hit and miss
        std::vector<long> m_random_states; // container of rndm num 
                                               // generator states

        std::map< int,HepMC::GenVertex*,std::greater<int> >   m_vertex_barcodes;
        std::map< int,HepMC::GenParticle*,std::less<int> >    m_particle_barcodes;
        GenCrossSection*         m_cross_section;             // undefined by default
        HeavyIon*             m_heavy_ion;            // undefined by default
        PdfInfo*              m_pdf_info;             // undefined by default
        Units::MomentumUnit   m_momentum_unit;    // default value set by configure switch
        Units::LengthUnit     m_position_unit;    // default value set by configure switch

    };


    // IO Free Functions     //
  
    std::ostream & operator << (std::ostream &, GenEvent &);
    std::istream & operator >> (std::istream &, GenEvent &);
    std::istream & set_input_units(std::istream &, 
                                   Units::MomentumUnit, Units::LengthUnit);
    std::ostream & write_HepMC_IO_block_begin(std::ostream & );
    std::ostream & write_HepMC_IO_block_end(std::ostream & );


    // INLINE Free Functions //

    // Implemented in terms of GenEvent::use_...
    inline GenEvent& convert_units(GenEvent & evt, Units::MomentumUnit m, Units::LengthUnit l)
    {
      evt.use_units(m, l);
      return evt;
    }

    // INLINE Access Methods //

    inline int GenEvent::signal_process_id() const 
    { return m_signal_process_id; }

    inline int GenEvent::event_number() const { return m_event_number; }

    inline int GenEvent::mpi() const { return m_mpi; }

    inline double GenEvent::event_scale() const { return m_event_scale; }

    inline double GenEvent::alphaQCD() const { return m_alphaQCD; }

    inline double GenEvent::alphaQED() const { return m_alphaQED; }
 
    inline GenVertex* GenEvent::signal_process_vertex() const {
        return m_signal_process_vertex;
    }  

    inline WeightContainer& GenEvent::weights() { return m_weights; }

    inline const WeightContainer& GenEvent::weights() const 
    { return m_weights; }

    inline GenCrossSection const * GenEvent::cross_section() const 
    { return m_cross_section; }

    inline GenCrossSection*  GenEvent::cross_section()  
    { return m_cross_section; }

    inline HeavyIon const * GenEvent::heavy_ion() const 
    { return m_heavy_ion; }

    inline HeavyIon*  GenEvent::heavy_ion()  
    { return m_heavy_ion; }

    inline PdfInfo const * GenEvent::pdf_info() const 
    { return m_pdf_info; }

    inline PdfInfo*  GenEvent::pdf_info()  
    { return m_pdf_info; }

    inline const std::vector<long>& GenEvent::random_states() const 
    { return m_random_states; }

    inline void GenEvent::set_signal_process_id( int id )
    { m_signal_process_id = id; }

    inline void GenEvent::set_event_number( int eventno )
    { m_event_number = eventno; }

    inline void GenEvent::set_mpi( int nmpi )
    { m_mpi = nmpi; }


    inline void GenEvent::set_event_scale( double sc ) { m_event_scale = sc; }

    inline void GenEvent::set_alphaQCD( double a ) { m_alphaQCD = a; }

    inline void GenEvent::set_alphaQED( double a ) { m_alphaQED = a; }

    inline void GenEvent::set_signal_process_vertex( GenVertex* vtx ) {
        m_signal_process_vertex = vtx;
        if ( m_signal_process_vertex ) add_vertex( m_signal_process_vertex );
    }

    inline void GenEvent::set_cross_section( const GenCrossSection& xs )
    { 
        delete m_cross_section;
        m_cross_section = new GenCrossSection(xs); 
    }

    inline void GenEvent::set_heavy_ion( const HeavyIon& ion )
    { 
        delete m_heavy_ion;
        m_heavy_ion = new HeavyIon(ion); 
    }

    inline void GenEvent::set_pdf_info( const PdfInfo& p )
    { 
        delete m_pdf_info;
        m_pdf_info = new PdfInfo(p); 
    }

    inline void GenEvent::set_random_states( const std::vector<long>&
                                             randomstates )
    { m_random_states = randomstates; }

    inline void GenEvent::remove_barcode( GenParticle* p )
    { m_particle_barcodes.erase( p->barcode() ); }

    inline void GenEvent::remove_barcode( GenVertex* v )
    { m_vertex_barcodes.erase( v->barcode() ); }

    inline GenParticle* GenEvent::barcode_to_particle( int barCode ) const
    { 
        std::map<int,HepMC::GenParticle*>::const_iterator i 
            = m_particle_barcodes.find(barCode);
        return ( i != m_particle_barcodes.end() ) ? (*i).second : 0;
    }

    inline GenVertex* GenEvent::barcode_to_vertex( int barCode ) const
    {
        std::map<int,GenVertex*,std::greater<int> >::const_iterator i 
            = m_vertex_barcodes.find(barCode);
        return ( i != m_vertex_barcodes.end() ) ? (*i).second : 0;
    }

    inline int GenEvent::particles_size() const {
        return (int)m_particle_barcodes.size();
    }
    inline bool GenEvent::particles_empty() const {
        return (bool)m_particle_barcodes.empty();
    }
    inline int GenEvent::vertices_size() const {
        return (int)m_vertex_barcodes.size();
    }
    inline bool GenEvent::vertices_empty() const {
        return (bool)m_vertex_barcodes.empty();
    }
    
    // beam particles
    inline std::pair<HepMC::GenParticle *,HepMC::GenParticle *> GenEvent::beam_particles() const {
        return std::pair<GenParticle *,GenParticle *> (m_beam_particle_1, m_beam_particle_2);
    }

    // units
    inline Units::MomentumUnit GenEvent::momentum_unit() const {
        return m_momentum_unit; 
    }
    inline Units::LengthUnit   GenEvent::length_unit()   const {
        return m_position_unit; 
    }
    
    inline void GenEvent::use_units( Units::MomentumUnit new_m, Units::LengthUnit new_l ) { 
       use_momentum_unit( new_m );
       use_length_unit( new_l );
    }
    
    inline void GenEvent::use_units( std::string& new_m, std::string& new_l ) { 
       use_momentum_unit( new_m );
       use_length_unit( new_l );
    }

} // HepMC

#endif  // HEPMC_GEN_EVENT_H

//--------------------------------------------------------------------------

---