Rev 124 | Rev 168 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
| Rev | Author | Line No. | Line |
|---|---|---|---|
| 2 | garren | 1 | ////////////////////////////////////////////////////////////////////////// |
| 2 | // Matt.Dobbs@Cern.CH, September 1999 |
||
| 3 | // GenVertex within an event |
||
| 4 | ////////////////////////////////////////////////////////////////////////// |
||
| 5 | |||
| 6 | #include "HepMC/GenParticle.h" |
||
| 7 | #include "HepMC/GenVertex.h" |
||
| 8 | #include "HepMC/GenEvent.h" |
||
| 25 | garren | 9 | #include <iomanip> // needed for formatted output |
| 2 | garren | 10 | |
| 11 | namespace HepMC { |
||
| 12 | |||
| 43 | garren | 13 | GenVertex::GenVertex( const FourVector& position, |
| 2 | garren | 14 | int id, const WeightContainer& weights ) |
| 15 | : m_position(position), m_id(id), m_weights(weights), m_event(0), |
||
| 16 | m_barcode(0) |
||
| 17 | { |
||
| 18 | s_counter++; |
||
| 19 | } |
||
| 20 | |||
| 143 | garren | 21 | GenVertex::GenVertex( const GenVertex& invertex ) |
| 22 | : m_position( invertex.position() ), |
||
| 23 | m_particles_in(), |
||
| 24 | m_particles_out(), |
||
| 25 | m_id( invertex.id() ), |
||
| 26 | m_weights( invertex.weights() ), |
||
| 27 | m_event(0), |
||
| 28 | m_barcode(0) |
||
| 29 | { |
||
| 65 | garren | 30 | /// Shallow copy: does not copy the FULL list of particle pointers. |
| 31 | /// Creates a copy of - invertex |
||
| 32 | /// - outgoing particles of invertex, but sets the |
||
| 33 | /// decay vertex of these particles to NULL |
||
| 34 | /// - all incoming particles which do not have a |
||
| 35 | /// creation vertex. |
||
| 36 | /// (i.e. it creates copies of all particles which it owns) |
||
| 37 | /// (note - impossible to copy the FULL list of particle pointers |
||
| 38 | /// while having the vertex |
||
| 39 | /// and particles in/out point-back to one another -- unless you |
||
| 40 | /// copy the entire tree -- which we don't want to do) |
||
| 143 | garren | 41 | // |
| 42 | for ( particles_in_const_iterator |
||
| 43 | part1 = invertex.particles_in_const_begin(); |
||
| 44 | part1 != invertex.particles_in_const_end(); ++part1 ) { |
||
| 45 | if ( !(*part1)->production_vertex() ) { |
||
| 46 | GenParticle* pin = new GenParticle(**part1); |
||
| 47 | add_particle_in(pin); |
||
| 48 | } |
||
| 49 | } |
||
| 50 | for ( particles_out_const_iterator |
||
| 51 | part2 = invertex.particles_out_const_begin(); |
||
| 52 | part2 != invertex.particles_out_const_end(); part2++ ) { |
||
| 53 | GenParticle* pin = new GenParticle(**part2); |
||
| 54 | add_particle_out(pin); |
||
| 55 | } |
||
| 56 | suggest_barcode( invertex.barcode() ); |
||
| 57 | // |
||
| 2 | garren | 58 | s_counter++; |
| 59 | } |
||
| 60 | |||
| 61 | GenVertex::~GenVertex() { |
||
| 62 | // |
||
| 63 | // need to delete any particles previously owned by this vertex |
||
| 64 | if ( parent_event() ) parent_event()->remove_barcode(this); |
||
| 65 | delete_adopted_particles(); |
||
| 66 | s_counter--; |
||
| 67 | } |
||
| 68 | |||
| 143 | garren | 69 | void GenVertex::swap( GenVertex & other) |
| 70 | { |
||
| 71 | m_position.swap( other.m_position ); |
||
| 72 | m_particles_in.swap( other.m_particles_in ); |
||
| 73 | m_particles_out.swap( other.m_particles_out ); |
||
| 74 | std::swap( m_id, other.m_id ); |
||
| 75 | m_weights.swap( other.m_weights ); |
||
| 76 | std::swap( m_event, other.m_event ); |
||
| 77 | std::swap( m_barcode, other.m_barcode ); |
||
| 78 | } |
||
| 79 | |||
| 2 | garren | 80 | GenVertex& GenVertex::operator=( const GenVertex& invertex ) { |
| 65 | garren | 81 | /// Shallow: does not copy the FULL list of particle pointers. |
| 82 | /// Creates a copy of - invertex |
||
| 83 | /// - outgoing particles of invertex, but sets the |
||
| 84 | /// decay vertex of these particles to NULL |
||
| 85 | /// - all incoming particles which do not have a |
||
| 86 | /// creation vertex. |
||
| 87 | /// - it does not alter *this's m_event (!) |
||
| 88 | /// (i.e. it creates copies of all particles which it owns) |
||
| 89 | /// (note - impossible to copy the FULL list of particle pointers |
||
| 90 | /// while having the vertex |
||
| 91 | /// and particles in/out point-back to one another -- unless you |
||
| 92 | /// copy the entire tree -- which we don't want to do) |
||
| 93 | /// |
||
| 143 | garren | 94 | |
| 95 | // best practices implementation |
||
| 96 | GenVertex tmp( invertex ); |
||
| 97 | swap( tmp ); |
||
| 2 | garren | 98 | return *this; |
| 99 | } |
||
| 100 | |||
| 101 | bool GenVertex::operator==( const GenVertex& a ) const { |
||
| 65 | garren | 102 | /// Returns true if the positions and the particles in the lists of a |
| 103 | /// and this are identical. Does not compare barcodes. |
||
| 104 | /// Note that it is impossible for two vertices to point to the same |
||
| 105 | /// particle's address, so we need to do more than just compare the |
||
| 106 | /// particle pointers |
||
| 2 | garren | 107 | // |
| 108 | if ( a.position() != this->position() ) return 0; |
||
| 109 | // if the size of the inlist differs, return false. |
||
| 110 | if ( a.particles_in_size() != this->particles_in_size() ) return 0; |
||
| 111 | // if the size of the outlist differs, return false. |
||
| 112 | if ( a.particles_out_size() != this->particles_out_size() ) return 0; |
||
| 113 | // loop over the inlist and ensure particles are identical |
||
| 114 | // (only do this if the lists aren't empty - we already know |
||
| 115 | // if one isn't, then other isn't either!) |
||
| 116 | if ( a.particles_in_const_begin() != a.particles_in_const_end() ) { |
||
| 117 | for ( GenVertex::particles_in_const_iterator |
||
| 118 | ia = a.particles_in_const_begin(), |
||
| 119 | ib = this->particles_in_const_begin(); |
||
| 120 | ia != a.particles_in_const_end(); ia++, ib++ ){ |
||
| 121 | if ( **ia != **ib ) return 0; |
||
| 122 | } |
||
| 123 | } |
||
| 124 | // loop over the outlist and ensure particles are identical |
||
| 125 | // (only do this if the lists aren't empty - we already know |
||
| 126 | // if one isn't, then other isn't either!) |
||
| 127 | if ( a.particles_out_const_begin() != a.particles_out_const_end() ) { |
||
| 128 | for ( GenVertex::particles_out_const_iterator |
||
| 129 | ia = a.particles_out_const_begin(), |
||
| 130 | ib = this->particles_out_const_begin(); |
||
| 131 | ia != a.particles_out_const_end(); ia++, ib++ ){ |
||
| 132 | if ( **ia != **ib ) return 0; |
||
| 133 | } |
||
| 134 | } |
||
| 135 | return 1; |
||
| 136 | } |
||
| 137 | |||
| 138 | bool GenVertex::operator!=( const GenVertex& a ) const { |
||
| 139 | // Returns true if the positions and lists of a and this are not equal. |
||
| 140 | return !( a == *this ); |
||
| 141 | } |
||
| 142 | |||
| 143 | void GenVertex::print( std::ostream& ostr ) const { |
||
| 144 | if ( barcode()!=0 ) { |
||
| 43 | garren | 145 | if ( position() != FourVector(0,0,0,0) ) { |
| 25 | garren | 146 | ostr << "Vertex:"; |
| 147 | ostr.width(9); |
||
| 148 | ostr << barcode(); |
||
| 149 | ostr << " ID:"; |
||
| 150 | ostr.width(5); |
||
| 151 | ostr << id(); |
||
| 152 | ostr << " (X,cT)="; |
||
| 153 | ostr.width(9); |
||
| 154 | ostr.precision(2); |
||
| 155 | ostr.setf(std::ios::scientific, std::ios::floatfield); |
||
| 156 | ostr.setf(std::ios_base::showpos); |
||
| 157 | ostr << position().x() << ","; |
||
| 158 | ostr.width(9); |
||
| 159 | ostr.precision(2); |
||
| 160 | ostr << position().y() << ","; |
||
| 161 | ostr.width(9); |
||
| 162 | ostr.precision(2); |
||
| 163 | ostr << position().z() << ","; |
||
| 164 | ostr.width(9); |
||
| 165 | ostr.precision(2); |
||
| 166 | ostr << position().t(); |
||
| 167 | ostr.setf(std::ios::fmtflags(0), std::ios::floatfield); |
||
| 168 | ostr.unsetf(std::ios_base::showpos); |
||
| 169 | ostr << std::endl; |
||
| 2 | garren | 170 | } else { |
| 25 | garren | 171 | ostr << "GenVertex:"; |
| 172 | ostr.width(9); |
||
| 173 | ostr << barcode(); |
||
| 174 | ostr << " ID:"; |
||
| 175 | ostr.width(5); |
||
| 176 | ostr << id(); |
||
| 177 | ostr << " (X,cT):0"; |
||
| 178 | ostr << std::endl; |
||
| 2 | garren | 179 | } |
| 180 | } else { |
||
| 181 | // If the vertex doesn't have a unique barcode assigned, then |
||
| 182 | // we print its memory address instead... so that the |
||
| 183 | // print out gives us a unique tag for the particle. |
||
| 43 | garren | 184 | if ( position() != FourVector(0,0,0,0) ) { |
| 25 | garren | 185 | ostr << "Vertex:"; |
| 186 | ostr.width(9); |
||
| 187 | ostr << (void*)this; |
||
| 188 | ostr << " ID:"; |
||
| 189 | ostr.width(5); |
||
| 190 | ostr << id(); |
||
| 191 | ostr << " (X,cT)="; |
||
| 192 | ostr.width(9); |
||
| 193 | ostr.precision(2); |
||
| 194 | ostr.setf(std::ios::scientific, std::ios::floatfield); |
||
| 195 | ostr.setf(std::ios_base::showpos); |
||
| 196 | ostr << position().x(); |
||
| 197 | ostr.width(9); |
||
| 198 | ostr.precision(2); |
||
| 199 | ostr << position().y(); |
||
| 200 | ostr.width(9); |
||
| 201 | ostr.precision(2); |
||
| 202 | ostr << position().z(); |
||
| 203 | ostr.width(9); |
||
| 204 | ostr.precision(2); |
||
| 205 | ostr << position().t(); |
||
| 206 | ostr.setf(std::ios::fmtflags(0), std::ios::floatfield); |
||
| 207 | ostr.unsetf(std::ios_base::showpos); |
||
| 208 | ostr << std::endl; |
||
| 2 | garren | 209 | } else { |
| 25 | garren | 210 | ostr << "GenVertex:"; |
| 211 | ostr.width(9); |
||
| 212 | ostr << (void*)this; |
||
| 213 | ostr << " ID:"; |
||
| 214 | ostr.width(5); |
||
| 215 | ostr << id(); |
||
| 216 | ostr << " (X,cT):0"; |
||
| 217 | ostr << std::endl; |
||
| 2 | garren | 218 | } |
| 219 | } |
||
| 220 | |||
| 221 | // print the weights if there are any |
||
| 222 | if ( ! weights().empty() ) { |
||
| 223 | ostr << " Wgts(" << weights().size() << ")="; |
||
| 224 | for ( WeightContainer::const_iterator wgt = weights().begin(); |
||
| 225 | wgt != weights().end(); wgt++ ) { ostr << *wgt << " "; } |
||
| 226 | ostr << std::endl; |
||
| 227 | } |
||
| 228 | // print out all the incoming, then outgoing particles |
||
| 229 | for ( particles_in_const_iterator part1 = particles_in_const_begin(); |
||
| 230 | part1 != particles_in_const_end(); part1++ ) { |
||
| 231 | if ( part1 == particles_in_const_begin() ) { |
||
| 25 | garren | 232 | ostr << " I:"; |
| 233 | ostr.width(2); |
||
| 234 | ostr << m_particles_in.size(); |
||
| 2 | garren | 235 | } else { ostr << " "; } |
| 236 | //(*part1)->print( ostr ); //uncomment for long debugging printout |
||
| 237 | ostr << **part1 << std::endl; |
||
| 238 | } |
||
| 239 | for ( particles_out_const_iterator part2 = particles_out_const_begin(); |
||
| 240 | part2 != particles_out_const_end(); part2++ ) { |
||
| 241 | if ( part2 == particles_out_const_begin() ) { |
||
| 25 | garren | 242 | ostr << " O:"; |
| 243 | ostr.width(2); |
||
| 244 | ostr << m_particles_out.size(); |
||
| 2 | garren | 245 | } else { ostr << " "; } |
| 246 | //(*part2)->print( ostr ); // uncomment for long debugging printout |
||
| 247 | ostr << **part2 << std::endl; |
||
| 248 | } |
||
| 249 | } |
||
| 250 | |||
| 251 | double GenVertex::check_momentum_conservation() const { |
||
| 65 | garren | 252 | /// finds the difference between the total momentum out and the total |
| 253 | /// momentum in vectors, and returns the magnitude of this vector |
||
| 254 | /// i.e. returns | \vec{p_in} - \vec{p_out} | |
||
| 2 | garren | 255 | double sumpx = 0, sumpy = 0, sumpz = 0; |
| 256 | for ( particles_in_const_iterator part1 = particles_in_const_begin(); |
||
| 257 | part1 != particles_in_const_end(); part1++ ) { |
||
| 258 | sumpx += (*part1)->momentum().px(); |
||
| 259 | sumpy += (*part1)->momentum().py(); |
||
| 260 | sumpz += (*part1)->momentum().pz(); |
||
| 261 | } |
||
| 262 | for ( particles_out_const_iterator part2 = particles_out_const_begin(); |
||
| 263 | part2 != particles_out_const_end(); part2++ ) { |
||
| 264 | sumpx -= (*part2)->momentum().px(); |
||
| 265 | sumpy -= (*part2)->momentum().py(); |
||
| 266 | sumpz -= (*part2)->momentum().pz(); |
||
| 267 | } |
||
| 268 | return sqrt( sumpx*sumpx + sumpy*sumpy + sumpz*sumpz ); |
||
| 269 | } |
||
| 270 | |||
| 271 | void GenVertex::add_particle_in( GenParticle* inparticle ) { |
||
| 272 | if ( !inparticle ) return; |
||
| 273 | // if inparticle previously had a decay vertex, remove it from that |
||
| 274 | // vertex's list |
||
| 275 | if ( inparticle->end_vertex() ) { |
||
| 276 | inparticle->end_vertex()->m_particles_in.erase( inparticle ); |
||
| 277 | } |
||
| 278 | m_particles_in.insert( inparticle ); |
||
| 279 | inparticle->set_end_vertex_( this ); |
||
| 280 | } |
||
| 281 | |||
| 282 | void GenVertex::add_particle_out( GenParticle* outparticle ) { |
||
| 283 | if ( !outparticle ) return; |
||
| 284 | // if outparticle previously had a production vertex, |
||
| 285 | // remove it from that vertex's list |
||
| 286 | if ( outparticle->production_vertex() ) { |
||
| 287 | outparticle->production_vertex()->m_particles_out.erase( |
||
| 288 | outparticle ); |
||
| 289 | } |
||
| 290 | m_particles_out.insert( outparticle ); |
||
| 291 | outparticle->set_production_vertex_( this ); |
||
| 292 | } |
||
| 293 | |||
| 294 | GenParticle* GenVertex::remove_particle( GenParticle* particle ) { |
||
| 65 | garren | 295 | /// this finds *particle in the in and/or out list and removes it from |
| 296 | /// these lists ... it DOES NOT DELETE THE PARTICLE or its relations. |
||
| 297 | /// you could delete the particle too as follows: |
||
| 298 | /// delete vtx->remove_particle( particle ); |
||
| 299 | /// or if the particle has an end vertex, you could: |
||
| 300 | /// delete vtx->remove_particle( particle )->end_vertex(); |
||
| 301 | /// which would delete the particle's end vertex, and thus would |
||
| 302 | /// also delete the particle, since the particle would be |
||
| 303 | /// owned by the end vertex. |
||
| 2 | garren | 304 | if ( !particle ) return 0; |
| 305 | if ( particle->end_vertex() == this ) { |
||
| 306 | particle->set_end_vertex_( 0 ); |
||
| 307 | m_particles_in.erase(particle); |
||
| 308 | } |
||
| 309 | if ( particle->production_vertex() == this ) { |
||
| 310 | particle->set_production_vertex_(0); |
||
| 311 | m_particles_out.erase(particle); |
||
| 312 | } |
||
| 313 | return particle; |
||
| 314 | } |
||
| 315 | |||
| 316 | void GenVertex::delete_adopted_particles() { |
||
| 65 | garren | 317 | /// deletes all particles which this vertex owns |
| 318 | /// to be used by the vertex destructor and operator= |
||
| 2 | garren | 319 | // |
| 320 | if ( m_particles_out.empty() && m_particles_in.empty() ) return; |
||
| 321 | // 1. delete all outgoing particles which don't have decay vertices. |
||
| 322 | // those that do become the responsibility of the decay vertex |
||
| 323 | // and have their productionvertex pointer set to NULL |
||
| 95 | garren | 324 | for ( std::set<GenParticle*,GenParticleComparison>::iterator part1 = m_particles_out.begin(); |
| 2 | garren | 325 | part1 != m_particles_out.end(); ) { |
| 326 | if ( !(*part1)->end_vertex() ) { |
||
| 327 | delete *(part1++); |
||
| 328 | } else { |
||
| 329 | (*part1)->set_production_vertex_(0); |
||
| 330 | ++part1; |
||
| 331 | } |
||
| 332 | } |
||
| 333 | m_particles_out.clear(); |
||
| 334 | // |
||
| 335 | // 2. delete all incoming particles which don't have production |
||
| 336 | // vertices. those that do become the responsibility of the |
||
| 337 | // production vertex and have their decayvertex pointer set to NULL |
||
| 95 | garren | 338 | for ( std::set<GenParticle*,GenParticleComparison>::iterator part2 = m_particles_in.begin(); |
| 2 | garren | 339 | part2 != m_particles_in.end(); ) { |
| 340 | if ( !(*part2)->production_vertex() ) { |
||
| 341 | delete *(part2++); |
||
| 342 | } else { |
||
| 343 | (*part2)->set_end_vertex_(0); |
||
| 344 | ++part2; |
||
| 345 | } |
||
| 346 | } |
||
| 347 | m_particles_in.clear(); |
||
| 348 | } |
||
| 349 | |||
| 350 | bool GenVertex::suggest_barcode( int the_bar_code ) |
||
| 351 | { |
||
| 65 | garren | 352 | /// allows a barcode to be suggested for this vertex. |
| 353 | /// In general it is better to let the event pick the barcode for |
||
| 354 | /// you, which is automatic. |
||
| 355 | /// Returns TRUE if the suggested barcode has been accepted (i.e. the |
||
| 356 | /// suggested barcode has not already been used in the event, |
||
| 357 | /// and so it was used). |
||
| 358 | /// Returns FALSE if the suggested barcode was rejected, or if the |
||
| 359 | /// vertex is not yet part of an event, such that it is not yet |
||
| 360 | /// possible to know if the suggested barcode will be accepted). |
||
| 2 | garren | 361 | if ( the_bar_code >0 ) { |
| 362 | std::cerr << "GenVertex::suggest_barcode WARNING, vertex bar codes" |
||
| 363 | << "\n MUST be negative integers. Positive integers " |
||
| 364 | << "\n are reserved for particles only. Your suggestion " |
||
| 365 | << "\n has been rejected." << std::endl; |
||
| 366 | return false; |
||
| 367 | } |
||
| 368 | bool success = false; |
||
| 369 | if ( parent_event() ) { |
||
| 370 | success = parent_event()->set_barcode( this, the_bar_code ); |
||
| 371 | } else { set_barcode_( the_bar_code ); } |
||
| 372 | return success; |
||
| 373 | } |
||
| 374 | |||
| 375 | void GenVertex::set_parent_event_( GenEvent* new_evt ) |
||
| 376 | { |
||
| 377 | GenEvent* orig_evt = m_event; |
||
| 378 | m_event = new_evt; |
||
| 379 | // |
||
| 380 | // every time a vertex's parent event changes, the map of barcodes |
||
| 381 | // in the new and old parent event needs to be modified to |
||
| 382 | // reflect this |
||
| 383 | if ( orig_evt != new_evt ) { |
||
| 384 | if (new_evt) new_evt->set_barcode( this, barcode() ); |
||
| 385 | if (orig_evt) orig_evt->remove_barcode( this ); |
||
| 386 | // we also need to loop over all the particles which are owned by |
||
| 387 | // this vertex, and remove their barcodes from the old event. |
||
| 388 | for ( particles_in_const_iterator part1=particles_in_const_begin(); |
||
| 389 | part1 != particles_in_const_end(); part1++ ) { |
||
| 390 | if ( !(*part1)->production_vertex() ) { |
||
| 391 | if ( orig_evt ) orig_evt->remove_barcode( *part1 ); |
||
| 392 | if ( new_evt ) new_evt->set_barcode( *part1, |
||
| 393 | (*part1)->barcode() ); |
||
| 394 | } |
||
| 395 | } |
||
| 396 | for ( particles_out_const_iterator |
||
| 397 | part2 = particles_out_const_begin(); |
||
| 398 | part2 != particles_out_const_end(); part2++ ) { |
||
| 399 | if ( orig_evt ) orig_evt->remove_barcode( *part2 ); |
||
| 400 | if ( new_evt ) new_evt->set_barcode( *part2, |
||
| 401 | (*part2)->barcode() ); |
||
| 402 | } |
||
| 403 | } |
||
| 404 | } |
||
| 405 | |||
| 406 | ///////////// |
||
| 407 | // Static // |
||
| 408 | ///////////// |
||
| 409 | unsigned int GenVertex::counter() { return s_counter; } |
||
| 410 | unsigned int GenVertex::s_counter = 0; |
||
| 411 | |||
| 412 | ///////////// |
||
| 413 | // Friends // |
||
| 414 | ///////////// |
||
| 415 | |||
| 69 | garren | 416 | /// send vertex information to ostr for printing |
| 2 | garren | 417 | std::ostream& operator<<( std::ostream& ostr, const GenVertex& vtx ) { |
| 418 | if ( vtx.barcode()!=0 ) ostr << "BarCode " << vtx.barcode(); |
||
| 419 | else ostr << "Address " << &vtx; |
||
| 420 | ostr << " (X,cT)="; |
||
| 43 | garren | 421 | if ( vtx.position() != FourVector(0,0,0,0)) { |
| 2 | garren | 422 | ostr << vtx.position().x() << "," |
| 423 | << vtx.position().y() << "," |
||
| 424 | << vtx.position().z() << "," |
||
| 425 | << vtx.position().t(); |
||
| 426 | } else { ostr << 0; } |
||
| 427 | ostr << " #in:" << vtx.particles_in_size() |
||
| 428 | << " #out:" << vtx.particles_out_size(); |
||
| 429 | return ostr; |
||
| 430 | } |
||
| 431 | |||
| 432 | ///////////////////////////// |
||
| 433 | // edge_iterator // (protected - for internal use only) |
||
| 434 | ///////////////////////////// |
||
| 435 | // If the user wants the functionality of the edge_iterator, he should |
||
| 436 | // use particle_iterator with IteratorRange = family, parents, or children |
||
| 437 | // |
||
| 438 | |||
| 439 | GenVertex::edge_iterator::edge_iterator() : m_vertex(0), m_range(family), |
||
| 440 | m_is_inparticle_iter(0), m_is_past_end(1) |
||
| 441 | {} |
||
| 442 | |||
| 443 | GenVertex::edge_iterator::edge_iterator( const GenVertex& vtx, |
||
| 444 | IteratorRange range ) : |
||
| 445 | m_vertex(&vtx), m_range(family) |
||
| 446 | { |
||
| 447 | // Note: (26.1.2000) the original version of edge_iterator inheritted |
||
| 448 | // from set<GenParticle*>::const_iterator() rather than using |
||
| 449 | // composition as it does now. |
||
| 450 | // The inheritted version suffered from a strange bug, which |
||
| 451 | // I have not fully understood --- it only occurred after many |
||
| 452 | // events were processed and only when I called the delete |
||
| 453 | // function on past events. I believe it had something to do with |
||
| 454 | // the past the end values, which are now robustly coded in this |
||
| 455 | // version as boolean members. |
||
| 456 | // |
||
| 457 | // default range is family, only other choices are children/parents |
||
| 458 | // descendants/ancestors not allowed & recasted ot children/parents |
||
| 459 | if ( range == descendants || range == children ) m_range = children; |
||
| 460 | if ( range == ancestors || range == parents ) m_range = parents; |
||
| 461 | // |
||
| 462 | if ( m_vertex->m_particles_in.empty() && |
||
| 463 | m_vertex->m_particles_out.empty() ) { |
||
| 464 | // Case: particles_in and particles_out is empty. |
||
| 465 | m_is_inparticle_iter = 0; |
||
| 466 | m_is_past_end = 1; |
||
| 467 | } else if ( m_range == parents && m_vertex->m_particles_in.empty() ){ |
||
| 468 | // Case: particles in is empty and parents is requested. |
||
| 469 | m_is_inparticle_iter = 1; |
||
| 470 | m_is_past_end = 1; |
||
| 471 | } else if ( m_range == children && m_vertex->m_particles_out.empty() ){ |
||
| 472 | // Case: particles out is empty and children is requested. |
||
| 473 | m_is_inparticle_iter = 0; |
||
| 474 | m_is_past_end = 1; |
||
| 475 | } else if ( m_range == children ) { |
||
| 476 | // Case: particles out is NOT empty, and children is requested |
||
| 477 | m_set_iter = m_vertex->m_particles_out.begin(); |
||
| 478 | m_is_inparticle_iter = 0; |
||
| 479 | m_is_past_end = 0; |
||
| 480 | } else if ( m_range == family && m_vertex->m_particles_in.empty() ) { |
||
| 481 | // Case: particles in is empty, particles out is NOT empty, |
||
| 482 | // and family is requested. Then skip ahead to partilces out. |
||
| 483 | m_set_iter = m_vertex->m_particles_out.begin(); |
||
| 484 | m_is_inparticle_iter = 0; |
||
| 485 | m_is_past_end = 0; |
||
| 486 | } else { |
||
| 487 | // Normal scenario: start with the first incoming particle |
||
| 488 | m_set_iter = m_vertex->m_particles_in.begin(); |
||
| 489 | m_is_inparticle_iter = 1; |
||
| 490 | m_is_past_end = 0; |
||
| 491 | } |
||
| 492 | } |
||
| 493 | |||
| 494 | GenVertex::edge_iterator::edge_iterator( const edge_iterator& p ) { |
||
| 495 | *this = p; |
||
| 496 | } |
||
| 497 | |||
| 498 | GenVertex::edge_iterator::~edge_iterator() {} |
||
| 499 | |||
| 500 | GenVertex::edge_iterator& GenVertex::edge_iterator::operator=( |
||
| 501 | const edge_iterator& p ) { |
||
| 502 | m_vertex = p.m_vertex; |
||
| 503 | m_range = p.m_range; |
||
| 504 | m_set_iter = p.m_set_iter; |
||
| 505 | m_is_inparticle_iter = p.m_is_inparticle_iter; |
||
| 506 | m_is_past_end = p.m_is_past_end; |
||
| 507 | return *this; |
||
| 508 | } |
||
| 509 | |||
| 510 | GenParticle* GenVertex::edge_iterator::operator*(void) const { |
||
| 511 | if ( !m_vertex || m_is_past_end ) return 0; |
||
| 512 | return *m_set_iter; |
||
| 513 | } |
||
| 514 | |||
| 515 | GenVertex::edge_iterator& GenVertex::edge_iterator::operator++(void){ |
||
| 516 | // Pre-fix increment |
||
| 517 | // |
||
| 518 | // increment the set iterator (unless we're past the end value) |
||
| 519 | if ( m_is_past_end ) return *this; |
||
| 520 | ++m_set_iter; |
||
| 521 | // handle cases where m_set_iter points past the end |
||
| 522 | if ( m_range == family && m_is_inparticle_iter && |
||
| 523 | m_set_iter == m_vertex->m_particles_in.end() ) { |
||
| 524 | // at the end on in particle set, and range is family, so move to |
||
| 525 | // out particle set |
||
| 526 | m_set_iter = m_vertex->m_particles_out.begin(); |
||
| 527 | m_is_inparticle_iter = 0; |
||
| 528 | } else if ( m_range == parents && |
||
| 529 | m_set_iter == m_vertex->m_particles_in.end() ) { |
||
| 530 | // at the end on in particle set, and range is parents only, so |
||
| 531 | // move into past the end state |
||
| 532 | m_is_past_end = 1; |
||
| 533 | } |
||
| 534 | // the following is not else if because we might have range=family |
||
| 535 | // with an empty particles_out set. |
||
| 536 | if ( m_set_iter == m_vertex->m_particles_out.end() ) { |
||
| 537 | //whenever out particles end is reached, go into past the end state |
||
| 538 | m_is_past_end = 1; |
||
| 539 | } |
||
| 540 | return *this; |
||
| 541 | } |
||
| 542 | |||
| 543 | GenVertex::edge_iterator GenVertex::edge_iterator::operator++(int){ |
||
| 544 | // Post-fix increment |
||
| 545 | edge_iterator returnvalue = *this; |
||
| 546 | ++*this; |
||
| 547 | return returnvalue; |
||
| 548 | } |
||
| 549 | |||
| 550 | bool GenVertex::edge_iterator::is_parent() const { |
||
| 551 | if ( **this && (**this)->end_vertex() == m_vertex ) return 1; |
||
| 552 | return 0; |
||
| 553 | } |
||
| 554 | |||
| 555 | bool GenVertex::edge_iterator::is_child() const { |
||
| 556 | if ( **this && (**this)->production_vertex() == m_vertex ) return 1; |
||
| 557 | return 0; |
||
| 558 | } |
||
| 559 | |||
| 560 | int GenVertex::edges_size( IteratorRange range ) const { |
||
| 561 | if ( range == children ) return m_particles_out.size(); |
||
| 562 | if ( range == parents ) return m_particles_in.size(); |
||
| 563 | if ( range == family ) return m_particles_out.size() |
||
| 564 | + m_particles_in.size(); |
||
| 565 | return 0; |
||
| 566 | } |
||
| 567 | |||
| 568 | ///////////////////// |
||
| 569 | // vertex_iterator // |
||
| 570 | ///////////////////// |
||
| 571 | |||
| 572 | GenVertex::vertex_iterator::vertex_iterator() |
||
| 84 | garren | 573 | : m_vertex(0), m_range(), m_visited_vertices(0), m_it_owns_set(0), |
| 574 | m_recursive_iterator(0) |
||
| 2 | garren | 575 | {} |
| 576 | |||
| 577 | GenVertex::vertex_iterator::vertex_iterator( GenVertex& vtx_root, |
||
| 578 | IteratorRange range ) |
||
| 579 | : m_vertex(&vtx_root), m_range(range) |
||
| 580 | { |
||
| 581 | // standard public constructor |
||
| 582 | // |
||
| 583 | m_visited_vertices = new std::set<const GenVertex*>; |
||
| 584 | m_it_owns_set = 1; |
||
| 585 | m_visited_vertices->insert( m_vertex ); |
||
| 586 | m_recursive_iterator = 0; |
||
| 587 | m_edge = m_vertex->edges_begin( m_range ); |
||
| 588 | // advance to the first good return value |
||
| 589 | if ( !follow_edge_() && |
||
| 590 | m_edge != m_vertex->edges_end( m_range )) ++*this; |
||
| 591 | } |
||
| 592 | |||
| 593 | GenVertex::vertex_iterator::vertex_iterator( GenVertex& vtx_root, |
||
| 594 | IteratorRange range, std::set<const GenVertex*>& visited_vertices ) : |
||
| 595 | m_vertex(&vtx_root), m_range(range), |
||
| 596 | m_visited_vertices(&visited_vertices), m_it_owns_set(0), |
||
| 597 | m_recursive_iterator(0) |
||
| 598 | { |
||
| 599 | // This constuctor is only to be called internally by this class |
||
| 600 | // for use with its recursive member pointer (m_recursive_iterator). |
||
| 601 | // Note: we do not need to insert m_vertex_root in the vertex - that is |
||
| 602 | // the responsibility of this iterator's mother, which is normally |
||
| 603 | // done just before calling this protected constructor. |
||
| 604 | m_edge = m_vertex->edges_begin( m_range ); |
||
| 605 | // advance to the first good return value |
||
| 606 | if ( !follow_edge_() && |
||
| 607 | m_edge != m_vertex->edges_end( m_range )) ++*this; |
||
| 608 | } |
||
| 609 | |||
| 610 | GenVertex::vertex_iterator::vertex_iterator( const vertex_iterator& v_iter) |
||
| 611 | : m_vertex(0), m_visited_vertices(0), m_it_owns_set(0), |
||
| 612 | m_recursive_iterator(0) |
||
| 613 | { |
||
| 614 | *this = v_iter; |
||
| 615 | } |
||
| 616 | |||
| 617 | GenVertex::vertex_iterator::~vertex_iterator() { |
||
| 618 | if ( m_recursive_iterator ) delete m_recursive_iterator; |
||
| 619 | if ( m_it_owns_set ) delete m_visited_vertices; |
||
| 620 | } |
||
| 621 | |||
| 622 | GenVertex::vertex_iterator& GenVertex::vertex_iterator::operator=( |
||
| 623 | const vertex_iterator& v_iter ) |
||
| 624 | { |
||
| 625 | // Note: when copying a vertex_iterator that is NOT the owner |
||
| 626 | // of its set container, the pointer to the set is copied. Beware! |
||
| 627 | // (see copy_with_own_set() if you want a different set pointed to) |
||
| 628 | // In practise the user never needs to worry |
||
| 629 | // since such iterators are only intended to be used internally. |
||
| 630 | // |
||
| 631 | // destruct data member pointers |
||
| 632 | if ( m_recursive_iterator ) delete m_recursive_iterator; |
||
| 633 | m_recursive_iterator = 0; |
||
| 634 | if ( m_it_owns_set ) delete m_visited_vertices; |
||
| 635 | m_visited_vertices = 0; |
||
| 636 | m_it_owns_set = 0; |
||
| 637 | // copy the target vertex_iterator to this iterator |
||
| 638 | m_vertex = v_iter.m_vertex; |
||
| 639 | m_range = v_iter.m_range; |
||
| 640 | if ( v_iter.m_it_owns_set ) { |
||
| 641 | // i.e. this vertex will own its set if v_iter points to any |
||
| 642 | // vertex set regardless of whether v_iter owns the set or not! |
||
| 643 | m_visited_vertices = |
||
| 644 | new std::set<const GenVertex*>(*v_iter.m_visited_vertices); |
||
| 645 | m_it_owns_set = 1; |
||
| 646 | } else { |
||
| 647 | m_visited_vertices = v_iter.m_visited_vertices; |
||
| 648 | m_it_owns_set = 0; |
||
| 649 | } |
||
| 650 | // |
||
| 651 | // Note: m_vertex_root is already included in the set of |
||
| 652 | // tv_iter.m_visited_vertices, we do not need to insert it. |
||
| 653 | // |
||
| 654 | m_edge = v_iter.m_edge; |
||
| 655 | copy_recursive_iterator_( v_iter.m_recursive_iterator ); |
||
| 656 | return *this; |
||
| 657 | } |
||
| 658 | |||
| 659 | GenVertex* GenVertex::vertex_iterator::operator*(void) const { |
||
| 660 | // de-reference operator |
||
| 661 | // |
||
| 662 | // if this iterator has an iterator_node, then we return the iterator |
||
| 663 | // node. |
||
| 664 | if ( m_recursive_iterator ) return **m_recursive_iterator; |
||
| 665 | // |
||
| 666 | // an iterator can only return its m_vertex -- any other vertex |
||
| 667 | // is returned by means of a recursive iterator_node |
||
| 668 | // (so this is the only place in the iterator code that a vertex |
||
| 669 | // is returned!) |
||
| 670 | if ( m_vertex ) return m_vertex; |
||
| 671 | return 0; |
||
| 672 | } |
||
| 673 | |||
| 674 | GenVertex::vertex_iterator& GenVertex::vertex_iterator::operator++(void) { |
||
| 675 | // Pre-fix incremental operator |
||
| 676 | // |
||
| 677 | // check for "past the end condition" denoted by m_vertex=0 |
||
| 678 | if ( !m_vertex ) return *this; |
||
| 679 | // if at the last edge, move into the "past the end condition" |
||
| 680 | if ( m_edge == m_vertex->edges_end( m_range ) ) { |
||
| 681 | m_vertex = 0; |
||
| 682 | return *this; |
||
| 683 | } |
||
| 684 | // check to see if we need to create a new recursive iterator by |
||
| 685 | // following the current edge only if a recursive iterator doesn't |
||
| 686 | // already exist. If a new recursive_iterator is created, return it. |
||
| 687 | if ( follow_edge_() ) { |
||
| 688 | return *this; |
||
| 689 | } |
||
| 690 | // |
||
| 691 | // if a recursive iterator already exists, increment it, and return its |
||
| 692 | // value (unless the recursive iterator has null root_vertex [its |
||
| 693 | // root vertex is set to null if it has already returned its root] |
||
| 694 | // - in which case we delete it) |
||
| 695 | // and return the vertex pointed to by the edge. |
||
| 696 | if ( m_recursive_iterator ) { |
||
| 697 | ++(*m_recursive_iterator); |
||
| 698 | if ( **m_recursive_iterator ) { |
||
| 699 | return *this; |
||
| 700 | } else { |
||
| 701 | delete m_recursive_iterator; |
||
| 702 | m_recursive_iterator = 0; |
||
| 703 | } |
||
| 704 | } |
||
| 705 | // |
||
| 706 | // increment to the next particle edge |
||
| 707 | ++m_edge; |
||
| 708 | // if m_edge is at the end, then we have incremented through all |
||
| 709 | // edges, and it is time to return m_vertex, which we accomplish |
||
| 710 | // by returning *this |
||
| 711 | if ( m_edge == m_vertex->edges_end( m_range ) ) return *this; |
||
| 712 | // otherwise we follow the current edge by recursively ++ing. |
||
| 713 | return ++(*this); |
||
| 714 | } |
||
| 715 | |||
| 716 | GenVertex::vertex_iterator GenVertex::vertex_iterator::operator++(int) { |
||
| 717 | // Post-fix increment |
||
| 718 | vertex_iterator returnvalue(*this); |
||
| 719 | ++(*this); |
||
| 720 | return returnvalue; |
||
| 721 | } |
||
| 722 | |||
| 723 | void GenVertex::vertex_iterator::copy_with_own_set( |
||
| 724 | const vertex_iterator& v_iter, |
||
| 725 | std::set<const GenVertex*>& visited_vertices ) { |
||
| 65 | garren | 726 | /// intended for internal use only. (use with care!) |
| 727 | /// this is the same as the operator= method, but it allows the |
||
| 728 | /// user to specify which set container m_visited_vertices points to. |
||
| 729 | /// in all cases, this vertex will NOT own its set. |
||
| 2 | garren | 730 | // |
| 731 | // destruct data member pointers |
||
| 732 | if ( m_recursive_iterator ) delete m_recursive_iterator; |
||
| 733 | m_recursive_iterator = 0; |
||
| 734 | if ( m_it_owns_set ) delete m_visited_vertices; |
||
| 735 | m_visited_vertices = 0; |
||
| 736 | m_it_owns_set = 0; |
||
| 737 | // copy the target vertex_iterator to this iterator |
||
| 738 | m_vertex = v_iter.m_vertex; |
||
| 739 | m_range = v_iter.m_range; |
||
| 740 | m_visited_vertices = &visited_vertices; |
||
| 741 | m_it_owns_set = 0; |
||
| 742 | m_edge = v_iter.m_edge; |
||
| 743 | copy_recursive_iterator_( v_iter.m_recursive_iterator ); |
||
| 744 | } |
||
| 745 | |||
| 746 | GenVertex* GenVertex::vertex_iterator::follow_edge_() { |
||
| 747 | // follows the edge pointed to by m_edge by creating a |
||
| 748 | // recursive iterator for it. |
||
| 749 | // |
||
| 750 | // if a m_recursive_iterator already exists, |
||
| 751 | // this routine has nothing to do, |
||
| 752 | // if there's no m_vertex, there's no point following anything, |
||
| 753 | // also there's no point trying to follow a null edge. |
||
| 754 | if ( m_recursive_iterator || !m_vertex || !*m_edge ) return 0; |
||
| 755 | // |
||
| 756 | // if the range is parents, children, or family (i.e. <= family) |
||
| 757 | // then only the iterator which owns the set is allowed to create |
||
| 758 | // recursive iterators (i.e. recursivity is only allowed to go one |
||
| 759 | // layer deep) |
||
| 760 | if ( m_range <= family && m_it_owns_set == 0 ) return 0; |
||
| 761 | // |
||
| 762 | // M.Dobbs 2001-07-16 |
||
| 763 | // Take care of the very special-rare case where a particle might |
||
| 764 | // point to the same vertex for both production and end |
||
| 765 | if ( (*m_edge)->production_vertex() == |
||
| 766 | (*m_edge)->end_vertex() ) return 0; |
||
| 767 | // |
||
| 768 | // figure out which vertex m_edge is pointing to |
||
| 769 | GenVertex* vtx = ( m_edge.is_parent() ? |
||
| 770 | (*m_edge)->production_vertex() : |
||
| 771 | (*m_edge)->end_vertex() ); |
||
| 772 | // if the pointed to vertex doesn't exist or has already been visited, |
||
| 773 | // then return null |
||
| 774 | if ( !vtx || !(m_visited_vertices->insert(vtx).second) ) return 0; |
||
| 775 | // follow that edge by creating a recursive iterator |
||
| 776 | m_recursive_iterator = new vertex_iterator( *vtx, m_range, |
||
| 777 | *m_visited_vertices); |
||
| 778 | // and return the vertex pointed to by m_recursive_iterator |
||
| 779 | return **m_recursive_iterator; |
||
| 780 | } |
||
| 781 | |||
| 782 | void GenVertex::vertex_iterator::copy_recursive_iterator_( |
||
| 783 | const vertex_iterator* recursive_v_iter ) { |
||
| 784 | // to properly copy the recursive iterator, we need to ensure |
||
| 785 | // the proper set container is transfered ... then do this |
||
| 786 | // operation .... you guessed it .... recursively! |
||
| 787 | // |
||
| 788 | if ( !recursive_v_iter ) return; |
||
| 789 | m_recursive_iterator = new vertex_iterator(); |
||
| 790 | m_recursive_iterator->m_vertex = recursive_v_iter->m_vertex; |
||
| 791 | m_recursive_iterator->m_range = recursive_v_iter->m_range; |
||
| 792 | m_recursive_iterator->m_visited_vertices = m_visited_vertices; |
||
| 793 | m_recursive_iterator->m_it_owns_set = 0; |
||
| 794 | m_recursive_iterator->m_edge = recursive_v_iter->m_edge; |
||
| 795 | m_recursive_iterator->copy_recursive_iterator_( |
||
| 796 | recursive_v_iter->m_recursive_iterator ); |
||
| 797 | } |
||
| 798 | |||
| 799 | /////////////////////////////// |
||
| 800 | // particle_iterator // |
||
| 801 | /////////////////////////////// |
||
| 802 | |||
| 803 | GenVertex::particle_iterator::particle_iterator() {} |
||
| 804 | |||
| 805 | GenVertex::particle_iterator::particle_iterator( GenVertex& vertex_root, |
||
| 806 | IteratorRange range ) { |
||
| 807 | // General Purpose Constructor |
||
| 808 | // |
||
| 809 | if ( range <= family ) { |
||
| 810 | m_edge = GenVertex::edge_iterator( vertex_root, range ); |
||
| 811 | } else { |
||
| 812 | m_vertex_iterator = GenVertex::vertex_iterator(vertex_root, range); |
||
| 813 | m_edge = GenVertex::edge_iterator( **m_vertex_iterator, |
||
| 814 | m_vertex_iterator.range() ); |
||
| 815 | } |
||
| 816 | advance_to_first_(); |
||
| 817 | } |
||
| 818 | |||
| 819 | GenVertex::particle_iterator::particle_iterator( |
||
| 820 | const particle_iterator& p_iter ){ |
||
| 821 | *this = p_iter; |
||
| 822 | } |
||
| 823 | |||
| 824 | GenVertex::particle_iterator::~particle_iterator() {} |
||
| 825 | |||
| 826 | GenVertex::particle_iterator& |
||
| 827 | GenVertex::particle_iterator::operator=( const particle_iterator& p_iter ) |
||
| 828 | { |
||
| 829 | m_vertex_iterator = p_iter.m_vertex_iterator; |
||
| 830 | m_edge = p_iter.m_edge; |
||
| 831 | return *this; |
||
| 832 | } |
||
| 833 | |||
| 834 | GenParticle* GenVertex::particle_iterator::operator*(void) const { |
||
| 835 | return *m_edge; |
||
| 836 | } |
||
| 837 | |||
| 838 | GenVertex::particle_iterator& |
||
| 839 | GenVertex::particle_iterator::operator++(void) { |
||
| 840 | //Pre-fix increment |
||
| 841 | // |
||
| 842 | if ( !*m_edge && !*m_vertex_iterator ) { |
||
| 843 | // past the end condition: do nothing |
||
| 844 | return *this; |
||
| 845 | } else if ( !*m_edge && *m_vertex_iterator ) { |
||
| 846 | // past end of edge, but still have more vertices to visit |
||
| 847 | // increment the vertex, checking that the result is valid |
||
| 848 | if ( !*(++m_vertex_iterator) ) return *this; |
||
| 849 | m_edge = GenVertex::edge_iterator( **m_vertex_iterator, |
||
| 850 | m_vertex_iterator.range() ); |
||
| 851 | } else { |
||
| 852 | ++m_edge; |
||
| 853 | } |
||
| 854 | advance_to_first_(); |
||
| 855 | return *this; |
||
| 856 | } |
||
| 857 | |||
| 858 | GenVertex::particle_iterator GenVertex::particle_iterator::operator++(int){ |
||
| 859 | //Post-fix increment |
||
| 860 | particle_iterator returnvalue(*this); |
||
| 861 | ++(*this); |
||
| 862 | return returnvalue; |
||
| 863 | } |
||
| 864 | |||
| 865 | GenParticle* GenVertex::particle_iterator::advance_to_first_() { |
||
| 65 | garren | 866 | /// if the current edge is not a suitable return value ( because |
| 867 | /// it is a parent of the vertex root that itself belongs to a |
||
| 868 | /// different vertex ) it advances to the first suitable return value |
||
| 2 | garren | 869 | if ( !*m_edge ) return *(++*this); |
| 870 | // if the range is relatives, we need to uniquely assign each particle |
||
| 871 | // to a single vertex so as to guarantee particles are returned |
||
| 872 | // exactly once. |
||
| 873 | if ( m_vertex_iterator.range() == relatives && |
||
| 874 | m_edge.is_parent() && |
||
| 875 | (*m_edge)->production_vertex() ) return *(++*this); |
||
| 876 | return *m_edge; |
||
| 877 | } |
||
| 878 | |||
| 879 | } // HepMC |