Subversion Repositories hepmc

Rev

Rev 521 | Details | Compare with Previous | Last modification | View Log | RSS feed

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() << ")=";
432 garren 308
        weights().print(ostr);
2 garren 309
        ostr << " EventScale " << event_scale()
310
             << " [energy] \t alphaQCD=" << alphaQCD()
311
             << "\t alphaQED=" << alphaQED() << std::endl;
312
        // print a legend to describe the particle info
25 garren 313
        ostr << "                                    GenParticle Legend\n";
314
        ostr  << "        Barcode   PDG ID      "
315
              << "( Px,       Py,       Pz,     E )"
316
              << " Stat  DecayVtx\n";
2 garren 317
        ostr << "________________________________________"
318
             << "________________________________________\n";
319
        // Print all Vertices
320
        for ( GenEvent::vertex_const_iterator vtx = this->vertices_begin();
321
              vtx != this->vertices_end(); ++vtx ) {
322
            (*vtx)->print(ostr);
323
        }
324
        ostr << "________________________________________"
325
             << "________________________________________" << std::endl;
326
    }
327
 
173 garren 328
    void GenEvent::print_version( std::ostream& ostr ) const {
329
        ostr << "---------------------------------------------" << std::endl;
330
        writeVersion( ostr );
331
        ostr << "---------------------------------------------" << std::endl;
332
    }
333
 
2 garren 334
    bool GenEvent::add_vertex( GenVertex* vtx ) {
65 garren 335
        /// returns true if successful - generally will only return false
336
        /// if the inserted vertex is already included in the event.
147 garren 337
        if ( !vtx ) return false;
2 garren 338
        // if vtx previously pointed to another GenEvent, remove it from that
339
        // GenEvent's list
340
        if ( vtx->parent_event() && vtx->parent_event() != this ) {
341
            bool remove_status = vtx->parent_event()->remove_vertex( vtx );
342
            if ( !remove_status ) {            
343
                std::cerr << "GenEvent::add_vertex ERROR "
344
                          << "GenVertex::parent_event points to \n"
345
                          << "an event that does not point back to the "
346
                          << "GenVertex. \n This probably indicates a deeper "
347
                          << "problem. " << std::endl;
348
            }
349
        }
350
        //
351
        // setting the vertex parent also inserts the vertex into this
352
        // event
353
        vtx->set_parent_event_( this );
354
        return ( m_vertex_barcodes.count(vtx->barcode()) ? true : false );
355
    }
356
 
357
    bool GenEvent::remove_vertex( GenVertex* vtx ) {
65 garren 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
2 garren 360
        if ( m_signal_process_vertex == vtx ) m_signal_process_vertex = 0;
361
        if ( vtx->parent_event() == this ) vtx->set_parent_event_( 0 );
362
        return ( m_vertex_barcodes.count(vtx->barcode()) ? false : true );
363
    }
96 garren 364
 
365
    void GenEvent::clear()
366
    {
367
        /// remove all information from the event
368
        /// deletes all vertices/particles in this evt
369
        ///
370
        delete_all_vertices();
390 garren 371
        // remove existing objects and set pointers to null
372
        delete m_cross_section;
373
        m_cross_section = 0;
96 garren 374
        delete m_heavy_ion;
390 garren 375
        m_heavy_ion = 0;
96 garren 376
        delete m_pdf_info;
390 garren 377
        m_pdf_info = 0;
96 garren 378
        m_signal_process_id = 0;
105 garren 379
        m_beam_particle_1 = 0;
380
        m_beam_particle_2 = 0;
96 garren 381
        m_event_number = 0;
390 garren 382
        m_mpi = -1;
96 garren 383
        m_event_scale = -1;
384
        m_alphaQCD = -1;
385
        m_alphaQED = -1;
386
        m_weights = std::vector<double>();
210 garren 387
        m_random_states = std::vector<long>();
390 garren 388
        // resetting unit information
389
        m_momentum_unit = Units::default_momentum_unit();
390
        m_position_unit = Units::default_length_unit();
96 garren 391
        // error check just to be safe
392
        if ( m_vertex_barcodes.size() != 0
393
             || m_particle_barcodes.size() != 0 ) {
394
            std::cerr << "GenEvent::clear() strange result ... \n"
395
                      << "either the particle and/or the vertex map isn't empty" << std::endl;
396
            std::cerr << "Number vtx,particle the event after deleting = "
397
                      << m_vertex_barcodes.size() << "  "
398
                      << m_particle_barcodes.size() << std::endl;
399
        }
400
        return;
401
    }
