23 : m_event_number(0), m_weights(std::vector<double>()),
24 m_momentum_unit(mu), m_length_unit(lu),
25 m_rootvertex(std::make_shared<
GenVertex>()) {}
31 : m_event_number(0), m_weights(std::vector<double>()),
32 m_momentum_unit(mu), m_length_unit(lu),
33 m_rootvertex(std::make_shared<
GenVertex>()),
35 if ( run && !run->weight_names().empty() )
36 m_weights = std::vector<double>(run->weight_names().size(), 1.0);
40 return *(
reinterpret_cast<const std::vector<ConstGenParticlePtr>*
>(&
m_particles));
44 return *(
reinterpret_cast<const std::vector<ConstGenVertexPtr>*
>(&
m_vertices));
50 if( !p|| p->in_event() )
return;
58 if( !p->production_vertex() )
68 std::lock_guard<std::recursive_mutex> rhs_lk(e.
m_lock_attributes, std::adopt_lock);
76 for ( std::map< std::string, std::map<
int, std::shared_ptr<Attribute> > >::iterator attm=
m_attributes.begin(); attm!=
m_attributes.end(); ++attm)
77 for ( std::map<
int, std::shared_ptr<Attribute> >::iterator att=attm->second.begin(); att!=attm->second.end(); ++att)
if (att->second) att->second->m_event =
nullptr;
79 for ( std::vector<GenVertexPtr>::iterator v=
m_vertices.begin(); v!=
m_vertices.end(); ++v )
if (*v)
if ((*v)->m_event==
this) (*v)->m_event=
nullptr;
80 for ( std::vector<GenParticlePtr>::iterator p=
m_particles.begin(); p!=
m_particles.end(); ++p )
if (*p)
if ((*p)->m_event==
this) (*p)->m_event=
nullptr;
88 std::lock_guard<std::recursive_mutex> rhs_lk(e.
m_lock_attributes, std::adopt_lock);
98 if( !v|| v->in_event() )
return;
105 for(
auto p: v->particles_in() ) {
107 p->m_end_vertex = v->shared_from_this();
110 for(
auto p: v->particles_out() ) {
112 p->m_production_vertex = v;
118 if( !p || p->parent_event() !=
this )
return;
120 HEPMC3_DEBUG( 30,
"GenEvent::remove_particle - called with particle: "<<p->id() );
121 GenVertexPtr end_vtx = p->end_vertex();
123 end_vtx->remove_particle_in(p);
126 if( end_vtx->particles_in().size() == 0 )
remove_vertex(end_vtx);
129 GenVertexPtr prod_vtx = p->production_vertex();
131 prod_vtx->remove_particle_out(p);
134 if( prod_vtx->particles_out().size() == 0 )
remove_vertex(prod_vtx);
137 HEPMC3_DEBUG( 30,
"GenEvent::remove_particle - erasing particle: " << p->id() )
144 std::vector<std::string> atts = p->attribute_names();
145 for(
const std::string &s: atts) {
146 p->remove_attribute(s);
152 std::vector< std::pair< int, std::shared_ptr<Attribute> > > changed_attributes;
155 changed_attributes.clear();
157 for ( std::map<
int, std::shared_ptr<Attribute> >::iterator vt2=vt1.second.begin(); vt2!=vt1.second.end(); ++vt2) {
158 if( (*vt2).first > p->id() ) {
159 changed_attributes.push_back(*vt2);
163 for(
att_val_t val: changed_attributes ) {
164 vt1.second.erase(val.first);
165 vt1.second[val.first-1] = val.second;
174 p->m_event =
nullptr;
180 inline bool operator()(
const GenParticlePtr& p1,
const GenParticlePtr& p2) {
181 return (p1->id() > p2->id());
189 for (std::vector<GenParticlePtr>::iterator p=v.begin(); p!=v.end(); ++p) {
195 if( !v || v->parent_event() !=
this )
return;
197 HEPMC3_DEBUG( 30,
"GenEvent::remove_vertex - called with vertex: "<<v->id() );
198 std::shared_ptr<GenVertex> null_vtx;
200 for(
auto p: v->particles_in() ) {
201 p->m_end_vertex = std::weak_ptr<GenVertex>();
204 for(
auto p: v->particles_out() ) {
205 p->m_production_vertex = std::weak_ptr<GenVertex>();
212 HEPMC3_DEBUG( 30,
"GenEvent::remove_vertex - erasing vertex: " << v->id() )
218 std::vector<std::string> atts = v->attribute_names();
219 for(std::string s: atts) {
220 v->remove_attribute(s);
227 std::vector< std::pair< int, std::shared_ptr<Attribute> > > changed_attributes;
230 changed_attributes.clear();
232 for ( std::map<
int, std::shared_ptr<Attribute> >::iterator vt2=vt1.second.begin(); vt2!=vt1.second.end(); ++vt2) {
233 if( (*vt2).first < v->id() ) {
234 changed_attributes.push_back(*vt2);
238 for(
att_val_t val: changed_attributes ) {
239 vt1.second.erase(val.first);
240 vt1.second[val.first+1] = val.second;
250 v->m_event =
nullptr;
255static bool visit_children(std::map<ConstGenVertexPtr,int> &a, ConstGenVertexPtr v)
257 for ( ConstGenParticlePtr p: v->particles_out())
260 if (a[p->end_vertex()]!=0)
return true;
261 else a[p->end_vertex()]++;
269 std::shared_ptr<IntAttribute> existing_hc=attribute<IntAttribute>(
"cycles");
270 bool has_cycles=
false;
271 std::map<GenVertexPtr,int> sortingv;
272 std::vector<GenVertexPtr> noinv;
273 if (existing_hc)
if (existing_hc->value()!=0) has_cycles=
true;
276 for(GenParticlePtr p: parts ) {
277 GenVertexPtr v = p->production_vertex();
279 if( !v || v->particles_in().size()==0 ) {
280 GenVertexPtr v2 = p->end_vertex();
281 if(v2) {noinv.push_back(v2); sortingv[v2]=0;}
284 for (GenVertexPtr v: noinv) {
285 std::map<ConstGenVertexPtr,int> sorting_temp(sortingv.begin(), sortingv.end());
297 std::deque<GenVertexPtr> sorting;
300 for(
auto p: parts ) {
301 const GenVertexPtr &v = p->production_vertex();
302 if( !v || v->particles_in().size()==0 ) {
303 const GenVertexPtr &v2 = p->end_vertex();
304 if(v2) sorting.push_back(v2);
309 unsigned int sorting_loop_count = 0;
310 unsigned int max_deque_size = 0;
314 while( !sorting.empty() ) {
316 if( sorting.size() > max_deque_size ) max_deque_size = sorting.size();
317 ++sorting_loop_count;
320 GenVertexPtr &v = sorting.front();
325 for(
auto p: v->particles_in() ) {
326 GenVertexPtr v2 = p->production_vertex();
327 if( v2 && !v2->in_event() && find(sorting.begin(),sorting.end(),v2)==sorting.end()) {
328 sorting.push_front(v2);
335 if( added )
continue;
338 if( !v->in_event() ) {
343 for(
auto p: v->particles_out() ) {
344 GenVertexPtr v2 = p->end_vertex();
345 if( v2 && !v2->in_event()&& find(sorting.begin(),sorting.end(),v2)==sorting.end() ) {
346 sorting.push_back(v2);
362 std::vector< std::pair< int, std::shared_ptr<Attribute> > > changed_attributes;
363 for (
auto vt2 : vt1.second )
364 if( vt2.first <= rootid )
365 changed_attributes.push_back(vt2);
366 for(
auto val : changed_attributes ) {
367 vt1.second.erase(val.first);
368 vt1.second[val.first == rootid? 0: val.first + 1] = val.second;
376 HEPMC3_WARNING(
"ReaderAsciiHepMC2: Suspicious looking rootvertex found. Will try to cope." )
382 <<this->
particles().size()<<
", max deque size: "
383 <<max_deque_size<<
", iterations: "<<sorting_loop_count )
421 return std::const_pointer_cast<const GenVertex>(
m_rootvertex)->particles_out();
433 if ( v->has_set_position() )
434 v->set_position( v->position() + delta );
444 long double tempX=mom.
x();
445 long double tempY=mom.
y();
446 long double tempZ=mom.
z();
453 long double cosa=cos(delta.
x());
454 long double sina=sin(delta.
x());
456 tempY_= cosa*tempY+sina*tempZ;
457 tempZ_=-sina*tempY+cosa*tempZ;
462 long double cosb=cos(delta.
y());
463 long double sinb=sin(delta.
y());
465 tempX_= cosb*tempX-sinb*tempZ;
466 tempZ_= sinb*tempX+cosb*tempZ;
470 long double cosg=cos(delta.
z());
471 long double sing=sin(delta.
z());
473 tempX_= cosg*tempX+sing*tempY;
474 tempY_=-sing*tempX+cosg*tempY;
479 p->set_momentum(temp);
484 long double tempX=pos.
x();
485 long double tempY=pos.
y();
486 long double tempZ=pos.
z();
493 long double cosa=cos(delta.
x());
494 long double sina=sin(delta.
x());
496 tempY_= cosa*tempY+sina*tempZ;
497 tempZ_=-sina*tempY+cosa*tempZ;
502 long double cosb=cos(delta.
y());
503 long double sinb=sin(delta.
y());
505 tempX_= cosb*tempX-sinb*tempZ;
506 tempZ_= sinb*tempX+cosb*tempZ;
510 long double cosg=cos(delta.
z());
511 long double sing=sin(delta.
z());
513 tempX_= cosg*tempX+sing*tempY;
514 tempY_=-sing*tempX+cosg*tempY;
519 v->set_position(temp);
561 double deltalength2d=delta.
length2();
562 if (deltalength2d>1.0)
564 HEPMC3_WARNING(
"GenEvent::boost: wrong large boost vector. Will leave event as is." )
567 if (
std::abs(deltalength2d-1.0)<std::numeric_limits<double>::epsilon())
569 HEPMC3_WARNING(
"GenEvent::boost: too large gamma. Will leave event as is." )
572 if (
std::abs(deltalength2d)<std::numeric_limits<double>::epsilon())
574 HEPMC3_WARNING(
"GenEvent::boost: wrong small boost vector. Will leave event as is." )
577 long double deltaX=delta.
x();
578 long double deltaY=delta.
y();
579 long double deltaZ=delta.
z();
580 long double deltalength2=deltaX*deltaX+deltaY*deltaY+deltaZ*deltaZ;
581 long double deltalength=std::sqrt(deltalength2 );
582 long double gamma=1.0/std::sqrt(1.0-deltalength2);
588 long double tempX=mom.
x();
589 long double tempY=mom.
y();
590 long double tempZ=mom.
z();
591 long double tempE=mom.
e();
592 long double nr=(deltaX*tempX+deltaY*tempY+deltaZ*tempZ)/deltalength;
594 tempX+=(deltaX*((gamma-1)*nr/deltalength)-deltaX*(tempE*gamma));
595 tempY+=(deltaY*((gamma-1)*nr/deltalength)-deltaY*(tempE*gamma));
596 tempZ+=(deltaZ*((gamma-1)*nr/deltalength)-deltaZ*(tempE*gamma));
597 tempE=gamma*(tempE-deltalength*nr);
599 p->set_momentum(temp);
623 std:: map< std::string, std::map<int, std::shared_ptr<Attribute> > >::iterator i1 =
627 std::map<int, std::shared_ptr<Attribute> >::iterator i2 = i1->second.find(
id);
628 if( i2 == i1->second.end() )
return;
630 i1->second.erase(i2);
634 std::vector<std::string> results;
637 if( vt1.second.count(
id ) == 1 ) {
638 results.push_back( vt1.first );
664 for( ConstGenParticlePtr p: this->
particles() ) {
668 for(ConstGenVertexPtr v: this->
vertices() ) {
669 data.
vertices.push_back( v->data() );
672 for(ConstGenParticlePtr p: v->particles_in() ) {
673 data.
links1.push_back( p->id() );
674 data.
links2.push_back( v_id );
677 for(ConstGenParticlePtr p: v->particles_out() ) {
678 data.
links1.push_back( v_id );
679 data.
links2.push_back( p->id() );
684 for(
const att_val_t& vt2: vt1.second ) {
688 bool status = vt2.second->to_string(st);
691 HEPMC3_WARNING(
"GenEvent::write_data: problem serializing attribute: "<<vt1.first )
716 GenParticlePtr p = std::make_shared<GenParticle>(pd);
726 GenVertexPtr v = std::make_shared<GenVertex>(vd);
735 for(
unsigned int i=0; i<data.
links1.size(); ++i) {
743 if ( (id1<0&&id2<0)|| (id1>0&&id2>0) ) {
HEPMC3_WARNING(
"GenEvent::read_data: wrong link: "<<id1<<
" "<<id2 );
continue;}
751 for(
unsigned int i=0; i<data.
attribute_id.size(); ++i) {
782 HEPMC3_WARNING(
"Attempting to add an empty particle as beam particle. Ignored.")
785 if( p1->in_event())
if (p1->parent_event()!=
this)
787 HEPMC3_WARNING(
"Attempting to add particle from another event. Ignored.")
790 if (p1->production_vertex()) p1->production_vertex()->remove_particle_out(p1);
800 std::map< std::string, std::map<int, std::shared_ptr<Attribute> > >::iterator i1 =
804 return run_info()->attribute_as_string(name);
806 return std::string();
809 std::map<int, std::shared_ptr<Attribute> >::iterator i2 = i1->second.find(
id);
810 if (i2 == i1->second.end() )
return std::string();
812 if( !i2->second )
return std::string();
815 i2->second->to_string(ret);
#define HEPMC3_DEBUG_CODE_BLOCK(x)
Macro for storing code useful for debugging.
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
#define HEPMC3_DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
Definition of struct GenEventData.
Definition of class GenEvent.
Definition of class GenParticle.
Definition of class GenVertex.
double e() const
Energy component of momentum.
double t() const
Time component of position/displacement.
bool is_zero() const
Check if the length of this vertex is zero.
double x() const
x-component of position/displacement
double length2() const
Squared magnitude of (x, y, z) 3-vector.
double y() const
y-component of position/displacement
double z() const
z-component of position/displacement
Stores event-related information.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
std::recursive_mutex m_lock_attributes
Mutex lock for the m_attibutes map.
void add_vertex(GenVertexPtr v)
Add vertex.
void shift_position_to(const FourVector &newpos)
Shift position of all vertices in the event to op.
int event_number() const
Get event number.
std::vector< std::string > attribute_names(const int &id=0) const
Get list of attribute names.
void write_data(GenEventData &data) const
Fill GenEventData object.
std::shared_ptr< GenRunInfo > run_info() const
Get a pointer to the the GenRunInfo object.
void remove_particle(GenParticlePtr v)
Remove particle from the event.
void remove_particles(std::vector< GenParticlePtr > v)
Remove a set of particles.
void add_particle(GenParticlePtr p)
Add particle.
void set_units(Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit)
Change event units Converts event from current units to new ones.
std::vector< double > m_weights
Event weights.
std::map< int, std::shared_ptr< Attribute > >::value_type att_val_t
Attribute map value type.
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > >::value_type att_key_t
Attribute map key type.
void shift_position_by(const FourVector &delta)
Shift position of all vertices in the event by delta.
void remove_vertex(GenVertexPtr v)
Remove vertex from the event.
void set_event_number(const int &num)
Set event number.
std::vector< ConstGenParticlePtr > beams() const
Vector of beam particles.
int m_event_number
Event number.
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
const Units::LengthUnit & length_unit() const
Get length unit.
bool reflect(const int axis)
Change sign of axis.
GenEvent(Units::MomentumUnit momentum_unit=Units::GEV, Units::LengthUnit length_unit=Units::MM)
Event constructor without a run.
const FourVector & event_pos() const
Vertex representing the overall event position.
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > > attributes() const
Get a copy of the list of attributes.
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
void read_data(const GenEventData &data)
Fill GenEvent based on GenEventData.
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > > m_attributes
Map of event, particle and vertex attributes.
Units::MomentumUnit m_momentum_unit
Momentum unit.
void remove_attribute(const std::string &name, const int &id=0)
Remove attribute.
std::vector< GenVertexPtr > m_vertices
List of vertices.
GenVertexPtr m_rootvertex
The root vertex is stored outside the normal vertices list to block user access to it.
std::string attribute_as_string(const std::string &name, const int &id=0) const
Get attribute of any type as string.
const std::vector< double > & weights() const
Get event weight values as a vector.
void clear()
Remove contents of this event.
Units::LengthUnit m_length_unit
Length unit.
bool boost(const FourVector &v)
Boost event using x,y,z components of v as velocities.
std::vector< GenParticlePtr > m_particles
List of particles.
void add_beam_particle(GenParticlePtr p1)
Add particle to root vertex.
void set_beam_particles(GenParticlePtr p1, GenParticlePtr p2)
Set incoming beam particles.
bool rotate(const FourVector &v)
Rotate event using x,y,z components of v as rotation angles.
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
GenEvent & operator=(const GenEvent &)
Assignment operator.
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
Stores particle-related information.
Stores vertex-related information.
static void convert(T &m, MomentumUnit from, MomentumUnit to)
Convert FourVector to different momentum unit.
LengthUnit
Position units.
MomentumUnit
Momentum units.
static bool visit_children(std::map< ConstGenVertexPtr, int > &a, ConstGenVertexPtr v)
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,...
Stores serializable event information.
std::vector< GenVertexData > vertices
Vertices.
std::vector< int > links2
Second id of the vertex links.
int event_number
Event number.
std::vector< std::string > attribute_string
Attribute serialized as string.
std::vector< GenParticleData > particles
Particles.
std::vector< int > links1
First id of the vertex links.
std::vector< std::string > attribute_name
Attribute name.
Units::LengthUnit length_unit
Length unit.
std::vector< int > attribute_id
Attribute owner id.
FourVector event_pos
Event position.
std::vector< double > weights
Weights.
Units::MomentumUnit momentum_unit
Momentum unit.
Stores serializable particle information.
Stores serializable vertex information.
Comparison of two particle by id.
bool operator()(const GenParticlePtr &p1, const GenParticlePtr &p2)
Comparison of two particle by id.