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;
253static 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)
417 return std::const_pointer_cast<const GenVertex>(
m_rootvertex)->particles_out();
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 )
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];
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;
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];
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();
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;
890 { att.second->m_particle =
m_particles[att.first - 1]; }
893 att.second->m_vertex =
m_vertices[-att.first - 1];
#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 vertices_size() const
Vertices size, HepMC2 compatiility.
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)
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.
int particles_size() const
Particles size, HepMC2 compatiility.
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 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.
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)
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.