Line data Source code
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 "HepMC/GenEvent.h"
16 : #include "HepMC/GenCrossSection.h"
17 : #include "HepMC/Version.h"
18 : #include "HepMC/StreamHelpers.h"
19 :
20 : namespace 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 0 : m_signal_process_id(signal_process_id),
30 0 : m_event_number(event_number),
31 0 : m_mpi(-1),
32 0 : m_event_scale(-1),
33 0 : m_alphaQCD(-1),
34 0 : m_alphaQED(-1),
35 0 : m_signal_process_vertex(signal_vertex),
36 0 : m_beam_particle_1(0),
37 0 : m_beam_particle_2(0),
38 0 : m_weights(weights),
39 0 : m_random_states(random_states),
40 0 : m_vertex_barcodes(),
41 0 : m_particle_barcodes(),
42 0 : m_cross_section(0),
43 0 : m_heavy_ion(0),
44 0 : m_pdf_info(0),
45 0 : m_momentum_unit(mom),
46 0 : m_position_unit(len)
47 0 : {
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 0 : }
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 0 : m_signal_process_id(signal_process_id),
64 0 : m_event_number(event_number),
65 0 : m_mpi(-1),
66 0 : m_event_scale(-1),
67 0 : m_alphaQCD(-1),
68 0 : m_alphaQED(-1),
69 0 : m_signal_process_vertex(signal_vertex),
70 0 : m_beam_particle_1(0),
71 0 : m_beam_particle_2(0),
72 0 : m_weights(weights),
73 0 : m_random_states(random_states),
74 0 : m_vertex_barcodes(),
75 0 : m_particle_barcodes(),
76 0 : m_cross_section(0),
77 0 : m_heavy_ion( new HeavyIon(ion) ),
78 0 : m_pdf_info( new PdfInfo(pdf) ),
79 0 : m_momentum_unit(mom),
80 0 : m_position_unit(len)
81 0 : {
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 0 : }
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 0 : m_signal_process_id(signal_process_id),
96 0 : m_event_number(event_number),
97 0 : m_mpi(-1),
98 0 : m_event_scale(-1),
99 0 : m_alphaQCD(-1),
100 0 : m_alphaQED(-1),
101 0 : m_signal_process_vertex(signal_vertex),
102 0 : m_beam_particle_1(0),
103 0 : m_beam_particle_2(0),
104 0 : m_weights(weights),
105 0 : m_random_states(random_states),
106 0 : m_vertex_barcodes(),
107 0 : m_particle_barcodes(),
108 0 : m_cross_section(0),
109 0 : m_heavy_ion(0),
110 0 : m_pdf_info(0),
111 0 : m_momentum_unit(mom),
112 0 : m_position_unit(len)
113 0 : {
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 0 : }
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 0 : m_signal_process_id(signal_process_id),
131 0 : m_event_number(event_number),
132 0 : m_mpi(-1),
133 0 : m_event_scale(-1),
134 0 : m_alphaQCD(-1),
135 0 : m_alphaQED(-1),
136 0 : m_signal_process_vertex(signal_vertex),
137 0 : m_beam_particle_1(0),
138 0 : m_beam_particle_2(0),
139 0 : m_weights(weights),
140 0 : m_random_states(random_states),
141 0 : m_vertex_barcodes(),
142 0 : m_particle_barcodes(),
143 0 : m_cross_section(0),
144 0 : m_heavy_ion( new HeavyIon(ion) ),
145 0 : m_pdf_info( new PdfInfo(pdf) ),
146 0 : m_momentum_unit(mom),
147 0 : m_position_unit(len)
148 0 : {
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 0 : }
155 :
156 : GenEvent::GenEvent( const GenEvent& inevent )
157 0 : : m_signal_process_id ( inevent.signal_process_id() ),
158 0 : m_event_number ( inevent.event_number() ),
159 0 : m_mpi ( inevent.mpi() ),
160 0 : m_event_scale ( inevent.event_scale() ),
161 0 : m_alphaQCD ( inevent.alphaQCD() ),
162 0 : m_alphaQED ( inevent.alphaQED() ),
163 0 : m_signal_process_vertex( /* inevent.m_signal_process_vertex */ ),
164 0 : m_beam_particle_1 ( /* inevent.m_beam_particle_1 */ ),
165 0 : m_beam_particle_2 ( /* inevent.m_beam_particle_2 */ ),
166 0 : m_weights ( /* inevent.m_weights */ ),
167 0 : m_random_states ( /* inevent.m_random_states */ ),
168 0 : m_vertex_barcodes ( /* inevent.m_vertex_barcodes */ ),
169 0 : m_particle_barcodes ( /* inevent.m_particle_barcodes */ ),
170 0 : m_cross_section ( inevent.cross_section() ? new GenCrossSection(*inevent.cross_section()) : 0 ),
171 0 : m_heavy_ion ( inevent.heavy_ion() ? new HeavyIon(*inevent.heavy_ion()) : 0 ),
172 0 : m_pdf_info ( inevent.pdf_info() ? new PdfInfo(*inevent.pdf_info()) : 0 ),
173 0 : m_momentum_unit ( inevent.momentum_unit() ),
174 0 : m_position_unit ( inevent.length_unit() )
175 0 : {
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 0 : std::map<const GenVertex*,GenVertex*> map_in_to_new;
185 0 : for ( GenEvent::vertex_const_iterator v = inevent.vertices_begin();
186 0 : v != inevent.vertices_end(); ++v ) {
187 0 : GenVertex* newvertex = new GenVertex(
188 0 : (*v)->position(), (*v)->id(), (*v)->weights() );
189 0 : newvertex->suggest_barcode( (*v)->barcode() );
190 0 : map_in_to_new[*v] = newvertex;
191 0 : add_vertex( newvertex );
192 : }
193 : // 2. copy the signal process vertex info.
194 0 : if ( inevent.signal_process_vertex() ) {
195 0 : set_signal_process_vertex(
196 0 : map_in_to_new[inevent.signal_process_vertex()] );
197 0 : } 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 0 : for ( GenEvent::particle_const_iterator p = inevent.particles_begin();
204 0 : p != inevent.particles_end(); ++p )
205 : {
206 0 : GenParticle* oldparticle = *p;
207 0 : GenParticle* newparticle = new GenParticle(*oldparticle);
208 0 : if ( oldparticle->end_vertex() ) {
209 0 : map_in_to_new[ oldparticle->end_vertex() ]->
210 : add_particle_in(newparticle);
211 : }
212 0 : if ( oldparticle->production_vertex() ) {
213 0 : map_in_to_new[ oldparticle->production_vertex() ]->
214 : add_particle_out(newparticle);
215 : }
216 0 : if ( oldparticle == inevent.beam_particles().first ) beam1 = newparticle;
217 0 : if ( oldparticle == inevent.beam_particles().second ) beam2 = newparticle;
218 : }
219 0 : set_beam_particles( beam1, beam2 );
220 : //
221 : // 4. now that vtx/particles are copied, copy weights and random states
222 0 : set_random_states( inevent.random_states() );
223 0 : weights() = inevent.weights();
224 0 : }
225 :
226 : void GenEvent::swap( GenEvent & other )
227 : {
228 : // if a container has a swap method, use that for improved performance
229 0 : std::swap(m_signal_process_id , other.m_signal_process_id );
230 0 : std::swap(m_event_number , other.m_event_number );
231 0 : std::swap(m_mpi , other.m_mpi );
232 0 : std::swap(m_event_scale , other.m_event_scale );
233 0 : std::swap(m_alphaQCD , other.m_alphaQCD );
234 0 : std::swap(m_alphaQED , other.m_alphaQED );
235 0 : std::swap(m_signal_process_vertex, other.m_signal_process_vertex);
236 0 : std::swap(m_beam_particle_1 , other.m_beam_particle_1 );
237 0 : std::swap(m_beam_particle_2 , other.m_beam_particle_2 );
238 0 : m_weights.swap( other.m_weights );
239 0 : m_random_states.swap( other.m_random_states );
240 0 : m_vertex_barcodes.swap( other.m_vertex_barcodes );
241 0 : m_particle_barcodes.swap( other.m_particle_barcodes );
242 0 : std::swap(m_cross_section , other.m_cross_section );
243 0 : std::swap(m_heavy_ion , other.m_heavy_ion );
244 0 : std::swap(m_pdf_info , other.m_pdf_info );
245 0 : std::swap(m_momentum_unit , other.m_momentum_unit );
246 0 : std::swap(m_position_unit , other.m_position_unit );
247 : // must now adjust GenVertex back pointers
248 0 : for ( GenEvent::vertex_const_iterator vthis = vertices_begin();
249 0 : vthis != vertices_end(); ++vthis ) {
250 0 : (*vthis)->change_parent_event_( this );
251 : }
252 0 : for ( GenEvent::vertex_const_iterator voth = other.vertices_begin();
253 0 : voth != other.vertices_end(); ++voth ) {
254 0 : (*voth)->change_parent_event_( &other );
255 : }
256 0 : }
257 :
258 : GenEvent::~GenEvent()
259 0 : {
260 : /// Deep destructor.
261 : /// deletes all vertices/particles in this GenEvent
262 : /// deletes the associated HeavyIon and PdfInfo
263 0 : delete_all_vertices();
264 0 : delete m_cross_section;
265 0 : delete m_heavy_ion;
266 0 : delete m_pdf_info;
267 0 : }
268 :
269 : GenEvent& GenEvent::operator=( const GenEvent& inevent )
270 : {
271 : /// best practices implementation
272 0 : GenEvent tmp( inevent );
273 0 : swap( tmp );
274 : return *this;
275 0 : }
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 0 : ostr << "________________________________________"
283 0 : << "________________________________________\n";
284 0 : ostr << "GenEvent: #" << event_number()
285 0 : << " ID=" << signal_process_id()
286 0 : << " SignalProcessGenVertex Barcode: "
287 0 : << ( signal_process_vertex() ? signal_process_vertex()->barcode()
288 : : 0 )
289 0 : << "\n";
290 0 : write_units( ostr );
291 0 : write_cross_section(ostr);
292 0 : ostr << " Entries this event: " << vertices_size() << " vertices, "
293 0 : << particles_size() << " particles.\n";
294 0 : if( m_beam_particle_1 && m_beam_particle_2 ) {
295 0 : ostr << " Beam Particle barcodes: " << beam_particles().first->barcode() << " "
296 0 : << beam_particles().second->barcode() << " \n";
297 0 : } else {
298 0 : ostr << " Beam Particles are not defined.\n";
299 : }
300 : // Random State
301 0 : ostr << " RndmState(" << m_random_states.size() << ")=";
302 0 : for ( std::vector<long>::const_iterator rs
303 0 : = m_random_states.begin();
304 0 : rs != m_random_states.end(); ++rs ) { ostr << *rs << " "; }
305 0 : ostr << "\n";
306 : // Weights
307 0 : ostr << " Wgts(" << weights().size() << ")=";
308 0 : weights().print(ostr);
309 0 : ostr << " EventScale " << event_scale()
310 0 : << " [energy] \t alphaQCD=" << alphaQCD()
311 0 : << "\t alphaQED=" << alphaQED() << std::endl;
312 : // print a legend to describe the particle info
313 0 : ostr << " GenParticle Legend\n";
314 0 : ostr << " Barcode PDG ID "
315 0 : << "( Px, Py, Pz, E )"
316 0 : << " Stat DecayVtx\n";
317 0 : ostr << "________________________________________"
318 0 : << "________________________________________\n";
319 : // Print all Vertices
320 0 : for ( GenEvent::vertex_const_iterator vtx = this->vertices_begin();
321 0 : vtx != this->vertices_end(); ++vtx ) {
322 0 : (*vtx)->print(ostr);
323 : }
324 0 : ostr << "________________________________________"
325 0 : << "________________________________________" << std::endl;
326 0 : }
327 :
328 : void GenEvent::print_version( std::ostream& ostr ) const {
329 0 : ostr << "---------------------------------------------" << std::endl;
330 0 : writeVersion( ostr );
331 0 : ostr << "---------------------------------------------" << std::endl;
332 0 : }
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 0 : if ( !vtx ) return false;
338 : // if vtx previously pointed to another GenEvent, remove it from that
339 : // GenEvent's list
340 0 : if ( vtx->parent_event() && vtx->parent_event() != this ) {
341 0 : bool remove_status = vtx->parent_event()->remove_vertex( vtx );
342 0 : if ( !remove_status ) {
343 0 : std::cerr << "GenEvent::add_vertex ERROR "
344 0 : << "GenVertex::parent_event points to \n"
345 0 : << "an event that does not point back to the "
346 0 : << "GenVertex. \n This probably indicates a deeper "
347 0 : << "problem. " << std::endl;
348 0 : }
349 0 : }
350 : //
351 : // setting the vertex parent also inserts the vertex into this
352 : // event
353 0 : vtx->set_parent_event_( this );
354 0 : return ( m_vertex_barcodes.count(vtx->barcode()) ? true : false );
355 0 : }
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 0 : if ( m_signal_process_vertex == vtx ) m_signal_process_vertex = 0;
361 0 : if ( vtx->parent_event() == this ) vtx->set_parent_event_( 0 );
362 0 : 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 0 : delete_all_vertices();
371 : // remove existing objects and set pointers to null
372 0 : delete m_cross_section;
373 0 : m_cross_section = 0;
374 0 : delete m_heavy_ion;
375 0 : m_heavy_ion = 0;
376 0 : delete m_pdf_info;
377 0 : m_pdf_info = 0;
378 0 : m_signal_process_id = 0;
379 0 : m_beam_particle_1 = 0;
380 0 : m_beam_particle_2 = 0;
381 0 : m_event_number = 0;
382 0 : m_mpi = -1;
383 0 : m_event_scale = -1;
384 0 : m_alphaQCD = -1;
385 0 : m_alphaQED = -1;
386 0 : m_weights = std::vector<double>();
387 0 : m_random_states = std::vector<long>();
388 : // resetting unit information
389 0 : m_momentum_unit = Units::default_momentum_unit();
390 0 : m_position_unit = Units::default_length_unit();
391 : // error check just to be safe
392 0 : if ( m_vertex_barcodes.size() != 0
393 0 : || m_particle_barcodes.size() != 0 ) {
394 0 : std::cerr << "GenEvent::clear() strange result ... \n"
395 0 : << "either the particle and/or the vertex map isn't empty" << std::endl;
396 0 : std::cerr << "Number vtx,particle the event after deleting = "
397 0 : << m_vertex_barcodes.size() << " "
398 0 : << m_particle_barcodes.size() << std::endl;
399 0 : }
400 0 : return;
401 0 : }
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 0 : while ( !vertices_empty() ) {
412 0 : GenVertex* vtx = ( m_vertex_barcodes.begin() )->second;
413 0 : m_vertex_barcodes.erase( m_vertex_barcodes.begin() );
414 0 : delete vtx;
415 : }
416 : //
417 : // Error checking:
418 0 : if ( !vertices_empty() || ! particles_empty() ) {
419 0 : std::cerr << "GenEvent::delete_all_vertices strange result ... "
420 0 : << "after deleting all vertices, \nthe particle and "
421 0 : << "vertex maps aren't empty.\n This probably "
422 0 : << "indicates deeper problems or memory leak in the "
423 0 : << "code." << std::endl;
424 0 : std::cerr << "Number vtx,particle the event after deleting = "
425 0 : << m_vertex_barcodes.size() << " "
426 0 : << m_particle_barcodes.size() << std::endl;
427 0 : }
428 0 : }
429 :
430 : bool GenEvent::set_barcode( GenParticle* p, int suggested_barcode )
431 : {
432 0 : if ( p->parent_event() != this ) {
433 0 : std::cerr << "GenEvent::set_barcode attempted, but the argument's"
434 0 : << "\n parent_event is not this ... request rejected."
435 0 : << std::endl;
436 0 : 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 0 : if ( p->barcode() != 0 && p->barcode() != suggested_barcode ) {
443 0 : if ( m_particle_barcodes.count(p->barcode()) &&
444 0 : m_particle_barcodes[p->barcode()] == p ) {
445 0 : m_particle_barcodes.erase( p->barcode() );
446 0 : }
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 0 : if ( suggested_barcode > 0 ) {
456 0 : if ( m_particle_barcodes.count(suggested_barcode) ) {
457 : // the suggested_barcode is already used.
458 0 : if ( m_particle_barcodes[suggested_barcode] == p ) {
459 : // but it was used for this particle ... so everythings ok
460 0 : p->set_barcode_( suggested_barcode );
461 0 : return true;
462 : }
463 : insert_success = false;
464 0 : suggested_barcode = 0;
465 : } else { // suggested barcode is OK, proceed to insert
466 0 : m_particle_barcodes[suggested_barcode] = p;
467 0 : p->set_barcode_( suggested_barcode );
468 0 : return true;
469 : }
470 0 : }
471 : //
472 : // Other possibility -- a valid barcode has not been suggested,
473 : // so one is automatically generated
474 0 : if ( suggested_barcode < 0 ) insert_success = false;
475 0 : if ( suggested_barcode <= 0 ) {
476 0 : if ( !m_particle_barcodes.empty() ) {
477 : // in this case we find the highest barcode that was used,
478 : // and increment it by 1
479 0 : suggested_barcode = m_particle_barcodes.rbegin()->first;
480 0 : ++suggested_barcode;
481 0 : }
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 0 : if ( suggested_barcode <= 10000 ) suggested_barcode = 10001;
489 : }
490 : // At this point we should have a valid barcode
491 0 : if ( m_particle_barcodes.count(suggested_barcode) ) {
492 0 : std::cerr << "GenEvent::set_barcode ERROR, this should never "
493 0 : << "happen \n report bug to matt.dobbs@cern.ch"
494 0 : << std::endl;
495 0 : }
496 0 : m_particle_barcodes[suggested_barcode] = p;
497 0 : p->set_barcode_( suggested_barcode );
498 0 : return insert_success;
499 0 : }
500 :
501 : bool GenEvent::set_barcode( GenVertex* v, int suggested_barcode )
502 : {
503 0 : if ( v->parent_event() != this ) {
504 0 : std::cerr << "GenEvent::set_barcode attempted, but the argument's"
505 0 : << "\n parent_event is not this ... request rejected."
506 0 : << std::endl;
507 0 : 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 0 : if ( v->barcode() != 0 && v->barcode() != suggested_barcode ) {
514 0 : if ( m_vertex_barcodes.count(v->barcode()) &&
515 0 : m_vertex_barcodes[v->barcode()] == v ) {
516 0 : m_vertex_barcodes.erase( v->barcode() );
517 0 : }
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 0 : if ( suggested_barcode < 0 ) {
528 0 : if ( m_vertex_barcodes.count(suggested_barcode) ) {
529 : // the suggested_barcode is already used.
530 0 : if ( m_vertex_barcodes[suggested_barcode] == v ) {
531 : // but it was used for this vertex ... so everythings ok
532 0 : v->set_barcode_( suggested_barcode );
533 0 : return true;
534 : }
535 : insert_success = false;
536 0 : suggested_barcode = 0;
537 : } else { // suggested barcode is OK, proceed to insert
538 0 : m_vertex_barcodes[suggested_barcode] = v;
539 0 : v->set_barcode_( suggested_barcode );
540 0 : return true;
541 : }
542 0 : }
543 : //
544 : // Other possibility -- a valid barcode has not been suggested,
545 : // so one is automatically generated
546 0 : if ( suggested_barcode > 0 ) insert_success = false;
547 0 : if ( suggested_barcode >= 0 ) {
548 0 : 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 0 : suggested_barcode = m_vertex_barcodes.rbegin()->first;
552 0 : --suggested_barcode;
553 0 : }
554 0 : if ( suggested_barcode >= 0 ) suggested_barcode = -1;
555 : }
556 : // At this point we should have a valid barcode
557 0 : if ( m_vertex_barcodes.count(suggested_barcode) ) {
558 0 : std::cerr << "GenEvent::set_barcode ERROR, this should never "
559 0 : << "happen \n report bug to matt.dobbs@cern.ch"
560 0 : << std::endl;
561 0 : }
562 0 : m_vertex_barcodes[suggested_barcode] = v;
563 0 : v->set_barcode_( suggested_barcode );
564 0 : return insert_success;
565 0 : }
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 0 : if(m_beam_particle_1 && m_beam_particle_2) {
573 : // now look for a match with the particle "list"
574 0 : for ( particle_const_iterator p = particles_begin();
575 0 : p != particles_end(); ++p ) {
576 0 : if( m_beam_particle_1 == *p ) have1 = true;
577 0 : if( m_beam_particle_2 == *p ) have2 = true;
578 : }
579 0 : }
580 0 : if( have1 && have2 ) return true;
581 0 : return false;
582 0 : }
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 0 : m_beam_particle_1 = bp1;
588 0 : m_beam_particle_2 = bp2;
589 0 : if( m_beam_particle_1 && m_beam_particle_2 ) return true;
590 0 : return false;
591 0 : }
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 0 : return set_beam_particles(bp.first,bp.second);
597 : }
598 :
599 : void GenEvent::write_units( std::ostream & os ) const {
600 0 : os << " Momenutm units:" << std::setw(8) << name(momentum_unit());
601 0 : os << " Position units:" << std::setw(8) << name(length_unit());
602 0 : os << std::endl;
603 0 : }
604 :
605 : void GenEvent::write_cross_section( std::ostream& os ) const
606 : {
607 : // write the GenCrossSection information if the cross section was set
608 0 : if( !cross_section() ) return;
609 0 : if( cross_section()->is_set() ) {
610 0 : os << " Cross Section: " << cross_section()->cross_section() ;
611 0 : os << " +/- " << cross_section()->cross_section_error() ;
612 0 : os << std::endl;
613 0 : }
614 0 : }
615 :
616 : bool GenEvent::use_momentum_unit( Units::MomentumUnit newunit ) {
617 : // currently not exception-safe.
618 : // Easy to fix, though, if needed.
619 0 : if ( m_momentum_unit != newunit ) {
620 0 : 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 0 : for ( GenEvent::particle_iterator p = particles_begin();
624 0 : p != particles_end(); ++p )
625 : {
626 0 : (*p)->convert_momentum(factor);
627 : }
628 : // ...
629 0 : m_momentum_unit = newunit;
630 0 : }
631 0 : return true;
632 0 : }
633 :
634 : bool GenEvent::use_length_unit( Units::LengthUnit newunit ) {
635 : // currently not exception-safe.
636 : // Easy to fix, though, if needed.
637 0 : if ( m_position_unit != newunit ) {
638 0 : 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 0 : for ( GenEvent::vertex_iterator vtx = vertices_begin();
642 0 : vtx != vertices_end(); ++vtx ) {
643 0 : (*vtx)->convert_position(factor);
644 : }
645 : // ...
646 0 : m_position_unit = newunit;
647 0 : }
648 0 : return true;
649 0 : }
650 :
651 : bool GenEvent::use_momentum_unit( std::string& newunit ) {
652 0 : if ( newunit == "MEV" ) return use_momentum_unit( Units::MEV );
653 0 : else if( newunit == "GEV" ) return use_momentum_unit( Units::GEV );
654 0 : else std::cerr << "GenEvent::use_momentum_unit ERROR: use either MEV or GEV\n";
655 0 : return false;
656 0 : }
657 :
658 : bool GenEvent::use_length_unit( std::string& newunit ) {
659 0 : if ( newunit == "MM" ) return use_length_unit( Units::MM );
660 0 : else if( newunit == "CM" ) return use_length_unit( Units::CM );
661 0 : else std::cerr << "GenEvent::use_length_unit ERROR: use either MM or CM\n";
662 0 : return false;
663 0 : }
664 :
665 : void GenEvent::define_units( std::string& new_m, std::string& new_l ) {
666 :
667 0 : if ( new_m == "MEV" ) m_momentum_unit = Units::MEV ;
668 0 : else if( new_m == "GEV" ) m_momentum_unit = Units::GEV ;
669 0 : else std::cerr << "GenEvent::define_units ERROR: use either MEV or GEV\n";
670 :
671 0 : if ( new_l == "MM" ) m_position_unit = Units::MM ;
672 0 : else if( new_l == "CM" ) m_position_unit = Units::CM ;
673 0 : else std::cerr << "GenEvent::define_units ERROR: use either MM or CM\n";
674 :
675 0 : }
676 :
677 : bool GenEvent::is_valid() const {
678 : /// A GenEvent is presumed valid if it has both associated
679 : /// particles and vertices. No other information is checked.
680 0 : if ( vertices_empty() ) return false;
681 0 : if ( particles_empty() ) return false;
682 0 : return true;
683 0 : }
684 :
685 : std::ostream & GenEvent::write_beam_particles(std::ostream & os,
686 : std::pair<HepMC::GenParticle *,HepMC::GenParticle *> pr )
687 : {
688 : GenParticle* p = pr.first;
689 0 : if(!p) {
690 0 : detail::output( os, 0 );
691 0 : } else {
692 0 : detail::output( os, p->barcode() );
693 : }
694 : p = pr.second;
695 0 : if(!p) {
696 0 : detail::output( os, 0 );
697 0 : } else {
698 0 : detail::output( os, p->barcode() );
699 : }
700 :
701 0 : return os;
702 : }
703 :
704 : std::ostream & GenEvent::write_vertex(std::ostream & os, GenVertex const * v)
705 : {
706 0 : if ( !v || !os ) {
707 0 : std::cerr << "GenEvent::write_vertex !v||!os, "
708 0 : << "v="<< v << " setting badbit" << std::endl;
709 0 : os.clear(std::ios::badbit);
710 0 : return os;
711 : }
712 : // First collect info we need
713 : // count the number of orphan particles going into v
714 0 : int num_orphans_in = 0;
715 0 : for ( GenVertex::particles_in_const_iterator p1
716 0 : = v->particles_in_const_begin();
717 0 : p1 != v->particles_in_const_end(); ++p1 ) {
718 0 : if ( !(*p1)->production_vertex() ) ++num_orphans_in;
719 : }
720 : //
721 0 : os << 'V';
722 0 : detail::output( os, v->barcode() ); // v's unique identifier
723 0 : detail::output( os, v->id() );
724 0 : detail::output( os, v->position().x() );
725 0 : detail::output( os, v->position().y() );
726 0 : detail::output( os, v->position().z() );
727 0 : detail::output( os, v->position().t() );
728 0 : detail::output( os, num_orphans_in );
729 0 : detail::output( os, (int)v->particles_out_size() );
730 0 : detail::output( os, (int)v->weights().size() );
731 0 : for ( WeightContainer::const_iterator w = v->weights().begin();
732 0 : w != v->weights().end(); ++w ) {
733 0 : detail::output( os, *w );
734 : }
735 0 : detail::output( os,'\n');
736 : // incoming particles
737 0 : for ( GenVertex::particles_in_const_iterator p2
738 0 : = v->particles_in_const_begin();
739 0 : p2 != v->particles_in_const_end(); ++p2 ) {
740 0 : if ( !(*p2)->production_vertex() ) {
741 0 : write_particle( os, *p2 );
742 0 : }
743 : }
744 : // outgoing particles
745 0 : for ( GenVertex::particles_out_const_iterator p3
746 0 : = v->particles_out_const_begin();
747 0 : p3 != v->particles_out_const_end(); ++p3 ) {
748 0 : write_particle( os, *p3 );
749 : }
750 : return os;
751 0 : }
752 :
753 : std::ostream & GenEvent::write_particle( std::ostream & os, GenParticle const * p )
754 : {
755 0 : if ( !p || !os ) {
756 0 : std::cerr << "GenEvent::write_particle !p||!os, "
757 0 : << "p="<< p << " setting badbit" << std::endl;
758 0 : os.clear(std::ios::badbit);
759 0 : return os;
760 : }
761 0 : os << 'P';
762 0 : detail::output( os, p->barcode() );
763 0 : detail::output( os, p->pdg_id() );
764 0 : detail::output( os, p->momentum().px() );
765 0 : detail::output( os, p->momentum().py() );
766 0 : detail::output( os, p->momentum().pz() );
767 0 : detail::output( os, p->momentum().e() );
768 0 : detail::output( os, p->generated_mass() );
769 0 : detail::output( os, p->status() );
770 0 : detail::output( os, p->polarization().theta() );
771 0 : detail::output( os, p->polarization().phi() );
772 : // since end_vertex is oftentimes null, this CREATES a null vertex
773 : // in the map
774 0 : detail::output( os, ( p->end_vertex() ? p->end_vertex()->barcode() : 0 ) );
775 0 : os << ' ' << p->flow() << "\n";
776 :
777 0 : return os;
778 0 : }
779 :
780 : } // HepMC
|
