hepmc - Blame information for rev 65

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