TEvtGen/HepMC/GenVertex.cc

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