2 garren 402
 
403
    void GenEvent::delete_all_vertices() {
65 garren 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.
147 garren 409
 
2 garren 410
        // delete each vertex individually (this deletes particles as well)
147 garren 411
        while ( !vertices_empty() ) {
2 garren 412
            GenVertex* vtx = ( m_vertex_barcodes.begin() )->second;
413
            m_vertex_barcodes.erase( m_vertex_barcodes.begin() );
414
            delete vtx;
415
        }
416
        //
417
        // Error checking:
418
        if ( !vertices_empty() || ! particles_empty() ) {
419
            std::cerr << "GenEvent::delete_all_vertices strange result ... "
420
                      << "after deleting all vertices, \nthe particle and "
421
                      << "vertex maps aren't empty.\n  This probably "
422
                      << "indicates deeper problems or memory leak in the "
423
                      << "code." << std::endl;
424
            std::cerr << "Number vtx,particle the event after deleting = "
425
                      << m_vertex_barcodes.size() << "  "
426
                      << m_particle_barcodes.size() << std::endl;
427
        }
428
    }
429
 
430
    bool GenEvent::set_barcode( GenParticle* p, int suggested_barcode )
431
    {
432
        if ( p->parent_event() != this ) {
433
            std::cerr << "GenEvent::set_barcode attempted, but the argument's"
434
                      << "\n parent_event is not this ... request rejected."
435
                      << std::endl;
436
            return false;
437
        }
438
        // M.Dobbs  Nov 4, 2002
439
        // First we must check to see if the particle already has a
440
        // barcode which is different from the suggestion. If yes, we
441
        // remove it from the particle map.
442
        if ( p->barcode() != 0 && p->barcode() != suggested_barcode ) {
443
            if ( m_particle_barcodes.count(p->barcode()) &&
444
                 m_particle_barcodes[p->barcode()] == p ) {
445
                m_particle_barcodes.erase( p->barcode() );
446
            }
447
            // At this point either the particle is NOT in
448
            // m_particle_barcodes, or else it is in the map, but
449
            // already with the suggested barcode.
450
        }
451
        //
452
        // First case --- a valid barcode has been suggested
453
        //     (valid barcodes are numbers greater than zero)
454
        bool insert_success = true;
455
        if ( suggested_barcode > 0 ) {
456
            if ( m_particle_barcodes.count(suggested_barcode) ) {
457
                // the suggested_barcode is already used.
458
                if ( m_particle_barcodes[suggested_barcode] == p ) {
459
                    // but it was used for this particle ... so everythings ok
460
                    p->set_barcode_( suggested_barcode );
461
                    return true;
462
                }
463
                insert_success = false;
464
                suggested_barcode = 0;
465
            } else { // suggested barcode is OK, proceed to insert
466
                m_particle_barcodes[suggested_barcode] = p;
467
                p->set_barcode_( suggested_barcode );
468
                return true;
469
            }
470
        }
471
        //
472
        // Other possibility -- a valid barcode has not been suggested,
473
        //    so one is automatically generated
474
        if ( suggested_barcode < 0 ) insert_success = false;
475
        if ( suggested_barcode <= 0 ) {
476
            if ( !m_particle_barcodes.empty() ) {
477
                // in this case we find the highest barcode that was used,
478
                // and increment it by 1
479
                suggested_barcode = m_particle_barcodes.rbegin()->first;
480
                ++suggested_barcode;
481
            }
482
            // For the automatically assigned barcodes, the first one
483
            //   we use is 10001 ... this is just because when we read 
484
            //   events from HEPEVT, we will suggest their index as the
485
            //   barcode, and that index has maximum value 10000.
486
            //  This way we will immediately be able to recognize the hepevt
487
            //   particles from the auto-assigned ones.
488
            if ( suggested_barcode <= 10000 ) suggested_barcode = 10001;
489
        }
490
        // At this point we should have a valid barcode
491
        if ( m_particle_barcodes.count(suggested_barcode) ) {
492
            std::cerr << "GenEvent::set_barcode ERROR, this should never "
493
                      << "happen \n report bug to matt.dobbs@cern.ch"
494
                      << std::endl;
495
        }
496
        m_particle_barcodes[suggested_barcode] = p;
497
        p->set_barcode_( suggested_barcode );
498
        return insert_success;
499
    }
