HepMC event record
GenEvent.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2015 The HepMC collaboration (see AUTHORS for details)
5 //
6 /**
7  * @file GenEvent.cc
8  * @brief Implementation of \b class GenEvent
9  *
10  */
11 #include "HepMC/GenEvent.h"
12 #include "HepMC/GenParticle.h"
13 #include "HepMC/GenVertex.h"
14 
15 #include "HepMC/Data/GenEventData.h"
16 #include "HepMC/Search/FindParticles.h"
17 
18 #include <deque>
19 #include <algorithm> // sort
20 using namespace std;
21 
22 namespace HepMC {
23 
24 
25 GenEvent::GenEvent(Units::MomentumUnit mu,
27  : m_event_number(0), m_momentum_unit(mu),
28  m_length_unit(lu),
29  m_rootvertex(make_shared<GenVertex>()) {}
30 
31 
32 GenEvent::GenEvent(shared_ptr<GenRunInfo> run,
35  : m_event_number(0), m_momentum_unit(mu),
36  m_length_unit(lu),
37  m_rootvertex(make_shared<GenVertex>()),
38  m_run_info(run) {}
39 
40 
41 // void GenEvent::add_particle( const GenParticlePtr &p ) {
43  if( p->in_event() ) return;
44 
45  m_particles.push_back(p);
46 
47  p->m_event = this;
48  p->m_id = particles().size();
49 
50  // Particles without production vertex are added to the root vertex
51  if( !p->production_vertex() )
52  m_rootvertex->add_particle_out(p);
53 }
54 
55 
56 // void GenEvent::add_vertex( const GenVertexPtr &v ) {
58  if( v->in_event() ) return;
59 
60  m_vertices.push_back(v);
61 
62  v->m_event = this;
63  v->m_id = -(int)vertices().size();
64 
65  // Add all incoming and outgoing particles and restore their production/end vertices
66  FOREACH( GenParticlePtr &p, v->m_particles_in ) {
67  if(!p->in_event()) add_particle(p);
68  p->m_end_vertex = v->m_this.lock();
69  }
70 
71  FOREACH( GenParticlePtr &p, v->m_particles_out ) {
72  if(!p->in_event()) add_particle(p);
73  p->m_production_vertex = v->m_this.lock();
74  }
75 }
76 
77 
79  if( !p || p->parent_event() != this ) return;
80 
81  DEBUG( 30, "GenEvent::remove_particle - called with particle: "<<p->id() );
82  GenVertexPtr end_vtx = p->end_vertex();
83  if( end_vtx ) {
84  end_vtx->remove_particle_in(p);
85 
86  // If that was the only incoming particle, remove vertex from the event
87  if( end_vtx->particles_in().size() == 0 ) remove_vertex(end_vtx);
88  }
89 
90  GenVertexPtr prod_vtx = p->production_vertex();
91  if( prod_vtx ) {
92  prod_vtx->remove_particle_out(p);
93 
94  // If that was the only outgoing particle, remove vertex from the event
95  if( prod_vtx->particles_out().size() == 0 ) remove_vertex(prod_vtx);
96  }
97 
98  DEBUG( 30, "GenEvent::remove_particle - erasing particle: " << p->id() )
99 
100  int idx = p->id();
101  vector<GenParticlePtr>::iterator it = m_particles.erase(m_particles.begin() + idx-1 );
102 
103  // Remove attributes of this particle
104  vector<string> atts = p->attribute_names();
105  FOREACH( string s, atts) {
106  p->remove_attribute(s);
107  }
108 
109  //
110  // Reassign id of attributes with id above this one
111  //
112 #ifndef HEPMC_HAS_CXX0X_GCC_ONLY
113  vector<att_val_t> changed_attributes;
114 
115  FOREACH( att_key_t& vt1, m_attributes ) {
116  changed_attributes.clear();
117 
118  FOREACH( const att_val_t& vt2, vt1.second ) {
119  if( vt2.first > p->id() ) {
120  changed_attributes.push_back(vt2);
121  }
122  }
123 
124  FOREACH( att_val_t val, changed_attributes ) {
125  vt1.second.erase(val.first);
126  vt1.second[val.first-1] = val.second;
127  }
128  }
129 #else
130  vector< pair< int, shared_ptr<Attribute> > > changed_attributes;
131 
132  FOREACH( att_key_t& vt1, m_attributes ) {
133  changed_attributes.clear();
134 
135  for ( std::map<int, shared_ptr<Attribute> >::iterator vt2=vt1.second.begin();vt2!=vt1.second.end();vt2++){
136  if( (*vt2).first > p->id() ) {
137  changed_attributes.push_back(*vt2);
138  }
139  }
140 
141  FOREACH( att_val_t val, changed_attributes ) {
142  vt1.second.erase(val.first);
143  vt1.second[val.first-1] = val.second;
144  }
145  }
146 #endif
147  // Reassign id of particles with id above this one
148  for(;it != m_particles.end(); ++it) {
149  --((*it)->m_id);
150  }
151 
152  // Finally - set parent event and id of this particle to 0
153  p->m_event = NULL;
154  p->m_id = 0;
155 }
156 
157  struct sort_by_id_asc {
158  inline bool operator()(const GenParticlePtr& p1, const GenParticlePtr& p2) {
159  return (p1->id() > p2->id());
160  }
161  };
162 
163 void GenEvent::remove_particles( vector<GenParticlePtr> v ) {
164  /// @todo Currently the only optimization is sort by id in ascending order.
165  /// Needs better optimization!
166 
167  sort( v.begin(), v.end(), sort_by_id_asc() );
168 
169  FOREACH( GenParticlePtr &p, v ) {
170  remove_particle(p);
171  }
172 }
173 
175  if( !v || v->parent_event() != this ) return;
176 
177  DEBUG( 30, "GenEvent::remove_vertex - called with vertex: "<<v->id() );
178  shared_ptr<GenVertex> null_vtx;
179 
180  FOREACH( GenParticlePtr &p, v->m_particles_in ) {
181  p->m_end_vertex.reset();
182  }
183 
184  FOREACH( GenParticlePtr &p, v->m_particles_out ) {
185  p->m_production_vertex.reset();
186 
187  // recursive delete rest of the tree
188  remove_particle(p);
189  }
190 
191  // Erase this vertex from vertices list
192  DEBUG( 30, "GenEvent::remove_vertex - erasing vertex: " << v->id() )
193 
194  int idx = -v->id();
195  vector<GenVertexPtr>::iterator it = m_vertices.erase(m_vertices.begin() + idx-1 );
196 
197  // Remove attributes of this vertex
198  vector<string> atts = v->attribute_names();
199  FOREACH( string s, atts) {
200  v->remove_attribute(s);
201  }
202 
203  //
204  // Reassign id of attributes with id below this one
205  //
206 #ifndef HEPMC_HAS_CXX0X_GCC_ONLY
207  vector<att_val_t> changed_attributes;
208 
209  FOREACH( att_key_t& vt1, m_attributes ) {
210  changed_attributes.clear();
211 
212  FOREACH( const att_val_t& vt2, vt1.second ) {
213  if( vt2.first < v->id() ) {
214  changed_attributes.push_back(vt2);
215  }
216  }
217 
218  FOREACH( att_val_t val, changed_attributes ) {
219  vt1.second.erase(val.first);
220  vt1.second[val.first+1] = val.second;
221  }
222  }
223 
224  // Reassign id of particles with id above this one
225  for(;it != m_vertices.end(); ++it) {
226  ++((*it)->m_id);
227  }
228 
229  // Finally - set parent event and id of this vertex to 0
230  v->m_event = NULL;
231  v->m_id = 0;
232 }
233 #else
234  vector< pair< int, shared_ptr<Attribute> > > changed_attributes;
235 
236  FOREACH( att_key_t& vt1, m_attributes ) {
237  changed_attributes.clear();
238 
239  for ( std::map<int, shared_ptr<Attribute> >::iterator vt2=vt1.second.begin();vt2!=vt1.second.end();vt2++){
240  if( (*vt2).first < v->id() ) {
241  changed_attributes.push_back(*vt2);
242  }
243  }
244 
245  FOREACH( att_val_t val, changed_attributes ) {
246  vt1.second.erase(val.first);
247  vt1.second[val.first+1] = val.second;
248  }
249  }
250 
251  // Reassign id of particles with id above this one
252  for(;it != m_vertices.end(); ++it) {
253  ++((*it)->m_id);
254  }
255 
256  // Finally - set parent event and id of this vertex to 0
257  v->m_event = NULL;
258  v->m_id = 0;
259 }
260 #endif
261 
262 void GenEvent::add_tree( const vector<GenParticlePtr> &parts ) {
263 
264  deque<GenVertexPtr> sorting;
265 
266  // Find all starting vertices (end vertex of particles that have no production vertex)
267  FOREACH( const GenParticlePtr &p, parts ) {
268  const GenVertexPtr &v = p->production_vertex();
269  if( !v || v->particles_in().size()==0 ) {
270  const GenVertexPtr &v2 = p->end_vertex();
271  if(v2) sorting.push_back(v2);
272  }
273  }
274 
275  DEBUG_CODE_BLOCK(
276  unsigned int sorting_loop_count = 0;
277  unsigned int max_deque_size = 0;
278  )
279 
280  // Add vertices to the event in topological order
281  while( !sorting.empty() ) {
282  DEBUG_CODE_BLOCK(
283  if( sorting.size() > max_deque_size ) max_deque_size = sorting.size();
284  ++sorting_loop_count;
285  )
286 
287  GenVertexPtr &v = sorting.front();
288 
289  bool added = false;
290 
291  // Add all mothers to the front of the list
292  FOREACH( const GenParticlePtr &p, v->particles_in() ) {
293  GenVertexPtr v2 = p->production_vertex();
294  if( v2 && !v2->in_event() ) {
295  sorting.push_front(v2);
296  added = true;
297  }
298  }
299 
300  // If we have added at least one production vertex,
301  // our vertex is not the first one on the list
302  if( added ) continue;
303 
304  // If vertex not yet added
305  if( !v->in_event() ) {
306 
307  add_vertex(v);
308 
309  // Add all end vertices to the end of the list
310  FOREACH( const GenParticlePtr &p, v->particles_out() ) {
311  GenVertexPtr v2 = p->end_vertex();
312  if( v2 && !v2->in_event() ) {
313  sorting.push_back(v2);
314  }
315  }
316  }
317 
318  sorting.pop_front();
319  }
320 
321  DEBUG_CODE_BLOCK(
322  DEBUG( 6, "GenEvent - particles sorted: "
323  <<this->particles().size()<<", max deque size: "
324  <<max_deque_size<<", iterations: "<<sorting_loop_count )
325  )
326 }
327 
328 
329 void GenEvent::reserve(unsigned int parts, unsigned int verts) {
330  m_particles.reserve(parts);
331  m_vertices.reserve(verts);
332 }
333 
334 
335 void GenEvent::set_units( Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit) {
336  if( new_momentum_unit != m_momentum_unit ) {
337  FOREACH( GenParticlePtr &p, m_particles ) {
338  Units::convert( p->m_data.momentum, m_momentum_unit, new_momentum_unit );
339  }
340 
341  m_momentum_unit = new_momentum_unit;
342  }
343 
344  if( new_length_unit != m_length_unit ) {
345  FOREACH( GenVertexPtr &v, m_vertices ) {
346  FourVector &fv = v->m_data.position;
347  if( !fv.is_zero() ) Units::convert( fv, m_length_unit, new_length_unit );
348  }
349 
350  m_length_unit = new_length_unit;
351  }
352 }
353 
354 
355 const FourVector& GenEvent::event_pos() const {
356  return m_rootvertex->data().position;
357 }
358 
359 const vector<GenParticlePtr>& GenEvent::beams() const {
360  return m_rootvertex->particles_out();
361 }
362 
363 void GenEvent::shift_position_by( const FourVector & delta ) {
364  m_rootvertex->set_position( event_pos() + delta );
365 
366  // Offset all vertices
367  FOREACH ( GenVertexPtr &v, m_vertices ) {
368  if ( v->has_set_position() )
369  v->set_position( v->position() + delta );
370  }
371 }
372 
373 
374 void GenEvent::clear() {
375  m_event_number = 0;
376  m_rootvertex = make_shared<GenVertex>();
377  m_weights.clear();
378  m_attributes.clear();
379  m_particles.clear();
380  m_vertices.clear();
381 }
382 
383 
384 //
385 // Deprecated functions
386 //
387 
388 void GenEvent::add_particle( GenParticle *p ) {
389  add_particle( GenParticlePtr(p) );
390 }
391 
392 
393 void GenEvent::add_vertex( GenVertex *v ) {
394  add_vertex( GenVertexPtr(v) );
395 }
396 
397 
398 void GenEvent::remove_attribute(const string &name, int id) {
399  map< string, map<int, shared_ptr<Attribute> > >::iterator i1 = m_attributes.find(name);
400  if( i1 == m_attributes.end() ) return;
401 
402  map<int, shared_ptr<Attribute> >::iterator i2 = i1->second.find(id);
403  if( i2 == i1->second.end() ) return;
404 
405  i1->second.erase(i2);
406 }
407 
408 vector<string> GenEvent::attribute_names(int id) const {
409  vector<string> results;
410 
411  FOREACH( const att_key_t& vt1, this->attributes() ) {
412  FOREACH( const att_val_t& vt2, vt1.second ) {
413  if( vt2.first == id ) results.push_back( vt1.first );
414  }
415  }
416 
417  return results;
418 }
419 
420 void GenEvent::write_data(GenEventData& data) const {
421  // Reserve memory for containers
422  data.particles.reserve( this->particles().size() );
423  data.vertices.reserve( this->vertices().size() );
424  data.links1.reserve( this->particles().size()*2 );
425  data.links2.reserve( this->particles().size()*2 );
426  data.attribute_id.reserve( this->attributes().size() );
427  data.attribute_name.reserve( this->attributes().size() );
428  data.attribute_string.reserve( this->attributes().size() );
429 
430  // Fill event data
431  data.event_number = this->event_number();
432  data.momentum_unit = this->momentum_unit();
433  data.length_unit = this->length_unit();
434  data.event_pos = this->event_pos();
435 
436  // Fill containers
437  data.weights = this->weights();
438 
439  FOREACH( const GenParticlePtr &p, this->particles() ) {
440  data.particles.push_back( p->data() );
441  }
442 
443  FOREACH( const GenVertexPtr &v, this->vertices() ) {
444  data.vertices.push_back( v->data() );
445  int v_id = v->id();
446 
447  FOREACH( const GenParticlePtr &p, v->particles_in() ) {
448  data.links1.push_back( p->id() );
449  data.links2.push_back( v_id );
450  }
451 
452  FOREACH( const GenParticlePtr &p, v->particles_out() ) {
453  data.links1.push_back( v_id );
454  data.links2.push_back( p->id() );
455  }
456  }
457 
458  FOREACH( const att_key_t& vt1, this->attributes() ) {
459  FOREACH( const att_val_t& vt2, vt1.second ) {
460 
461  string st;
462 
463  bool status = vt2.second->to_string(st);
464 
465  if( !status ) {
466  WARNING( "GenEvent::write_data: problem serializing attribute: "<<vt1.first )
467  }
468  else {
469  data.attribute_id.push_back(vt2.first);
470  data.attribute_name.push_back(vt1.first);
471  data.attribute_string.push_back(st);
472  }
473  }
474  }
475 }
476 
477 
478 void GenEvent::read_data(const GenEventData &data) {
479  this->clear();
480  this->set_event_number( data.event_number );
481  this->set_units( data.momentum_unit, data.length_unit );
482  this->shift_position_to( data.event_pos );
483 
484  // Fill weights
485  this->weights() = data.weights;
486 
487  // Fill particle information
488  FOREACH( const GenParticleData &pd, data.particles ) {
489  GenParticlePtr p = make_shared<GenParticle>(pd);
490 
491  m_particles.push_back(p);
492 
493  p->m_event = this;
494  p->m_id = particles().size();
495  }
496 
497  // Fill vertex information
498  FOREACH( const GenVertexData &vd, data.vertices ) {
499  GenVertexPtr v = make_shared<GenVertex>(vd);
500 
501  m_vertices.push_back(v);
502 
503  v->m_event = this;
504  v->m_id = -(int)vertices().size();
505  }
506 
507  // Restore links
508  for( unsigned int i=0; i<data.links1.size(); ++i) {
509  int id1 = data.links1[i];
510  int id2 = data.links2[i];
511 
512  if( id1 > 0 ) m_vertices[ (-id2)-1 ]->add_particle_in ( m_particles[ id1-1 ] );
513  else m_vertices[ (-id1)-1 ]->add_particle_out( m_particles[ id2-1 ] );
514  }
515 
516  // Read attributes
517  for( unsigned int i=0; i<data.attribute_id.size(); ++i) {
518  add_attribute( data.attribute_name[i],
519  make_shared<StringAttribute>(data.attribute_string[i]),
520  data.attribute_id[i] );
521  }
522 }
523 
524 
525 #ifndef HEPMC_NO_DEPRECATED
526 
527 bool GenEvent::valid_beam_particles() const {
528  /// @todo Change this definition to require status = 4... and in principle there don't have to be two of them
529  return (m_rootvertex->particles_out().size()==2);
530 }
531 
532 pair<GenParticlePtr,GenParticlePtr> GenEvent::beam_particles() const {
533  /// @todo Change this definition to require status = 4... and in principle there don't have to be two of them
534  switch( m_rootvertex->particles_out().size() ) {
535  case 0: return make_pair(GenParticlePtr(), GenParticlePtr());
536  case 1: return make_pair(m_rootvertex->particles_out()[0], GenParticlePtr());
537  default: return make_pair(m_rootvertex->particles_out()[0], m_rootvertex->particles_out()[1]);
538  }
539 }
540 
541 void GenEvent::set_beam_particles(const GenParticlePtr& p1, const GenParticlePtr& p2) {
542  /// @todo Require/set status = 4
543  m_rootvertex->add_particle_out(p1);
544  m_rootvertex->add_particle_out(p2);
545 }
546 
547 void GenEvent::set_beam_particles(const pair<GenParticlePtr,GenParticlePtr>& p) {
548  /// @todo Require/set status = 4
549  m_rootvertex->add_particle_out(p.first);
550  m_rootvertex->add_particle_out(p.second);
551 }
552 
553 #endif
554 
555 string GenEvent::attribute_as_string(const string &name, int id) const {
556 
557  std::map< string, std::map<int, shared_ptr<Attribute> > >::iterator i1 = m_attributes.find(name);
558  if( i1 == m_attributes.end() ) {
559  if ( id == 0 && run_info() ) {
560  return run_info()->attribute_as_string(name);
561  }
562  return string();
563  }
564 
565  std::map<int, shared_ptr<Attribute> >::iterator i2 = i1->second.find(id);
566  if (i2 == i1->second.end() ) return string();
567 
568  if( !i2->second ) return string();
569 
570  string ret;
571  i2->second->to_string(ret);
572 
573  return ret;
574 }
575 
576 } // namespace HepMC
const std::vector< GenParticlePtr > & particles() const
Get list of particles (const)
Stores serializable particle information.
std::vector< std::string > attribute_string
Attribute serialized as string.
void remove_particles(std::vector< GenParticlePtr > v)
Remove a set of particles.
Definition: GenEvent.cc:163
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:57
std::vector< int > links1
First id of the vertex links.
std::map< string, std::map< int, shared_ptr< Attribute > > > m_attributes
Map of event, particle and vertex attributes.
STL namespace.
void remove_particle(GenParticlePtr v)
Remove particle from the event.
Definition: GenEvent.cc:78
std::vector< GenVertexData > vertices
Vertices.
void remove_vertex(GenVertexPtr v)
Remove vertex from the event.
Definition: GenEvent.cc:174
bool is_zero() const
Check if the length of this vertex is zero.
Stores vertex-related information.
std::map< int, shared_ptr< Attribute > >::value_type att_val_t
Attribute map value type.
std::vector< GenParticlePtr > m_particles
List of particles.
Stores serializable vertex information.
GenVertexPtr m_rootvertex
The root vertex is stored outside the normal vertices list to block user access to it...
std::vector< int > attribute_id
Attribute owner id.
std::vector< int > links2
Second id of the vertex links.
void add_particle(GenParticlePtr p)
Add particle.
Definition: GenEvent.cc:42
std::vector< GenVertexPtr > m_vertices
List of vertices.
std::map< string, std::map< int, shared_ptr< Attribute > > >::value_type att_key_t
Attribute map key type.
GenEvent(Units::MomentumUnit momentum_unit=Units::GEV, Units::LengthUnit length_unit=Units::MM)
Event constructor without a run.
Definition: GenEvent.cc:25
const std::vector< GenVertexPtr > & vertices() const
Get list of vertices (const)
std::vector< std::string > attribute_name
Attribute name.
Definition of template class SmartPointer.
std::vector< GenParticleData > particles
Particles.
Units::MomentumUnit momentum_unit
Momentum unit.
Stores serializable event information.
Stores particle-related information.