Fork me on GitHub

source: svn/trunk/Utilities/HepMC/src/GenEvent.cc@ 987

Last change on this file since 987 was 571, checked in by cp3-support, 13 years ago

upgrade HepMC to version 2.06.05
File size: 27.6 KB
RevLine 
[349]1//////////////////////////////////////////////////////////////////////////
2// Matt.Dobbs@Cern.CH, September 1999
3// Updated: 7.1.2000 iterators complete and working!
4// Updated: 10.1.2000 GenEvent::vertex, particle iterators are made
5// constant WRT this event ... note that
6// GenVertex::***_iterator is not const, since it must
7// be able to return a mutable pointer to itself.
8// Updated: 08.2.2000 the event now holds a set of all attached vertices
9// rather than just the roots of the graph
10// Event record for MC generators (for use at any stage of generation)
11//////////////////////////////////////////////////////////////////////////
12
13#include <iomanip>
14
15#include "GenEvent.h"
[571]16#include "GenCrossSection.h"
[349]17#include "Version.h"
[571]18#include "StreamHelpers.h"
[349]19
20namespace HepMC {
21
22 GenEvent::GenEvent( int signal_process_id,
23 int event_number,
24 GenVertex* signal_vertex,
25 const WeightContainer& weights,
26 const std::vector<long>& random_states,
27 Units::MomentumUnit mom,
28 Units::LengthUnit len ) :
29 m_signal_process_id(signal_process_id),
30 m_event_number(event_number),
31 m_mpi(-1),
32 m_event_scale(-1),
33 m_alphaQCD(-1),
34 m_alphaQED(-1),
35 m_signal_process_vertex(signal_vertex),
36 m_beam_particle_1(0),
37 m_beam_particle_2(0),
38 m_weights(weights),
39 m_random_states(random_states),
40 m_vertex_barcodes(),
41 m_particle_barcodes(),
[571]42 m_cross_section(0),
[349]43 m_heavy_ion(0),
44 m_pdf_info(0),
45 m_momentum_unit(mom),
46 m_position_unit(len)
47 {
48 /// This constructor only allows null pointers to HeavyIon and PdfInfo
49 ///
50 /// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
51 /// are as suggested in hep-ph/0109068, "Generic Interface..."
52 ///
53 }
54
55 GenEvent::GenEvent( int signal_process_id, int event_number,
56 GenVertex* signal_vertex,
57 const WeightContainer& weights,
58 const std::vector<long>& random_states,
59 const HeavyIon& ion,
60 const PdfInfo& pdf,
61 Units::MomentumUnit mom,
62 Units::LengthUnit len ) :
63 m_signal_process_id(signal_process_id),
64 m_event_number(event_number),
65 m_mpi(-1),
66 m_event_scale(-1),
67 m_alphaQCD(-1),
68 m_alphaQED(-1),
69 m_signal_process_vertex(signal_vertex),
70 m_beam_particle_1(0),
71 m_beam_particle_2(0),
72 m_weights(weights),
73 m_random_states(random_states),
74 m_vertex_barcodes(),
75 m_particle_barcodes(),
[571]76 m_cross_section(0),
[349]77 m_heavy_ion( new HeavyIon(ion) ),
78 m_pdf_info( new PdfInfo(pdf) ),
79 m_momentum_unit(mom),
80 m_position_unit(len)
81 {
82 /// GenEvent makes its own copy of HeavyIon and PdfInfo
83 ///
84 /// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
85 /// are as suggested in hep-ph/0109068, "Generic Interface..."
86 }
87
88 GenEvent::GenEvent( Units::MomentumUnit mom,
89 Units::LengthUnit len,
90 int signal_process_id,
91 int event_number,
92 GenVertex* signal_vertex,
93 const WeightContainer& weights,
94 const std::vector<long>& random_states ) :
95 m_signal_process_id(signal_process_id),
96 m_event_number(event_number),
97 m_mpi(-1),
98 m_event_scale(-1),
99 m_alphaQCD(-1),
100 m_alphaQED(-1),
101 m_signal_process_vertex(signal_vertex),
102 m_beam_particle_1(0),
103 m_beam_particle_2(0),
104 m_weights(weights),
105 m_random_states(random_states),
106 m_vertex_barcodes(),
107 m_particle_barcodes(),
[571]108 m_cross_section(0),
[349]109 m_heavy_ion(0),
110 m_pdf_info(0),
111 m_momentum_unit(mom),
112 m_position_unit(len)
113 {
114 /// constructor requiring units - all else is default
115 /// This constructor only allows null pointers to HeavyIon and PdfInfo
116 ///
117 /// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
118 /// are as suggested in hep-ph/0109068, "Generic Interface..."
119 ///
120 }
121
122 GenEvent::GenEvent( Units::MomentumUnit mom,
123 Units::LengthUnit len,
124 int signal_process_id, int event_number,
125 GenVertex* signal_vertex,
126 const WeightContainer& weights,
127 const std::vector<long>& random_states,
128 const HeavyIon& ion,
129 const PdfInfo& pdf ) :
130 m_signal_process_id(signal_process_id),
131 m_event_number(event_number),
132 m_mpi(-1),
133 m_event_scale(-1),
134 m_alphaQCD(-1),
135 m_alphaQED(-1),
136 m_signal_process_vertex(signal_vertex),
137 m_beam_particle_1(0),
138 m_beam_particle_2(0),
139 m_weights(weights),
140 m_random_states(random_states),
141 m_vertex_barcodes(),
142 m_particle_barcodes(),
[571]143 m_cross_section(0),
[349]144 m_heavy_ion( new HeavyIon(ion) ),
145 m_pdf_info( new PdfInfo(pdf) ),
146 m_momentum_unit(mom),
147 m_position_unit(len)
148 {
149 /// explicit constructor with units first that takes HeavyIon and PdfInfo
150 /// GenEvent makes its own copy of HeavyIon and PdfInfo
151 ///
152 /// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
153 /// are as suggested in hep-ph/0109068, "Generic Interface..."
154 }
155
156 GenEvent::GenEvent( const GenEvent& inevent )
157 : m_signal_process_id ( inevent.signal_process_id() ),
158 m_event_number ( inevent.event_number() ),
159 m_mpi ( inevent.mpi() ),
160 m_event_scale ( inevent.event_scale() ),
161 m_alphaQCD ( inevent.alphaQCD() ),
162 m_alphaQED ( inevent.alphaQED() ),
163 m_signal_process_vertex( /* inevent.m_signal_process_vertex */ ),
164 m_beam_particle_1 ( /* inevent.m_beam_particle_1 */ ),
165 m_beam_particle_2 ( /* inevent.m_beam_particle_2 */ ),
166 m_weights ( /* inevent.m_weights */ ),
167 m_random_states ( /* inevent.m_random_states */ ),
168 m_vertex_barcodes ( /* inevent.m_vertex_barcodes */ ),
169 m_particle_barcodes ( /* inevent.m_particle_barcodes */ ),
[571]170 m_cross_section ( inevent.cross_section() ? new GenCrossSection(*inevent.cross_section()) : 0 ),
[349]171 m_heavy_ion ( inevent.heavy_ion() ? new HeavyIon(*inevent.heavy_ion()) : 0 ),
172 m_pdf_info ( inevent.pdf_info() ? new PdfInfo(*inevent.pdf_info()) : 0 ),
173 m_momentum_unit ( inevent.momentum_unit() ),
174 m_position_unit ( inevent.length_unit() )
175 {
176 /// deep copy - makes a copy of all vertices!
177 //
178
179 // 1. create a NEW copy of all vertices from inevent
180 // taking care to map new vertices onto the vertices being copied
181 // and add these new vertices to this event.
182 // We do not use GenVertex::operator= because that would copy
183 // the attached particles as well.
184 std::map<const GenVertex*,GenVertex*> map_in_to_new;
185 for ( GenEvent::vertex_const_iterator v = inevent.vertices_begin();
186 v != inevent.vertices_end(); ++v ) {
187 GenVertex* newvertex = new GenVertex(
188 (*v)->position(), (*v)->id(), (*v)->weights() );
189 newvertex->suggest_barcode( (*v)->barcode() );
190 map_in_to_new[*v] = newvertex;
191 add_vertex( newvertex );
192 }
193 // 2. copy the signal process vertex info.
194 if ( inevent.signal_process_vertex() ) {
195 set_signal_process_vertex(
196 map_in_to_new[inevent.signal_process_vertex()] );
197 } else set_signal_process_vertex( 0 );
198 //
199 // 3. create a NEW copy of all particles from inevent
200 // taking care to attach them to the appropriate vertex
201 GenParticle* beam1(0);
202 GenParticle* beam2(0);
203 for ( GenEvent::particle_const_iterator p = inevent.particles_begin();
204 p != inevent.particles_end(); ++p )
205 {
206 GenParticle* oldparticle = *p;
207 GenParticle* newparticle = new GenParticle(*oldparticle);
208 if ( oldparticle->end_vertex() ) {
209 map_in_to_new[ oldparticle->end_vertex() ]->
210 add_particle_in(newparticle);
211 }
212 if ( oldparticle->production_vertex() ) {
213 map_in_to_new[ oldparticle->production_vertex() ]->
214 add_particle_out(newparticle);
215 }
216 if ( oldparticle == inevent.beam_particles().first ) beam1 = newparticle;
217 if ( oldparticle == inevent.beam_particles().second ) beam2 = newparticle;
218 }
219 set_beam_particles( beam1, beam2 );
220 //
221 // 4. now that vtx/particles are copied, copy weights and random states
222 set_random_states( inevent.random_states() );
223 weights() = inevent.weights();
224 }
225
226 void GenEvent::swap( GenEvent & other )
227 {
228 // if a container has a swap method, use that for improved performance
229 std::swap(m_signal_process_id , other.m_signal_process_id );
230 std::swap(m_event_number , other.m_event_number );
231 std::swap(m_mpi , other.m_mpi );
232 std::swap(m_event_scale , other.m_event_scale );
233 std::swap(m_alphaQCD , other.m_alphaQCD );
234 std::swap(m_alphaQED , other.m_alphaQED );
235 std::swap(m_signal_process_vertex, other.m_signal_process_vertex);
236 std::swap(m_beam_particle_1 , other.m_beam_particle_1 );
237 std::swap(m_beam_particle_2 , other.m_beam_particle_2 );
238 m_weights.swap( other.m_weights );
239 m_random_states.swap( other.m_random_states );
240 m_vertex_barcodes.swap( other.m_vertex_barcodes );
241 m_particle_barcodes.swap( other.m_particle_barcodes );
[571]242 std::swap(m_cross_section , other.m_cross_section );
[349]243 std::swap(m_heavy_ion , other.m_heavy_ion );
244 std::swap(m_pdf_info , other.m_pdf_info );
245 std::swap(m_momentum_unit , other.m_momentum_unit );
246 std::swap(m_position_unit , other.m_position_unit );
247 // must now adjust GenVertex back pointers
248 for ( GenEvent::vertex_const_iterator vthis = vertices_begin();
249 vthis != vertices_end(); ++vthis ) {
250 (*vthis)->change_parent_event_( this );
251 }
252 for ( GenEvent::vertex_const_iterator voth = other.vertices_begin();
253 voth != other.vertices_end(); ++voth ) {
254 (*voth)->change_parent_event_( &other );
255 }
256 }
257
258 GenEvent::~GenEvent()
259 {
260 /// Deep destructor.
[571]261 /// deletes all vertices/particles in this GenEvent
262 /// deletes the associated HeavyIon and PdfInfo
[349]263 delete_all_vertices();
[571]264 delete m_cross_section;
[349]265 delete m_heavy_ion;
266 delete m_pdf_info;
267 }
268
269 GenEvent& GenEvent::operator=( const GenEvent& inevent )
270 {
271 /// best practices implementation
272 GenEvent tmp( inevent );
273 swap( tmp );
274 return *this;
275 }
276
277 void GenEvent::print( std::ostream& ostr ) const {
278 /// dumps the content of this event to ostr
279 /// to dump to cout use: event.print();
280 /// if you want to write this event to file outfile.txt you could use:
281 /// std::ofstream outfile("outfile.txt"); event.print( outfile );
282 ostr << "________________________________________"
283 << "________________________________________\n";
284 ostr << "GenEvent: #" << event_number()
285 << " ID=" << signal_process_id()
286 << " SignalProcessGenVertex Barcode: "
287 << ( signal_process_vertex() ? signal_process_vertex()->barcode()
288 : 0 )
289 << "\n";
290 write_units( ostr );
[571]291 write_cross_section(ostr);
[349]292 ostr << " Entries this event: " << vertices_size() << " vertices, "
293 << particles_size() << " particles.\n";
294 if( m_beam_particle_1 && m_beam_particle_2 ) {
295 ostr << " Beam Particle barcodes: " << beam_particles().first->barcode() << " "
296 << beam_particles().second->barcode() << " \n";
297 } else {
298 ostr << " Beam Particles are not defined.\n";
299 }
300 // Random State
301 ostr << " RndmState(" << m_random_states.size() << ")=";
302 for ( std::vector<long>::const_iterator rs
303 = m_random_states.begin();
304 rs != m_random_states.end(); ++rs ) { ostr << *rs << " "; }
305 ostr << "\n";
306 // Weights
307 ostr << " Wgts(" << weights().size() << ")=";
[571]308 weights().print(ostr);
[349]309 ostr << " EventScale " << event_scale()
310 << " [energy] \t alphaQCD=" << alphaQCD()
311 << "\t alphaQED=" << alphaQED() << std::endl;
312 // print a legend to describe the particle info
313 ostr << " GenParticle Legend\n";
314 ostr << " Barcode PDG ID "
315 << "( Px, Py, Pz, E )"
316 << " Stat DecayVtx\n";
317 ostr << "________________________________________"
318 << "________________________________________\n";
319 // Print all Vertices
320 for ( GenEvent::vertex_const_iterator vtx = this->vertices_begin();
321 vtx != this->vertices_end(); ++vtx ) {
322 (*vtx)->print(ostr);
323 }
324 ostr << "________________________________________"
325 << "________________________________________" << std::endl;
326 }
327
328 void GenEvent::print_version( std::ostream& ostr ) const {
329 ostr << "---------------------------------------------" << std::endl;
330 writeVersion( ostr );
331 ostr << "---------------------------------------------" << std::endl;
332 }
333
334 bool GenEvent::add_vertex( GenVertex* vtx ) {
335 /// returns true if successful - generally will only return false
336 /// if the inserted vertex is already included in the event.
337 if ( !vtx ) return false;
338 // if vtx previously pointed to another GenEvent, remove it from that
339 // GenEvent's list
340 if ( vtx->parent_event() && vtx->parent_event() != this ) {
341 bool remove_status = vtx->parent_event()->remove_vertex( vtx );
342 if ( !remove_status ) {
343 std::cerr << "GenEvent::add_vertex ERROR "
344 << "GenVertex::parent_event points to \n"
345 << "an event that does not point back to the "
346 << "GenVertex. \n This probably indicates a deeper "
347 << "problem. " << std::endl;
348 }
349 }
350 //
351 // setting the vertex parent also inserts the vertex into this
352 // event
353 vtx->set_parent_event_( this );
354 return ( m_vertex_barcodes.count(vtx->barcode()) ? true : false );
355 }
356
357 bool GenEvent::remove_vertex( GenVertex* vtx ) {
358 /// this removes vtx from the event but does NOT delete it.
359 /// returns True if an entry vtx existed in the table and was erased
360 if ( m_signal_process_vertex == vtx ) m_signal_process_vertex = 0;
361 if ( vtx->parent_event() == this ) vtx->set_parent_event_( 0 );
362 return ( m_vertex_barcodes.count(vtx->barcode()) ? false : true );
363 }
364
365 void GenEvent::clear()
366 {
367 /// remove all information from the event
368 /// deletes all vertices/particles in this evt
369 ///
370 delete_all_vertices();
[571]371 // remove existing objects and set pointers to null
372 delete m_cross_section;
373 m_cross_section = 0;
[349]374 delete m_heavy_ion;
[571]375 m_heavy_ion = 0;
[349]376 delete m_pdf_info;
[571]377 m_pdf_info = 0;
[349]378 m_signal_process_id = 0;
379 m_beam_particle_1 = 0;
380 m_beam_particle_2 = 0;
381 m_event_number = 0;
[571]382 m_mpi = -1;
[349]383 m_event_scale = -1;
384 m_alphaQCD = -1;
385 m_alphaQED = -1;
386 m_weights = std::vector<double>();
387 m_random_states = std::vector<long>();
[571]388 // resetting unit information
389 m_momentum_unit = Units::default_momentum_unit();
390 m_position_unit = Units::default_length_unit();
[349]391 // error check just to be safe
392 if ( m_vertex_barcodes.size() != 0
393 || m_particle_barcodes.size() != 0 ) {
394 std::cerr << "GenEvent::clear() strange result ... \n"
395 << "either the particle and/or the vertex map isn't empty" << std::endl;
396 std::cerr << "Number vtx,particle the event after deleting = "
397 << m_vertex_barcodes.size() << " "
398 << m_particle_barcodes.size() << std::endl;
399 }
400 return;
401 }
402
403 void GenEvent::delete_all_vertices() {
404 /// deletes all vertices in the vertex container
405 /// (i.e. all vertices owned by this event)
406 /// The vertices are the "owners" of the particles, so as we delete
407 /// the vertices, the vertex desctructors are automatically
408 /// deleting their particles.
409
410 // delete each vertex individually (this deletes particles as well)
411 while ( !vertices_empty() ) {
412 GenVertex* vtx = ( m_vertex_barcodes.begin() )->second;
413 m_vertex_barcodes.erase( m_vertex_barcodes.begin() );
414 delete vtx;
415 }
416 //
417 // Error checking:
418 if ( !vertices_empty() || ! particles_empty() ) {
419 std::cerr << "GenEvent::delete_all_vertices strange result ... "
420 << "after deleting all vertices, \nthe particle and "
421 << "vertex maps aren't empty.\n This probably "
422 << "indicates deeper problems or memory leak in the "
423 << "code." << std::endl;
424 std::cerr << "Number vtx,particle the event after deleting = "
425 << m_vertex_barcodes.size() << " "
426 << m_particle_barcodes.size() << std::endl;
427 }
428 }
429
430 bool GenEvent::set_barcode( GenParticle* p, int suggested_barcode )
431 {
432 if ( p->parent_event() != this ) {
433 std::cerr << "GenEvent::set_barcode attempted, but the argument's"
434 << "\n parent_event is not this ... request rejected."
435 << std::endl;
436 return false;
437 }
438 // M.Dobbs Nov 4, 2002
439 // First we must check to see if the particle already has a
440 // barcode which is different from the suggestion. If yes, we
441 // remove it from the particle map.
442 if ( p->barcode() != 0 && p->barcode() != suggested_barcode ) {
443 if ( m_particle_barcodes.count(p->barcode()) &&
444 m_particle_barcodes[p->barcode()] == p ) {
445 m_particle_barcodes.erase( p->barcode() );
446 }
447 // At this point either the particle is NOT in
448 // m_particle_barcodes, or else it is in the map, but
449 // already with the suggested barcode.
450 }
451 //
452 // First case --- a valid barcode has been suggested
453 // (valid barcodes are numbers greater than zero)
454 bool insert_success = true;
455 if ( suggested_barcode > 0 ) {
456 if ( m_particle_barcodes.count(suggested_barcode) ) {
457 // the suggested_barcode is already used.
458 if ( m_particle_barcodes[suggested_barcode] == p ) {
459 // but it was used for this particle ... so everythings ok
460 p->set_barcode_( suggested_barcode );
461 return true;
462 }
463 insert_success = false;
464 suggested_barcode = 0;
465 } else { // suggested barcode is OK, proceed to insert
466 m_particle_barcodes[suggested_barcode] = p;
467 p->set_barcode_( suggested_barcode );
468 return true;
469 }
470 }
471 //
472 // Other possibility -- a valid barcode has not been suggested,
473 // so one is automatically generated
474 if ( suggested_barcode < 0 ) insert_success = false;
475 if ( suggested_barcode <= 0 ) {
476 if ( !m_particle_barcodes.empty() ) {
477 // in this case we find the highest barcode that was used,
478 // and increment it by 1
479 suggested_barcode = m_particle_barcodes.rbegin()->first;
480 ++suggested_barcode;
481 }
482 // For the automatically assigned barcodes, the first one
483 // we use is 10001 ... this is just because when we read
484 // events from HEPEVT, we will suggest their index as the
485 // barcode, and that index has maximum value 10000.
486 // This way we will immediately be able to recognize the hepevt
487 // particles from the auto-assigned ones.
488 if ( suggested_barcode <= 10000 ) suggested_barcode = 10001;
489 }
490 // At this point we should have a valid barcode
491 if ( m_particle_barcodes.count(suggested_barcode) ) {
492 std::cerr << "GenEvent::set_barcode ERROR, this should never "
493 << "happen \n report bug to matt.dobbs@cern.ch"
494 << std::endl;
495 }
496 m_particle_barcodes[suggested_barcode] = p;
497 p->set_barcode_( suggested_barcode );
498 return insert_success;
499 }
500
501 bool GenEvent::set_barcode( GenVertex* v, int suggested_barcode )
502 {
503 if ( v->parent_event() != this ) {
504 std::cerr << "GenEvent::set_barcode attempted, but the argument's"
505 << "\n parent_event is not this ... request rejected."
506 << std::endl;
507 return false;
508 }
509 // M.Dobbs Nov 4, 2002
510 // First we must check to see if the vertex already has a
511 // barcode which is different from the suggestion. If yes, we
512 // remove it from the vertex map.
513 if ( v->barcode() != 0 && v->barcode() != suggested_barcode ) {
514 if ( m_vertex_barcodes.count(v->barcode()) &&
515 m_vertex_barcodes[v->barcode()] == v ) {
516 m_vertex_barcodes.erase( v->barcode() );
517 }
518 // At this point either the vertex is NOT in
519 // m_vertex_barcodes, or else it is in the map, but
520 // already with the suggested barcode.
521 }
522
523 //
524 // First case --- a valid barcode has been suggested
525 // (valid barcodes are numbers greater than zero)
526 bool insert_success = true;
527 if ( suggested_barcode < 0 ) {
528 if ( m_vertex_barcodes.count(suggested_barcode) ) {
529 // the suggested_barcode is already used.
530 if ( m_vertex_barcodes[suggested_barcode] == v ) {
531 // but it was used for this vertex ... so everythings ok
532 v->set_barcode_( suggested_barcode );
533 return true;
534 }
535 insert_success = false;
536 suggested_barcode = 0;
537 } else { // suggested barcode is OK, proceed to insert
538 m_vertex_barcodes[suggested_barcode] = v;
539 v->set_barcode_( suggested_barcode );
540 return true;
541 }
542 }
543 //
544 // Other possibility -- a valid barcode has not been suggested,
545 // so one is automatically generated
546 if ( suggested_barcode > 0 ) insert_success = false;
547 if ( suggested_barcode >= 0 ) {
548 if ( !m_vertex_barcodes.empty() ) {
549 // in this case we find the highest barcode that was used,
550 // and increment it by 1, (vertex barcodes are negative)
551 suggested_barcode = m_vertex_barcodes.rbegin()->first;
552 --suggested_barcode;
553 }
554 if ( suggested_barcode >= 0 ) suggested_barcode = -1;
555 }
556 // At this point we should have a valid barcode
557 if ( m_vertex_barcodes.count(suggested_barcode) ) {
558 std::cerr << "GenEvent::set_barcode ERROR, this should never "
559 << "happen \n report bug to matt.dobbs@cern.ch"
560 << std::endl;
561 }
562 m_vertex_barcodes[suggested_barcode] = v;
563 v->set_barcode_( suggested_barcode );
564 return insert_success;
565 }
566
567 /// test to see if we have two valid beam particles
568 bool GenEvent::valid_beam_particles() const {
569 bool have1 = false;
570 bool have2 = false;
571 // first check that both are defined
572 if(m_beam_particle_1 && m_beam_particle_2) {
573 // now look for a match with the particle "list"
574 for ( particle_const_iterator p = particles_begin();
575 p != particles_end(); ++p ) {
576 if( m_beam_particle_1 == *p ) have1 = true;
577 if( m_beam_particle_2 == *p ) have2 = true;
578 }
579 }
580 if( have1 && have2 ) return true;
581 return false;
582 }
583
584 /// construct the beam particle information using pointers to GenParticle
585 /// returns false if either GenParticle* is null
586 bool GenEvent::set_beam_particles(GenParticle* bp1, GenParticle* bp2) {
587 m_beam_particle_1 = bp1;
588 m_beam_particle_2 = bp2;
589 if( m_beam_particle_1 && m_beam_particle_2 ) return true;
590 return false;
591 }
592
593 /// construct the beam particle information using a std::pair of pointers to GenParticle
594 /// returns false if either GenParticle* is null
595 bool GenEvent::set_beam_particles(std::pair<HepMC::GenParticle*, HepMC::GenParticle*> const & bp) {
596 return set_beam_particles(bp.first,bp.second);
597 }
598
599 void GenEvent::write_units( std::ostream & os ) const {
600 os << " Momenutm units:" << std::setw(8) << name(momentum_unit());
601 os << " Position units:" << std::setw(8) << name(length_unit());
602 os << std::endl;
603 }
604
[571]605 void GenEvent::write_cross_section( std::ostream& os ) const
606 {
607 // write the GenCrossSection information if the cross section was set
608 if( !cross_section() ) return;
609 if( cross_section()->is_set() ) {
610 os << " Cross Section: " << cross_section()->cross_section() ;
611 os << " +/- " << cross_section()->cross_section_error() ;
612 os << std::endl;
613 }
614 }
615
[349]616 bool GenEvent::use_momentum_unit( Units::MomentumUnit newunit ) {
617 // currently not exception-safe.
618 // Easy to fix, though, if needed.
619 if ( m_momentum_unit != newunit ) {
620 const double factor = Units::conversion_factor( m_momentum_unit, newunit );
621 // multiply all momenta by 'factor',
622 // loop is entered only if particle list is not empty
623 for ( GenEvent::particle_iterator p = particles_begin();
624 p != particles_end(); ++p )
625 {
626 (*p)->convert_momentum(factor);
627 }
628 // ...
629 m_momentum_unit = newunit;
630 }
631 return true;
632 }
633
634 bool GenEvent::use_length_unit( Units::LengthUnit newunit ) {
635 // currently not exception-safe.
636 // Easy to fix, though, if needed.
637 if ( m_position_unit != newunit ) {
638 const double factor = Units::conversion_factor( m_position_unit, newunit );
639 // multiply all lengths by 'factor',
640 // loop is entered only if vertex list is not empty
641 for ( GenEvent::vertex_iterator vtx = vertices_begin();
642 vtx != vertices_end(); ++vtx ) {
643 (*vtx)->convert_position(factor);
644 }
645 // ...
646 m_position_unit = newunit;
647 }
648 return true;
649 }
650
651 bool GenEvent::use_momentum_unit( std::string& newunit ) {
652 if ( newunit == "MEV" ) return use_momentum_unit( Units::MEV );
653 else if( newunit == "GEV" ) return use_momentum_unit( Units::GEV );
654 else std::cerr << "GenEvent::use_momentum_unit ERROR: use either MEV or GEV\n";
655 return false;
656 }
657
658 bool GenEvent::use_length_unit( std::string& newunit ) {
659 if ( newunit == "MM" ) return use_length_unit( Units::MM );
660 else if( newunit == "CM" ) return use_length_unit( Units::CM );
661 else std::cerr << "GenEvent::use_length_unit ERROR: use either MEV or GEV\n";
662 return false;
663 }
664
[571]665 bool GenEvent::is_valid() const {
666 /// A GenEvent is presumed valid if it has both associated
667 /// particles and vertices. No other information is checked.
668 if ( vertices_empty() ) return false;
669 if ( particles_empty() ) return false;
670 return true;
671 }
672
673 std::ostream & GenEvent::write_beam_particles(std::ostream & os,
674 std::pair<HepMC::GenParticle *,HepMC::GenParticle *> pr )
675 {
676 GenParticle* p = pr.first;
677 if(!p) {
678 detail::output( os, 0 );
679 } else {
680 detail::output( os, p->barcode() );
681 }
682 p = pr.second;
683 if(!p) {
684 detail::output( os, 0 );
685 } else {
686 detail::output( os, p->barcode() );
687 }
688
689 return os;
690 }
691
692 std::ostream & GenEvent::write_vertex(std::ostream & os, GenVertex const * v)
693 {
694 if ( !v || !os ) {
695 std::cerr << "GenEvent::write_vertex !v||!os, "
696 << "v="<< v << " setting badbit" << std::endl;
697 os.clear(std::ios::badbit);
698 return os;
699 }
700 // First collect info we need
701 // count the number of orphan particles going into v
702 int num_orphans_in = 0;
703 for ( GenVertex::particles_in_const_iterator p1
704 = v->particles_in_const_begin();
705 p1 != v->particles_in_const_end(); ++p1 ) {
706 if ( !(*p1)->production_vertex() ) ++num_orphans_in;
707 }
708 //
709 os << 'V';
710 detail::output( os, v->barcode() ); // v's unique identifier
711 detail::output( os, v->id() );
712 detail::output( os, v->position().x() );
713 detail::output( os, v->position().y() );
714 detail::output( os, v->position().z() );
715 detail::output( os, v->position().t() );
716 detail::output( os, num_orphans_in );
717 detail::output( os, (int)v->particles_out_size() );
718 detail::output( os, (int)v->weights().size() );
719 for ( WeightContainer::const_iterator w = v->weights().begin();
720 w != v->weights().end(); ++w ) {
721 detail::output( os, *w );
722 }
723 detail::output( os,'\n');
724 // incoming particles
725 for ( GenVertex::particles_in_const_iterator p2
726 = v->particles_in_const_begin();
727 p2 != v->particles_in_const_end(); ++p2 ) {
728 if ( !(*p2)->production_vertex() ) {
729 write_particle( os, *p2 );
730 }
731 }
732 // outgoing particles
733 for ( GenVertex::particles_out_const_iterator p3
734 = v->particles_out_const_begin();
735 p3 != v->particles_out_const_end(); ++p3 ) {
736 write_particle( os, *p3 );
737 }
738 return os;
739 }
740
741 std::ostream & GenEvent::write_particle( std::ostream & os, GenParticle const * p )
742 {
743 if ( !p || !os ) {
744 std::cerr << "GenEvent::write_particle !p||!os, "
745 << "p="<< p << " setting badbit" << std::endl;
746 os.clear(std::ios::badbit);
747 return os;
748 }
749 os << 'P';
750 detail::output( os, p->barcode() );
751 detail::output( os, p->pdg_id() );
752 detail::output( os, p->momentum().px() );
753 detail::output( os, p->momentum().py() );
754 detail::output( os, p->momentum().pz() );
755 detail::output( os, p->momentum().e() );
756 detail::output( os, p->generated_mass() );
757 detail::output( os, p->status() );
758 detail::output( os, p->polarization().theta() );
759 detail::output( os, p->polarization().phi() );
760 // since end_vertex is oftentimes null, this CREATES a null vertex
761 // in the map
762 detail::output( os, ( p->end_vertex() ? p->end_vertex()->barcode() : 0 ) );
763 os << ' ' << p->flow() << "\n";
764
765 return os;
766 }
767
[349]768} // HepMC
Note: See TracBrowser for help on using the repository browser.