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-2021 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
19namespace 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
28GenEvent::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
39const std::vector<ConstGenParticlePtr>& GenEvent::particles() const {
40 return *(reinterpret_cast<const std::vector<ConstGenParticlePtr>*>(&m_particles));
41}
42
43const std::vector<ConstGenVertexPtr>& GenEvent::vertices() const {
44 return *(reinterpret_cast<const std::vector<ConstGenVertexPtr>*>(&m_vertices));
45}
46
47
48void 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
96void 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
116void 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
184void 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
192void 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
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
253static 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
265void 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
385void GenEvent::reserve(const size_t& parts, const size_t& verts) {
386 m_particles.reserve(parts);
387 m_vertices.reserve(verts);
388}
389
390
391void 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
416std::vector<ConstGenParticlePtr> GenEvent::beams() const {
417 return std::const_pointer_cast<const GenVertex>(m_rootvertex)->particles_out();
418}
419
420const 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
434bool 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
521bool 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
553bool 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
609void 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
621std::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.
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
760void GenEvent::set_beam_particles(GenParticlePtr p1, GenParticlePtr p2) {
761 m_rootvertex->add_particle_out(p1);
762 m_rootvertex->add_particle_out(p2);
763}
764
765void 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
784std::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
806void GenEvent::add_attribute(const std::string &name, const std::shared_ptr<Attribute> &att, const int& id ) {
807 ///Disallow empty strings
808 if (name.length() == 0) return;
809 if (!att) return;
810 std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
811 if (m_attributes.find(name) == m_attributes.end()) m_attributes[name] = std::map<int, std::shared_ptr<Attribute> >();
812 m_attributes[name][id] = att;
813 att->m_event = this;
814 if ( id > 0 && id <= int(particles().size()) )
815 att->m_particle = particles()[id - 1];
816 if ( id < 0 && -id <= int(vertices().size()) )
817 att->m_vertex = vertices()[-id - 1];
818}
819
820
821void GenEvent::add_attributes(const std::vector<std::string> &names, const std::vector<std::shared_ptr<Attribute> > &atts, const std::vector<int>& ids) {
822 size_t N = names.size();
823 if ( N == 0 ) return;
824 if (N != atts.size()) return;
825 if (N != ids.size()) return;
826
827 std::vector<std::string> unames = names;
828 vector<std::string>::iterator ip;
829 ip = std::unique(unames.begin(), unames.end());
830 unames.resize(std::distance(unames.begin(), ip));
831 std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
832 for (auto name: unames)
833 if (m_attributes.count(name) == 0) m_attributes[name] = std::map<int, std::shared_ptr<Attribute> >();
834
835 const int particles_size = int(m_particles.size());
836 const int vertices_size = int(m_vertices.size());
837 for (size_t i = 0; i < N; i++) {
838 ///Disallow empty strings
839 if (names.at(i).length() == 0) continue;
840 if (!atts[i]) continue;
841 m_attributes[names.at(i)][ids.at(i)] = atts[i];
842 atts[i]->m_event = this;
843 if ( ids.at(i) > 0 && ids.at(i) <= particles_size )
844 { atts[i]->m_particle = m_particles[ids.at(i) - 1]; }
845 else {
846 if ( ids.at(i) < 0 && -ids.at(i) <= vertices_size )
847 atts[i]->m_vertex = m_vertices[-ids.at(i) - 1];
848 }
849 }
850}
851
852void GenEvent::add_attributes(const std::string& name, const std::vector<std::shared_ptr<Attribute> > &atts, const std::vector<int>& ids) {
853 if (name.length() == 0) return;
854 size_t N = ids.size();
855 if(!N) return;
856 if ( N != atts.size()) return;
857
858 std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
859 if (m_attributes.find(name) == m_attributes.end()) m_attributes[name] = std::map<int, std::shared_ptr<Attribute> >();
860 auto& tmap = m_attributes[name];
861 const int particles_size = int(m_particles.size());
862 const int vertices_size = int(m_vertices.size());
863 for (size_t i = 0; i < N; i++) {
864 ///Disallow empty strings
865 if (!atts[i]) continue;
866 tmap[ids.at(i)] = atts[i];
867 atts[i]->m_event = this;
868 if ( ids.at(i) > 0 && ids.at(i) <= particles_size )
869 { atts[i]->m_particle = m_particles[ids.at(i) - 1]; }
870 else {
871 if ( ids.at(i) < 0 && -ids.at(i) <= vertices_size )
872 atts[i]->m_vertex = m_vertices[-ids.at(i) - 1];
873 }
874 }
875}
876void GenEvent::add_attributes(const std::string& name, const std::vector<std::pair<int, std::shared_ptr<Attribute> > > &atts) {
877 if (name.length() == 0) return;
878 if (!atts.size()) return;
879 std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
880 if (m_attributes.find(name) == m_attributes.end()) m_attributes[name] = std::map<int, std::shared_ptr<Attribute> >();
881 auto& tmap = m_attributes[name];
882 const int particles_size = int(m_particles.size());
883 const int vertices_size = int(m_vertices.size());
884 for (auto att: atts) {
885 ///Disallow empty strings
886 if (!att.second) continue;
887 tmap.insert(att);
888 att.second->m_event = this;
889 if ( att.first > 0 && att.first <= particles_size )
890 { att.second->m_particle = m_particles[att.first - 1]; }
891 else {
892 if ( att.first < 0 && -att.first <= vertices_size )
893 att.second->m_vertex = m_vertices[-att.first - 1];
894 }
895 }
896}
897
898} // namespace HepMC3
#define HEPMC3_DEBUG_CODE_BLOCK(x)
Macro for storing code useful for debugging.
Definition: Errors.h:35
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
Definition: Errors.h:27
#define HEPMC3_DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
Definition: Errors.h:33
Definition of struct GenEventData.
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
Generic 4-vector.
Definition: FourVector.h:36
double e() const
Energy component of momentum.
Definition: FourVector.h:131
void setT(double tt)
Definition: FourVector.h:106
double t() const
Time component of position/displacement.
Definition: FourVector.h:102
bool is_zero() const
Check if the length of this vertex is zero.
Definition: FourVector.h:193
void setY(double yy)
Definition: FourVector.h:92
void setX(double xx)
Definition: FourVector.h:85
double x() const
x-component of position/displacement
Definition: FourVector.h:81
double length2() const
Squared magnitude of (x, y, z) 3-vector.
Definition: FourVector.h:144
double y() const
y-component of position/displacement
Definition: FourVector.h:88
double z() const
z-component of position/displacement
Definition: FourVector.h:95
void setZ(double zz)
Definition: FourVector.h:99
Stores event-related information.
Definition: GenEvent.h:41
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
Definition: GenEvent.cc:265
std::recursive_mutex m_lock_attributes
Mutex lock for the m_attibutes map.
Definition: GenEvent.h:397
void add_vertex(GenVertexPtr v)
Add vertex.
Definition: GenEvent.cc:96
void shift_position_to(const FourVector &newpos)
Shift position of all vertices in the event to op.
Definition: GenEvent.h:200
int vertices_size() const
Vertices size, HepMC2 compatiility.
Definition: GenEvent.h:89
int event_number() const
Get event number.
Definition: GenEvent.h:148
std::vector< std::string > attribute_names(const int &id=0) const
Get list of attribute names.
Definition: GenEvent.cc:621
void write_data(GenEventData &data) const
Fill GenEventData object.
Definition: GenEvent.cc:633
std::shared_ptr< GenRunInfo > run_info() const
Get a pointer to the the GenRunInfo object.
Definition: GenEvent.h:137
void remove_particle(GenParticlePtr v)
Remove particle from the event.
Definition: GenEvent.cc:116
void remove_particles(std::vector< GenParticlePtr > v)
Remove a set of particles.
Definition: GenEvent.cc:184
void add_particle(GenParticlePtr p)
Add particle.
Definition: GenEvent.cc:48
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::vector< double > m_weights
Event weights.
Definition: GenEvent.h:372
std::map< int, std::shared_ptr< Attribute > >::value_type att_val_t
Attribute map value type.
Definition: GenEvent.h:394
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > >::value_type att_key_t
Attribute map key type.
Definition: GenEvent.h:391
void shift_position_by(const FourVector &delta)
Shift position of all vertices in the event by delta.
Definition: GenEvent.cc:424
void remove_vertex(GenVertexPtr v)
Remove vertex from the event.
Definition: GenEvent.cc:192
void set_event_number(const int &num)
Set event number.
Definition: GenEvent.h:150
std::vector< ConstGenParticlePtr > beams() const
Vector of beam particles.
Definition: GenEvent.cc:416
int m_event_number
Event number.
Definition: GenEvent.h:369
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:43
~GenEvent()
Destructor.
Definition: GenEvent.cc:74
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
Definition: GenEvent.h:153
const Units::LengthUnit & length_unit() const
Get length unit.
Definition: GenEvent.h:155
bool reflect(const int axis)
Change sign of axis.
Definition: GenEvent.cc:521
GenEvent(Units::MomentumUnit momentum_unit=Units::GEV, Units::LengthUnit length_unit=Units::MM)
Event constructor without a run.
Definition: GenEvent.cc:21
const FourVector & event_pos() const
Vertex representing the overall event position.
Definition: GenEvent.cc:412
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > > attributes() const
Get a copy of the list of attributes.
Definition: GenEvent.h:257
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Definition: GenEvent.cc:806
void read_data(const GenEventData &data)
Fill GenEvent based on GenEventData.
Definition: GenEvent.cc:690
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > > m_attributes
Map of event, particle and vertex attributes.
Definition: GenEvent.h:388
Units::MomentumUnit m_momentum_unit
Momentum unit.
Definition: GenEvent.h:375
void remove_attribute(const std::string &name, const int &id=0)
Remove attribute.
Definition: GenEvent.cc:609
std::vector< GenVertexPtr > m_vertices
List of vertices.
Definition: GenEvent.h:366
GenVertexPtr m_rootvertex
The root vertex is stored outside the normal vertices list to block user access to it.
Definition: GenEvent.h:380
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
const std::vector< double > & weights() const
Get event weight values as a vector.
Definition: GenEvent.h:98
void clear()
Remove contents of this event.
Definition: GenEvent.cc:599
int particles_size() const
Particles size, HepMC2 compatiility.
Definition: GenEvent.h:85
Units::LengthUnit m_length_unit
Length unit.
Definition: GenEvent.h:377
bool boost(const FourVector &v)
Boost event using x,y,z components of v as velocities.
Definition: GenEvent.cc:553
std::vector< GenParticlePtr > m_particles
List of particles.
Definition: GenEvent.h:364
void add_beam_particle(GenParticlePtr p1)
Add particle to root vertex.
Definition: GenEvent.cc:765
void add_attributes(const std::vector< std::string > &names, const std::vector< std::shared_ptr< Attribute > > &atts, const std::vector< int > &ids)
Add multiple attributes to event.
Definition: GenEvent.cc:821
void set_beam_particles(GenParticlePtr p1, GenParticlePtr p2)
Set incoming beam particles.
Definition: GenEvent.cc:760
bool rotate(const FourVector &v)
Rotate event using x,y,z components of v as rotation angles.
Definition: GenEvent.cc:434
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
Definition: GenEvent.cc:39
GenEvent & operator=(const GenEvent &)
Assignment operator.
Definition: GenEvent.cc:82
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
Definition: GenEvent.cc:385
Stores particle-related information.
Definition: GenParticle.h:31
Stores vertex-related information.
Definition: GenVertex.h:26
static void convert(T &m, MomentumUnit from, MomentumUnit to)
Convert FourVector to different momentum unit.
Definition: Units.h:81
LengthUnit
Position units.
Definition: Units.h:32
MomentumUnit
Momentum units.
Definition: Units.h:29
HepMC3 main namespace.
static bool visit_children(std::map< ConstGenVertexPtr, int > &a, ConstGenVertexPtr v)
Definition: GenEvent.cc:253
Stores serializable event information.
Definition: GenEventData.h:26
std::vector< GenVertexData > vertices
Vertices.
Definition: GenEventData.h:32
std::vector< int > links2
Second id of the vertex links.
Definition: GenEventData.h:52
int event_number
Event number.
Definition: GenEventData.h:27
std::vector< std::string > attribute_string
Attribute serialized as string.
Definition: GenEventData.h:56
std::vector< GenParticleData > particles
Particles.
Definition: GenEventData.h:31
std::vector< int > links1
First id of the vertex links.
Definition: GenEventData.h:51
std::vector< std::string > attribute_name
Attribute name.
Definition: GenEventData.h:55
Units::LengthUnit length_unit
Length unit.
Definition: GenEventData.h:29
std::vector< int > attribute_id
Attribute owner id.
Definition: GenEventData.h:54
FourVector event_pos
Event position.
Definition: GenEventData.h:35
std::vector< double > weights
Weights.
Definition: GenEventData.h:33
Units::MomentumUnit momentum_unit
Momentum unit.
Definition: GenEventData.h:28
Stores serializable particle information.
Stores serializable vertex information.
Definition: GenVertexData.h:22
Comparison of two particle by id.
Definition: GenEvent.cc:177
bool operator()(const GenParticlePtr &p1, const GenParticlePtr &p2)
Comparison of two particle by id.
Definition: GenEvent.cc:179