HepMC3 event record library
GenEvent.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details)
5 //
6 /**
7  * @file GenEvent.cc
8  * @brief Implementation of \b class GenEvent
9  *
10  */
11 #include "HepMC3/GenEvent.h"
12 #include "HepMC3/GenParticle.h"
13 #include "HepMC3/GenVertex.h"
15 
16 #include <deque>
17 #include <algorithm> // sort
18 using namespace std;
19 
20 namespace HepMC3 {
21 
22 GenEvent::GenEvent(Units::MomentumUnit mu,
24  : m_event_number(0), m_weights(std::vector<double>()), //m_weights(std::vector<double>(1, 1.0)),//Prevent from different number of weights and names
25  m_momentum_unit(mu), m_length_unit(lu),
26  m_rootvertex(make_shared<GenVertex>()) {}
27 
28 
29 GenEvent::GenEvent(shared_ptr<GenRunInfo> run,
32  : m_event_number(0), m_weights(std::vector<double>()), //m_weights(std::vector<double>(1, 1.0)),//Prevent from different number of weights and names
33  m_momentum_unit(mu), m_length_unit(lu),
34  m_rootvertex(make_shared<GenVertex>()),
35  m_run_info(run) {
36  if ( run && !run->weight_names().empty() )
37  m_weights = std::vector<double>(run->weight_names().size(), 1.0);
38 }
39 
40 const std::vector<ConstGenParticlePtr>& GenEvent::particles() const {
41  return *(reinterpret_cast<const std::vector<ConstGenParticlePtr>*>(&m_particles));
42 }
43 
44 const std::vector<ConstGenVertexPtr>& GenEvent::vertices() const {
45  return *(reinterpret_cast<const std::vector<ConstGenVertexPtr>*>(&m_vertices));
46 }
47 
48 
49 // void GenEvent::add_particle( const GenParticlePtr &p ) {
50 void GenEvent::add_particle( GenParticlePtr p ) {
51  if( !p|| p->in_event() ) return;
52 
53  m_particles.push_back(p);
54 
55  p->m_event = this;
56  p->m_id = particles().size();
57 
58  // Particles without production vertex are added to the root vertex
59  if( !p->production_vertex() )
60  m_rootvertex->add_particle_out(p);
61 }
62 
63 
65  if (this != &e)
66  {
68  std::lock_guard<std::recursive_mutex> lhs_lk(m_lock_attributes, std::adopt_lock);
69  std::lock_guard<std::recursive_mutex> rhs_lk(e.m_lock_attributes, std::adopt_lock);
70  GenEventData tdata;
71  e.write_data(tdata);
72  read_data(tdata);
73  }
74 }
75 
77  for ( std::map< string, std::map<int, shared_ptr<Attribute> > >::iterator attm=m_attributes.begin(); attm!=m_attributes.end(); ++attm)
78  for ( std::map<int, shared_ptr<Attribute> >::iterator att=attm->second.begin(); att!=attm->second.end(); ++att) att->second->m_event = nullptr;
79 
80  for ( std::vector<GenVertexPtr>::iterator v=m_vertices.begin(); v!=m_vertices.end(); ++v ) if ((*v)->m_event==this) (*v)->m_event=nullptr;
81  for ( std::vector<GenParticlePtr>::iterator p=m_particles.begin(); p!=m_particles.end(); ++p ) if ((*p)->m_event==this) (*p)->m_event=nullptr;
82 }
83 
85  if (this != &e)
86  {
88  std::lock_guard<std::recursive_mutex> lhs_lk(m_lock_attributes, std::adopt_lock);
89  std::lock_guard<std::recursive_mutex> rhs_lk(e.m_lock_attributes, std::adopt_lock);
90  GenEventData tdata;
91  e.write_data(tdata);
92  read_data(tdata);
93  }
94  return *this;
95 }
96 
97 
98 void GenEvent::add_vertex( GenVertexPtr v ) {
99  if( !v|| v->in_event() ) return;
100  m_vertices.push_back(v);
101 
102  v->m_event = this;
103  v->m_id = -(int)vertices().size();
104 
105  // Add all incoming and outgoing particles and restore their production/end vertices
106  for(auto p: v->particles_in() ) {
107  if(!p->in_event()) add_particle(p);
108  p->m_end_vertex = v->shared_from_this();
109  }
110 
111  for(auto p: v->particles_out() ) {
112  if(!p->in_event()) add_particle(p);
113  p->m_production_vertex = v;
114  }
115 }
116 
117 
118 void GenEvent::remove_particle( GenParticlePtr p ) {
119  if( !p || p->parent_event() != this ) return;
120 
121  HEPMC3_DEBUG( 30, "GenEvent::remove_particle - called with particle: "<<p->id() );
122  GenVertexPtr end_vtx = p->end_vertex();
123  if( end_vtx ) {
124  end_vtx->remove_particle_in(p);
125 
126  // If that was the only incoming particle, remove vertex from the event
127  if( end_vtx->particles_in().size() == 0 ) remove_vertex(end_vtx);
128  }
129 
130  GenVertexPtr prod_vtx = p->production_vertex();
131  if( prod_vtx ) {
132  prod_vtx->remove_particle_out(p);
133 
134  // If that was the only outgoing particle, remove vertex from the event
135  if( prod_vtx->particles_out().size() == 0 ) remove_vertex(prod_vtx);
136  }
137 
138  HEPMC3_DEBUG( 30, "GenEvent::remove_particle - erasing particle: " << p->id() )
139 
140  int idx = p->id();
141  vector<GenParticlePtr>::iterator it = m_particles.erase(m_particles.begin() + idx-1 );
142 
143  // Remove attributes of this particle
144  std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
145  vector<string> atts = p->attribute_names();
146  for(const string &s: atts) {
147  p->remove_attribute(s);
148  }
149 
150  //
151  // Reassign id of attributes with id above this one
152  //
153  vector< pair< int, shared_ptr<Attribute> > > changed_attributes;
154 
155  for(att_key_t& vt1: m_attributes ) {
156  changed_attributes.clear();
157 
158  for ( std::map<int, shared_ptr<Attribute> >::iterator vt2=vt1.second.begin(); vt2!=vt1.second.end(); ++vt2) {
159  if( (*vt2).first > p->id() ) {
160  changed_attributes.push_back(*vt2);
161  }
162  }
163 
164  for( att_val_t val: changed_attributes ) {
165  vt1.second.erase(val.first);
166  vt1.second[val.first-1] = val.second;
167  }
168  }
169  // Reassign id of particles with id above this one
170  for(; it != m_particles.end(); ++it) {
171  --((*it)->m_id);
172  }
173 
174  // Finally - set parent event and id of this particle to 0
175  p->m_event = nullptr;
176  p->m_id = 0;
177 }
178 /** @brief Comparison of two particle by id */
180  /** @brief Comparison of two particle by id */
181  inline bool operator()(const GenParticlePtr& p1, const GenParticlePtr& p2) {
182  return (p1->id() > p2->id());
183  }
184 };
185 
186 void GenEvent::remove_particles( vector<GenParticlePtr> v ) {
187 
188  sort( v.begin(), v.end(), sort_by_id_asc() );
189 
190  for (std::vector<GenParticlePtr>::iterator p=v.begin(); p!=v.end(); ++p) {
191  remove_particle(*p);
192  }
193 }
194 
195 void GenEvent::remove_vertex( GenVertexPtr v ) {
196  if( !v || v->parent_event() != this ) return;
197 
198  HEPMC3_DEBUG( 30, "GenEvent::remove_vertex - called with vertex: "<<v->id() );
199  shared_ptr<GenVertex> null_vtx;
200 
201  for(auto p: v->particles_in() ) {
202  p->m_end_vertex = std::weak_ptr<GenVertex>();
203  }
204 
205  for(auto p: v->particles_out() ) {
206  p->m_production_vertex = std::weak_ptr<GenVertex>();
207 
208  // recursive delete rest of the tree
209  remove_particle(p);
210  }
211 
212  // Erase this vertex from vertices list
213  HEPMC3_DEBUG( 30, "GenEvent::remove_vertex - erasing vertex: " << v->id() )
214 
215  int idx = -v->id();
216  vector<GenVertexPtr>::iterator it = m_vertices.erase(m_vertices.begin() + idx-1 );
217  // Remove attributes of this vertex
218  std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
219  vector<string> atts = v->attribute_names();
220  for(string s: atts) {
221  v->remove_attribute(s);
222  }
223 
224  //
225  // Reassign id of attributes with id below this one
226  //
227 
228  vector< pair< int, shared_ptr<Attribute> > > changed_attributes;
229 
230  for( att_key_t& vt1: m_attributes ) {
231  changed_attributes.clear();
232 
233  for ( std::map<int, shared_ptr<Attribute> >::iterator vt2=vt1.second.begin(); vt2!=vt1.second.end(); ++vt2) {
234  if( (*vt2).first < v->id() ) {
235  changed_attributes.push_back(*vt2);
236  }
237  }
238 
239  for( att_val_t val: changed_attributes ) {
240  vt1.second.erase(val.first);
241  vt1.second[val.first+1] = val.second;
242  }
243  }
244 
245  // Reassign id of particles with id above this one
246  for(; it != m_vertices.end(); ++it) {
247  ++((*it)->m_id);
248  }
249 
250  // Finally - set parent event and id of this vertex to 0
251  v->m_event = nullptr;
252  v->m_id = 0;
253 }
254 /// @todo This looks dangerously similar to the recusive event traversel that we forbade in the
255 /// Core library due to wories about generator dependence
256 static bool visit_children(std::map<ConstGenVertexPtr,int> &a, ConstGenVertexPtr v)
257 {
258  for ( ConstGenParticlePtr p: v->particles_out())
259  if (p->end_vertex())
260  {
261  if (a[p->end_vertex()]!=0) return true;
262  else a[p->end_vertex()]++;
263  if (visit_children(a, p->end_vertex())) return true;
264  }
265  return false;
266 }
267 
268 void GenEvent::add_tree( const vector<GenParticlePtr> &parts ) {
269 
270  shared_ptr<IntAttribute> existing_hc=attribute<IntAttribute>("cycles");
271  bool has_cycles=false;
272  std::map<GenVertexPtr,int> sortingv;
273  std::vector<GenVertexPtr> noinv;
274  if (existing_hc) if (existing_hc->value()!=0) has_cycles=true;
275  if(!existing_hc)
276  {
277  for(GenParticlePtr p: parts ) {
278  GenVertexPtr v = p->production_vertex();
279  if(v) sortingv[v]=0;
280  if( !v || v->particles_in().size()==0 ) {
281  GenVertexPtr v2 = p->end_vertex();
282  if(v2) {noinv.push_back(v2); sortingv[v2]=0;}
283  }
284  }
285  for (GenVertexPtr v: noinv) {
286  std::map<ConstGenVertexPtr,int> sorting_temp(sortingv.begin(), sortingv.end());
287  has_cycles=(has_cycles||visit_children(sorting_temp, v));
288  }
289  }
290  if (has_cycles) {
291  add_attribute("cycles", std::make_shared<IntAttribute>(1));
292  /* Commented out as improvemnts allow us to do sorting in other way.
293  for( std::map<GenVertexPtr,int>::iterator vi=sortingv.begin();vi!=sortingv.end();++vi) if( !vi->first->in_event() ) add_vertex(vi->first);
294  return;
295  */
296  }
297 
298  deque<GenVertexPtr> sorting;
299 
300  // Find all starting vertices (end vertex of particles that have no production vertex)
301  for(auto p: parts ) {
302  const GenVertexPtr &v = p->production_vertex();
303  if( !v || v->particles_in().size()==0 ) {
304  const GenVertexPtr &v2 = p->end_vertex();
305  if(v2) sorting.push_back(v2);
306  }
307  }
308 
310  unsigned int sorting_loop_count = 0;
311  unsigned int max_deque_size = 0;
312  )
313 
314  // Add vertices to the event in topological order
315  while( !sorting.empty() ) {
317  if( sorting.size() > max_deque_size ) max_deque_size = sorting.size();
318  ++sorting_loop_count;
319  )
320 
321  GenVertexPtr &v = sorting.front();
322 
323  bool added = false;
324 
325  // Add all mothers to the front of the list
326  for( auto p: v->particles_in() ) {
327  GenVertexPtr v2 = p->production_vertex();
328  if( v2 && !v2->in_event() && find(sorting.begin(),sorting.end(),v2)==sorting.end()) {
329  sorting.push_front(v2);
330  added = true;
331  }
332  }
333 
334  // If we have added at least one production vertex,
335  // our vertex is not the first one on the list
336  if( added ) continue;
337 
338  // If vertex not yet added
339  if( !v->in_event() ) {
340 
341  add_vertex(v);
342 
343  // Add all end vertices to the end of the list
344  for(auto p: v->particles_out() ) {
345  GenVertexPtr v2 = p->end_vertex();
346  if( v2 && !v2->in_event()&& find(sorting.begin(),sorting.end(),v2)==sorting.end() ) {
347  sorting.push_back(v2);
348  }
349  }
350  }
351 
352  sorting.pop_front();
353  }
354 
355  // LL: Make sure root vertex has index zero and is not written out
356  if ( m_rootvertex->id() != 0 ) {
357  const int vx = -1 - m_rootvertex->id();
358  const int rootid = m_rootvertex->id();
359  if ( vx >= 0 && vx < m_vertices.size() && m_vertices[vx] == m_rootvertex ) {
360  auto next = m_vertices.erase(m_vertices.begin() + vx);
361  std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
362  for(auto & vt1: m_attributes ) {
363  vector< pair< int, shared_ptr<Attribute> > > changed_attributes;
364  for ( auto vt2 : vt1.second )
365  if( vt2.first <= rootid )
366  changed_attributes.push_back(vt2);
367  for( auto val : changed_attributes ) {
368  vt1.second.erase(val.first);
369  vt1.second[val.first == rootid? 0: val.first + 1] = val.second;
370  }
371  }
372  m_rootvertex->set_id(0);
373  while ( next != m_vertices.end() ) {
374  ++((*next++)->m_id);
375  }
376  } else {
377  HEPMC3_WARNING( "ReaderAsciiHepMC2: Suspicious looking rootvertex found. Will try to cope." )
378  }
379  }
380 
382  HEPMC3_DEBUG( 6, "GenEvent - particles sorted: "
383  <<this->particles().size()<<", max deque size: "
384  <<max_deque_size<<", iterations: "<<sorting_loop_count )
385  )
386  return;
387 }
388 
389 
390 void GenEvent::reserve(const size_t& parts, const size_t& verts) {
391  m_particles.reserve(parts);
392  m_vertices.reserve(verts);
393 }
394 
395 
396 void GenEvent::set_units( Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit) {
397  if( new_momentum_unit != m_momentum_unit ) {
398  for( GenParticlePtr p: m_particles ) {
399  Units::convert( p->m_data.momentum, m_momentum_unit, new_momentum_unit );
400  Units::convert( p->m_data.mass, m_momentum_unit, new_momentum_unit );
401  }
402 
403  m_momentum_unit = new_momentum_unit;
404  }
405 
406  if( new_length_unit != m_length_unit ) {
407  for(GenVertexPtr &v: m_vertices ) {
408  FourVector &fv = v->m_data.position;
409  if( !fv.is_zero() ) Units::convert( fv, m_length_unit, new_length_unit );
410  }
411 
412  m_length_unit = new_length_unit;
413  }
414 }
415 
416 
418  return m_rootvertex->data().position;
419 }
420 
421 vector<ConstGenParticlePtr> GenEvent::beams() const {
422  return std::const_pointer_cast<const GenVertex>(m_rootvertex)->particles_out();
423 }
424 
425 const vector<GenParticlePtr> & GenEvent::beams() {
426  return m_rootvertex->particles_out();
427 }
428 
430  m_rootvertex->set_position( event_pos() + delta );
431 
432  // Offset all vertices
433  for ( GenVertexPtr v: m_vertices ) {
434  if ( v->has_set_position() )
435  v->set_position( v->position() + delta );
436  }
437 }
438 
439 bool GenEvent::rotate( const FourVector& delta )
440 {
441 
442  for ( auto p: m_particles)
443  {
444  FourVector mom=p->momentum();
445  long double tempX=mom.x();
446  long double tempY=mom.y();
447  long double tempZ=mom.z();
448 
449  long double tempX_;
450  long double tempY_;
451  long double tempZ_;
452 
453 
454  long double cosa=cos(delta.x());
455  long double sina=sin(delta.x());
456 
457  tempY_= cosa*tempY+sina*tempZ;
458  tempZ_=-sina*tempY+cosa*tempZ;
459  tempY=tempY_;
460  tempZ=tempZ_;
461 
462 
463  long double cosb=cos(delta.y());
464  long double sinb=sin(delta.y());
465 
466  tempX_= cosb*tempX-sinb*tempZ;
467  tempZ_= sinb*tempX+cosb*tempZ;
468  tempX=tempX_;
469  tempZ=tempZ_;
470 
471  long double cosg=cos(delta.z());
472  long double sing=sin(delta.z());
473 
474  tempX_= cosg*tempX+sing*tempY;
475  tempY_=-sing*tempX+cosg*tempY;
476  tempX=tempX_;
477  tempY=tempY_;
478 
479  FourVector temp(tempX,tempY,tempZ,mom.e());
480  p->set_momentum(temp);
481  }
482  for ( auto v: m_vertices)
483  {
484  FourVector pos=v->position();
485  long double tempX=pos.x();
486  long double tempY=pos.y();
487  long double tempZ=pos.z();
488 
489  long double tempX_;
490  long double tempY_;
491  long double tempZ_;
492 
493 
494  long double cosa=cos(delta.x());
495  long double sina=sin(delta.x());
496 
497  tempY_= cosa*tempY+sina*tempZ;
498  tempZ_=-sina*tempY+cosa*tempZ;
499  tempY=tempY_;
500  tempZ=tempZ_;
501 
502 
503  long double cosb=cos(delta.y());
504  long double sinb=sin(delta.y());
505 
506  tempX_= cosb*tempX-sinb*tempZ;
507  tempZ_= sinb*tempX+cosb*tempZ;
508  tempX=tempX_;
509  tempZ=tempZ_;
510 
511  long double cosg=cos(delta.z());
512  long double sing=sin(delta.z());
513 
514  tempX_= cosg*tempX+sing*tempY;
515  tempY_=-sing*tempX+cosg*tempY;
516  tempX=tempX_;
517  tempY=tempY_;
518 
519  FourVector temp(tempX,tempY,tempZ,pos.t());
520  v->set_position(temp);
521  }
522 
523 
524  return true;
525 }
526 
527 bool GenEvent::reflect(const int axis)
528 {
529  if (axis>3||axis<0)
530  {
531  HEPMC3_WARNING( "GenEvent::reflect: wrong axis" )
532  return false;
533  }
534  switch (axis)
535  {
536  case 0:
537  for ( auto p: m_particles) { FourVector temp=p->momentum(); temp.setX(-p->momentum().x()); p->set_momentum(temp);}
538  for ( auto v: m_vertices) { FourVector temp=v->position(); temp.setX(-v->position().x()); v->set_position(temp);}
539  break;
540  case 1:
541  for ( auto p: m_particles) { FourVector temp=p->momentum(); temp.setY(-p->momentum().y()); p->set_momentum(temp);}
542  for ( auto v: m_vertices) { FourVector temp=v->position(); temp.setY(-v->position().y()); v->set_position(temp);}
543  break;
544  case 2:
545  for ( auto p: m_particles) { FourVector temp=p->momentum(); temp.setZ(-p->momentum().z()); p->set_momentum(temp);}
546  for ( auto v: m_vertices) { FourVector temp=v->position(); temp.setZ(-v->position().z()); v->set_position(temp);}
547  break;
548  case 3:
549  for ( auto p: m_particles) { FourVector temp=p->momentum(); temp.setT(-p->momentum().e()); p->set_momentum(temp);}
550  for ( auto v: m_vertices) { FourVector temp=v->position(); temp.setT(-v->position().t()); v->set_position(temp);}
551  break;
552  default:
553  return false;
554  }
555 
556  return true;
557 }
558 
559 bool GenEvent::boost( const FourVector& delta )
560 {
561 
562  double deltalength2d=delta.length2();
563  if (deltalength2d>1.0)
564  {
565  HEPMC3_WARNING( "GenEvent::boost: wrong large boost vector. Will leave event as is." )
566  return false;
567  }
568  if (std::abs(deltalength2d-1.0)<std::numeric_limits<double>::epsilon())
569  {
570  HEPMC3_WARNING( "GenEvent::boost: too large gamma. Will leave event as is." )
571  return false;
572  }
573  if (std::abs(deltalength2d)<std::numeric_limits<double>::epsilon())
574  {
575  HEPMC3_WARNING( "GenEvent::boost: wrong small boost vector. Will leave event as is." )
576  return true;
577  }
578  long double deltaX=delta.x();
579  long double deltaY=delta.y();
580  long double deltaZ=delta.z();
581  long double deltalength2=deltaX*deltaX+deltaY*deltaY+deltaZ*deltaZ;
582  long double deltalength=std::sqrt(deltalength2 );
583  long double gamma=1.0/std::sqrt(1.0-deltalength2);
584 
585  for ( auto p: m_particles)
586  {
587  FourVector mom=p->momentum();
588 
589  long double tempX=mom.x();
590  long double tempY=mom.y();
591  long double tempZ=mom.z();
592  long double tempE=mom.e();
593  long double nr=(deltaX*tempX+deltaY*tempY+deltaZ*tempZ)/deltalength;
594 
595  tempX+=(deltaX*((gamma-1)*nr/deltalength)-deltaX*(tempE*gamma));
596  tempY+=(deltaY*((gamma-1)*nr/deltalength)-deltaY*(tempE*gamma));
597  tempZ+=(deltaZ*((gamma-1)*nr/deltalength)-deltaZ*(tempE*gamma));
598  tempE=gamma*(tempE-deltalength*nr);
599  FourVector temp(tempX,tempY,tempZ,tempE);
600  p->set_momentum(temp);
601  }
602 
603  return true;
604 }
605 
606 
607 
608 
610  std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
611  m_event_number = 0;
612  m_rootvertex = make_shared<GenVertex>();
613  m_weights.clear();
614  m_attributes.clear();
615  m_particles.clear();
616  m_vertices.clear();
617 }
618 
619 
620 
621 
622 void GenEvent::remove_attribute(const string &name, const int& id) {
623  std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
624  map< string, map<int, shared_ptr<Attribute> > >::iterator i1 =
625  m_attributes.find(name);
626  if( i1 == m_attributes.end() ) return;
627 
628  map<int, shared_ptr<Attribute> >::iterator i2 = i1->second.find(id);
629  if( i2 == i1->second.end() ) return;
630 
631  i1->second.erase(i2);
632 }
633 
634 vector<string> GenEvent::attribute_names( const int& id) const {
635  vector<string> results;
636 
637  for(const att_key_t& vt1: m_attributes ) {
638  for( const att_val_t& vt2: vt1.second ) {
639  if( vt2.first == id ) results.push_back( vt1.first );
640  }
641  }
642 
643  return results;
644 }
645 
647  // Reserve memory for containers
648  data.particles.reserve( this->particles().size() );
649  data.vertices.reserve( this->vertices().size() );
650  data.links1.reserve( this->particles().size()*2 );
651  data.links2.reserve( this->particles().size()*2 );
652  data.attribute_id.reserve( m_attributes.size() );
653  data.attribute_name.reserve( m_attributes.size() );
654  data.attribute_string.reserve( m_attributes.size() );
655 
656  // Fill event data
657  data.event_number = this->event_number();
658  data.momentum_unit = this->momentum_unit();
659  data.length_unit = this->length_unit();
660  data.event_pos = this->event_pos();
661 
662  // Fill containers
663  data.weights = this->weights();
664 
665  for( ConstGenParticlePtr p: this->particles() ) {
666  data.particles.push_back( p->data() );
667  }
668 
669  for(ConstGenVertexPtr v: this->vertices() ) {
670  data.vertices.push_back( v->data() );
671  int v_id = v->id();
672 
673  for(ConstGenParticlePtr p: v->particles_in() ) {
674  data.links1.push_back( p->id() );
675  data.links2.push_back( v_id );
676  }
677 
678  for(ConstGenParticlePtr p: v->particles_out() ) {
679  data.links1.push_back( v_id );
680  data.links2.push_back( p->id() );
681  }
682  }
683 
684  for( const att_key_t& vt1: this->attributes() ) {
685  for( const att_val_t& vt2: vt1.second ) {
686 
687  string st;
688 
689  bool status = vt2.second->to_string(st);
690 
691  if( !status ) {
692  HEPMC3_WARNING( "GenEvent::write_data: problem serializing attribute: "<<vt1.first )
693  }
694  else {
695  data.attribute_id.push_back(vt2.first);
696  data.attribute_name.push_back(vt1.first);
697  data.attribute_string.push_back(st);
698  }
699  }
700  }
701 }
702 
703 
705  this->clear();
706  this->set_event_number( data.event_number );
707 //Note: set_units checks the current unit of event, i.e. applicable only for fully constructed event.
710  this->shift_position_to( data.event_pos );
711 
712  // Fill weights
713  this->weights() = data.weights;
714 
715  // Fill particle information
716  for( const GenParticleData &pd: data.particles ) {
717  GenParticlePtr p = make_shared<GenParticle>(pd);
718 
719  m_particles.push_back(p);
720 
721  p->m_event = this;
722  p->m_id = m_particles.size();
723  }
724 
725  // Fill vertex information
726  for( const GenVertexData &vd: data.vertices ) {
727  GenVertexPtr v = make_shared<GenVertex>(vd);
728 
729  m_vertices.push_back(v);
730 
731  v->m_event = this;
732  v->m_id = -(int)m_vertices.size();
733  }
734 
735  // Restore links
736  for( unsigned int i=0; i<data.links1.size(); ++i) {
737  int id1 = data.links1[i];
738  int id2 = data.links2[i];
739  /* @note:
740  The meaningfull combinations for (id1,id2) are:
741  (+-) -- particle has end vertex
742  (-+) -- particle has production vertex
743  */
744  if ( (id1<0&&id2<0)|| (id1>0&&id2>0) ) { HEPMC3_WARNING( "GenEvent::read_data: wrong link: "<<id1<<" "<<id2 ); continue;}
745 
746  if ( id1 > 0 ) { m_vertices[ (-id2)-1 ]->add_particle_in ( m_particles[ id1-1 ] ); continue; }
747  if ( id1 < 0 ) { m_vertices[ (-id1)-1 ]->add_particle_out( m_particles[ id2-1 ] ); continue; }
748  }
749  for (auto p: m_particles) if (!p->production_vertex()) m_rootvertex->add_particle_out(p);
750 
751  // Read attributes
752  for( unsigned int i=0; i<data.attribute_id.size(); ++i) {
753  add_attribute( data.attribute_name[i],
754  make_shared<StringAttribute>(data.attribute_string[i]),
755  data.attribute_id[i] );
756  }
757 }
758 
759 
760 
761 //
762 // Deprecated functions
763 //
764 
766  add_particle( GenParticlePtr(p) );
767 }
768 
769 
771  add_vertex( GenVertexPtr(v) );
772 }
773 
774 
775 void GenEvent::set_beam_particles(GenParticlePtr p1, GenParticlePtr p2) {
776  m_rootvertex->add_particle_out(p1);
777  m_rootvertex->add_particle_out(p2);
778 }
779 
780 void GenEvent::add_beam_particle(GenParticlePtr p1) {
781  if (!p1)
782  {
783  HEPMC3_WARNING("Attempting to add an empty particle as beam particle. Ignored.")
784  return;
785  }
786  if( p1->in_event()) if (p1->parent_event()!=this)
787  {
788  HEPMC3_WARNING("Attempting to add particle from another event. Ignored.")
789  return;
790  }
791  if (p1->production_vertex()) p1->production_vertex()->remove_particle_out(p1);
792 //Particle w/o production vertex is added to root vertex.
793  add_particle(p1);
794  p1->set_status(4);
795  return;
796 }
797 
798 
799 string GenEvent::attribute_as_string(const string &name, const int& id) const {
800  std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
801  std::map< string, std::map<int, shared_ptr<Attribute> > >::iterator i1 =
802  m_attributes.find(name);
803  if( i1 == m_attributes.end() ) {
804  if ( id == 0 && run_info() ) {
805  return run_info()->attribute_as_string(name);
806  }
807  return string();
808  }
809 
810  std::map<int, shared_ptr<Attribute> >::iterator i2 = i1->second.find(id);
811  if (i2 == i1->second.end() ) return string();
812 
813  if( !i2->second ) return string();
814 
815  string ret;
816  i2->second->to_string(ret);
817 
818  return ret;
819 }
820 
821 } // namespace HepMC3
Units::MomentumUnit m_momentum_unit
Momentum unit.
Definition: GenEvent.h:353
int event_number() const
Get event number.
Definition: GenEvent.h:135
void remove_particle(GenParticlePtr v)
Remove particle from the event.
Definition: GenEvent.cc:118
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
Definition: GenEvent.h:140
HepMC3 main namespace.
void add_beam_particle(GenParticlePtr p1)
Add particle to root vertex.
Definition: GenEvent.cc:780
#define HEPMC3_DEBUG_CODE_BLOCK(x)
Macro for storing code useful for debugging.
Definition: Errors.h:34
const Units::LengthUnit & length_unit() const
Get length unit.
Definition: GenEvent.h:142
void remove_particles(std::vector< GenParticlePtr > v)
Remove a set of particles.
Definition: GenEvent.cc:186
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
Definition: Errors.h:26
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:98
std::vector< int > attribute_id
Attribute owner id.
Definition: GenEventData.h:54
bool boost(const FourVector &v)
Boost event using x,y,z components of v as velocities.
Definition: GenEvent.cc:559
void write_data(GenEventData &data) const
Fill GenEventData object.
Definition: GenEvent.cc:646
Definition of class GenParticle.
void remove_vertex(GenVertexPtr v)
Remove vertex from the event.
Definition: GenEvent.cc:195
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:44
void shift_position_by(const FourVector &delta)
Shift position of all vertices in the event by delta.
Definition: GenEvent.cc:429
GenEvent(Units::MomentumUnit momentum_unit=Units::GEV, Units::LengthUnit length_unit=Units::MM)
Event constructor without a run.
Definition: GenEvent.cc:22
Stores vertex-related information.
Definition: GenVertex.h:26
std::vector< int > links2
Second id of the vertex links.
Definition: GenEventData.h:52
std::vector< GenVertexPtr > m_vertices
List of vertices.
Definition: GenEvent.h:344
static void convert(T &m, MomentumUnit from, MomentumUnit to)
Convert FourVector to different momentum unit.
Definition: Units.h:81
std::vector< string > attribute_names(const int &id=0) const
Get list of attribute names.
Definition: GenEvent.cc:634
Definition of class GenVertex.
#define HEPMC3_DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
Definition: Errors.h:32
Stores particle-related information.
Definition: GenParticle.h:31
int event_number
Event number.
Definition: GenEventData.h:27
bool reflect(const int axis)
Change sign of axis.
Definition: GenEvent.cc:527
void add_particle(GenParticlePtr p)
Add particle.
Definition: GenEvent.cc:50
shared_ptr< GenRunInfo > run_info() const
Get a pointer to the the GenRunInfo object.
Definition: GenEvent.h:124
Definition of struct GenEventData.
bool operator()(const GenParticlePtr &p1, const GenParticlePtr &p2)
Comparison of two particle by id.
Definition: GenEvent.cc:181
Units::MomentumUnit momentum_unit
Momentum unit.
Definition: GenEventData.h:28
bool is_zero() const
Check if the length of this vertex is zero.
Definition: FourVector.h:192
std::vector< int > links1
First id of the vertex links.
Definition: GenEventData.h:51
double x() const
x-component of position/displacement
Definition: FourVector.h:80
FourVector event_pos
Event position.
Definition: GenEventData.h:35
std::vector< std::string > attribute_string
Attribute serialized as string.
Definition: GenEventData.h:56
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:86
const FourVector & event_pos() const
Vertex representing the overall event position.
Definition: GenEvent.cc:417
LengthUnit
Position units.
Definition: Units.h:32
void setY(double yy)
Definition: FourVector.h:91
MomentumUnit
Momentum units.
Definition: Units.h:29
void setZ(double zz)
Definition: FourVector.h:98
Stores event-related information.
Definition: GenEvent.h:41
std::map< string, std::map< int, shared_ptr< Attribute > > > attributes() const
Get a copy of the list of attributes.
Definition: GenEvent.h:235
Stores serializable event information.
Definition: GenEventData.h:26
Generic 4-vector.
Definition: FourVector.h:35
std::recursive_mutex m_lock_attributes
Mutex lock for the m_attibutes map.
Definition: GenEvent.h:375
std::vector< std::string > attribute_name
Attribute name.
Definition: GenEventData.h:55
Stores serializable vertex information.
Definition: GenVertexData.h:22
void set_beam_particles(GenParticlePtr p1, GenParticlePtr p2)
Set incoming beam particles.
Definition: GenEvent.cc:775
void remove_attribute(const string &name, const int &id=0)
Remove attribute.
Definition: GenEvent.cc:622
void setT(double tt)
Definition: FourVector.h:105
bool rotate(const FourVector &v)
Rotate event using x,y,z components of v as rotation angles.
Definition: GenEvent.cc:439
Stores serializable particle information.
string attribute_as_string(const string &name, const int &id=0) const
Get attribute of any type as string.
Definition: GenEvent.cc:799
double y() const
y-component of position/displacement
Definition: FourVector.h:87
static bool visit_children(std::map< ConstGenVertexPtr, int > &a, ConstGenVertexPtr v)
Definition: GenEvent.cc:256
double t() const
Time component of position/displacement.
Definition: FourVector.h:101
void read_data(const GenEventData &data)
Fill GenEvent based on GenEventData.
Definition: GenEvent.cc:704
std::map< int, shared_ptr< Attribute > >::value_type att_val_t
Attribute map value type.
Definition: GenEvent.h:372
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
Definition: GenEvent.cc:40
std::vector< GenParticleData > particles
Particles.
Definition: GenEventData.h:31
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition: GenEvent.cc:268
Units::LengthUnit m_length_unit
Length unit.
Definition: GenEvent.h:355
std::vector< GenParticlePtr > m_particles
List of particles.
Definition: GenEvent.h:342
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
Definition: GenEvent.cc:390
std::map< string, std::map< int, shared_ptr< Attribute > > > m_attributes
Map of event, particle and vertex attributes.
Definition: GenEvent.h:366
void set_units(Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit)
Change event units Converts event from current units to new ones.
Definition: GenEvent.cc:396
void add_attribute(const string &name, const shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
Definition: GenEvent.h:208
int m_event_number
Event number.
Definition: GenEvent.h:347
std::vector< double > weights
Weights.
Definition: GenEventData.h:33
std::vector< GenVertexData > vertices
Vertices.
Definition: GenEventData.h:32
double e() const
Energy component of momentum.
Definition: FourVector.h:130
GenVertexPtr m_rootvertex
The root vertex is stored outside the normal vertices list to block user access to it.
Definition: GenEvent.h:358
std::vector< ConstGenParticlePtr > beams() const
Vector of beam particles.
Definition: GenEvent.cc:421
Definition of class GenEvent.
GenEvent & operator=(const GenEvent &)
Assignment operator.
Definition: GenEvent.cc:84
void setX(double xx)
Definition: FourVector.h:84
Comparison of two particle by id.
Definition: GenEvent.cc:179
std::map< string, std::map< int, shared_ptr< Attribute > > >::value_type att_key_t
Attribute map key type.
Definition: GenEvent.h:369
~GenEvent()
Destructor.
Definition: GenEvent.cc:76
void clear()
Remove contents of this event.
Definition: GenEvent.cc:609
std::vector< double > m_weights
Event weights.
Definition: GenEvent.h:350
double length2() const
Squared magnitude of (x, y, z) 3-vector.
Definition: FourVector.h:143
Units::LengthUnit length_unit
Length unit.
Definition: GenEventData.h:29
Feature< Feature_type > abs(const Feature< Feature_type > &input)
Obtain the absolute value of a Feature. This works as you'd expect. If foo is a valid Feature,...
Definition: Feature.h:316
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:137
void shift_position_to(const FourVector &newpos)
Shift position of all vertices in the event to op.
Definition: GenEvent.h:187
double z() const
z-component of position/displacement
Definition: FourVector.h:94