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