CLHEP/HepMC/GenEvent.h

00001 //--------------------------------------------------------------------------
00002 #ifndef HEPMC_GEN_EVENT_H
00003 #define HEPMC_GEN_EVENT_H
00004 
00006 // Matt.Dobbs@Cern.CH, September 1999, refer to:
00007 // M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
00008 // High Energy Physics", Computer Physics Communications (to be published).
00009 //
00010 // Event record for MC generators (for use at any stage of generation)
00012 //
00013 // This class is intended as both a "container class" ( to store a MC
00014 //  event for interface between MC generators and detector simulation )
00015 //  and also as a "work in progress class" ( that could be used inside
00016 //  a generator and modified as the event is built ).
00017 //
00018 // Iterators are provided which allow the user to easily obtain a
00019 //  list of particles or vertices in an event --- this list can be filled
00020 //  subject to some sort of selection criteria. Examples are given below
00021 //  ( see HepMC::copy_if and std::copy )
00023 //
00024 // L. Garren
00025 // replace std::forward_iterator with the standards compliant 
00026 // std::iterator<std::forward_iterator_tag, ...> using typedefs
00027 // 
00029 
00030 namespace HepMC {
00031     // To create a list from an iterator, use: (i.e. for a list of particles);
00032     // #include <algorithm>
00033     //     list<GenParticle*> thelist;
00034     //     copy( evt->particles_begin(), evt->particles_end(), 
00035     //           back_inserter(thelist) );
00036     // to create a list subject to a condition (predicate) use:
00037     //     list<GenParticle*> thelist;
00038     //     HepMC::copy_if( evt->particles_begin(), evt->particles_end(), 
00039     //                     back_inserter(thelist), is_photon() );
00040     // where is_photon() is a predicate like:
00041     //     class is_photon {
00042     //       public:
00043     //         bool operator() ( const GenParticle* p ) {
00044     //             if ( p && p->pdg_id() == 22 ) return 1;
00045     //             return 0;
00046     //         }
00047     //     };
00048     // which the user defines herself.
00049     template <class InputIterator, class OutputIterator, class Predicate>
00050     void copy_if( InputIterator first, InputIterator last, OutputIterator out,
00051                   Predicate pred ) {
00052         for ( ; first != last; ++first ) { if ( pred(*first) ) out = *first; }
00053     }
00054 } // HepMC
00055 
00056 // Since a container of all vertices in the event is maintained, the time
00057 //  required to loop over all vertices (or particles) is very fast -- and 
00058 //  the user does not gain much by first making his own list.
00059 //  (this is not true for the GenVertex:: versions of these iterators, which
00060 //   allow you to specify the vertex starting point and range)
00061 
00062 // Data Members:
00063 // signal_process_id()   The integer ID that uniquely specifies this signal
00064 //                       process, i.e. MSUB in Pythia. It is necessary to
00065 //                       package this with each event rather than with the run
00066 //                       because many processes may be generated within one
00067 //                       run.
00068 // event_number()        Strictly speaking we cannot think of any reason that
00069 //                       an event would need to know its own event number, it
00070 //                       is more likely something that would be assigned by
00071 //                       a database. It is included anyway (tradition?) since
00072 //                       we expect it may be useful for debugging. It can
00073 //                       be reset later by a database.
00074 // signal_process_vertex() pointer to the vertex containing the signal process
00075 // weights()             Vector of doubles which specify th weight of the evnt,
00076 //                       the first entry will be the "event weight" used for
00077 //                       hit and miss etc., but a general vector is used to
00078 //                       allow for reweighting etc. We envision a list of
00079 //                       WeightTags to be included with a run class which
00080 //                       would specify the meaning of the Weights .
00081 // random_states()       Vector of doubles which specify the random number 
00082 //                       generator's state for this event. It is left to the
00083 //                       generator to make use of this. We envision a vector of
00084 //                       RndmStatesTags to be included with a run class which
00085 //                       would specify the meaning of the random_states.
00086 //
00088 // Memory allocation //
00090 // -When a vertex (particle) is added to a event (vertex), it is "adopted" 
00091 //  and becomes the responsibility of the event (vertex) to delete that 
00092 //  particle. 
00093 // -objects responsible for deleting memory:
00094 //    -events delete included vertices
00095 //    -each vertex deletes its outgoing particles which do not have decay
00096 //     vertices
00097 //    -each vertex deletes its incoming particles which do not
00098 //     have creation vertices 
00099 //
00101 // About the Barcodes //
00103 // - each vertex or particle has a barcode, which is just an integer which
00104 //   uniquely identifies it inside the event (i.e. there is a one to one
00105 //   mapping between particle memory addresses and particle barcodes... and 
00106 //   the same applied for vertices)
00107 // - The value of a barcode has NO MEANING and NO ORDER!
00108 //   For the user's convenience, when an event is read in via an IO_method
00109 //   from an indexed list (like the HEPEVT common block), then the index will
00110 //   become the barcode for that particle.
00111 // - particle barcodes are always positive integers
00112 //   vertex barcodes are always negative integers
00113 //   The barcodes are chosen and set automatically when a vertex or particle
00114 //   comes under the ownership of an event (i.e. it is contained in an event).
00115 // - You can tell when a particle or vertex is owned, because its 
00116 //   parent_event() return value will return a pointer to the event which owns
00117 //   it (or null if its an orphan).
00118 // 
00119 
00120 #include "CLHEP/HepMC/GenVertex.h"
00121 #include "CLHEP/HepMC/GenParticle.h"
00122 #include "CLHEP/HepMC/WeightContainer.h"
00123 #include <map>
00124 #include <vector>
00125 #include <algorithm>
00126 #include "CLHEP/config/CLHEP.h"
00127 #include "CLHEP/config/iostream.h"
00128 
00129 namespace HepMC {
00130 
00135     class GenEvent {
00136         friend class GenParticle;
00137         friend class GenVertex;  
00138     public:
00139 
00140 #if defined __GNUG__ && ( __GNUG__ < 3 )
00141         typedef std::forward_iterator<GenVertex*,ptrdiff_t> forwardVertexIterType;
00142         typedef std::forward_iterator<GenParticle*,ptrdiff_t> forwardParticleIterType;
00143 #else
00144         typedef std::iterator<std::forward_iterator_tag,GenVertex*,ptrdiff_t> forwardVertexIterType;
00145         typedef std::iterator<std::forward_iterator_tag,GenParticle*,ptrdiff_t> forwardParticleIterType;
00146 #endif
00147 
00148         GenEvent( int signal_process_id = 0, int event_number = 0,
00149                   GenVertex* signal_vertex = 0,
00150                   const WeightContainer& weights = std::vector<double>(),
00151                   const std::vector<double>& randomstates
00152                   = std::vector<double>() );
00153         GenEvent( const GenEvent& inevent );          // deep copy
00154         GenEvent& operator=( const GenEvent& inevent ); // deep.
00155         virtual ~GenEvent(); //deletes all vertices/particles in this evt
00156     
00157         void print( std::ostream& ostr = std::cout ) const; // dumps to ostr
00158 
00159         GenParticle* barcode_to_particle( int barCode ) const;
00160         GenVertex*   barcode_to_vertex(   int barCode ) const;
00161 
00163         // access methods //
00165 
00166         int signal_process_id() const;
00167         int event_number() const;
00168         double event_scale() const;
00169         double alphaQCD() const;
00170         double alphaQED() const;
00171         GenVertex* signal_process_vertex() const;
00172 
00173         // direct access to the weights container is allowed. 
00174         // Thus you can user myevt.weights()[2];
00175         // to access element 2 of the weights.
00176         // or use myevt.weights().push_back( mywgt ); to add an element.
00177         // and you can set the weights with myevt.weights() = myvector;
00178         WeightContainer&        weights();
00179         const WeightContainer&  weights() const;
00180 
00181         std::vector<double> random_states() const;
00182         double              random_states( int i ) const;
00183 
00184         void set_signal_process_id( int id );
00185         void set_event_number( int eventno );
00186         void set_event_scale( double scale );
00187         void set_alphaQCD( double a );
00188         void set_alphaQED( double a );
00189         void set_signal_process_vertex( GenVertex* );
00190         void set_random_states( const std::vector<double>& randomstates );
00191 
00192         int     particles_size() const;
00193         bool    particles_empty() const;
00194         int     vertices_size() const;
00195         bool    vertices_empty() const;
00196 
00197         bool    add_vertex( GenVertex* vtx );    // adds to evt and adopts
00198         bool    remove_vertex( GenVertex* vtx ); // erases vtx from evt, 
00199 
00200         GenParticle* particle( int i ) const { return barcode_to_particle(i); }
00201 
00202     public:
00204         // vertex_iterators          //
00206         // Note:  the XXX_iterator is "resolvable" as XXX_const_iterator, but 
00207         //  not the reverse, which is consistent with STL, 
00208         //  see Musser, Derge, Saini 2ndEd. p. 69,70.
00209         class vertex_const_iterator :
00210             public forwardVertexIterType {
00211             // Iterates over all vertices in this event
00212         public:
00213             vertex_const_iterator(
00214                 const 
00215                 std::map<int,GenVertex*,std::greater<int> >::const_iterator& i)
00216                 : m_map_iterator(i) {}
00217             vertex_const_iterator() {}
00218             vertex_const_iterator( const vertex_const_iterator& i )
00219                 { *this = i; }
00220             virtual ~vertex_const_iterator() {}
00221             vertex_const_iterator&  operator=( const vertex_const_iterator& i )
00222                 { m_map_iterator = i.m_map_iterator; return *this; }
00223             GenVertex* operator*(void) const { return m_map_iterator->second; }
00224             vertex_const_iterator&  operator++(void)  //Pre-fix increment 
00225                 { ++m_map_iterator; return *this; }
00226             vertex_const_iterator   operator++(int)   //Post-fix increment
00227                 { vertex_const_iterator out(*this); ++(*this); return out; }
00228             bool  operator==( const vertex_const_iterator& a ) const
00229                 { return m_map_iterator == a.m_map_iterator; }
00230             bool  operator!=( const vertex_const_iterator& a ) const
00231                 { return !(m_map_iterator == a.m_map_iterator); }
00232         protected:
00233             std::map<int,GenVertex*,std::greater<int> >::const_iterator 
00234                                                                 m_map_iterator;
00235         };
00236         friend class vertex_const_iterator;
00237         vertex_const_iterator      vertices_begin() const
00238             { return GenEvent::vertex_const_iterator( 
00239                 m_vertex_barcodes.begin() ); }
00240         vertex_const_iterator      vertices_end() const
00241             { return GenEvent::vertex_const_iterator(
00242                 m_vertex_barcodes.end() ); }
00243 
00244         class vertex_iterator :
00245             public forwardVertexIterType {
00246             // Iterates over all vertices in this event
00247         public:
00248             vertex_iterator( 
00249                 const 
00250                 std::map<int,GenVertex*,std::greater<int> >::iterator& i )
00251                 : m_map_iterator( i ) {}
00252             vertex_iterator() {}
00253             vertex_iterator( const vertex_iterator& i ) { *this = i; }
00254             virtual ~vertex_iterator() {}
00255             vertex_iterator&  operator=( const vertex_iterator& i ) {
00256                 m_map_iterator = i.m_map_iterator;
00257                 return *this;
00258             }
00259             operator vertex_const_iterator() const
00260                 { return vertex_const_iterator(m_map_iterator); }
00261             GenVertex*        operator*(void) const
00262                 { return m_map_iterator->second; }
00263             vertex_iterator&  operator++(void)  //Pre-fix increment 
00264                 { ++m_map_iterator;     return *this; }
00265             vertex_iterator   operator++(int)   //Post-fix increment
00266                 { vertex_iterator out(*this); ++(*this); return out; }
00267             bool              operator==( const vertex_iterator& a ) const
00268                 { return m_map_iterator == a.m_map_iterator; }
00269             bool              operator!=( const vertex_iterator& a ) const
00270                 { return !(m_map_iterator == a.m_map_iterator); }
00271         protected:
00272             std::map<int,GenVertex*,std::greater<int> >::iterator 
00273                                                                m_map_iterator;
00274         };
00275         friend class vertex_iterator;
00276         vertex_iterator            vertices_begin() 
00277             { return GenEvent::vertex_iterator( 
00278                 m_vertex_barcodes.begin() ); }
00279         vertex_iterator            vertices_end()
00280             { return GenEvent::vertex_iterator(
00281                 m_vertex_barcodes.end() ); }
00282 
00283     public:
00285         // particle_iterator         //
00287         // Example of iterating over all particles in the event:
00288         //      for ( GenEvent::particle_const_iterator p = particles_begin();
00289         //            p != particles_end(); ++p ) {
00290         //         (*p)->print();
00291         //      }
00292         //
00293         class particle_const_iterator :
00294             public forwardParticleIterType {
00295             // Iterates over all vertices in this event
00296         public:
00297             particle_const_iterator(
00298                 const std::map<int,GenParticle*>::const_iterator& i )
00299                 : m_map_iterator(i) {}
00300             particle_const_iterator() {}
00301             particle_const_iterator( const particle_const_iterator& i )
00302                 { *this = i; }
00303             virtual ~particle_const_iterator() {}
00304             particle_const_iterator& operator=(
00305                 const particle_const_iterator& i )
00306                 { m_map_iterator = i.m_map_iterator; return *this; }
00307             GenParticle*        operator*(void) const
00308                 { return m_map_iterator->second; }
00309             particle_const_iterator&  operator++(void)  //Pre-fix increment 
00310                 { ++m_map_iterator; return *this; }
00311             particle_const_iterator   operator++(int)   //Post-fix increment
00312                 { particle_const_iterator out(*this); ++(*this); return out; }
00313             bool  operator==( const particle_const_iterator& a ) const
00314                 { return m_map_iterator == a.m_map_iterator; }
00315             bool  operator!=( const particle_const_iterator& a ) const
00316                 { return !(m_map_iterator == a.m_map_iterator); }
00317         protected:
00318             std::map<int,GenParticle*>::const_iterator m_map_iterator;
00319         };      
00320         friend class particle_const_iterator;
00321         particle_const_iterator      particles_begin() const
00322             { return GenEvent::particle_const_iterator( 
00323                 m_particle_barcodes.begin() ); }
00324         particle_const_iterator      particles_end() const
00325             { return GenEvent::particle_const_iterator(
00326                 m_particle_barcodes.end() ); }
00327 
00328         class particle_iterator :
00329             public forwardParticleIterType {
00330             // Iterates over all vertices in this event
00331         public:
00332             particle_iterator( const std::map<int,GenParticle*>::iterator& i )
00333                 : m_map_iterator( i ) {}
00334             particle_iterator() {}
00335             particle_iterator( const particle_iterator& i ) { *this = i; }
00336             virtual ~particle_iterator() {}
00337             particle_iterator&  operator=( const particle_iterator& i ) {
00338                 m_map_iterator = i.m_map_iterator;
00339                 return *this;
00340             }
00341             operator particle_const_iterator() const
00342                 { return particle_const_iterator(m_map_iterator); }
00343             GenParticle*        operator*(void) const
00344                 { return m_map_iterator->second; }
00345             particle_iterator&  operator++(void)  //Pre-fix increment 
00346                 { ++m_map_iterator;     return *this; }
00347             particle_iterator   operator++(int)   //Post-fix increment
00348                 { particle_iterator out(*this); ++(*this); return out; }
00349             bool              operator==( const particle_iterator& a ) const
00350                 { return m_map_iterator == a.m_map_iterator; }
00351             bool              operator!=( const particle_iterator& a ) const
00352                 { return !(m_map_iterator == a.m_map_iterator); }
00353         protected:
00354             std::map<int,GenParticle*>::iterator m_map_iterator;
00355         };
00356         friend class particle_iterator;
00357         particle_iterator particles_begin() 
00358             { return GenEvent::particle_iterator(
00359                 m_particle_barcodes.begin() ); }
00360         particle_iterator particles_end()
00361             { return GenEvent::particle_iterator(
00362                 m_particle_barcodes.end() ); }
00363 
00365     protected:
00366         //
00367         // Following methods intended for use by GenParticle/Vertex classes:
00368         // In general there is no reason they should be used elsewhere.
00369         bool         set_barcode( GenParticle* p, int suggested_barcode =0 );
00370         bool         set_barcode( GenVertex*   v, int suggested_barcode =0 );
00371         void         remove_barcode( GenParticle* p );
00372         void         remove_barcode( GenVertex*   v );
00373 
00374         static unsigned int counter(); //num GenEvent objects in memory
00375         void delete_all_vertices();
00376 
00377     private: // data members
00378         int                   m_signal_process_id;
00379         int                   m_event_number;  
00380         double                m_event_scale;// GeV, see hep-ph/0109068
00381         double                m_alphaQCD;   // QCD coupling, see hep-ph/0109068
00382         double                m_alphaQED;   // QED coupling, see hep-ph/0109068
00383         GenVertex*            m_signal_process_vertex;
00384         WeightContainer       m_weights; // weights for this event first weight
00385                                          // is used by default for hit and miss
00386         std::vector<double>   m_random_states; // container of rndm num 
00387                                                // generator states
00388 
00389         std::map< int,GenVertex*,std::greater<int> >   m_vertex_barcodes;
00390         std::map< int,GenParticle*,std::less<int> >    m_particle_barcodes;
00391 
00392         static unsigned int   s_counter;
00393     };
00394 
00396     // INLINE Access Methods //
00398 
00399     inline int GenEvent::signal_process_id() const 
00400     { return m_signal_process_id; }
00401 
00402     inline int GenEvent::event_number() const { return m_event_number; }
00403 
00404     inline double GenEvent::event_scale() const { return m_event_scale; }
00405 
00406     inline double GenEvent::alphaQCD() const { return m_alphaQCD; }
00407 
00408     inline double GenEvent::alphaQED() const { return m_alphaQED; }
00409  
00410     inline GenVertex* GenEvent::signal_process_vertex() const {
00411         // returns a (mutable) pointer to the signal process vertex
00412         return m_signal_process_vertex;
00413     }  
00414 
00415     inline WeightContainer& GenEvent::weights() { return m_weights; }
00416 
00417     inline const WeightContainer& GenEvent::weights() const 
00418     { return m_weights; }
00419 
00420     inline std::vector<double> GenEvent::random_states() const 
00421     { return m_random_states; }
00422 
00423     inline double GenEvent::random_states( int i ) const 
00424     { return m_random_states[i]; }
00425 
00426     inline void GenEvent::set_signal_process_id( int id )
00427     { m_signal_process_id = id; }
00428 
00429     inline void GenEvent::set_event_number( int eventno )
00430     { m_event_number = eventno; }
00431 
00432 
00433     inline void GenEvent::set_event_scale( double sc ) { m_event_scale = sc; }
00434 
00435     inline void GenEvent::set_alphaQCD( double a ) { m_alphaQCD = a; }
00436 
00437     inline void GenEvent::set_alphaQED( double a ) { m_alphaQED = a; }
00438 
00439     inline void GenEvent::set_signal_process_vertex( GenVertex* vtx ) {
00440         m_signal_process_vertex = vtx;
00441         if ( m_signal_process_vertex ) add_vertex( m_signal_process_vertex );
00442     }
00443 
00444     inline void GenEvent::set_random_states( const std::vector<double>&
00445                                              randomstates )
00446     { m_random_states = randomstates; }
00447 
00448     inline void GenEvent::remove_barcode( GenParticle* p )
00449     { m_particle_barcodes.erase( p->barcode() ); }
00450 
00451     inline void GenEvent::remove_barcode( GenVertex* v )
00452     { m_vertex_barcodes.erase( v->barcode() ); }
00453 
00454     inline GenParticle* GenEvent::barcode_to_particle( int barCode ) const
00455     { 
00456         std::map<int,GenParticle*>::const_iterator i 
00457             = m_particle_barcodes.find(barCode);
00458         return ( i != m_particle_barcodes.end() ) ? (*i).second : 0;
00459     }
00460 
00461     inline GenVertex* GenEvent::barcode_to_vertex( int barCode ) const
00462     {
00463         std::map<int,GenVertex*,std::greater<int> >::const_iterator i 
00464             = m_vertex_barcodes.find(barCode);
00465         return ( i != m_vertex_barcodes.end() ) ? (*i).second : 0;
00466     }
00467 
00468     inline int GenEvent::particles_size() const {
00469         return (int)m_particle_barcodes.size();
00470     }
00471     inline bool GenEvent::particles_empty() const {
00472         return (bool)m_particle_barcodes.empty();
00473     }
00474     inline int GenEvent::vertices_size() const {
00475         return (int)m_vertex_barcodes.size();
00476     }
00477     inline bool GenEvent::vertices_empty() const {
00478         return (bool)m_vertex_barcodes.empty();
00479     }
00480 
00481 } // HepMC
00482 
00483 #endif  // HEPMC_GEN_EVENT_H
00484 
00485 //--------------------------------------------------------------------------
00486 
00487 

Class Library for High Energy Physics (version 1.8)