hepmc - Blame information for rev 179

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