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