SSO Logout

Subversion Repositories hepmc

Rev

Rev 124 | Rev 168 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

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