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