500
 
501
    bool GenEvent::set_barcode( GenVertex* v, int suggested_barcode )
502
    {
503
        if ( v->parent_event() != this ) {
504
            std::cerr << "GenEvent::set_barcode attempted, but the argument's"
505
                      << "\n parent_event is not this ... request rejected."
506
                      << std::endl;
507
            return false;
508
        }
509
        // M.Dobbs Nov 4, 2002
510
        // First we must check to see if the vertex already has a
511
        // barcode which is different from the suggestion. If yes, we
512
        // remove it from the vertex map.
513
        if ( v->barcode() != 0 && v->barcode() != suggested_barcode ) {
514
            if ( m_vertex_barcodes.count(v->barcode()) &&
515
                 m_vertex_barcodes[v->barcode()] == v ) {
516
                m_vertex_barcodes.erase( v->barcode() );
517
            }
518
            // At this point either the vertex is NOT in
519
            // m_vertex_barcodes, or else it is in the map, but
520
            // already with the suggested barcode.
521
        }
522
 
523
        //
524
        // First case --- a valid barcode has been suggested
525
        //     (valid barcodes are numbers greater than zero)
526
        bool insert_success = true;
527
        if ( suggested_barcode < 0 ) {
528
            if ( m_vertex_barcodes.count(suggested_barcode) ) {
529
                // the suggested_barcode is already used.
530
                if ( m_vertex_barcodes[suggested_barcode] == v ) {
531
                    // but it was used for this vertex ... so everythings ok
532
                    v->set_barcode_( suggested_barcode );
533
                    return true;
534
                }
535
                insert_success = false;
536
                suggested_barcode = 0;
537
            } else { // suggested barcode is OK, proceed to insert
538
                m_vertex_barcodes[suggested_barcode] = v;
539
                v->set_barcode_( suggested_barcode );
540
                return true;
541
            }
542
        }
543
        //
544
        // Other possibility -- a valid barcode has not been suggested,
545
        //    so one is automatically generated
546
        if ( suggested_barcode > 0 ) insert_success = false;
547
        if ( suggested_barcode >= 0 ) {
548
            if ( !m_vertex_barcodes.empty() ) {
549
                // in this case we find the highest barcode that was used,
550
                // and increment it by 1, (vertex barcodes are negative)
551
                suggested_barcode = m_vertex_barcodes.rbegin()->first;
552
                --suggested_barcode;
553
            }
554
            if ( suggested_barcode >= 0 ) suggested_barcode = -1;
555
        }
556
        // At this point we should have a valid barcode
557
        if ( m_vertex_barcodes.count(suggested_barcode) ) {
558
            std::cerr << "GenEvent::set_barcode ERROR, this should never "
559
                      << "happen \n report bug to matt.dobbs@cern.ch"
560
                      << std::endl;
561
        }
562
        m_vertex_barcodes[suggested_barcode] = v;
563
        v->set_barcode_( suggested_barcode );
564
        return insert_success;
565
    }
566
 
105 garren 567
    /// test to see if we have two valid beam particles
568
    bool  GenEvent::valid_beam_particles() const {
569
        bool have1 = false;
570
        bool have2 = false;
571
        // first check that both are defined
572
        if(m_beam_particle_1 && m_beam_particle_2) {
573
            // now look for a match with the particle "list"
574
            for ( particle_const_iterator p = particles_begin();
147 garren 575
                  p != particles_end(); ++p ) {
105 garren 576
                if( m_beam_particle_1 == *p ) have1 = true;
577
                if( m_beam_particle_2 == *p ) have2 = true;
578
            }
579
        }
580
        if( have1 && have2 ) return true;
581
        return false;
582
    }
583
 
584
    /// construct the beam particle information using pointers to GenParticle
