hepmc - Blame information for rev 537

Subversion Repositories:
Rev:
Rev Author Line No. Line
2 garren 1
  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