hepmc - Blame information for rev 238

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 )
238 garren 76       : m_signal_process_id    ( inevent.signal_process_id() ),
77         m_event_number         ( inevent.event_number() ),
78         m_mpi                  ( inevent.mpi() ),
79         m_event_scale          ( inevent.event_scale() ),
80         m_alphaQCD             ( inevent.alphaQCD() ),
81         m_alphaQED             ( inevent.alphaQED() ),
147 garren 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
238 garren 118         //    taking care to attach them to the appropriate vertex
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         //
238 garren 139         // 4. now that vtx/particles are copied, copy weights and random states
2 garren 140         set_random_states( inevent.random_states() );
141         weights() = inevent.weights();
147 garren 142     }
143  
144     void GenEvent::swap( GenEvent & other )
145     {
146         // if a container has a swap method, use that for improved performance
147         std::swap(m_signal_process_id    , other.m_signal_process_id    );
148         std::swap(m_event_number         , other.m_event_number         );
149         std::swap(m_mpi                  , other.m_mpi                  );
150         std::swap(m_event_scale          , other.m_event_scale          );
151         std::swap(m_alphaQCD             , other.m_alphaQCD             );
152         std::swap(m_alphaQED             , other.m_alphaQED             );
153         std::swap(m_signal_process_vertex, other.m_signal_process_vertex);
154         std::swap(m_beam_particle_1      , other.m_beam_particle_1      );
155         std::swap(m_beam_particle_2      , other.m_beam_particle_2      );
156         m_weights.swap(           other.m_weights  );
157         m_random_states.swap(     other.m_random_states  );
158         m_vertex_barcodes.swap(   other.m_vertex_barcodes );
159         m_particle_barcodes.swap( other.m_particle_barcodes );
160         std::swap(m_heavy_ion            , other.m_heavy_ion            );
161         std::swap(m_pdf_info             , other.m_pdf_info             );
162     }
163  
164     GenEvent::~GenEvent()
165     {
166         /// Deep destructor.
167         /// deletes all vertices/particles in this evt
168         ///
169         delete_all_vertices();
170         delete m_heavy_ion;
171         delete m_pdf_info;
172         --s_counter;
173     }
174  
175     GenEvent& GenEvent::operator=( const GenEvent& inevent )
176     {
177         /// best practices implementation
178         GenEvent tmp( inevent );
179         swap( tmp );
2 garren 180         return *this;
181     }
182  
183     void GenEvent::print( std::ostream& ostr ) const {
65 garren 184         /// dumps the content of this event to ostr
185         ///   to dump to cout use: event.print();
186         ///   if you want to write this event to file outfile.txt you could use:
187         ///      std::ofstream outfile("outfile.txt"); event.print( outfile );
2 garren 188         ostr << "________________________________________"
189              << "________________________________________\n";
190         ostr << "GenEvent: #" << event_number()
191              << " ID=" << signal_process_id()
192              << " SignalProcessGenVertex Barcode: "
193              << ( signal_process_vertex() ? signal_process_vertex()->barcode()
194                   : 0 )
195              << "\n";
196         ostr << " Current Memory Usage: "
197              << GenEvent::counter() << " events, "
198              << GenVertex::counter() << " vertices, "
199              << GenParticle::counter() << " particles.\n";
200         ostr << " Entries this event: " << vertices_size() << " vertices, "
201              << particles_size() << " particles.\n";
105 garren 202         if( m_beam_particle_1 && m_beam_particle_2 ) {
203           ostr << " Beam Particle barcodes: " << beam_particles().first->barcode() << " "
204                << beam_particles().second->barcode() << " \n";
205         } else {
206           ostr << " Beam Particles are not defined.\n";
207         }
2 garren 208         // Random State
209         ostr << " RndmState(" << m_random_states.size() << ")=";
210 garren 210         for ( std::vector<long>::const_iterator rs
2 garren 211                   = m_random_states.begin();
147 garren 212               rs != m_random_states.end(); ++rs ) { ostr << *rs << " "; }
2 garren 213         ostr << "\n";
214         // Weights
215         ostr << " Wgts(" << weights().size() << ")=";
216         for ( WeightContainer::const_iterator wgt = weights().begin();
147 garren 217               wgt != weights().end(); ++wgt ) { ostr << *wgt << " "; }
2 garren 218         ostr << "\n";
219         ostr << " EventScale " << event_scale()
220              << " [energy] \t alphaQCD=" << alphaQCD()
221              << "\t alphaQED=" << alphaQED() << std::endl;
222         // print a legend to describe the particle info
25 garren 223         ostr << "                                    GenParticle Legend\n";
224         ostr  << "        Barcode   PDG ID      "
225               << "( Px,       Py,       Pz,     E )"
226               << " Stat  DecayVtx\n";
2 garren 227         ostr << "________________________________________"
228              << "________________________________________\n";
229         // Print all Vertices
230         for ( GenEvent::vertex_const_iterator vtx = this->vertices_begin();
231               vtx != this->vertices_end(); ++vtx ) {
232             (*vtx)->print(ostr);
233         }
234         ostr << "________________________________________"
235              << "________________________________________" << std::endl;
236     }
237  
173 garren 238     void GenEvent::print_version( std::ostream& ostr ) const {
239         ostr << "---------------------------------------------" << std::endl;
240         writeVersion( ostr );
241         ostr << "---------------------------------------------" << std::endl;
242     }
243  
2 garren 244     bool GenEvent::add_vertex( GenVertex* vtx ) {
65 garren 245         /// returns true if successful - generally will only return false
246         /// if the inserted vertex is already included in the event.
147 garren 247         if ( !vtx ) return false;
2 garren 248         // if vtx previously pointed to another GenEvent, remove it from that
249         // GenEvent's list
250         if ( vtx->parent_event() && vtx->parent_event() != this ) {
251             bool remove_status = vtx->parent_event()->remove_vertex( vtx );
252             if ( !remove_status ) {            
253                 std::cerr << "GenEvent::add_vertex ERROR "
254                           << "GenVertex::parent_event points to \n"
255                           << "an event that does not point back to the "
256                           << "GenVertex. \n This probably indicates a deeper "
257                           << "problem. " << std::endl;
258             }
259         }
260         //
261         // setting the vertex parent also inserts the vertex into this
262         // event
263         vtx->set_parent_event_( this );
264         return ( m_vertex_barcodes.count(vtx->barcode()) ? true : false );
265     }
266  
267     bool GenEvent::remove_vertex( GenVertex* vtx ) {
65 garren 268         /// this removes vtx from the event but does NOT delete it.
269         /// returns True if an entry vtx existed in the table and was erased
2 garren 270         if ( m_signal_process_vertex == vtx ) m_signal_process_vertex = 0;
271         if ( vtx->parent_event() == this ) vtx->set_parent_event_( 0 );
272         return ( m_vertex_barcodes.count(vtx->barcode()) ? false : true );
273     }
96 garren 274  
275     void GenEvent::clear()
276     {
277         /// remove all information from the event
278         /// deletes all vertices/particles in this evt
279         ///
280         delete_all_vertices();
281         delete m_heavy_ion;
282         delete m_pdf_info;
283         m_signal_process_id = 0;
105 garren 284         m_beam_particle_1 = 0;
285         m_beam_particle_2 = 0;
96 garren 286         m_event_number = 0;
287         m_event_scale = -1;
288         m_alphaQCD = -1;
289         m_alphaQED = -1;
290         m_weights = std::vector<double>();
210 garren 291         m_random_states = std::vector<long>();
96 garren 292         // error check just to be safe
293         if ( m_vertex_barcodes.size() != 0
294              || m_particle_barcodes.size() != 0 ) {
295             std::cerr << "GenEvent::clear() strange result ... \n"
296                       << "either the particle and/or the vertex map isn't empty" << std::endl;
297             std::cerr << "Number vtx,particle the event after deleting = "
298                       << m_vertex_barcodes.size() << "  "
299                       << m_particle_barcodes.size() << std::endl;
300             std::cerr << "Total Number vtx,particle in memory "
301                       << "after method called = "
302                       << GenVertex::counter() << "\t"
303                       << GenParticle::counter() << std::endl;
304         }
305         return;
306     }
2 garren 307  
308     void GenEvent::delete_all_vertices() {
65 garren 309         /// deletes all vertices in the vertex container
310         /// (i.e. all vertices owned by this event)
311         /// The vertices are the "owners" of the particles, so as we delete
312         ///   the vertices, the vertex desctructors are automatically
313         ///   deleting their particles.
147 garren 314  
2 garren 315         // delete each vertex individually (this deletes particles as well)
147 garren 316         while ( !vertices_empty() ) {
2 garren 317             GenVertex* vtx = ( m_vertex_barcodes.begin() )->second;
318             m_vertex_barcodes.erase( m_vertex_barcodes.begin() );
319             delete vtx;
320         }
321         //
322         // Error checking:
323         if ( !vertices_empty() || ! particles_empty() ) {
324             std::cerr << "GenEvent::delete_all_vertices strange result ... "
325                       << "after deleting all vertices, \nthe particle and "
326                       << "vertex maps aren't empty.\n  This probably "
327                       << "indicates deeper problems or memory leak in the "
328                       << "code." << std::endl;
329             std::cerr << "Number vtx,particle the event after deleting = "
330                       << m_vertex_barcodes.size() << "  "
331                       << m_particle_barcodes.size() << std::endl;
332             std::cerr << "Total Number vtx,particle in memory "
333                       << "after method called = "
334                       << GenVertex::counter() << "\t"
335                       << GenParticle::counter() << std::endl;
336         }
337     }
338  
339     bool GenEvent::set_barcode( GenParticle* p, int suggested_barcode )
340     {
341         if ( p->parent_event() != this ) {
342             std::cerr << "GenEvent::set_barcode attempted, but the argument's"
343                       << "\n parent_event is not this ... request rejected."
344                       << std::endl;
345             return false;
346         }
347         // M.Dobbs  Nov 4, 2002
348         // First we must check to see if the particle already has a
349         // barcode which is different from the suggestion. If yes, we
350         // remove it from the particle map.
351         if ( p->barcode() != 0 && p->barcode() != suggested_barcode ) {
352             if ( m_particle_barcodes.count(p->barcode()) &&
353                  m_particle_barcodes[p->barcode()] == p ) {
354                 m_particle_barcodes.erase( p->barcode() );
355             }
356             // At this point either the particle is NOT in
357             // m_particle_barcodes, or else it is in the map, but
358             // already with the suggested barcode.
359         }
360         //
361         // First case --- a valid barcode has been suggested
362         //     (valid barcodes are numbers greater than zero)
363         bool insert_success = true;
364         if ( suggested_barcode > 0 ) {
365             if ( m_particle_barcodes.count(suggested_barcode) ) {
366                 // the suggested_barcode is already used.
367                 if ( m_particle_barcodes[suggested_barcode] == p ) {
368                     // but it was used for this particle ... so everythings ok
369                     p->set_barcode_( suggested_barcode );
370                     return true;
371                 }
372                 insert_success = false;
373                 suggested_barcode = 0;
374             } else { // suggested barcode is OK, proceed to insert
375                 m_particle_barcodes[suggested_barcode] = p;
376                 p->set_barcode_( suggested_barcode );
377                 return true;
378             }
379         }
380         //
381         // Other possibility -- a valid barcode has not been suggested,
382         //    so one is automatically generated
383         if ( suggested_barcode < 0 ) insert_success = false;
384         if ( suggested_barcode <= 0 ) {
385             if ( !m_particle_barcodes.empty() ) {
386                 // in this case we find the highest barcode that was used,
387                 // and increment it by 1
388                 suggested_barcode = m_particle_barcodes.rbegin()->first;
389                 ++suggested_barcode;
390             }
391             // For the automatically assigned barcodes, the first one
392             //   we use is 10001 ... this is just because when we read
393             //   events from HEPEVT, we will suggest their index as the
394             //   barcode, and that index has maximum value 10000.
395             //  This way we will immediately be able to recognize the hepevt
396             //   particles from the auto-assigned ones.
397             if ( suggested_barcode <= 10000 ) suggested_barcode = 10001;
398         }
399         // At this point we should have a valid barcode
400         if ( m_particle_barcodes.count(suggested_barcode) ) {
401             std::cerr << "GenEvent::set_barcode ERROR, this should never "
402                       << "happen \n report bug to matt.dobbs@cern.ch"
403                       << std::endl;
404         }
405         m_particle_barcodes[suggested_barcode] = p;
406         p->set_barcode_( suggested_barcode );
407         return insert_success;
408     }
409  
410     bool GenEvent::set_barcode( GenVertex* v, int suggested_barcode )
411     {
412         if ( v->parent_event() != this ) {
413             std::cerr << "GenEvent::set_barcode attempted, but the argument's"
414                       << "\n parent_event is not this ... request rejected."
415                       << std::endl;
416             return false;
417         }
418         // M.Dobbs Nov 4, 2002
419         // First we must check to see if the vertex already has a
420         // barcode which is different from the suggestion. If yes, we
421         // remove it from the vertex map.
422         if ( v->barcode() != 0 && v->barcode() != suggested_barcode ) {
423             if ( m_vertex_barcodes.count(v->barcode()) &&
424                  m_vertex_barcodes[v->barcode()] == v ) {
425                 m_vertex_barcodes.erase( v->barcode() );
426             }
427             // At this point either the vertex is NOT in
428             // m_vertex_barcodes, or else it is in the map, but
429             // already with the suggested barcode.
430         }
431  
432         //
433         // First case --- a valid barcode has been suggested
434         //     (valid barcodes are numbers greater than zero)
435         bool insert_success = true;
436         if ( suggested_barcode < 0 ) {
437             if ( m_vertex_barcodes.count(suggested_barcode) ) {
438                 // the suggested_barcode is already used.
439                 if ( m_vertex_barcodes[suggested_barcode] == v ) {
440                     // but it was used for this vertex ... so everythings ok
441                     v->set_barcode_( suggested_barcode );
442                     return true;
443                 }
444                 insert_success = false;
445                 suggested_barcode = 0;
446             } else { // suggested barcode is OK, proceed to insert
447                 m_vertex_barcodes[suggested_barcode] = v;
448                 v->set_barcode_( suggested_barcode );
449                 return true;
450             }
451         }
452         //
453         // Other possibility -- a valid barcode has not been suggested,
454         //    so one is automatically generated
455         if ( suggested_barcode > 0 ) insert_success = false;
456         if ( suggested_barcode >= 0 ) {
457             if ( !m_vertex_barcodes.empty() ) {
458                 // in this case we find the highest barcode that was used,
459                 // and increment it by 1, (vertex barcodes are negative)
460                 suggested_barcode = m_vertex_barcodes.rbegin()->first;
461                 --suggested_barcode;
462             }
463             if ( suggested_barcode >= 0 ) suggested_barcode = -1;
464         }
465         // At this point we should have a valid barcode
466         if ( m_vertex_barcodes.count(suggested_barcode) ) {
467             std::cerr << "GenEvent::set_barcode ERROR, this should never "
468                       << "happen \n report bug to matt.dobbs@cern.ch"
469                       << std::endl;
470         }
471         m_vertex_barcodes[suggested_barcode] = v;
472         v->set_barcode_( suggested_barcode );
473         return insert_success;
474     }
475  
105 garren 476     /// test to see if we have two valid beam particles
477     bool  GenEvent::valid_beam_particles() const {
478         bool have1 = false;
479         bool have2 = false;
480         // first check that both are defined
481         if(m_beam_particle_1 && m_beam_particle_2) {
482             // now look for a match with the particle "list"
483             for ( particle_const_iterator p = particles_begin();
147 garren 484                   p != particles_end(); ++p ) {
105 garren 485                 if( m_beam_particle_1 == *p ) have1 = true;
486                 if( m_beam_particle_2 == *p ) have2 = true;
487             }
488         }
489         if( have1 && have2 ) return true;
490         return false;
491     }
492  
493     /// construct the beam particle information using pointers to GenParticle
494     /// returns false if either GenParticle* is null
495     bool  GenEvent::set_beam_particles(GenParticle* bp1, GenParticle* bp2) {
496         m_beam_particle_1 = bp1;
497         m_beam_particle_2 = bp2;
498         if( m_beam_particle_1 && m_beam_particle_2 ) return true;
499         return false;
500     }
501  
502     /// construct the beam particle information using a std::pair of pointers to GenParticle
503     /// returns false if either GenParticle* is null
179 garren 504     bool  GenEvent::set_beam_particles(std::pair<HepMC::GenParticle*, HepMC::GenParticle*> const & bp) {
105 garren 505         return set_beam_particles(bp.first,bp.second);
506     }
507  
2 garren 508     /////////////
509     // Static  //
510     /////////////
511     unsigned int GenEvent::counter() { return s_counter; }
512     unsigned int GenEvent::s_counter = 0;
513  
514 } // HepMC