hepmc - Blame information for rev 128

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