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);
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;
255 static 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 =
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);