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