hepmc - Blame information for rev 422

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