hepmc - Blame information for rev 96

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