585
    /// returns false if either GenParticle* is null
586
    bool  GenEvent::set_beam_particles(GenParticle* bp1, GenParticle* bp2) {
587
        m_beam_particle_1 = bp1;
588
        m_beam_particle_2 = bp2;
589
        if( m_beam_particle_1 && m_beam_particle_2 ) return true;
590
        return false;
591
    }
592
 
593
    /// construct the beam particle information using a std::pair of pointers to GenParticle
594
    /// returns false if either GenParticle* is null
179 garren 595
    bool  GenEvent::set_beam_particles(std::pair<HepMC::GenParticle*, HepMC::GenParticle*> const & bp) {
105 garren 596
        return set_beam_particles(bp.first,bp.second);
597
    }
598
 
274 garren 599
    void GenEvent::write_units( std::ostream & os ) const {
537 garren 600
        os << " Momentum units:" << std::setw(8) << name(momentum_unit());
335 garren 601
        os << "     Position units:" << std::setw(8) << name(length_unit());
274 garren 602
        os << std::endl;
603
    }
2 garren 604
 
421 garren 605
    void GenEvent::write_cross_section( std::ostream& os ) const
606
    {
607
        // write the GenCrossSection information if the cross section was set
608
        if( !cross_section() ) return;
609
        if( cross_section()->is_set() ) {
610
            os << " Cross Section: " << cross_section()->cross_section() ;
611
            os << " +/- " << cross_section()->cross_section_error() ;
612
            os << std::endl;
613
        }
614
    }
615
 
335 garren 616
   bool GenEvent::use_momentum_unit( Units::MomentumUnit newunit ) {
617
        // currently not exception-safe. 
618
        // Easy to fix, though, if needed.
619
        if ( m_momentum_unit != newunit ) {
620
            const double factor = Units::conversion_factor( m_momentum_unit, newunit );
621
            // multiply all momenta by 'factor',  
338 garren 622
            // loop is entered only if particle list is not empty
335 garren 623
            for ( GenEvent::particle_iterator p = particles_begin();
624
                                              p != particles_end(); ++p )
625
            {
626
                (*p)->convert_momentum(factor);
627
            }
628
            // ... 
629
            m_momentum_unit = newunit;
274 garren 630
        }
335 garren 631
        return true;
274 garren 632
    }
335 garren 633
 
634
    bool GenEvent::use_length_unit( Units::LengthUnit newunit ) {
635
        // currently not exception-safe. 
636
        // Easy to fix, though, if needed.
637
        if ( m_position_unit != newunit ) {
638
            const double factor = Units::conversion_factor( m_position_unit, newunit );
639
            // multiply all lengths by 'factor', 
338 garren 640
            // loop is entered only if vertex list is not empty
335 garren 641
            for ( GenEvent::vertex_iterator vtx = vertices_begin();
642
                                            vtx != vertices_end(); ++vtx ) {
643
                (*vtx)->convert_position(factor);
644
            }
645
            // ... 
646
            m_position_unit = newunit;
647
        }
648
        return true;
649
    }  
274 garren 650
 
335 garren 651
    bool GenEvent::use_momentum_unit( std::string& newunit ) {
652
        if     ( newunit == "MEV" ) return use_momentum_unit( Units::MEV );
653
        else if( newunit == "GEV" ) return use_momentum_unit( Units::GEV );
654
        else std::cerr << "GenEvent::use_momentum_unit ERROR: use either MEV or GEV\n";
655
        return false;
274 garren 656
    }
335 garren 657
 
658
    bool GenEvent::use_length_unit( std::string& newunit ) {
659
        if     ( newunit == "MM" ) return use_length_unit( Units::MM );
660
        else if( newunit == "CM" ) return use_length_unit( Units::CM );
508 garren 661
        else std::cerr << "GenEvent::use_length_unit ERROR: use either MM or CM\n";
335 garren 662
        return false;
663
    }  
521 garren 664
 
