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));
49 if ( !p || p->in_event() )
return;
57 if ( !p->production_vertex() )
67 std::lock_guard<std::recursive_mutex> rhs_lk(e.
m_lock_attributes, std::adopt_lock);
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;
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;
87 std::lock_guard<std::recursive_mutex> rhs_lk(e.
m_lock_attributes, std::adopt_lock);
97 if ( !v|| v->in_event() )
return;
104 for (
auto p: v->particles_in()) {
106 p->m_end_vertex = v->shared_from_this();
109 for (
auto p: v->particles_out()) {
111 p->m_production_vertex = v;
117 if ( !p || p->parent_event() != this )
return;
119 HEPMC3_DEBUG(30,
"GenEvent::remove_particle - called with particle: " << p->id());
120 GenVertexPtr end_vtx = p->end_vertex();
122 end_vtx->remove_particle_in(p);
125 if ( end_vtx->particles_in().size() == 0 )
remove_vertex(end_vtx);
128 GenVertexPtr prod_vtx = p->production_vertex();
130 prod_vtx->remove_particle_out(p);
133 if ( prod_vtx->particles_out().size() == 0 )
remove_vertex(prod_vtx);
136 HEPMC3_DEBUG(30,
"GenEvent::remove_particle - erasing particle: " << p->id())
143 std::vector<std::string> atts = p->attribute_names();
144 for (
const std::string &s: atts) {
145 p->remove_attribute(s);
151 std::vector< std::pair< int, std::shared_ptr<Attribute> > > changed_attributes;
154 changed_attributes.clear();
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);
162 for (
att_val_t val: changed_attributes ) {
163 vt1.second.erase(val.first);
164 vt1.second[val.first-1] = val.second;
173 p->m_event =
nullptr;
179 inline bool operator()(
const GenParticlePtr& p1,
const GenParticlePtr& p2) {
180 return (p1->id() > p2->id());
187 for (std::vector<GenParticlePtr>::iterator p = v.begin(); p != v.end(); ++p) {
193 if ( !v || v->parent_event() != this )
return;
195 HEPMC3_DEBUG(30,
"GenEvent::remove_vertex - called with vertex: " << v->id());
196 std::shared_ptr<GenVertex> null_vtx;
198 for (
auto p: v->particles_in()) {
199 p->m_end_vertex = std::weak_ptr<GenVertex>();
202 for (
auto p: v->particles_out()) {
203 p->m_production_vertex = std::weak_ptr<GenVertex>();
210 HEPMC3_DEBUG(30,
"GenEvent::remove_vertex - erasing vertex: " << v->id())
216 std::vector<std::string> atts = v->attribute_names();
217 for (std::string s: atts) {
218 v->remove_attribute(s);
225 std::vector< std::pair< int, std::shared_ptr<Attribute> > > changed_attributes;
228 changed_attributes.clear();
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);
236 for (
att_val_t val: changed_attributes ) {
237 vt1.second.erase(val.first);
238 vt1.second[val.first+1] = val.second;
248 v->m_event =
nullptr;
253 static bool visit_children(std::map<ConstGenVertexPtr, int> &a, ConstGenVertexPtr v)
255 for ( ConstGenParticlePtr p: v->particles_out())
258 if (a[p->end_vertex()] != 0)
return true;
259 else a[p->end_vertex()]++;
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;
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;}
281 for (GenVertexPtr v: noinv) {
282 std::map<ConstGenVertexPtr, int> sorting_temp(sortingv.begin(), sortingv.end());
294 std::deque<GenVertexPtr> sorting;
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);
306 unsigned int sorting_loop_count = 0;
307 unsigned int max_deque_size = 0;
311 while ( !sorting.empty() ) {
313 if ( sorting.size() > max_deque_size ) max_deque_size = sorting.size();
314 ++sorting_loop_count;
317 GenVertexPtr &v = sorting.front();
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);
332 if ( added )
continue;
335 if ( !v->in_event() ) {
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);
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;
372 HEPMC3_WARNING(
"ReaderAsciiHepMC2: Suspicious looking rootvertex found. Will try to cope.")
378 << this->
particles().size() <<
", max deque size: "
379 << max_deque_size <<
", iterations: " << sorting_loop_count)
429 if ( v->has_set_position() )
430 v->set_position(v->position() + delta);
439 long double tempX = mom.
x();
440 long double tempY = mom.
y();
441 long double tempZ = mom.
z();
448 long double cosa = std::cos(delta.
x());
449 long double sina = std::sin(delta.
x());
451 tempY_ = cosa*tempY+sina*tempZ;
452 tempZ_ = -sina*tempY+cosa*tempZ;
457 long double cosb = std::cos(delta.
y());
458 long double sinb = std::sin(delta.
y());
460 tempX_ = cosb*tempX-sinb*tempZ;
461 tempZ_ = sinb*tempX+cosb*tempZ;
465 long double cosg = std::cos(delta.
z());
466 long double sing = std::sin(delta.
z());
468 tempX_ = cosg*tempX+sing*tempY;
469 tempY_ = -sing*tempX+cosg*tempY;
474 p->set_momentum(temp);
479 long double tempX = pos.
x();
480 long double tempY = pos.
y();
481 long double tempZ = pos.
z();
488 long double cosa = std::cos(delta.
x());
489 long double sina = std::sin(delta.
x());
491 tempY_ = cosa*tempY+sina*tempZ;
492 tempZ_ = -sina*tempY+cosa*tempZ;
497 long double cosb = std::cos(delta.
y());
498 long double sinb = std::sin(delta.
y());
500 tempX_ = cosb*tempX-sinb*tempZ;
501 tempZ_ = sinb*tempX+cosb*tempZ;
505 long double cosg = std::cos(delta.
z());
506 long double sing = std::sin(delta.
z());
508 tempX_ = cosg*tempX+sing*tempY;
509 tempY_ = -sing*tempX+cosg*tempY;
514 v->set_position(temp);
523 if ( axis > 3 || axis < 0 )
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);}
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);}
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);}
555 double deltalength2d = delta.
length2();
556 if (deltalength2d > 1.0)
558 HEPMC3_WARNING(
"GenEvent::boost: wrong large boost vector. Will leave event as is.")
561 if (
std::abs(deltalength2d-1.0) < std::numeric_limits<double>::epsilon())
563 HEPMC3_WARNING(
"GenEvent::boost: too large gamma. Will leave event as is.")
566 if (
std::abs(deltalength2d) < std::numeric_limits<double>::epsilon())
568 HEPMC3_WARNING(
"GenEvent::boost: wrong small boost vector. Will leave event as is.")
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);
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;
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);
593 p->set_momentum(temp);
611 std:: map< std::string, std::map<int, std::shared_ptr<Attribute> > >::iterator i1 =
615 std::map<int, std::shared_ptr<Attribute> >::iterator i2 = i1->second.find(
id);
616 if ( i2 == i1->second.end() )
return;
618 i1->second.erase(i2);
622 std::vector<std::string> results;
625 if ( vt1.second.count(
id) == 1 ) {
626 results.push_back(vt1.first);
652 for (ConstGenParticlePtr p: this->
particles()) {
656 for (ConstGenVertexPtr v: this->
vertices()) {
660 for (ConstGenParticlePtr p: v->particles_in()) {
661 data.
links1.push_back(p->id());
662 data.
links2.push_back(v_id);
665 for (ConstGenParticlePtr p: v->particles_out()) {
666 data.
links1.push_back(v_id);
667 data.
links2.push_back(p->id());
675 bool status = vt2.second->to_string(st);
678 HEPMC3_WARNING(
"GenEvent::write_data: problem serializing attribute: " << vt1.first)
703 GenParticlePtr p = std::make_shared<GenParticle>(pd);
713 GenVertexPtr v = std::make_shared<GenVertex>(vd);
722 for (
unsigned int i = 0; i < data.
links1.size(); ++i) {
730 if ((id1 < 0 && id2 <0) || (id1 > 0 && id2 > 0)) {
HEPMC3_WARNING(
"GenEvent::read_data: wrong link: " << id1 <<
" " << id2);
continue;}
738 for (
unsigned int i = 0; i < data.
attribute_id.size(); ++i) {
768 HEPMC3_WARNING(
"Attempting to add an empty particle as beam particle. Ignored.")
771 if (p1->in_event())
if (p1->parent_event() !=
this)
773 HEPMC3_WARNING(
"Attempting to add particle from another event. Ignored.")
776 if (p1->production_vertex()) p1->production_vertex()->remove_particle_out(p1);
786 std::map< std::string, std::map<int, std::shared_ptr<Attribute> > >::iterator i1 =
790 return run_info()->attribute_as_string(name);
792 return std::string();
795 std::map<int, std::shared_ptr<Attribute> >::iterator i2 = i1->second.find(
id);
796 if (i2 == i1->second.end() )
return std::string();
798 if ( !i2->second )
return std::string();
801 i2->second->to_string(ret);
808 if (name.length() == 0)
return;
814 if (
id > 0 &&
id <=
int(
particles().size()) )
816 if (
id < 0 && -
id <=
int(
vertices().size()) )
817 att->m_vertex =
vertices()[-
id - 1];
821 void 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;
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));
832 for (
auto name: unames)
837 for (
size_t i = 0; i < N; i++) {
839 if (names.at(i).length() == 0)
continue;
840 if (!atts[i])
continue;
842 atts[i]->m_event =
this;
844 { atts[i]->m_particle =
m_particles[ids.at(i) - 1]; }
847 atts[i]->m_vertex =
m_vertices[-ids.at(i) - 1];
852 void 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();
856 if ( N != atts.size())
return;
863 for (
size_t i = 0; i < N; i++) {
865 if (!atts[i])
continue;
866 tmap[ids.at(i)] = atts[i];
867 atts[i]->m_event =
this;
869 { atts[i]->m_particle =
m_particles[ids.at(i) - 1]; }
872 atts[i]->m_vertex =
m_vertices[-ids.at(i) - 1];
877 if (name.length() == 0)
return;
878 if (!atts.size())
return;
884 for (
auto att: atts) {
886 if (!att.second)
continue;
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]; }
892 if ( att.first < 0 && -att.first <= vertices_size )
893 att.second->m_vertex =
m_vertices[-att.first - 1];
Units::MomentumUnit m_momentum_unit
Momentum unit.
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
int particles_size() const
Particles size, HepMC2 compatiility.
void remove_particle(GenParticlePtr v)
Remove particle from the event.
double length2() const
Squared magnitude of (x, y, z) 3-vector.
void add_beam_particle(GenParticlePtr p1)
Add particle to root vertex.
#define HEPMC3_DEBUG_CODE_BLOCK(x)
Macro for storing code useful for debugging.
void remove_particles(std::vector< GenParticlePtr > v)
Remove a set of particles.
#define HEPMC3_WARNING(MESSAGE)
Macro for printing HEPMC3_HEPMC3_WARNING messages.
void add_vertex(GenVertexPtr v)
Add vertex.
std::vector< int > attribute_id
Attribute owner id.
bool boost(const FourVector &v)
Boost event using x,y,z components of v as velocities.
double t() const
Time component of position/displacement.
Definition of class GenParticle.
void remove_vertex(GenVertexPtr v)
Remove vertex from the event.
void shift_position_by(const FourVector &delta)
Shift position of all vertices in the event by delta.
GenEvent(Units::MomentumUnit momentum_unit=Units::GEV, Units::LengthUnit length_unit=Units::MM)
Event constructor without a run.
Stores vertex-related information.
std::vector< int > links2
Second id of the vertex links.
std::vector< GenVertexPtr > m_vertices
List of vertices.
static void convert(T &m, MomentumUnit from, MomentumUnit to)
Convert FourVector to different momentum unit.
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 of class GenVertex.
bool is_zero() const
Check if the length of this vertex is zero.
#define HEPMC3_DEBUG(LEVEL, MESSAGE)
Macro for printing debug messages with appropriate debug level.
Stores particle-related information.
const Units::LengthUnit & length_unit() const
Get length unit.
double z() const
z-component of position/displacement
int event_number
Event number.
bool reflect(const int axis)
Change sign of axis.
double x() const
x-component of position/displacement
void add_particle(GenParticlePtr p)
Add particle.
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > >::value_type att_key_t
Attribute map key type.
Definition of struct GenEventData.
int event_number() const
Get event number.
void add_attribute(const std::string &name, const std::shared_ptr< Attribute > &att, const int &id=0)
bool operator()(const GenParticlePtr &p1, const GenParticlePtr &p2)
Comparison of two particle by id.
Units::MomentumUnit momentum_unit
Momentum unit.
const FourVector & event_pos() const
Vertex representing the overall event position.
std::vector< int > links1
First id of the vertex links.
FourVector event_pos
Event position.
std::vector< std::string > attribute_string
Attribute serialized as string.
std::map< int, std::shared_ptr< Attribute > >::value_type att_val_t
Attribute map value type.
LengthUnit
Position units.
std::vector< ConstGenParticlePtr > beams() const
Vector of beam particles.
std::string attribute_as_string(const std::string &name, const int &id=0) const
Get attribute of any type as string.
MomentumUnit
Momentum units.
const Units::MomentumUnit & momentum_unit() const
Get momentum unit.
double e() const
Energy component of momentum.
Stores event-related information.
Stores serializable event information.
void write_data(GenEventData &data) const
Fill GenEventData object.
std::recursive_mutex m_lock_attributes
Mutex lock for the m_attibutes map.
std::vector< std::string > attribute_name
Attribute name.
Stores serializable vertex information.
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.
Stores serializable particle information.
static bool visit_children(std::map< ConstGenVertexPtr, int > &a, ConstGenVertexPtr v)
void read_data(const GenEventData &data)
Fill GenEvent based on GenEventData.
std::vector< GenParticleData > particles
Particles.
void add_tree(const std::vector< GenParticlePtr > &particles)
Add whole tree in topological order.
void remove_attribute(const std::string &name, const int &id=0)
Remove attribute.
Units::LengthUnit m_length_unit
Length unit.
std::vector< GenParticlePtr > m_particles
List of particles.
void reserve(const size_t &particles, const size_t &vertices=0)
Reserve memory for particles and vertices.
void set_units(Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit)
Change event units Converts event from current units to new ones.
int m_event_number
Event number.
double y() const
y-component of position/displacement
std::vector< double > weights
Weights.
std::vector< GenVertexData > vertices
Vertices.
std::vector< std::string > attribute_names(const int &id=0) const
Get list of attribute names.
const std::vector< double > & weights() const
Get event weight values as a vector.
int vertices_size() const
Vertices size, HepMC2 compatiility.
GenVertexPtr m_rootvertex
The root vertex is stored outside the normal vertices list to block user access to it...
Definition of class GenEvent.
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > > attributes() const
Get a copy of the list of attributes.
GenEvent & operator=(const GenEvent &)
Assignment operator.
std::map< std::string, std::map< int, std::shared_ptr< Attribute > > > m_attributes
Map of event, particle and vertex attributes.
Comparison of two particle by id.
void clear()
Remove contents of this event.
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
std::vector< double > m_weights
Event weights.
std::shared_ptr< GenRunInfo > run_info() const
Get a pointer to the the GenRunInfo object.
Units::LengthUnit length_unit
Length unit.
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...
void set_event_number(const int &num)
Set event number.
void shift_position_to(const FourVector &newpos)
Shift position of all vertices in the event to op.