Fork me on GitHub

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

Last change on this file since 387 was 349, checked in by severine ovyn, 16 years ago

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