665
    void GenEvent::define_units( std::string& new_m, std::string& new_l ) {
274 garren 666
 
521 garren 667
        if     ( new_m == "MEV" ) m_momentum_unit = Units::MEV ;
668
        else if( new_m == "GEV" ) m_momentum_unit = Units::GEV ;
669
        else std::cerr << "GenEvent::define_units ERROR: use either MEV or GEV\n";
670
 
671
        if     ( new_l == "MM" ) m_position_unit = Units::MM ;
672
        else if( new_l == "CM" ) m_position_unit = Units::CM ;
673
        else std::cerr << "GenEvent::define_units ERROR: use either MM or CM\n";
674
 
675
    }
676
 
390 garren 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
        if ( vertices_empty() ) return false;
681
        if ( particles_empty() ) return false;
682
        return true;
683
    }
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
        if(!p) {
690
           detail::output( os, 0 );
691
        } else {
692
           detail::output( os, p->barcode() );
693
        }
694
        p = pr.second;
695
        if(!p) {
696
           detail::output( os, 0 );
697
        } else {
698
           detail::output( os, p->barcode() );
699
        }
700
 
701
        return os;
702
    }
703
 
704
    std::ostream & GenEvent::write_vertex(std::ostream & os, GenVertex const * v)
705
    {
706
        if ( !v || !os ) {
707
            std::cerr << "GenEvent::write_vertex !v||!os, "
708
                      << "v="<< v << " setting badbit" << std::endl;
709
            os.clear(std::ios::badbit);
710
            return os;
711
        }
712
        // First collect info we need
713
        // count the number of orphan particles going into v
714
        int num_orphans_in = 0;
715
        for ( GenVertex::particles_in_const_iterator p1
716
                  = v->particles_in_const_begin();
717
              p1 != v->particles_in_const_end(); ++p1 ) {
718
            if ( !(*p1)->production_vertex() ) ++num_orphans_in;
719
        }
720
        //
721
        os << 'V';
722
        detail::output( os, v->barcode() ); // v's unique identifier
723
        detail::output( os, v->id() );
724
        detail::output( os, v->position().x() );
725
        detail::output( os, v->position().y() );
726
        detail::output( os, v->position().z() );
727
        detail::output( os, v->position().t() );
728
        detail::output( os, num_orphans_in );
729
        detail::output( os, (int)v->particles_out_size() );
730
        detail::output( os, (int)v->weights().size() );
731
        for ( WeightContainer::const_iterator w = v->weights().begin();
732
              w != v->weights().end(); ++w ) {
733
            detail::output( os, *w );
734
        }
735
        detail::output( os,'\n');
736
        // incoming particles
737
        for ( GenVertex::particles_in_const_iterator p2
738
                  = v->particles_in_const_begin();
739
              p2 != v->particles_in_const_end(); ++p2 ) {
740
            if ( !(*p2)->production_vertex() ) {
741
                write_particle( os, *p2 );
742
            }
743
        }
744
        // outgoing particles
745
        for ( GenVertex::particles_out_const_iterator p3
746
                  = v->particles_out_const_begin();
747
              p3 != v->particles_out_const_end(); ++p3 ) {
748
            write_particle( os, *p3 );
749
        }
750
        return os;
751
    }
752
 
753
    std::ostream & GenEvent::write_particle( std::ostream & os, GenParticle const * p )
754
    {
755
        if ( !p || !os ) {
756
            std::cerr << "GenEvent::write_particle !p||!os, "
757
                      << "p="<< p << " setting badbit" << std::endl;
758
            os.clear(std::ios::badbit);
759
            return os;
760
        }
761
        os << 'P';
762
        detail::output( os, p->barcode() );
763
        detail::output( os, p->pdg_id() );
764
        detail::output( os, p->momentum().px() );
765
        detail::output( os, p->momentum().py() );
766
        detail::output( os, p->momentum().pz() );
767
        detail::output( os, p->momentum().e() );
768
        detail::output( os, p->generated_mass() );
769
        detail::output( os, p->status() );
770
        detail::output( os, p->polarization().theta() );
771
        detail::output( os, p->polarization().phi() );
772
        // since end_vertex is oftentimes null, this CREATES a null vertex
773
        // in the map
774
        detail::output( os,   ( p->end_vertex() ? p->end_vertex()->barcode() : 0 )  );
775
        os << ' ' << p->flow() << "\n";
776
 
777
        return os;
778
    }
779
 
2 garren 780
} // HepMC