24 : m_event_number(0), m_weights(std::vector<double>()),
25 m_momentum_unit(mu), m_length_unit(lu),
32 : m_event_number(0), m_weights(std::vector<double>()),
33 m_momentum_unit(mu), m_length_unit(lu),
36 if ( run && !run->weight_names().empty() )
37 m_weights = std::vector<double>(run->weight_names().size(), 1.0);
41 return *(
reinterpret_cast<const std::vector<ConstGenParticlePtr>*
>(&
m_particles));
45 return *(
reinterpret_cast<const std::vector<ConstGenVertexPtr>*
>(&
m_vertices));
51 if( !p|| p->in_event() )
return;
59 if( !p->production_vertex() )
69 std::lock_guard<std::recursive_mutex> rhs_lk(e.
m_lock_attributes, std::adopt_lock);
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;
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;
89 std::lock_guard<std::recursive_mutex> rhs_lk(e.
m_lock_attributes, std::adopt_lock);
99 if( !v|| v->in_event() )
return;
106 for(
auto p: v->particles_in() ) {
108 p->m_end_vertex = v->shared_from_this();
111 for(
auto p: v->particles_out() ) {
113 p->m_production_vertex = v;
119 if( !p || p->parent_event() !=
this )
return;
121 DEBUG( 30,
"GenEvent::remove_particle - called with particle: "<<p->id() );
122 GenVertexPtr end_vtx = p->end_vertex();
124 end_vtx->remove_particle_in(p);
127 if( end_vtx->particles_in().size() == 0 )
remove_vertex(end_vtx);
130 GenVertexPtr prod_vtx = p->production_vertex();
132 prod_vtx->remove_particle_out(p);
135 if( prod_vtx->particles_out().size() == 0 )
remove_vertex(prod_vtx);
138 DEBUG( 30,
"GenEvent::remove_particle - erasing particle: " << p->id() )
145 vector<string> atts = p->attribute_names();
146 for(
const string &s: atts) {
147 p->remove_attribute(s);
153 vector< pair< int, shared_ptr<Attribute> > > changed_attributes;
156 changed_attributes.clear();
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);
164 for(
att_val_t val: changed_attributes ) {
165 vt1.second.erase(val.first);
166 vt1.second[val.first-1] = val.second;
175 p->m_event =
nullptr;
181 inline bool operator()(
const GenParticlePtr& p1,
const GenParticlePtr& p2) {
182 return (p1->id() > p2->id());
190 for (std::vector<GenParticlePtr>::iterator p=v.begin(); p!=v.end(); ++p) {
196 if( !v || v->parent_event() !=
this )
return;
198 DEBUG( 30,
"GenEvent::remove_vertex - called with vertex: "<<v->id() );
199 shared_ptr<GenVertex> null_vtx;
201 for(
auto p: v->particles_in() ) {
202 p->m_end_vertex = std::weak_ptr<GenVertex>();
205 for(
auto p: v->particles_out() ) {
206 p->m_production_vertex = std::weak_ptr<GenVertex>();
213 DEBUG( 30,
"GenEvent::remove_vertex - erasing vertex: " << v->id() )
219 vector<string> atts = v->attribute_names();
220 for(
string s: atts) {
221 v->remove_attribute(s);
228 vector< pair< int, shared_ptr<Attribute> > > changed_attributes;
231 changed_attributes.clear();
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);
239 for(
att_val_t val: changed_attributes ) {
240 vt1.second.erase(val.first);
241 vt1.second[val.first+1] = val.second;
251 v->m_event =
nullptr;
256static bool visit_children(std::map<ConstGenVertexPtr,int> &a, ConstGenVertexPtr v)
258 for ( ConstGenParticlePtr p: v->particles_out())
261 if (a[p->end_vertex()]!=0)
return true;
262 else a[p->end_vertex()]++;
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;
277 for(GenParticlePtr p: parts ) {
278 GenVertexPtr v = p->production_vertex();
280 if( !v || v->particles_in().size()==0 ) {
281 GenVertexPtr v2 = p->end_vertex();
282 if(v2) {noinv.push_back(v2); sortingv[v2]=0;}
285 for (GenVertexPtr v: noinv) {
286 std::map<ConstGenVertexPtr,int> sorting_temp(sortingv.begin(), sortingv.end());
298 deque<GenVertexPtr> sorting;
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);
310 unsigned int sorting_loop_count = 0;
311 unsigned int max_deque_size = 0;
315 while( !sorting.empty() ) {
317 if( sorting.size() > max_deque_size ) max_deque_size = sorting.size();
318 ++sorting_loop_count;
321 GenVertexPtr &v = sorting.front();
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);
336 if( added )
continue;
339 if( !v->in_event() ) {
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);
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;
377 WARNING(
"ReaderAsciiHepMC2: Suspicious looking rootvertex found. Will try to cope." )
382 DEBUG( 6,
"GenEvent - particles sorted: "
383 <<this->
particles().size()<<
", max deque size: "
384 <<max_deque_size<<
", iterations: "<<sorting_loop_count )
422 return std::const_pointer_cast<const GenVertex>(
m_rootvertex)->particles_out();
434 if ( v->has_set_position() )
435 v->set_position( v->position() + delta );
445 long double tempX=mom.
x();
446 long double tempY=mom.
y();
447 long double tempZ=mom.
z();
454 long double cosa=cos(delta.
x());
455 long double sina=sin(delta.
x());
457 tempY_= cosa*tempY+sina*tempZ;
458 tempZ_=-sina*tempY+cosa*tempZ;
463 long double cosb=cos(delta.
y());
464 long double sinb=sin(delta.
y());
466 tempX_= cosb*tempX-sinb*tempZ;
467 tempZ_= sinb*tempX+cosb*tempZ;
471 long double cosg=cos(delta.
z());
472 long double sing=sin(delta.
z());
474 tempX_= cosg*tempX+sing*tempY;
475 tempY_=-sing*tempX+cosg*tempY;
480 p->set_momentum(temp);
485 long double tempX=pos.
x();
486 long double tempY=pos.
y();
487 long double tempZ=pos.
z();
494 long double cosa=cos(delta.
x());
495 long double sina=sin(delta.
x());
497 tempY_= cosa*tempY+sina*tempZ;
498 tempZ_=-sina*tempY+cosa*tempZ;
503 long double cosb=cos(delta.
y());
504 long double sinb=sin(delta.
y());
506 tempX_= cosb*tempX-sinb*tempZ;
507 tempZ_= sinb*tempX+cosb*tempZ;
511 long double cosg=cos(delta.
z());
512 long double sing=sin(delta.
z());
514 tempX_= cosg*tempX+sing*tempY;
515 tempY_=-sing*tempX+cosg*tempY;
520 v->set_position(temp);
531 WARNING(
"GenEvent::reflect: wrong axis" )
562 double deltalength2d=delta.
length2();
563 if (deltalength2d>1.0)
565 WARNING(
"GenEvent::boost: wrong large boost vector. Will leave event as is." )
568 if (
std::abs(deltalength2d-1.0)<std::numeric_limits<double>::epsilon())
570 WARNING(
"GenEvent::boost: too large gamma. Will leave event as is." )
573 if (
std::abs(deltalength2d)<std::numeric_limits<double>::epsilon())
575 WARNING(
"GenEvent::boost: wrong small boost vector. Will leave event as is." )
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);
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;
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);
601 p->set_momentum(temp);
625 map< string, map<int, shared_ptr<Attribute> > >::iterator i1 =
629 map<int, shared_ptr<Attribute> >::iterator i2 = i1->second.find(
id);
630 if( i2 == i1->second.end() )
return;
632 i1->second.erase(i2);
636 vector<string> results;
639 for(
const att_val_t& vt2: vt1.second ) {
640 if( vt2.first ==
id ) results.push_back( vt1.first );
666 for( ConstGenParticlePtr p: this->
particles() ) {
670 for(ConstGenVertexPtr v: this->
vertices() ) {
671 data.
vertices.push_back( v->data() );
674 for(ConstGenParticlePtr p: v->particles_in() ) {
675 data.
links1.push_back( p->id() );
676 data.
links2.push_back( v_id );
679 for(ConstGenParticlePtr p: v->particles_out() ) {
680 data.
links1.push_back( v_id );
681 data.
links2.push_back( p->id() );
686 for(
const att_val_t& vt2: vt1.second ) {
690 bool status = vt2.second->to_string(st);
693 WARNING(
"GenEvent::write_data: problem serializing attribute: "<<vt1.first )
718 GenParticlePtr p = make_shared<GenParticle>(pd);
728 GenVertexPtr v = make_shared<GenVertex>(vd);
737 for(
unsigned int i=0; i<data.
links1.size(); ++i) {
745 if ( (id1<0&&id2<0)|| (id1>0&&id2>0) ) {
WARNING(
"GenEvent::read_data: wrong link: "<<id1<<
" "<<id2 );
continue;}
753 for(
unsigned int i=0; i<data.
attribute_id.size(); ++i) {
784 WARNING(
"Attempting to add an empty particle as beam particle. Ignored.")
787 if( p1->in_event())
if (p1->parent_event()!=
this)
789 WARNING(
"Attempting to add particle from another event. Ignored.")
792 if (p1->production_vertex()) p1->production_vertex()->remove_particle_out(p1);
802 std::map< string, std::map<int, shared_ptr<Attribute> > >::iterator i1 =
806 return run_info()->attribute_as_string(name);
811 std::map<int, shared_ptr<Attribute> >::iterator i2 = i1->second.find(
id);
812 if (i2 == i1->second.end() )
return string();
814 if( !i2->second )
return string();
817 i2->second->to_string(ret);
#define DEBUG_CODE_BLOCK(x)
Macro for storing code useful for debugging.
#define WARNING(MESSAGE)
Macro for printing warning messages.
#define 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.
std::map< int, shared_ptr< Attribute > >::value_type att_val_t
Attribute map value type.
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.
void remove_attribute(const string &name, const int &id=0)
Remove attribute.
void write_data(GenEventData &data) const
Fill GenEventData 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< ConstGenParticlePtr > beams() const
Vector of beam particles.
std::vector< double > m_weights
Event weights.
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.
string attribute_as_string(const string &name, const int &id=0) const
Get attribute of any type as string.
void add_attribute(const string &name, const shared_ptr< Attribute > &att, const int &id=0)
Add event attribute to event.
void set_event_number(const int &num)
Set event number.
int m_event_number
Event number.
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
std::map< string, std::map< int, shared_ptr< Attribute > > > m_attributes
Map of event, particle and vertex attributes.
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.
shared_ptr< GenRunInfo > run_info() const
Get a pointer to the the GenRunInfo object.
const FourVector & event_pos() const
Vertex representing the overall event position.
void read_data(const GenEventData &data)
Fill GenEvent based on GenEventData.
Units::MomentumUnit m_momentum_unit
Momentum unit.
bool rotate(const FourVector v)
Rotate event using x,y,z components of v as rotation angles.
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.
const std::vector< double > & weights() const
Get event weight values as a vector.
bool boost(const FourVector v)
Boost event using x,y,z components of v as velocities.
void clear()
Remove contents of this event.
Units::LengthUnit m_length_unit
Length unit.
std::vector< GenParticlePtr > m_particles
List of particles.
std::vector< string > attribute_names(const int &id=0) const
Get list of attribute names.
void add_beam_particle(GenParticlePtr p1)
Add particle to root vertex.
void set_beam_particles(GenParticlePtr p1, GenParticlePtr p2)
Set incoming beam particles.
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
GenEvent & operator=(const GenEvent &)
Assignment operator.
std::map< string, std::map< int, shared_ptr< Attribute > > > attributes() const
Get a copy of the list of attributes.
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
std::map< string, std::map< int, shared_ptr< Attribute > > >::value_type att_key_t
Attribute map key type.
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.