hepmc - Blame information for rev 179

Subversion Repositories:
Rev:
Rev Author Line No. Line
2 garren 1 //////////////////////////////////////////////////////////////////////////
2 // Matt.Dobbs@Cern.CH, September 1999
3 // Updated: 7.1.2000 iterators complete and working!
4 // Updated: 10.1.2000 GenEvent::vertex, particle iterators are made
5 //                    constant WRT this event ... note that
6 //                    GenVertex::***_iterator is not const, since it must
7 //                    be able to return a mutable pointer to itself.
8 // Updated: 08.2.2000 the event now holds a set of all attached vertices
9 //                    rather than just the roots of the graph
10 // Event record for MC generators (for use at any stage of generation)
11 //////////////////////////////////////////////////////////////////////////
12  
13 #include "HepMC/GenEvent.h"
173 garren 14 #include "HepMC/Version.h"
2 garren 15  
16 namespace HepMC {
17  
92 garren 18     GenEvent::GenEvent( int signal_process_id,
19                         int event_number,
20                         GenVertex* signal_vertex,
21                         const WeightContainer& weights,
179 garren 22                         const std::vector<long>& random_states ) :
128 garren 23         m_signal_process_id(signal_process_id),
24         m_event_number(event_number),
103 garren 25         m_mpi(-1),
128 garren 26         m_event_scale(-1),
27         m_alphaQCD(-1),
28         m_alphaQED(-1),
105 garren 29         m_signal_process_vertex(signal_vertex),
128 garren 30         m_beam_particle_1(0),
31         m_beam_particle_2(0),
105 garren 32         m_weights(weights),
128 garren 33         m_random_states(random_states),
34         m_vertex_barcodes(),
35         m_particle_barcodes(),
92 garren 36         m_heavy_ion(0),
37         m_pdf_info(0)
38     {
39         /// This constructor only allows null pointers to HeavyIon and PdfInfo
40         ///
41         /// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
42         ///       are as suggested in hep-ph/0109068, "Generic Interface..."
128 garren 43         ++s_counter;
92 garren 44     }
45  
2 garren 46     GenEvent::GenEvent( int signal_process_id, int event_number,
92 garren 47                         GenVertex* signal_vertex,
48                         const WeightContainer& weights,
179 garren 49                         const std::vector<long>& random_states,
92 garren 50                         const HeavyIon& ion,
51                         const PdfInfo& pdf ) :
128 garren 52         m_signal_process_id(signal_process_id),
53         m_event_number(event_number),
103 garren 54         m_mpi(-1),
128 garren 55         m_event_scale(-1),
56         m_alphaQCD(-1),
57         m_alphaQED(-1),
105 garren 58         m_signal_process_vertex(signal_vertex),
128 garren 59         m_beam_particle_1(0),
60         m_beam_particle_2(0),
105 garren 61         m_weights(weights),
92 garren 62         m_random_states(random_states),
128 garren 63         m_vertex_barcodes(),
64         m_particle_barcodes(),
92 garren 65         m_heavy_ion( new HeavyIon(ion) ),
66         m_pdf_info( new PdfInfo(pdf) )
2 garren 67     {
92 garren 68         /// GenEvent makes its own copy of HeavyIon and PdfInfo
69         ///
65 garren 70         /// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
71         ///       are as suggested in hep-ph/0109068, "Generic Interface..."
128 garren 72         ++s_counter;
2 garren 73     }
74  
75     GenEvent::GenEvent( const GenEvent& inevent )
147 garren 76       : m_signal_process_id    ( /* inevent.m_signal_process_id */ ),
77         m_event_number         ( /* inevent.m_event_number */ ),
78         m_mpi                  ( /* inevent.m_mpi */ ),
79         m_event_scale          ( /* inevent.m_event_scale */ ),
80         m_alphaQCD             ( /* inevent.m_alphaQCD */ ),
81         m_alphaQED             ( /* inevent.m_alphaQED */ ),
82         m_signal_process_vertex( /* inevent.m_signal_process_vertex */ ),
83         m_beam_particle_1      ( /* inevent.m_beam_particle_1 */ ),
84         m_beam_particle_2      ( /* inevent.m_beam_particle_2 */ ),
85         m_weights              ( /* inevent.m_weights */ ),
86         m_random_states        ( /* inevent.m_random_states */ ),
87         m_vertex_barcodes      ( /* inevent.m_vertex_barcodes */ ),
88         m_particle_barcodes    ( /* inevent.m_particle_barcodes */ ),
89         m_heavy_ion            ( inevent.heavy_ion() ? new HeavyIon(*inevent.heavy_ion()) : 0 ),
90         m_pdf_info             ( inevent.pdf_info() ? new PdfInfo(*inevent.pdf_info()) : 0 )
2 garren 91     {
65 garren 92         /// deep copy
128 garren 93         ++s_counter;
65 garren 94         /// deep - makes a copy of all vertices!
2 garren 95         //
147 garren 96  
97         // 1. create a NEW copy of all vertices from inevent
2 garren 98         //    taking care to map new vertices onto the vertices being copied
99         //    and add these new vertices to this event.
100         //    We do not use GenVertex::operator= because that would copy
101         //    the attached particles as well.
102         std::map<const GenVertex*,GenVertex*> map_in_to_new;
103         for ( GenEvent::vertex_const_iterator v = inevent.vertices_begin();
147 garren 104               v != inevent.vertices_end(); ++v ) {
2 garren 105             GenVertex* newvertex = new GenVertex(
106                 (*v)->position(), (*v)->id(), (*v)->weights() );
107             newvertex->suggest_barcode( (*v)->barcode() );
108             map_in_to_new[*v] = newvertex;
109             add_vertex( newvertex );
110         }
147 garren 111         // 2. copy the signal process vertex info.
2 garren 112         if ( inevent.signal_process_vertex() ) {
113             set_signal_process_vertex(
114                 map_in_to_new[inevent.signal_process_vertex()] );
115         } else set_signal_process_vertex( 0 );
116         //
117         // 3. create a NEW copy of all particles from inevent
118         //    taking care to attach them to the appropriate
128 garren 119         GenParticle* beam1(0);
120         GenParticle* beam2(0);
2 garren 121         for ( GenEvent::particle_const_iterator p = inevent.particles_begin();
147 garren 122               p != inevent.particles_end(); ++p )
2 garren 123         {
124             GenParticle* oldparticle = *p;
125             GenParticle* newparticle = new GenParticle(*oldparticle);
126             if ( oldparticle->end_vertex() ) {
127                 map_in_to_new[ oldparticle->end_vertex() ]->
128                                          add_particle_in(newparticle);
129             }
130             if ( oldparticle->production_vertex() ) {
131                 map_in_to_new[ oldparticle->production_vertex() ]->
132                                          add_particle_out(newparticle);
133             }
128 garren 134             if ( oldparticle == inevent.beam_particles().first ) beam1 = newparticle;
135             if ( oldparticle == inevent.beam_particles().second ) beam2 = newparticle;
2 garren 136         }
128 garren 137         set_beam_particles( beam1, beam2 );
2 garren 138         //
139         // 4. now that vtx/particles are copied, do everything else
140         set_signal_process_id( inevent.signal_process_id() );
141         set_event_number( inevent.event_number() );
142         set_event_scale( inevent.event_scale() );
143         set_alphaQCD( inevent.alphaQCD() );
144         set_alphaQED( inevent.alphaQED() );
128 garren 145         set_mpi( inevent.mpi() );
2 garren 146         set_random_states( inevent.random_states() );
147         weights() = inevent.weights();
147 garren 148     }
149  
150     void GenEvent::swap( GenEvent & other )
151     {
152         // if a container has a swap method, use that for improved performance
153         std::swap(m_signal_process_id    , other.m_signal_process_id    );
154         std::swap(m_event_number         , other.m_event_number         );
155         std::swap(m_mpi                  , other.m_mpi                  );
156         std::swap(m_event_scale          , other.m_event_scale          );
157         std::swap(m_alphaQCD             , other.m_alphaQCD             );
158         std::swap(m_alphaQED             , other.m_alphaQED             );
159         std::swap(m_signal_process_vertex, other.m_signal_process_vertex);
160         std::swap(m_beam_particle_1      , other.m_beam_particle_1      );
161         std::swap(m_beam_particle_2      , other.m_beam_particle_2      );
162         m_weights.swap(           other.m_weights  );
163         m_random_states.swap(     other.m_random_states  );
164         m_vertex_barcodes.swap(   other.m_vertex_barcodes );
165         m_particle_barcodes.swap( other.m_particle_barcodes );
166         std::swap(m_heavy_ion            , other.m_heavy_ion            );
167         std::swap(m_pdf_info             , other.m_pdf_info             );
168     }
169  
170     GenEvent::~GenEvent()
171     {
172         /// Deep destructor.
173         /// deletes all vertices/particles in this evt
174         ///
175         delete_all_vertices();
176         delete m_heavy_ion;
177         delete m_pdf_info;
178         --s_counter;
179     }
180  
181     GenEvent& GenEvent::operator=( const GenEvent& inevent )
182     {
183         /// best practices implementation
184         GenEvent tmp( inevent );
185         swap( tmp );
2 garren 186         return *this;
187     }
188  
189     void GenEvent::print( std::ostream& ostr ) const {
65 garren 190         /// dumps the content of this event to ostr
191         ///   to dump to cout use: event.print();
192         ///   if you want to write this event to file outfile.txt you could use:
193         ///      std::ofstream outfile("outfile.txt"); event.print( outfile );
2 garren 194         ostr << "________________________________________"
195              << "________________________________________\n";
196         ostr << "GenEvent: #" << event_number()
197              << " ID=" << signal_process_id()
198              << " SignalProcessGenVertex Barcode: "
199              << ( signal_process_vertex() ? signal_process_vertex()->barcode()
200                   : 0 )
201              << "\n";
202         ostr << " Current Memory Usage: "
203              << GenEvent::counter() << " events, "
204              << GenVertex::counter() << " vertices, "
205              << GenParticle::counter() << " particles.\n";
206         ostr << " Entries this event: " << vertices_size() << " vertices, "
207              << particles_size() << " particles.\n";
105 garren 208         if( m_beam_particle_1 && m_beam_particle_2 ) {
209           ostr << " Beam Particle barcodes: " << beam_particles().first->barcode() << " "
210                << beam_particles().second->barcode() << " \n";
211         } else {
212           ostr << " Beam Particles are not defined.\n";
213         }
2 garren 214         // Random State
215         ostr << " RndmState(" << m_random_states.size() << ")=";
216         for ( std::vector<long int>::const_iterator rs
217                   = m_random_states.begin();
147 garren 218               rs != m_random_states.end(); ++rs ) { ostr << *rs << " "; }
2 garren 219         ostr << "\n";
220         // Weights
221         ostr << " Wgts(" << weights().size() << ")=";
222         for ( WeightContainer::const_iterator wgt = weights().begin();
147 garren 223               wgt != weights().end(); ++wgt ) { ostr << *wgt << " "; }
2 garren 224         ostr << "\n";
225         ostr << " EventScale " << event_scale()
226              << " [energy] \t alphaQCD=" << alphaQCD()
227              << "\t alphaQED=" << alphaQED() << std::endl;
228         // print a legend to describe the particle info
25 garren 229         ostr << "                                    GenParticle Legend\n";
230         ostr  << "        Barcode   PDG ID      "
231               << "( Px,       Py,       Pz,     E )"
232               << " Stat  DecayVtx\n";
2 garren 233         ostr << "________________________________________"
234              << "________________________________________\n";
235         // Print all Vertices
236         for ( GenEvent::vertex_const_iterator vtx = this->vertices_begin();
237               vtx != this->vertices_end(); ++vtx ) {
238             (*vtx)->print(ostr);
239         }
240         ostr << "________________________________________"
241              << "________________________________________" << std::endl;
242     }
243  
173 garren 244     void GenEvent::print_version( std::ostream& ostr ) const {
245         ostr << "---------------------------------------------" << std::endl;
246         writeVersion( ostr );
247         ostr << "---------------------------------------------" << std::endl;
248     }
249  
2 garren 250     bool GenEvent::add_vertex( GenVertex* vtx ) {
65 garren 251         /// returns true if successful - generally will only return false
252         /// if the inserted vertex is already included in the event.
147 garren 253         if ( !vtx ) return false;
2 garren 254         // if vtx previously pointed to another GenEvent, remove it from that
255         // GenEvent's list
256         if ( vtx->parent_event() && vtx->parent_event() != this ) {
257             bool remove_status = vtx->parent_event()->remove_vertex( vtx );
258             if ( !remove_status ) {            
259                 std::cerr << "GenEvent::add_vertex ERROR "
260                           << "GenVertex::parent_event points to \n"
261                           << "an event that does not point back to the "
262                           << "GenVertex. \n This probably indicates a deeper "
263                           << "problem. " << std::endl;
264             }
265         }
266         //
267         // setting the vertex parent also inserts the vertex into this
268         // event
269         vtx->set_parent_event_( this );
270         return ( m_vertex_barcodes.count(vtx->barcode()) ? true : false );
271     }
272  
273     bool GenEvent::remove_vertex( GenVertex* vtx ) {
65 garren 274         /// this removes vtx from the event but does NOT delete it.
275         /// returns True if an entry vtx existed in the table and was erased
2 garren 276         if ( m_signal_process_vertex == vtx ) m_signal_process_vertex = 0;
277         if ( vtx->parent_event() == this ) vtx->set_parent_event_( 0 );
278         return ( m_vertex_barcodes.count(vtx->barcode()) ? false : true );
279     }
96 garren 280  
281     void GenEvent::clear()
282     {
283         /// remove all information from the event
284         /// deletes all vertices/particles in this evt
285         ///
286         delete_all_vertices();
287         delete m_heavy_ion;
288         delete m_pdf_info;
289         m_signal_process_id = 0;
105 garren 290         m_beam_particle_1 = 0;
291         m_beam_particle_2 = 0;
96 garren 292         m_event_number = 0;
293         m_event_scale = -1;
294         m_alphaQCD = -1;
295         m_alphaQED = -1;
296         m_weights = std::vector<double>();
297         m_random_states = std::vector<long int>();
298         // error check just to be safe
299         if ( m_vertex_barcodes.size() != 0
300              || m_particle_barcodes.size() != 0 ) {
301             std::cerr << "GenEvent::clear() strange result ... \n"
302                       << "either the particle and/or the vertex map isn't empty" << std::endl;
303             std::cerr << "Number vtx,particle the event after deleting = "
304                       << m_vertex_barcodes.size() << "  "
305                       << m_particle_barcodes.size() << std::endl;
306             std::cerr << "Total Number vtx,particle in memory "
307                       << "after method called = "
308                       << GenVertex::counter() << "\t"
309                       << GenParticle::counter() << std::endl;
310         }
311         return;
312     }
2 garren 313  
314     void GenEvent::delete_all_vertices() {
65 garren 315         /// deletes all vertices in the vertex container
316         /// (i.e. all vertices owned by this event)
317         /// The vertices are the "owners" of the particles, so as we delete
318         ///   the vertices, the vertex desctructors are automatically
319         ///   deleting their particles.
147 garren 320  
2 garren 321         // delete each vertex individually (this deletes particles as well)
147 garren 322         while ( !vertices_empty() ) {
2 garren 323             GenVertex* vtx = ( m_vertex_barcodes.begin() )->second;
324             m_vertex_barcodes.erase( m_vertex_barcodes.begin() );
325             delete vtx;
326         }
327         //
328         // Error checking:
329         if ( !vertices_empty() || ! particles_empty() ) {
330             std::cerr << "GenEvent::delete_all_vertices strange result ... "
331                       << "after deleting all vertices, \nthe particle and "
332                       << "vertex maps aren't empty.\n  This probably "
333                       << "indicates deeper problems or memory leak in the "
334                       << "code." << std::endl;
335             std::cerr << "Number vtx,particle the event after deleting = "
336                       << m_vertex_barcodes.size() << "  "
337                       << m_particle_barcodes.size() << std::endl;
338             std::cerr << "Total Number vtx,particle in memory "
339                       << "after method called = "
340                       << GenVertex::counter() << "\t"
341                       << GenParticle::counter() << std::endl;
342         }
343     }
344  
345     bool GenEvent::set_barcode( GenParticle* p, int suggested_barcode )
346     {
347         if ( p->parent_event() != this ) {
348             std::cerr << "GenEvent::set_barcode attempted, but the argument's"
349                       << "\n parent_event is not this ... request rejected."
350                       << std::endl;
351             return false;
352         }
353         // M.Dobbs  Nov 4, 2002
354         // First we must check to see if the particle already has a
355         // barcode which is different from the suggestion. If yes, we
356         // remove it from the particle map.
357         if ( p->barcode() != 0 && p->barcode() != suggested_barcode ) {
358             if ( m_particle_barcodes.count(p->barcode()) &&
359                  m_particle_barcodes[p->barcode()] == p ) {
360                 m_particle_barcodes.erase( p->barcode() );
361             }
362             // At this point either the particle is NOT in
363             // m_particle_barcodes, or else it is in the map, but
364             // already with the suggested barcode.
365         }
366         //
367         // First case --- a valid barcode has been suggested
368         //     (valid barcodes are numbers greater than zero)
369         bool insert_success = true;
370         if ( suggested_barcode > 0 ) {
371             if ( m_particle_barcodes.count(suggested_barcode) ) {
372                 // the suggested_barcode is already used.
373                 if ( m_particle_barcodes[suggested_barcode] == p ) {
374                     // but it was used for this particle ... so everythings ok
375                     p->set_barcode_( suggested_barcode );
376                     return true;
377                 }
378                 insert_success = false;
379                 suggested_barcode = 0;
380             } else { // suggested barcode is OK, proceed to insert
381                 m_particle_barcodes[suggested_barcode] = p;
382                 p->set_barcode_( suggested_barcode );
383                 return true;
384             }
385         }
386         //
387         // Other possibility -- a valid barcode has not been suggested,
388         //    so one is automatically generated
389         if ( suggested_barcode < 0 ) insert_success = false;
390         if ( suggested_barcode <= 0 ) {
391             if ( !m_particle_barcodes.empty() ) {
392                 // in this case we find the highest barcode that was used,
393                 // and increment it by 1
394                 suggested_barcode = m_particle_barcodes.rbegin()->first;
395                 ++suggested_barcode;
396             }
397             // For the automatically assigned barcodes, the first one
398             //   we use is 10001 ... this is just because when we read
399             //   events from HEPEVT, we will suggest their index as the
400             //   barcode, and that index has maximum value 10000.
401             //  This way we will immediately be able to recognize the hepevt
402             //   particles from the auto-assigned ones.
403             if ( suggested_barcode <= 10000 ) suggested_barcode = 10001;
404         }
405         // At this point we should have a valid barcode
406         if ( m_particle_barcodes.count(suggested_barcode) ) {
407             std::cerr << "GenEvent::set_barcode ERROR, this should never "
408                       << "happen \n report bug to matt.dobbs@cern.ch"
409                       << std::endl;
410         }
411         m_particle_barcodes[suggested_barcode] = p;
412         p->set_barcode_( suggested_barcode );
413         return insert_success;
414     }
415  
416     bool GenEvent::set_barcode( GenVertex* v, int suggested_barcode )
417     {
418         if ( v->parent_event() != this ) {
419             std::cerr << "GenEvent::set_barcode attempted, but the argument's"
420                       << "\n parent_event is not this ... request rejected."
421                       << std::endl;
422             return false;
423         }
424         // M.Dobbs Nov 4, 2002
425         // First we must check to see if the vertex already has a
426         // barcode which is different from the suggestion. If yes, we
427         // remove it from the vertex map.
428         if ( v->barcode() != 0 && v->barcode() != suggested_barcode ) {
429             if ( m_vertex_barcodes.count(v->barcode()) &&
430                  m_vertex_barcodes[v->barcode()] == v ) {
431                 m_vertex_barcodes.erase( v->barcode() );
432             }
433             // At this point either the vertex is NOT in
434             // m_vertex_barcodes, or else it is in the map, but
435             // already with the suggested barcode.
436         }
437  
438         //
439         // First case --- a valid barcode has been suggested
440         //     (valid barcodes are numbers greater than zero)
441         bool insert_success = true;
442         if ( suggested_barcode < 0 ) {
443             if ( m_vertex_barcodes.count(suggested_barcode) ) {
444                 // the suggested_barcode is already used.
445                 if ( m_vertex_barcodes[suggested_barcode] == v ) {
446                     // but it was used for this vertex ... so everythings ok
447                     v->set_barcode_( suggested_barcode );
448                     return true;
449                 }
450                 insert_success = false;
451                 suggested_barcode = 0;
452             } else { // suggested barcode is OK, proceed to insert
453                 m_vertex_barcodes[suggested_barcode] = v;
454                 v->set_barcode_( suggested_barcode );
455                 return true;
456             }
457         }
458         //
459         // Other possibility -- a valid barcode has not been suggested,
460         //    so one is automatically generated
461         if ( suggested_barcode > 0 ) insert_success = false;
462         if ( suggested_barcode >= 0 ) {
463             if ( !m_vertex_barcodes.empty() ) {
464                 // in this case we find the highest barcode that was used,
465                 // and increment it by 1, (vertex barcodes are negative)
466                 suggested_barcode = m_vertex_barcodes.rbegin()->first;
467                 --suggested_barcode;
468             }
469             if ( suggested_barcode >= 0 ) suggested_barcode = -1;
470         }
471         // At this point we should have a valid barcode
472         if ( m_vertex_barcodes.count(suggested_barcode) ) {
473             std::cerr << "GenEvent::set_barcode ERROR, this should never "
474                       << "happen \n report bug to matt.dobbs@cern.ch"
475                       << std::endl;
476         }
477         m_vertex_barcodes[suggested_barcode] = v;
478         v->set_barcode_( suggested_barcode );
479         return insert_success;
480     }
481  
105 garren 482     /// test to see if we have two valid beam particles
483     bool  GenEvent::valid_beam_particles() const {
484         bool have1 = false;
485         bool have2 = false;
486         // first check that both are defined
487         if(m_beam_particle_1 && m_beam_particle_2) {
488             // now look for a match with the particle "list"
489             for ( particle_const_iterator p = particles_begin();
147 garren 490                   p != particles_end(); ++p ) {
105 garren 491                 if( m_beam_particle_1 == *p ) have1 = true;
492                 if( m_beam_particle_2 == *p ) have2 = true;
493             }
494         }
495         if( have1 && have2 ) return true;
496         return false;
497     }
498  
499     /// construct the beam particle information using pointers to GenParticle
500     /// returns false if either GenParticle* is null
501     bool  GenEvent::set_beam_particles(GenParticle* bp1, GenParticle* bp2) {
502         m_beam_particle_1 = bp1;
503         m_beam_particle_2 = bp2;
504         if( m_beam_particle_1 && m_beam_particle_2 ) return true;
505         return false;
506     }
507  
508     /// construct the beam particle information using a std::pair of pointers to GenParticle
509     /// returns false if either GenParticle* is null
179 garren 510     bool  GenEvent::set_beam_particles(std::pair<HepMC::GenParticle*, HepMC::GenParticle*> const & bp) {
105 garren 511         return set_beam_particles(bp.first,bp.second);
512     }
513  
2 garren 514     /////////////
515     // Static  //
516     /////////////
517     unsigned int GenEvent::counter() { return s_counter; }
518     unsigned int GenEvent::s_counter = 0;
519  
520 } // HepMC
521