[95a917c] | 1 | // -*- C++ -*-
|
---|
| 2 | //
|
---|
| 3 | // This file is part of HepMC
|
---|
| 4 | // Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details)
|
---|
| 5 | //
|
---|
| 6 | /**
|
---|
| 7 | * @file GenEvent.cc
|
---|
| 8 | * @brief Implementation of \b class GenEvent
|
---|
| 9 | *
|
---|
| 10 | */
|
---|
| 11 | #include "HepMC3/GenEvent.h"
|
---|
| 12 | #include "HepMC3/GenParticle.h"
|
---|
| 13 | #include "HepMC3/GenVertex.h"
|
---|
| 14 | #include "HepMC3/Data/GenEventData.h"
|
---|
| 15 |
|
---|
| 16 | #include <deque>
|
---|
| 17 | #include <algorithm> // sort
|
---|
| 18 |
|
---|
| 19 | namespace HepMC3 {
|
---|
| 20 |
|
---|
| 21 | GenEvent::GenEvent(Units::MomentumUnit mu,
|
---|
| 22 | Units::LengthUnit lu)
|
---|
| 23 | : m_event_number(0), m_weights(std::vector<double>()), //m_weights(std::vector<double>(1, 1.0)),//Prevent from different number of weights and names
|
---|
| 24 | m_momentum_unit(mu), m_length_unit(lu),
|
---|
| 25 | m_rootvertex(std::make_shared<GenVertex>()) {}
|
---|
| 26 |
|
---|
| 27 |
|
---|
| 28 | GenEvent::GenEvent(std::shared_ptr<GenRunInfo> run,
|
---|
| 29 | Units::MomentumUnit mu,
|
---|
| 30 | Units::LengthUnit lu)
|
---|
| 31 | : m_event_number(0), m_weights(std::vector<double>()), //m_weights(std::vector<double>(1, 1.0)),//Prevent from different number of weights and names
|
---|
| 32 | m_momentum_unit(mu), m_length_unit(lu),
|
---|
| 33 | m_rootvertex(std::make_shared<GenVertex>()),
|
---|
| 34 | m_run_info(run) {
|
---|
| 35 | if ( run && !run->weight_names().empty() )
|
---|
| 36 | m_weights = std::vector<double>(run->weight_names().size(), 1.0);
|
---|
| 37 | }
|
---|
| 38 |
|
---|
| 39 | const std::vector<ConstGenParticlePtr>& GenEvent::particles() const {
|
---|
| 40 | return *(reinterpret_cast<const std::vector<ConstGenParticlePtr>*>(&m_particles));
|
---|
| 41 | }
|
---|
| 42 |
|
---|
| 43 | const std::vector<ConstGenVertexPtr>& GenEvent::vertices() const {
|
---|
| 44 | return *(reinterpret_cast<const std::vector<ConstGenVertexPtr>*>(&m_vertices));
|
---|
| 45 | }
|
---|
| 46 |
|
---|
| 47 |
|
---|
| 48 | // void GenEvent::add_particle( const GenParticlePtr &p ) {
|
---|
| 49 | void GenEvent::add_particle( GenParticlePtr p ) {
|
---|
| 50 | if( !p|| p->in_event() ) return;
|
---|
| 51 |
|
---|
| 52 | m_particles.push_back(p);
|
---|
| 53 |
|
---|
| 54 | p->m_event = this;
|
---|
| 55 | p->m_id = particles().size();
|
---|
| 56 |
|
---|
| 57 | // Particles without production vertex are added to the root vertex
|
---|
| 58 | if( !p->production_vertex() )
|
---|
| 59 | m_rootvertex->add_particle_out(p);
|
---|
| 60 | }
|
---|
| 61 |
|
---|
| 62 |
|
---|
| 63 | GenEvent::GenEvent(const GenEvent&e) {
|
---|
| 64 | if (this != &e)
|
---|
| 65 | {
|
---|
| 66 | std::lock(m_lock_attributes, e.m_lock_attributes);
|
---|
| 67 | std::lock_guard<std::recursive_mutex> lhs_lk(m_lock_attributes, std::adopt_lock);
|
---|
| 68 | std::lock_guard<std::recursive_mutex> rhs_lk(e.m_lock_attributes, std::adopt_lock);
|
---|
| 69 | GenEventData tdata;
|
---|
| 70 | e.write_data(tdata);
|
---|
| 71 | read_data(tdata);
|
---|
| 72 | }
|
---|
| 73 | }
|
---|
| 74 |
|
---|
| 75 | GenEvent::~GenEvent() {
|
---|
| 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;
|
---|
| 78 |
|
---|
| 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;
|
---|
| 81 | }
|
---|
| 82 |
|
---|
| 83 | GenEvent& GenEvent::operator=(const GenEvent& e) {
|
---|
| 84 | if (this != &e)
|
---|
| 85 | {
|
---|
| 86 | std::lock(m_lock_attributes, e.m_lock_attributes);
|
---|
| 87 | std::lock_guard<std::recursive_mutex> lhs_lk(m_lock_attributes, std::adopt_lock);
|
---|
| 88 | std::lock_guard<std::recursive_mutex> rhs_lk(e.m_lock_attributes, std::adopt_lock);
|
---|
| 89 | GenEventData tdata;
|
---|
| 90 | e.write_data(tdata);
|
---|
| 91 | read_data(tdata);
|
---|
| 92 | }
|
---|
| 93 | return *this;
|
---|
| 94 | }
|
---|
| 95 |
|
---|
| 96 |
|
---|
| 97 | void GenEvent::add_vertex( GenVertexPtr v ) {
|
---|
| 98 | if( !v|| v->in_event() ) return;
|
---|
| 99 | m_vertices.push_back(v);
|
---|
| 100 |
|
---|
| 101 | v->m_event = this;
|
---|
| 102 | v->m_id = -(int)vertices().size();
|
---|
| 103 |
|
---|
| 104 | // Add all incoming and outgoing particles and restore their production/end vertices
|
---|
| 105 | for(auto p: v->particles_in() ) {
|
---|
| 106 | if(!p->in_event()) add_particle(p);
|
---|
| 107 | p->m_end_vertex = v->shared_from_this();
|
---|
| 108 | }
|
---|
| 109 |
|
---|
| 110 | for(auto p: v->particles_out() ) {
|
---|
| 111 | if(!p->in_event()) add_particle(p);
|
---|
| 112 | p->m_production_vertex = v;
|
---|
| 113 | }
|
---|
| 114 | }
|
---|
| 115 |
|
---|
| 116 |
|
---|
| 117 | void GenEvent::remove_particle( GenParticlePtr p ) {
|
---|
| 118 | if( !p || p->parent_event() != this ) return;
|
---|
| 119 |
|
---|
| 120 | HEPMC3_DEBUG( 30, "GenEvent::remove_particle - called with particle: "<<p->id() );
|
---|
| 121 | GenVertexPtr end_vtx = p->end_vertex();
|
---|
| 122 | if( end_vtx ) {
|
---|
| 123 | end_vtx->remove_particle_in(p);
|
---|
| 124 |
|
---|
| 125 | // If that was the only incoming particle, remove vertex from the event
|
---|
| 126 | if( end_vtx->particles_in().size() == 0 ) remove_vertex(end_vtx);
|
---|
| 127 | }
|
---|
| 128 |
|
---|
| 129 | GenVertexPtr prod_vtx = p->production_vertex();
|
---|
| 130 | if( prod_vtx ) {
|
---|
| 131 | prod_vtx->remove_particle_out(p);
|
---|
| 132 |
|
---|
| 133 | // If that was the only outgoing particle, remove vertex from the event
|
---|
| 134 | if( prod_vtx->particles_out().size() == 0 ) remove_vertex(prod_vtx);
|
---|
| 135 | }
|
---|
| 136 |
|
---|
| 137 | HEPMC3_DEBUG( 30, "GenEvent::remove_particle - erasing particle: " << p->id() )
|
---|
| 138 |
|
---|
| 139 | int idx = p->id();
|
---|
| 140 | std::vector<GenParticlePtr>::iterator it = m_particles.erase(m_particles.begin() + idx-1 );
|
---|
| 141 |
|
---|
| 142 | // Remove attributes of this particle
|
---|
| 143 | std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
|
---|
| 144 | std::vector<std::string> atts = p->attribute_names();
|
---|
| 145 | for(const std::string &s: atts) {
|
---|
| 146 | p->remove_attribute(s);
|
---|
| 147 | }
|
---|
| 148 |
|
---|
| 149 | //
|
---|
| 150 | // Reassign id of attributes with id above this one
|
---|
| 151 | //
|
---|
| 152 | std::vector< std::pair< int, std::shared_ptr<Attribute> > > changed_attributes;
|
---|
| 153 |
|
---|
| 154 | for(att_key_t& vt1: m_attributes ) {
|
---|
| 155 | changed_attributes.clear();
|
---|
| 156 |
|
---|
| 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);
|
---|
| 160 | }
|
---|
| 161 | }
|
---|
| 162 |
|
---|
| 163 | for( att_val_t val: changed_attributes ) {
|
---|
| 164 | vt1.second.erase(val.first);
|
---|
| 165 | vt1.second[val.first-1] = val.second;
|
---|
| 166 | }
|
---|
| 167 | }
|
---|
| 168 | // Reassign id of particles with id above this one
|
---|
| 169 | for(; it != m_particles.end(); ++it) {
|
---|
| 170 | --((*it)->m_id);
|
---|
| 171 | }
|
---|
| 172 |
|
---|
| 173 | // Finally - set parent event and id of this particle to 0
|
---|
| 174 | p->m_event = nullptr;
|
---|
| 175 | p->m_id = 0;
|
---|
| 176 | }
|
---|
| 177 | /** @brief Comparison of two particle by id */
|
---|
| 178 | struct sort_by_id_asc {
|
---|
| 179 | /** @brief Comparison of two particle by id */
|
---|
| 180 | inline bool operator()(const GenParticlePtr& p1, const GenParticlePtr& p2) {
|
---|
| 181 | return (p1->id() > p2->id());
|
---|
| 182 | }
|
---|
| 183 | };
|
---|
| 184 |
|
---|
| 185 | void GenEvent::remove_particles( std::vector<GenParticlePtr> v ) {
|
---|
| 186 |
|
---|
| 187 | std::sort( v.begin(), v.end(), sort_by_id_asc() );
|
---|
| 188 |
|
---|
| 189 | for (std::vector<GenParticlePtr>::iterator p=v.begin(); p!=v.end(); ++p) {
|
---|
| 190 | remove_particle(*p);
|
---|
| 191 | }
|
---|
| 192 | }
|
---|
| 193 |
|
---|
| 194 | void GenEvent::remove_vertex( GenVertexPtr v ) {
|
---|
| 195 | if( !v || v->parent_event() != this ) return;
|
---|
| 196 |
|
---|
| 197 | HEPMC3_DEBUG( 30, "GenEvent::remove_vertex - called with vertex: "<<v->id() );
|
---|
| 198 | std::shared_ptr<GenVertex> null_vtx;
|
---|
| 199 |
|
---|
| 200 | for(auto p: v->particles_in() ) {
|
---|
| 201 | p->m_end_vertex = std::weak_ptr<GenVertex>();
|
---|
| 202 | }
|
---|
| 203 |
|
---|
| 204 | for(auto p: v->particles_out() ) {
|
---|
| 205 | p->m_production_vertex = std::weak_ptr<GenVertex>();
|
---|
| 206 |
|
---|
| 207 | // recursive delete rest of the tree
|
---|
| 208 | remove_particle(p);
|
---|
| 209 | }
|
---|
| 210 |
|
---|
| 211 | // Erase this vertex from vertices list
|
---|
| 212 | HEPMC3_DEBUG( 30, "GenEvent::remove_vertex - erasing vertex: " << v->id() )
|
---|
| 213 |
|
---|
| 214 | int idx = -v->id();
|
---|
| 215 | std::vector<GenVertexPtr>::iterator it = m_vertices.erase(m_vertices.begin() + idx-1 );
|
---|
| 216 | // Remove attributes of this vertex
|
---|
| 217 | std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
|
---|
| 218 | std::vector<std::string> atts = v->attribute_names();
|
---|
| 219 | for(std::string s: atts) {
|
---|
| 220 | v->remove_attribute(s);
|
---|
| 221 | }
|
---|
| 222 |
|
---|
| 223 | //
|
---|
| 224 | // Reassign id of attributes with id below this one
|
---|
| 225 | //
|
---|
| 226 |
|
---|
| 227 | std::vector< std::pair< int, std::shared_ptr<Attribute> > > changed_attributes;
|
---|
| 228 |
|
---|
| 229 | for( att_key_t& vt1: m_attributes ) {
|
---|
| 230 | changed_attributes.clear();
|
---|
| 231 |
|
---|
| 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);
|
---|
| 235 | }
|
---|
| 236 | }
|
---|
| 237 |
|
---|
| 238 | for( att_val_t val: changed_attributes ) {
|
---|
| 239 | vt1.second.erase(val.first);
|
---|
| 240 | vt1.second[val.first+1] = val.second;
|
---|
| 241 | }
|
---|
| 242 | }
|
---|
| 243 |
|
---|
| 244 | // Reassign id of particles with id above this one
|
---|
| 245 | for(; it != m_vertices.end(); ++it) {
|
---|
| 246 | ++((*it)->m_id);
|
---|
| 247 | }
|
---|
| 248 |
|
---|
| 249 | // Finally - set parent event and id of this vertex to 0
|
---|
| 250 | v->m_event = nullptr;
|
---|
| 251 | v->m_id = 0;
|
---|
| 252 | }
|
---|
| 253 | /// @todo This looks dangerously similar to the recusive event traversel that we forbade in the
|
---|
| 254 | /// Core library due to wories about generator dependence
|
---|
| 255 | static bool visit_children(std::map<ConstGenVertexPtr,int> &a, ConstGenVertexPtr v)
|
---|
| 256 | {
|
---|
| 257 | for ( ConstGenParticlePtr p: v->particles_out())
|
---|
| 258 | if (p->end_vertex())
|
---|
| 259 | {
|
---|
| 260 | if (a[p->end_vertex()]!=0) return true;
|
---|
| 261 | else a[p->end_vertex()]++;
|
---|
| 262 | if (visit_children(a, p->end_vertex())) return true;
|
---|
| 263 | }
|
---|
| 264 | return false;
|
---|
| 265 | }
|
---|
| 266 |
|
---|
| 267 | void GenEvent::add_tree( const std::vector<GenParticlePtr> &parts ) {
|
---|
| 268 |
|
---|
| 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;
|
---|
| 274 | if(!existing_hc)
|
---|
| 275 | {
|
---|
| 276 | for(GenParticlePtr p: parts ) {
|
---|
| 277 | GenVertexPtr v = p->production_vertex();
|
---|
| 278 | if(v) sortingv[v]=0;
|
---|
| 279 | if( !v || v->particles_in().size()==0 ) {
|
---|
| 280 | GenVertexPtr v2 = p->end_vertex();
|
---|
| 281 | if(v2) {noinv.push_back(v2); sortingv[v2]=0;}
|
---|
| 282 | }
|
---|
| 283 | }
|
---|
| 284 | for (GenVertexPtr v: noinv) {
|
---|
| 285 | std::map<ConstGenVertexPtr,int> sorting_temp(sortingv.begin(), sortingv.end());
|
---|
| 286 | has_cycles=(has_cycles||visit_children(sorting_temp, v));
|
---|
| 287 | }
|
---|
| 288 | }
|
---|
| 289 | if (has_cycles) {
|
---|
| 290 | add_attribute("cycles", std::make_shared<IntAttribute>(1));
|
---|
| 291 | /* Commented out as improvemnts allow us to do sorting in other way.
|
---|
| 292 | for( std::map<GenVertexPtr,int>::iterator vi=sortingv.begin();vi!=sortingv.end();++vi) if( !vi->first->in_event() ) add_vertex(vi->first);
|
---|
| 293 | return;
|
---|
| 294 | */
|
---|
| 295 | }
|
---|
| 296 |
|
---|
| 297 | std::deque<GenVertexPtr> sorting;
|
---|
| 298 |
|
---|
| 299 | // Find all starting vertices (end vertex of particles that have no production vertex)
|
---|
| 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);
|
---|
| 305 | }
|
---|
| 306 | }
|
---|
| 307 |
|
---|
| 308 | HEPMC3_DEBUG_CODE_BLOCK(
|
---|
| 309 | unsigned int sorting_loop_count = 0;
|
---|
| 310 | unsigned int max_deque_size = 0;
|
---|
| 311 | )
|
---|
| 312 |
|
---|
| 313 | // Add vertices to the event in topological order
|
---|
| 314 | while( !sorting.empty() ) {
|
---|
| 315 | HEPMC3_DEBUG_CODE_BLOCK(
|
---|
| 316 | if( sorting.size() > max_deque_size ) max_deque_size = sorting.size();
|
---|
| 317 | ++sorting_loop_count;
|
---|
| 318 | )
|
---|
| 319 |
|
---|
| 320 | GenVertexPtr &v = sorting.front();
|
---|
| 321 |
|
---|
| 322 | bool added = false;
|
---|
| 323 |
|
---|
| 324 | // Add all mothers to the front of the list
|
---|
| 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);
|
---|
| 329 | added = true;
|
---|
| 330 | }
|
---|
| 331 | }
|
---|
| 332 |
|
---|
| 333 | // If we have added at least one production vertex,
|
---|
| 334 | // our vertex is not the first one on the list
|
---|
| 335 | if( added ) continue;
|
---|
| 336 |
|
---|
| 337 | // If vertex not yet added
|
---|
| 338 | if( !v->in_event() ) {
|
---|
| 339 |
|
---|
| 340 | add_vertex(v);
|
---|
| 341 |
|
---|
| 342 | // Add all end vertices to the end of the list
|
---|
| 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);
|
---|
| 347 | }
|
---|
| 348 | }
|
---|
| 349 | }
|
---|
| 350 |
|
---|
| 351 | sorting.pop_front();
|
---|
| 352 | }
|
---|
| 353 |
|
---|
| 354 | // LL: Make sure root vertex has index zero and is not written out
|
---|
| 355 | if ( m_rootvertex->id() != 0 ) {
|
---|
| 356 | const int vx = -1 - m_rootvertex->id();
|
---|
| 357 | const int rootid = m_rootvertex->id();
|
---|
| 358 | if ( vx >= 0 && vx < m_vertices.size() && m_vertices[vx] == m_rootvertex ) {
|
---|
| 359 | auto next = m_vertices.erase(m_vertices.begin() + vx);
|
---|
| 360 | std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
|
---|
| 361 | for(auto & vt1: m_attributes ) {
|
---|
| 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;
|
---|
| 369 | }
|
---|
| 370 | }
|
---|
| 371 | m_rootvertex->set_id(0);
|
---|
| 372 | while ( next != m_vertices.end() ) {
|
---|
| 373 | ++((*next++)->m_id);
|
---|
| 374 | }
|
---|
| 375 | } else {
|
---|
| 376 | HEPMC3_WARNING( "ReaderAsciiHepMC2: Suspicious looking rootvertex found. Will try to cope." )
|
---|
| 377 | }
|
---|
| 378 | }
|
---|
| 379 |
|
---|
| 380 | HEPMC3_DEBUG_CODE_BLOCK(
|
---|
| 381 | HEPMC3_DEBUG( 6, "GenEvent - particles sorted: "
|
---|
| 382 | <<this->particles().size()<<", max deque size: "
|
---|
| 383 | <<max_deque_size<<", iterations: "<<sorting_loop_count )
|
---|
| 384 | )
|
---|
| 385 | return;
|
---|
| 386 | }
|
---|
| 387 |
|
---|
| 388 |
|
---|
| 389 | void GenEvent::reserve(const size_t& parts, const size_t& verts) {
|
---|
| 390 | m_particles.reserve(parts);
|
---|
| 391 | m_vertices.reserve(verts);
|
---|
| 392 | }
|
---|
| 393 |
|
---|
| 394 |
|
---|
| 395 | void GenEvent::set_units( Units::MomentumUnit new_momentum_unit, Units::LengthUnit new_length_unit) {
|
---|
| 396 | if( new_momentum_unit != m_momentum_unit ) {
|
---|
| 397 | for( GenParticlePtr p: m_particles ) {
|
---|
| 398 | Units::convert( p->m_data.momentum, m_momentum_unit, new_momentum_unit );
|
---|
| 399 | Units::convert( p->m_data.mass, m_momentum_unit, new_momentum_unit );
|
---|
| 400 | }
|
---|
| 401 |
|
---|
| 402 | m_momentum_unit = new_momentum_unit;
|
---|
| 403 | }
|
---|
| 404 |
|
---|
| 405 | if( new_length_unit != m_length_unit ) {
|
---|
| 406 | for(GenVertexPtr &v: m_vertices ) {
|
---|
| 407 | FourVector &fv = v->m_data.position;
|
---|
| 408 | if( !fv.is_zero() ) Units::convert( fv, m_length_unit, new_length_unit );
|
---|
| 409 | }
|
---|
| 410 |
|
---|
| 411 | m_length_unit = new_length_unit;
|
---|
| 412 | }
|
---|
| 413 | }
|
---|
| 414 |
|
---|
| 415 |
|
---|
| 416 | const FourVector& GenEvent::event_pos() const {
|
---|
| 417 | return m_rootvertex->data().position;
|
---|
| 418 | }
|
---|
| 419 |
|
---|
| 420 | std::vector<ConstGenParticlePtr> GenEvent::beams() const {
|
---|
| 421 | return std::const_pointer_cast<const GenVertex>(m_rootvertex)->particles_out();
|
---|
| 422 | }
|
---|
| 423 |
|
---|
| 424 | const std::vector<GenParticlePtr> & GenEvent::beams() {
|
---|
| 425 | return m_rootvertex->particles_out();
|
---|
| 426 | }
|
---|
| 427 |
|
---|
| 428 | void GenEvent::shift_position_by( const FourVector & delta ) {
|
---|
| 429 | m_rootvertex->set_position( event_pos() + delta );
|
---|
| 430 |
|
---|
| 431 | // Offset all vertices
|
---|
| 432 | for ( GenVertexPtr v: m_vertices ) {
|
---|
| 433 | if ( v->has_set_position() )
|
---|
| 434 | v->set_position( v->position() + delta );
|
---|
| 435 | }
|
---|
| 436 | }
|
---|
| 437 |
|
---|
| 438 | bool GenEvent::rotate( const FourVector& delta )
|
---|
| 439 | {
|
---|
| 440 |
|
---|
| 441 | for ( auto p: m_particles)
|
---|
| 442 | {
|
---|
| 443 | FourVector mom=p->momentum();
|
---|
| 444 | long double tempX=mom.x();
|
---|
| 445 | long double tempY=mom.y();
|
---|
| 446 | long double tempZ=mom.z();
|
---|
| 447 |
|
---|
| 448 | long double tempX_;
|
---|
| 449 | long double tempY_;
|
---|
| 450 | long double tempZ_;
|
---|
| 451 |
|
---|
| 452 |
|
---|
| 453 | long double cosa=cos(delta.x());
|
---|
| 454 | long double sina=sin(delta.x());
|
---|
| 455 |
|
---|
| 456 | tempY_= cosa*tempY+sina*tempZ;
|
---|
| 457 | tempZ_=-sina*tempY+cosa*tempZ;
|
---|
| 458 | tempY=tempY_;
|
---|
| 459 | tempZ=tempZ_;
|
---|
| 460 |
|
---|
| 461 |
|
---|
| 462 | long double cosb=cos(delta.y());
|
---|
| 463 | long double sinb=sin(delta.y());
|
---|
| 464 |
|
---|
| 465 | tempX_= cosb*tempX-sinb*tempZ;
|
---|
| 466 | tempZ_= sinb*tempX+cosb*tempZ;
|
---|
| 467 | tempX=tempX_;
|
---|
| 468 | tempZ=tempZ_;
|
---|
| 469 |
|
---|
| 470 | long double cosg=cos(delta.z());
|
---|
| 471 | long double sing=sin(delta.z());
|
---|
| 472 |
|
---|
| 473 | tempX_= cosg*tempX+sing*tempY;
|
---|
| 474 | tempY_=-sing*tempX+cosg*tempY;
|
---|
| 475 | tempX=tempX_;
|
---|
| 476 | tempY=tempY_;
|
---|
| 477 |
|
---|
| 478 | FourVector temp(tempX,tempY,tempZ,mom.e());
|
---|
| 479 | p->set_momentum(temp);
|
---|
| 480 | }
|
---|
| 481 | for ( auto v: m_vertices)
|
---|
| 482 | {
|
---|
| 483 | FourVector pos=v->position();
|
---|
| 484 | long double tempX=pos.x();
|
---|
| 485 | long double tempY=pos.y();
|
---|
| 486 | long double tempZ=pos.z();
|
---|
| 487 |
|
---|
| 488 | long double tempX_;
|
---|
| 489 | long double tempY_;
|
---|
| 490 | long double tempZ_;
|
---|
| 491 |
|
---|
| 492 |
|
---|
| 493 | long double cosa=cos(delta.x());
|
---|
| 494 | long double sina=sin(delta.x());
|
---|
| 495 |
|
---|
| 496 | tempY_= cosa*tempY+sina*tempZ;
|
---|
| 497 | tempZ_=-sina*tempY+cosa*tempZ;
|
---|
| 498 | tempY=tempY_;
|
---|
| 499 | tempZ=tempZ_;
|
---|
| 500 |
|
---|
| 501 |
|
---|
| 502 | long double cosb=cos(delta.y());
|
---|
| 503 | long double sinb=sin(delta.y());
|
---|
| 504 |
|
---|
| 505 | tempX_= cosb*tempX-sinb*tempZ;
|
---|
| 506 | tempZ_= sinb*tempX+cosb*tempZ;
|
---|
| 507 | tempX=tempX_;
|
---|
| 508 | tempZ=tempZ_;
|
---|
| 509 |
|
---|
| 510 | long double cosg=cos(delta.z());
|
---|
| 511 | long double sing=sin(delta.z());
|
---|
| 512 |
|
---|
| 513 | tempX_= cosg*tempX+sing*tempY;
|
---|
| 514 | tempY_=-sing*tempX+cosg*tempY;
|
---|
| 515 | tempX=tempX_;
|
---|
| 516 | tempY=tempY_;
|
---|
| 517 |
|
---|
| 518 | FourVector temp(tempX,tempY,tempZ,pos.t());
|
---|
| 519 | v->set_position(temp);
|
---|
| 520 | }
|
---|
| 521 |
|
---|
| 522 |
|
---|
| 523 | return true;
|
---|
| 524 | }
|
---|
| 525 |
|
---|
| 526 | bool GenEvent::reflect(const int axis)
|
---|
| 527 | {
|
---|
| 528 | if (axis>3||axis<0)
|
---|
| 529 | {
|
---|
| 530 | HEPMC3_WARNING( "GenEvent::reflect: wrong axis" )
|
---|
| 531 | return false;
|
---|
| 532 | }
|
---|
| 533 | switch (axis)
|
---|
| 534 | {
|
---|
| 535 | case 0:
|
---|
| 536 | for ( auto p: m_particles) { FourVector temp=p->momentum(); temp.setX(-p->momentum().x()); p->set_momentum(temp);}
|
---|
| 537 | for ( auto v: m_vertices) { FourVector temp=v->position(); temp.setX(-v->position().x()); v->set_position(temp);}
|
---|
| 538 | break;
|
---|
| 539 | case 1:
|
---|
| 540 | for ( auto p: m_particles) { FourVector temp=p->momentum(); temp.setY(-p->momentum().y()); p->set_momentum(temp);}
|
---|
| 541 | for ( auto v: m_vertices) { FourVector temp=v->position(); temp.setY(-v->position().y()); v->set_position(temp);}
|
---|
| 542 | break;
|
---|
| 543 | case 2:
|
---|
| 544 | for ( auto p: m_particles) { FourVector temp=p->momentum(); temp.setZ(-p->momentum().z()); p->set_momentum(temp);}
|
---|
| 545 | for ( auto v: m_vertices) { FourVector temp=v->position(); temp.setZ(-v->position().z()); v->set_position(temp);}
|
---|
| 546 | break;
|
---|
| 547 | case 3:
|
---|
| 548 | for ( auto p: m_particles) { FourVector temp=p->momentum(); temp.setT(-p->momentum().e()); p->set_momentum(temp);}
|
---|
| 549 | for ( auto v: m_vertices) { FourVector temp=v->position(); temp.setT(-v->position().t()); v->set_position(temp);}
|
---|
| 550 | break;
|
---|
| 551 | default:
|
---|
| 552 | return false;
|
---|
| 553 | }
|
---|
| 554 |
|
---|
| 555 | return true;
|
---|
| 556 | }
|
---|
| 557 |
|
---|
| 558 | bool GenEvent::boost( const FourVector& delta )
|
---|
| 559 | {
|
---|
| 560 |
|
---|
| 561 | double deltalength2d=delta.length2();
|
---|
| 562 | if (deltalength2d>1.0)
|
---|
| 563 | {
|
---|
| 564 | HEPMC3_WARNING( "GenEvent::boost: wrong large boost vector. Will leave event as is." )
|
---|
| 565 | return false;
|
---|
| 566 | }
|
---|
| 567 | if (std::abs(deltalength2d-1.0)<std::numeric_limits<double>::epsilon())
|
---|
| 568 | {
|
---|
| 569 | HEPMC3_WARNING( "GenEvent::boost: too large gamma. Will leave event as is." )
|
---|
| 570 | return false;
|
---|
| 571 | }
|
---|
| 572 | if (std::abs(deltalength2d)<std::numeric_limits<double>::epsilon())
|
---|
| 573 | {
|
---|
| 574 | HEPMC3_WARNING( "GenEvent::boost: wrong small boost vector. Will leave event as is." )
|
---|
| 575 | return true;
|
---|
| 576 | }
|
---|
| 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);
|
---|
| 583 |
|
---|
| 584 | for ( auto p: m_particles)
|
---|
| 585 | {
|
---|
| 586 | FourVector mom=p->momentum();
|
---|
| 587 |
|
---|
| 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;
|
---|
| 593 |
|
---|
| 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);
|
---|
| 598 | FourVector temp(tempX,tempY,tempZ,tempE);
|
---|
| 599 | p->set_momentum(temp);
|
---|
| 600 | }
|
---|
| 601 |
|
---|
| 602 | return true;
|
---|
| 603 | }
|
---|
| 604 |
|
---|
| 605 |
|
---|
| 606 |
|
---|
| 607 |
|
---|
| 608 | void GenEvent::clear() {
|
---|
| 609 | std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
|
---|
| 610 | m_event_number = 0;
|
---|
| 611 | m_rootvertex = std::make_shared<GenVertex>();
|
---|
| 612 | m_weights.clear();
|
---|
| 613 | m_attributes.clear();
|
---|
| 614 | m_particles.clear();
|
---|
| 615 | m_vertices.clear();
|
---|
| 616 | }
|
---|
| 617 |
|
---|
| 618 |
|
---|
| 619 |
|
---|
| 620 |
|
---|
| 621 | void GenEvent::remove_attribute(const std::string &name, const int& id) {
|
---|
| 622 | std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
|
---|
| 623 | std:: map< std::string, std::map<int, std::shared_ptr<Attribute> > >::iterator i1 =
|
---|
| 624 | m_attributes.find(name);
|
---|
| 625 | if( i1 == m_attributes.end() ) return;
|
---|
| 626 |
|
---|
| 627 | std::map<int, std::shared_ptr<Attribute> >::iterator i2 = i1->second.find(id);
|
---|
| 628 | if( i2 == i1->second.end() ) return;
|
---|
| 629 |
|
---|
| 630 | i1->second.erase(i2);
|
---|
| 631 | }
|
---|
| 632 |
|
---|
| 633 | std::vector<std::string> GenEvent::attribute_names( const int& id) const {
|
---|
| 634 | std::vector<std::string> results;
|
---|
| 635 |
|
---|
| 636 | for(const att_key_t& vt1: m_attributes ) {
|
---|
| 637 | if( vt1.second.count( id ) == 1 ) {
|
---|
| 638 | results.push_back( vt1.first );
|
---|
| 639 | }
|
---|
| 640 | }
|
---|
| 641 |
|
---|
| 642 | return results;
|
---|
| 643 | }
|
---|
| 644 |
|
---|
| 645 | void GenEvent::write_data(GenEventData& data) const {
|
---|
| 646 | // Reserve memory for containers
|
---|
| 647 | data.particles.reserve( this->particles().size() );
|
---|
| 648 | data.vertices.reserve( this->vertices().size() );
|
---|
| 649 | data.links1.reserve( this->particles().size()*2 );
|
---|
| 650 | data.links2.reserve( this->particles().size()*2 );
|
---|
| 651 | data.attribute_id.reserve( m_attributes.size() );
|
---|
| 652 | data.attribute_name.reserve( m_attributes.size() );
|
---|
| 653 | data.attribute_string.reserve( m_attributes.size() );
|
---|
| 654 |
|
---|
| 655 | // Fill event data
|
---|
| 656 | data.event_number = this->event_number();
|
---|
| 657 | data.momentum_unit = this->momentum_unit();
|
---|
| 658 | data.length_unit = this->length_unit();
|
---|
| 659 | data.event_pos = this->event_pos();
|
---|
| 660 |
|
---|
| 661 | // Fill containers
|
---|
| 662 | data.weights = this->weights();
|
---|
| 663 |
|
---|
| 664 | for( ConstGenParticlePtr p: this->particles() ) {
|
---|
| 665 | data.particles.push_back( p->data() );
|
---|
| 666 | }
|
---|
| 667 |
|
---|
| 668 | for(ConstGenVertexPtr v: this->vertices() ) {
|
---|
| 669 | data.vertices.push_back( v->data() );
|
---|
| 670 | int v_id = v->id();
|
---|
| 671 |
|
---|
| 672 | for(ConstGenParticlePtr p: v->particles_in() ) {
|
---|
| 673 | data.links1.push_back( p->id() );
|
---|
| 674 | data.links2.push_back( v_id );
|
---|
| 675 | }
|
---|
| 676 |
|
---|
| 677 | for(ConstGenParticlePtr p: v->particles_out() ) {
|
---|
| 678 | data.links1.push_back( v_id );
|
---|
| 679 | data.links2.push_back( p->id() );
|
---|
| 680 | }
|
---|
| 681 | }
|
---|
| 682 |
|
---|
| 683 | for( const att_key_t& vt1: this->attributes() ) {
|
---|
| 684 | for( const att_val_t& vt2: vt1.second ) {
|
---|
| 685 |
|
---|
| 686 | std::string st;
|
---|
| 687 |
|
---|
| 688 | bool status = vt2.second->to_string(st);
|
---|
| 689 |
|
---|
| 690 | if( !status ) {
|
---|
| 691 | HEPMC3_WARNING( "GenEvent::write_data: problem serializing attribute: "<<vt1.first )
|
---|
| 692 | }
|
---|
| 693 | else {
|
---|
| 694 | data.attribute_id.push_back(vt2.first);
|
---|
| 695 | data.attribute_name.push_back(vt1.first);
|
---|
| 696 | data.attribute_string.push_back(st);
|
---|
| 697 | }
|
---|
| 698 | }
|
---|
| 699 | }
|
---|
| 700 | }
|
---|
| 701 |
|
---|
| 702 |
|
---|
| 703 | void GenEvent::read_data(const GenEventData &data) {
|
---|
| 704 | this->clear();
|
---|
| 705 | this->set_event_number( data.event_number );
|
---|
| 706 | //Note: set_units checks the current unit of event, i.e. applicable only for fully constructed event.
|
---|
| 707 | m_momentum_unit=data.momentum_unit;
|
---|
| 708 | m_length_unit=data.length_unit;
|
---|
| 709 | this->shift_position_to( data.event_pos );
|
---|
| 710 |
|
---|
| 711 | // Fill weights
|
---|
| 712 | this->weights() = data.weights;
|
---|
| 713 |
|
---|
| 714 | // Fill particle information
|
---|
| 715 | for( const GenParticleData &pd: data.particles ) {
|
---|
| 716 | GenParticlePtr p = std::make_shared<GenParticle>(pd);
|
---|
| 717 |
|
---|
| 718 | m_particles.push_back(p);
|
---|
| 719 |
|
---|
| 720 | p->m_event = this;
|
---|
| 721 | p->m_id = m_particles.size();
|
---|
| 722 | }
|
---|
| 723 |
|
---|
| 724 | // Fill vertex information
|
---|
| 725 | for( const GenVertexData &vd: data.vertices ) {
|
---|
| 726 | GenVertexPtr v = std::make_shared<GenVertex>(vd);
|
---|
| 727 |
|
---|
| 728 | m_vertices.push_back(v);
|
---|
| 729 |
|
---|
| 730 | v->m_event = this;
|
---|
| 731 | v->m_id = -(int)m_vertices.size();
|
---|
| 732 | }
|
---|
| 733 |
|
---|
| 734 | // Restore links
|
---|
| 735 | for( unsigned int i=0; i<data.links1.size(); ++i) {
|
---|
| 736 | int id1 = data.links1[i];
|
---|
| 737 | int id2 = data.links2[i];
|
---|
| 738 | /* @note:
|
---|
| 739 | The meaningfull combinations for (id1,id2) are:
|
---|
| 740 | (+-) -- particle has end vertex
|
---|
| 741 | (-+) -- particle has production vertex
|
---|
| 742 | */
|
---|
| 743 | if ( (id1<0&&id2<0)|| (id1>0&&id2>0) ) { HEPMC3_WARNING( "GenEvent::read_data: wrong link: "<<id1<<" "<<id2 ); continue;}
|
---|
| 744 |
|
---|
| 745 | if ( id1 > 0 ) { m_vertices[ (-id2)-1 ]->add_particle_in ( m_particles[ id1-1 ] ); continue; }
|
---|
| 746 | if ( id1 < 0 ) { m_vertices[ (-id1)-1 ]->add_particle_out( m_particles[ id2-1 ] ); continue; }
|
---|
| 747 | }
|
---|
| 748 | for (auto p: m_particles) if (!p->production_vertex()) m_rootvertex->add_particle_out(p);
|
---|
| 749 |
|
---|
| 750 | // Read attributes
|
---|
| 751 | for( unsigned int i=0; i<data.attribute_id.size(); ++i) {
|
---|
| 752 | add_attribute( data.attribute_name[i],
|
---|
| 753 | std::make_shared<StringAttribute>(data.attribute_string[i]),
|
---|
| 754 | data.attribute_id[i] );
|
---|
| 755 | }
|
---|
| 756 | }
|
---|
| 757 |
|
---|
| 758 |
|
---|
| 759 |
|
---|
| 760 | //
|
---|
| 761 | // Deprecated functions
|
---|
| 762 | //
|
---|
| 763 |
|
---|
| 764 | void GenEvent::add_particle( GenParticle *p ) {
|
---|
| 765 | add_particle( GenParticlePtr(p) );
|
---|
| 766 | }
|
---|
| 767 |
|
---|
| 768 |
|
---|
| 769 | void GenEvent::add_vertex( GenVertex *v ) {
|
---|
| 770 | add_vertex( GenVertexPtr(v) );
|
---|
| 771 | }
|
---|
| 772 |
|
---|
| 773 |
|
---|
| 774 | void GenEvent::set_beam_particles(GenParticlePtr p1, GenParticlePtr p2) {
|
---|
| 775 | m_rootvertex->add_particle_out(p1);
|
---|
| 776 | m_rootvertex->add_particle_out(p2);
|
---|
| 777 | }
|
---|
| 778 |
|
---|
| 779 | void GenEvent::add_beam_particle(GenParticlePtr p1) {
|
---|
| 780 | if (!p1)
|
---|
| 781 | {
|
---|
| 782 | HEPMC3_WARNING("Attempting to add an empty particle as beam particle. Ignored.")
|
---|
| 783 | return;
|
---|
| 784 | }
|
---|
| 785 | if( p1->in_event()) if (p1->parent_event()!=this)
|
---|
| 786 | {
|
---|
| 787 | HEPMC3_WARNING("Attempting to add particle from another event. Ignored.")
|
---|
| 788 | return;
|
---|
| 789 | }
|
---|
| 790 | if (p1->production_vertex()) p1->production_vertex()->remove_particle_out(p1);
|
---|
| 791 | //Particle w/o production vertex is added to root vertex.
|
---|
| 792 | add_particle(p1);
|
---|
| 793 | p1->set_status(4);
|
---|
| 794 | return;
|
---|
| 795 | }
|
---|
| 796 |
|
---|
| 797 |
|
---|
| 798 | std::string GenEvent::attribute_as_string(const std::string &name, const int& id) const {
|
---|
| 799 | std::lock_guard<std::recursive_mutex> lock(m_lock_attributes);
|
---|
| 800 | std::map< std::string, std::map<int, std::shared_ptr<Attribute> > >::iterator i1 =
|
---|
| 801 | m_attributes.find(name);
|
---|
| 802 | if( i1 == m_attributes.end() ) {
|
---|
| 803 | if ( id == 0 && run_info() ) {
|
---|
| 804 | return run_info()->attribute_as_string(name);
|
---|
| 805 | }
|
---|
| 806 | return std::string();
|
---|
| 807 | }
|
---|
| 808 |
|
---|
| 809 | std::map<int, std::shared_ptr<Attribute> >::iterator i2 = i1->second.find(id);
|
---|
| 810 | if (i2 == i1->second.end() ) return std::string();
|
---|
| 811 |
|
---|
| 812 | if( !i2->second ) return std::string();
|
---|
| 813 |
|
---|
| 814 | std::string ret;
|
---|
| 815 | i2->second->to_string(ret);
|
---|
| 816 |
|
---|
| 817 | return ret;
|
---|
| 818 | }
|
---|
| 819 |
|
---|
| 820 | } // namespace HepMC3
|
---|