hepmc - Rev 476
Subversion Repositories:
Rev 474 -
Rev 508 -
Go to most recent revision -
Blame -
Compare with Previous -
Last modification -
View Log
- Download
- 
//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, September 1999
// Updated: 7.1.2000 iterators complete and working!
// Updated: 10.1.2000 GenEvent::vertex, particle iterators are made
// constant WRT this event ... note that
// GenVertex::***_iterator is not const, since it must
// be able to return a mutable pointer to itself.
// Updated: 08.2.2000 the event now holds a set of all attached vertices
// rather than just the roots of the graph
// Event record for MC generators (for use at any stage of generation)
//////////////////////////////////////////////////////////////////////////
#include <iomanip>
#include "HepMC/GenEvent.h"
#include "HepMC/GenCrossSection.h"
#include "HepMC/Version.h"
#include "HepMC/StreamHelpers.h"
namespace HepMC {
GenEvent::GenEvent( int signal_process_id,
int event_number,
GenVertex* signal_vertex,
const WeightContainer& weights,
const std::vector<long>& random_states,
Units::MomentumUnit mom,
Units::LengthUnit len ) :
m_signal_process_id(signal_process_id),
m_event_number(event_number),
m_mpi(-1),
m_event_scale(-1),
m_alphaQCD(-1),
m_alphaQED(-1),
m_signal_process_vertex(signal_vertex),
m_beam_particle_1(0),
m_beam_particle_2(0),
m_weights(weights),
m_random_states(random_states),
m_vertex_barcodes(),
m_particle_barcodes(),
m_cross_section(0),
m_heavy_ion(0),
m_pdf_info(0),
m_momentum_unit(mom),
m_position_unit(len)
{
/// This constructor only allows null pointers to HeavyIon and PdfInfo
///
/// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
/// are as suggested in hep-ph/0109068, "Generic Interface..."
///
}
GenEvent::GenEvent( int signal_process_id, int event_number,
GenVertex* signal_vertex,
const WeightContainer& weights,
const std::vector<long>& random_states,
const HeavyIon& ion,
const PdfInfo& pdf,
Units::MomentumUnit mom,
Units::LengthUnit len ) :
m_signal_process_id(signal_process_id),
m_event_number(event_number),
m_mpi(-1),
m_event_scale(-1),
m_alphaQCD(-1),
m_alphaQED(-1),
m_signal_process_vertex(signal_vertex),
m_beam_particle_1(0),
m_beam_particle_2(0),
m_weights(weights),
m_random_states(random_states),
m_vertex_barcodes(),
m_particle_barcodes(),
m_cross_section(0),
m_heavy_ion( new HeavyIon(ion) ),
m_pdf_info( new PdfInfo(pdf) ),
m_momentum_unit(mom),
m_position_unit(len)
{
/// GenEvent makes its own copy of HeavyIon and PdfInfo
///
/// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
/// are as suggested in hep-ph/0109068, "Generic Interface..."
}
GenEvent::GenEvent( Units::MomentumUnit mom,
Units::LengthUnit len,
int signal_process_id,
int event_number,
GenVertex* signal_vertex,
const WeightContainer& weights,
const std::vector<long>& random_states ) :
m_signal_process_id(signal_process_id),
m_event_number(event_number),
m_mpi(-1),
m_event_scale(-1),
m_alphaQCD(-1),
m_alphaQED(-1),
m_signal_process_vertex(signal_vertex),
m_beam_particle_1(0),
m_beam_particle_2(0),
m_weights(weights),
m_random_states(random_states),
m_vertex_barcodes(),
m_particle_barcodes(),
m_cross_section(0),
m_heavy_ion(0),
m_pdf_info(0),
m_momentum_unit(mom),
m_position_unit(len)
{
/// constructor requiring units - all else is default
/// This constructor only allows null pointers to HeavyIon and PdfInfo
///
/// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
/// are as suggested in hep-ph/0109068, "Generic Interface..."
///
}
GenEvent::GenEvent( Units::MomentumUnit mom,
Units::LengthUnit len,
int signal_process_id, int event_number,
GenVertex* signal_vertex,
const WeightContainer& weights,
const std::vector<long>& random_states,
const HeavyIon& ion,
const PdfInfo& pdf ) :
m_signal_process_id(signal_process_id),
m_event_number(event_number),
m_mpi(-1),
m_event_scale(-1),
m_alphaQCD(-1),
m_alphaQED(-1),
m_signal_process_vertex(signal_vertex),
m_beam_particle_1(0),
m_beam_particle_2(0),
m_weights(weights),
m_random_states(random_states),
m_vertex_barcodes(),
m_particle_barcodes(),
m_cross_section(0),
m_heavy_ion( new HeavyIon(ion) ),
m_pdf_info( new PdfInfo(pdf) ),
m_momentum_unit(mom),
m_position_unit(len)
{
/// explicit constructor with units first that takes HeavyIon and PdfInfo
/// GenEvent makes its own copy of HeavyIon and PdfInfo
///
/// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
/// are as suggested in hep-ph/0109068, "Generic Interface..."
}
GenEvent::GenEvent( const GenEvent& inevent )
: m_signal_process_id ( inevent.signal_process_id() ),
m_event_number ( inevent.event_number() ),
m_mpi ( inevent.mpi() ),
m_event_scale ( inevent.event_scale() ),
m_alphaQCD ( inevent.alphaQCD() ),
m_alphaQED ( inevent.alphaQED() ),
m_signal_process_vertex( /* inevent.m_signal_process_vertex */ ),
m_beam_particle_1 ( /* inevent.m_beam_particle_1 */ ),
m_beam_particle_2 ( /* inevent.m_beam_particle_2 */ ),
m_weights ( /* inevent.m_weights */ ),
m_random_states ( /* inevent.m_random_states */ ),
m_vertex_barcodes ( /* inevent.m_vertex_barcodes */ ),
m_particle_barcodes ( /* inevent.m_particle_barcodes */ ),
m_cross_section ( inevent.cross_section() ? new GenCrossSection(*inevent.cross_section()) : 0 ),
m_heavy_ion ( inevent.heavy_ion() ? new HeavyIon(*inevent.heavy_ion()) : 0 ),
m_pdf_info ( inevent.pdf_info() ? new PdfInfo(*inevent.pdf_info()) : 0 ),
m_momentum_unit ( inevent.momentum_unit() ),
m_position_unit ( inevent.length_unit() )
{
/// deep copy - makes a copy of all vertices!
//
// 1. create a NEW copy of all vertices from inevent
// taking care to map new vertices onto the vertices being copied
// and add these new vertices to this event.
// We do not use GenVertex::operator= because that would copy
// the attached particles as well.
std::map<const GenVertex*,GenVertex*> map_in_to_new;
for ( GenEvent::vertex_const_iterator v = inevent.vertices_begin();
v != inevent.vertices_end(); ++v ) {
GenVertex* newvertex = new GenVertex(
(*v)->position(), (*v)->id(), (*v)->weights() );
newvertex->suggest_barcode( (*v)->barcode() );
map_in_to_new[*v] = newvertex;
add_vertex( newvertex );
}
// 2. copy the signal process vertex info.
if ( inevent.signal_process_vertex() ) {
set_signal_process_vertex(
map_in_to_new[inevent.signal_process_vertex()] );
} else set_signal_process_vertex( 0 );
//
// 3. create a NEW copy of all particles from inevent
// taking care to attach them to the appropriate vertex
GenParticle* beam1(0);
GenParticle* beam2(0);
for ( GenEvent::particle_const_iterator p = inevent.particles_begin();
p != inevent.particles_end(); ++p )
{
GenParticle* oldparticle = *p;
GenParticle* newparticle = new GenParticle(*oldparticle);
if ( oldparticle->end_vertex() ) {
map_in_to_new[ oldparticle->end_vertex() ]->
add_particle_in(newparticle);
}
if ( oldparticle->production_vertex() ) {
map_in_to_new[ oldparticle->production_vertex() ]->
add_particle_out(newparticle);
}
if ( oldparticle == inevent.beam_particles().first ) beam1 = newparticle;
if ( oldparticle == inevent.beam_particles().second ) beam2 = newparticle;
}
set_beam_particles( beam1, beam2 );
//
// 4. now that vtx/particles are copied, copy weights and random states
set_random_states( inevent.random_states() );
weights() = inevent.weights();
}
void GenEvent::swap( GenEvent & other )
{
// if a container has a swap method, use that for improved performance
std::swap(m_signal_process_id , other.m_signal_process_id );
std::swap(m_event_number , other.m_event_number );
std::swap(m_mpi , other.m_mpi );
std::swap(m_event_scale , other.m_event_scale );
std::swap(m_alphaQCD , other.m_alphaQCD );
std::swap(m_alphaQED , other.m_alphaQED );
std::swap(m_signal_process_vertex, other.m_signal_process_vertex);
std::swap(m_beam_particle_1 , other.m_beam_particle_1 );
std::swap(m_beam_particle_2 , other.m_beam_particle_2 );
m_weights.swap( other.m_weights );
m_random_states.swap( other.m_random_states );
m_vertex_barcodes.swap( other.m_vertex_barcodes );
m_particle_barcodes.swap( other.m_particle_barcodes );
std::swap(m_cross_section , other.m_cross_section );
std::swap(m_heavy_ion , other.m_heavy_ion );
std::swap(m_pdf_info , other.m_pdf_info );
std::swap(m_momentum_unit , other.m_momentum_unit );
std::swap(m_position_unit , other.m_position_unit );
// must now adjust GenVertex back pointers
for ( GenEvent::vertex_const_iterator vthis = vertices_begin();
vthis != vertices_end(); ++vthis ) {
(*vthis)->change_parent_event_( this );
}
for ( GenEvent::vertex_const_iterator voth = other.vertices_begin();
voth != other.vertices_end(); ++voth ) {
(*voth)->change_parent_event_( &other );
}
}
GenEvent::~GenEvent()
{
/// Deep destructor.
/// deletes all vertices/particles in this GenEvent
/// deletes the associated HeavyIon and PdfInfo
delete_all_vertices();
delete m_cross_section;
delete m_heavy_ion;
delete m_pdf_info;
}
GenEvent& GenEvent::operator=( const GenEvent& inevent )
{
/// best practices implementation
GenEvent tmp( inevent );
swap( tmp );
return *this;
}
void GenEvent::print( std::ostream& ostr ) const {
/// dumps the content of this event to ostr
/// to dump to cout use: event.print();
/// if you want to write this event to file outfile.txt you could use:
/// std::ofstream outfile("outfile.txt"); event.print( outfile );
ostr << "________________________________________"
<< "________________________________________\n";
ostr << "GenEvent: #" << event_number()
<< " ID=" << signal_process_id()
<< " SignalProcessGenVertex Barcode: "
<< ( signal_process_vertex() ? signal_process_vertex()->barcode()
: 0 )
<< "\n";
write_units( ostr );
write_cross_section(ostr);
ostr << " Entries this event: " << vertices_size() << " vertices, "
<< particles_size() << " particles.\n";
if( m_beam_particle_1 && m_beam_particle_2 ) {
ostr << " Beam Particle barcodes: " << beam_particles().first->barcode() << " "
<< beam_particles().second->barcode() << " \n";
} else {
ostr << " Beam Particles are not defined.\n";
}
// Random State
ostr << " RndmState(" << m_random_states.size() << ")=";
for ( std::vector<long>::const_iterator rs
= m_random_states.begin();
rs != m_random_states.end(); ++rs ) { ostr << *rs << " "; }
ostr << "\n";
// Weights
ostr << " Wgts(" << weights().size() << ")=";
weights().print(ostr);
ostr << " EventScale " << event_scale()
<< " [energy] \t alphaQCD=" << alphaQCD()
<< "\t alphaQED=" << alphaQED() << std::endl;
// print a legend to describe the particle info
ostr << " GenParticle Legend\n";
ostr << " Barcode PDG ID "
<< "( Px, Py, Pz, E )"
<< " Stat DecayVtx\n";
ostr << "________________________________________"
<< "________________________________________\n";
// Print all Vertices
for ( GenEvent::vertex_const_iterator vtx = this->vertices_begin();
vtx != this->vertices_end(); ++vtx ) {
(*vtx)->print(ostr);
}
ostr << "________________________________________"
<< "________________________________________" << std::endl;
}
void GenEvent::print_version( std::ostream& ostr ) const {
ostr << "---------------------------------------------" << std::endl;
writeVersion( ostr );
ostr << "---------------------------------------------" << std::endl;
}
bool GenEvent::add_vertex( GenVertex* vtx ) {
/// returns true if successful - generally will only return false
/// if the inserted vertex is already included in the event.
if ( !vtx ) return false;
// if vtx previously pointed to another GenEvent, remove it from that
// GenEvent's list
if ( vtx->parent_event() && vtx->parent_event() != this ) {
bool remove_status = vtx->parent_event()->remove_vertex( vtx );
if ( !remove_status ) {
std::cerr << "GenEvent::add_vertex ERROR "
<< "GenVertex::parent_event points to \n"
<< "an event that does not point back to the "
<< "GenVertex. \n This probably indicates a deeper "
<< "problem. " << std::endl;
}
}
//
// setting the vertex parent also inserts the vertex into this
// event
vtx->set_parent_event_( this );
return ( m_vertex_barcodes.count(vtx->barcode()) ? true : false );
}
bool GenEvent::remove_vertex( GenVertex* vtx ) {
/// this removes vtx from the event but does NOT delete it.
/// returns True if an entry vtx existed in the table and was erased
if ( m_signal_process_vertex == vtx ) m_signal_process_vertex = 0;
if ( vtx->parent_event() == this ) vtx->set_parent_event_( 0 );
return ( m_vertex_barcodes.count(vtx->barcode()) ? false : true );
}
void GenEvent::clear()
{
/// remove all information from the event
/// deletes all vertices/particles in this evt
///
delete_all_vertices();
// remove existing objects and set pointers to null
delete m_cross_section;
m_cross_section = 0;
delete m_heavy_ion;
m_heavy_ion = 0;
delete m_pdf_info;
m_pdf_info = 0;
m_signal_process_id = 0;
m_beam_particle_1 = 0;
m_beam_particle_2 = 0;
m_event_number = 0;
m_mpi = -1;
m_event_scale = -1;
m_alphaQCD = -1;
m_alphaQED = -1;
m_weights = std::vector<double>();
m_random_states = std::vector<long>();
// resetting unit information
m_momentum_unit = Units::default_momentum_unit();
m_position_unit = Units::default_length_unit();
// error check just to be safe
if ( m_vertex_barcodes.size() != 0
|| m_particle_barcodes.size() != 0 ) {
std::cerr << "GenEvent::clear() strange result ... \n"
<< "either the particle and/or the vertex map isn't empty" << std::endl;
std::cerr << "Number vtx,particle the event after deleting = "
<< m_vertex_barcodes.size() << " "
<< m_particle_barcodes.size() << std::endl;
}
return;
}
void GenEvent::delete_all_vertices() {
/// deletes all vertices in the vertex container
/// (i.e. all vertices owned by this event)
/// The vertices are the "owners" of the particles, so as we delete
/// the vertices, the vertex desctructors are automatically
/// deleting their particles.
// delete each vertex individually (this deletes particles as well)
while ( !vertices_empty() ) {
GenVertex* vtx = ( m_vertex_barcodes.begin() )->second;
m_vertex_barcodes.erase( m_vertex_barcodes.begin() );
delete vtx;
}
//
// Error checking:
if ( !vertices_empty() || ! particles_empty() ) {
std::cerr << "GenEvent::delete_all_vertices strange result ... "
<< "after deleting all vertices, \nthe particle and "
<< "vertex maps aren't empty.\n This probably "
<< "indicates deeper problems or memory leak in the "
<< "code." << std::endl;
std::cerr << "Number vtx,particle the event after deleting = "
<< m_vertex_barcodes.size() << " "
<< m_particle_barcodes.size() << std::endl;
}
}
bool GenEvent::set_barcode( GenParticle* p, int suggested_barcode )
{
if ( p->parent_event() != this ) {
std::cerr << "GenEvent::set_barcode attempted, but the argument's"
<< "\n parent_event is not this ... request rejected."
<< std::endl;
return false;
}
// M.Dobbs Nov 4, 2002
// First we must check to see if the particle already has a
// barcode which is different from the suggestion. If yes, we
// remove it from the particle map.
if ( p->barcode() != 0 && p->barcode() != suggested_barcode ) {
if ( m_particle_barcodes.count(p->barcode()) &&
m_particle_barcodes[p->barcode()] == p ) {
m_particle_barcodes.erase( p->barcode() );
}
// At this point either the particle is NOT in
// m_particle_barcodes, or else it is in the map, but
// already with the suggested barcode.
}
//
// First case --- a valid barcode has been suggested
// (valid barcodes are numbers greater than zero)
bool insert_success = true;
if ( suggested_barcode > 0 ) {
if ( m_particle_barcodes.count(suggested_barcode) ) {
// the suggested_barcode is already used.
if ( m_particle_barcodes[suggested_barcode] == p ) {
// but it was used for this particle ... so everythings ok
p->set_barcode_( suggested_barcode );
return true;
}
insert_success = false;
suggested_barcode = 0;
} else { // suggested barcode is OK, proceed to insert
m_particle_barcodes[suggested_barcode] = p;
p->set_barcode_( suggested_barcode );
return true;
}
}
//
// Other possibility -- a valid barcode has not been suggested,
// so one is automatically generated
if ( suggested_barcode < 0 ) insert_success = false;
if ( suggested_barcode <= 0 ) {
if ( !m_particle_barcodes.empty() ) {
// in this case we find the highest barcode that was used,
// and increment it by 1
suggested_barcode = m_particle_barcodes.rbegin()->first;
++suggested_barcode;
}
// For the automatically assigned barcodes, the first one
// we use is 10001 ... this is just because when we read
// events from HEPEVT, we will suggest their index as the
// barcode, and that index has maximum value 10000.
// This way we will immediately be able to recognize the hepevt
// particles from the auto-assigned ones.
if ( suggested_barcode <= 10000 ) suggested_barcode = 10001;
}
// At this point we should have a valid barcode
if ( m_particle_barcodes.count(suggested_barcode) ) {
std::cerr << "GenEvent::set_barcode ERROR, this should never "
<< "happen \n report bug to matt.dobbs@cern.ch"
<< std::endl;
}
m_particle_barcodes[suggested_barcode] = p;
p->set_barcode_( suggested_barcode );
return insert_success;
}
bool GenEvent::set_barcode( GenVertex* v, int suggested_barcode )
{
if ( v->parent_event() != this ) {
std::cerr << "GenEvent::set_barcode attempted, but the argument's"
<< "\n parent_event is not this ... request rejected."
<< std::endl;
return false;
}
// M.Dobbs Nov 4, 2002
// First we must check to see if the vertex already has a
// barcode which is different from the suggestion. If yes, we
// remove it from the vertex map.
if ( v->barcode() != 0 && v->barcode() != suggested_barcode ) {
if ( m_vertex_barcodes.count(v->barcode()) &&
m_vertex_barcodes[v->barcode()] == v ) {
m_vertex_barcodes.erase( v->barcode() );
}
// At this point either the vertex is NOT in
// m_vertex_barcodes, or else it is in the map, but
// already with the suggested barcode.
}
//
// First case --- a valid barcode has been suggested
// (valid barcodes are numbers greater than zero)
bool insert_success = true;
if ( suggested_barcode < 0 ) {
if ( m_vertex_barcodes.count(suggested_barcode) ) {
// the suggested_barcode is already used.
if ( m_vertex_barcodes[suggested_barcode] == v ) {
// but it was used for this vertex ... so everythings ok
v->set_barcode_( suggested_barcode );
return true;
}
insert_success = false;
suggested_barcode = 0;
} else { // suggested barcode is OK, proceed to insert
m_vertex_barcodes[suggested_barcode] = v;
v->set_barcode_( suggested_barcode );
return true;
}
}
//
// Other possibility -- a valid barcode has not been suggested,
// so one is automatically generated
if ( suggested_barcode > 0 ) insert_success = false;
if ( suggested_barcode >= 0 ) {
if ( !m_vertex_barcodes.empty() ) {
// in this case we find the highest barcode that was used,
// and increment it by 1, (vertex barcodes are negative)
suggested_barcode = m_vertex_barcodes.rbegin()->first;
--suggested_barcode;
}
if ( suggested_barcode >= 0 ) suggested_barcode = -1;
}
// At this point we should have a valid barcode
if ( m_vertex_barcodes.count(suggested_barcode) ) {
std::cerr << "GenEvent::set_barcode ERROR, this should never "
<< "happen \n report bug to matt.dobbs@cern.ch"
<< std::endl;
}
m_vertex_barcodes[suggested_barcode] = v;
v->set_barcode_( suggested_barcode );
return insert_success;
}
/// test to see if we have two valid beam particles
bool GenEvent::valid_beam_particles() const {
bool have1 = false;
bool have2 = false;
// first check that both are defined
if(m_beam_particle_1 && m_beam_particle_2) {
// now look for a match with the particle "list"
for ( particle_const_iterator p = particles_begin();
p != particles_end(); ++p ) {
if( m_beam_particle_1 == *p ) have1 = true;
if( m_beam_particle_2 == *p ) have2 = true;
}
}
if( have1 && have2 ) return true;
return false;
}
/// construct the beam particle information using pointers to GenParticle
/// returns false if either GenParticle* is null
bool GenEvent::set_beam_particles(GenParticle* bp1, GenParticle* bp2) {
m_beam_particle_1 = bp1;
m_beam_particle_2 = bp2;
if( m_beam_particle_1 && m_beam_particle_2 ) return true;
return false;
}
/// construct the beam particle information using a std::pair of pointers to GenParticle
/// returns false if either GenParticle* is null
bool GenEvent::set_beam_particles(std::pair<HepMC::GenParticle*, HepMC::GenParticle*> const & bp) {
return set_beam_particles(bp.first,bp.second);
}
void GenEvent::write_units( std::ostream & os ) const {
os << " Momenutm units:" << std::setw(8) << name(momentum_unit());
os << " Position units:" << std::setw(8) << name(length_unit());
os << std::endl;
}
void GenEvent::write_cross_section( std::ostream& os ) const
{
// write the GenCrossSection information if the cross section was set
if( !cross_section() ) return;
if( cross_section()->is_set() ) {
os << " Cross Section: " << cross_section()->cross_section() ;
os << " +/- " << cross_section()->cross_section_error() ;
os << std::endl;
}
}
bool GenEvent::use_momentum_unit( Units::MomentumUnit newunit ) {
// currently not exception-safe.
// Easy to fix, though, if needed.
if ( m_momentum_unit != newunit ) {
const double factor = Units::conversion_factor( m_momentum_unit, newunit );
// multiply all momenta by 'factor',
// loop is entered only if particle list is not empty
for ( GenEvent::particle_iterator p = particles_begin();
p != particles_end(); ++p )
{
(*p)->convert_momentum(factor);
}
// ...
m_momentum_unit = newunit;
}
return true;
}
bool GenEvent::use_length_unit( Units::LengthUnit newunit ) {
// currently not exception-safe.
// Easy to fix, though, if needed.
if ( m_position_unit != newunit ) {
const double factor = Units::conversion_factor( m_position_unit, newunit );
// multiply all lengths by 'factor',
// loop is entered only if vertex list is not empty
for ( GenEvent::vertex_iterator vtx = vertices_begin();
vtx != vertices_end(); ++vtx ) {
(*vtx)->convert_position(factor);
}
// ...
m_position_unit = newunit;
}
return true;
}
bool GenEvent::use_momentum_unit( std::string& newunit ) {
if ( newunit == "MEV" ) return use_momentum_unit( Units::MEV );
else if( newunit == "GEV" ) return use_momentum_unit( Units::GEV );
else std::cerr << "GenEvent::use_momentum_unit ERROR: use either MEV or GEV\n";
return false;
}
bool GenEvent::use_length_unit( std::string& newunit ) {
if ( newunit == "MM" ) return use_length_unit( Units::MM );
else if( newunit == "CM" ) return use_length_unit( Units::CM );
else std::cerr << "GenEvent::use_length_unit ERROR: use either MEV or GEV\n";
return false;
}
bool GenEvent::is_valid() const {
/// A GenEvent is presumed valid if it has both associated
/// particles and vertices. No other information is checked.
if ( vertices_empty() ) return false;
if ( particles_empty() ) return false;
return true;
}
std::ostream & GenEvent::write_beam_particles(std::ostream & os,
std::pair<HepMC::GenParticle *,HepMC::GenParticle *> pr )
{
GenParticle* p = pr.first;
if(!p) {
detail::output( os, 0 );
} else {
detail::output( os, p->barcode() );
}
p = pr.second;
if(!p) {
detail::output( os, 0 );
} else {
detail::output( os, p->barcode() );
}
return os;
}
std::ostream & GenEvent::write_vertex(std::ostream & os, GenVertex const * v)
{
if ( !v || !os ) {
std::cerr << "GenEvent::write_vertex !v||!os, "
<< "v="<< v << " setting badbit" << std::endl;
os.clear(std::ios::badbit);
return os;
}
// First collect info we need
// count the number of orphan particles going into v
int num_orphans_in = 0;
for ( GenVertex::particles_in_const_iterator p1
= v->particles_in_const_begin();
p1 != v->particles_in_const_end(); ++p1 ) {
if ( !(*p1)->production_vertex() ) ++num_orphans_in;
}
//
os << 'V';
detail::output( os, v->barcode() ); // v's unique identifier
detail::output( os, v->id() );
detail::output( os, v->position().x() );
detail::output( os, v->position().y() );
detail::output( os, v->position().z() );
detail::output( os, v->position().t() );
detail::output( os, num_orphans_in );
detail::output( os, (int)v->particles_out_size() );
detail::output( os, (int)v->weights().size() );
for ( WeightContainer::const_iterator w = v->weights().begin();
w != v->weights().end(); ++w ) {
detail::output( os, *w );
}
detail::output( os,'\n');
// incoming particles
for ( GenVertex::particles_in_const_iterator p2
= v->particles_in_const_begin();
p2 != v->particles_in_const_end(); ++p2 ) {
if ( !(*p2)->production_vertex() ) {
write_particle( os, *p2 );
}
}
// outgoing particles
for ( GenVertex::particles_out_const_iterator p3
= v->particles_out_const_begin();
p3 != v->particles_out_const_end(); ++p3 ) {
write_particle( os, *p3 );
}
return os;
}
std::ostream & GenEvent::write_particle( std::ostream & os, GenParticle const * p )
{
if ( !p || !os ) {
std::cerr << "GenEvent::write_particle !p||!os, "
<< "p="<< p << " setting badbit" << std::endl;
os.clear(std::ios::badbit);
return os;
}
os << 'P';
detail::output( os, p->barcode() );
detail::output( os, p->pdg_id() );
detail::output( os, p->momentum().px() );
detail::output( os, p->momentum().py() );
detail::output( os, p->momentum().pz() );
detail::output( os, p->momentum().e() );
detail::output( os, p->generated_mass() );
detail::output( os, p->status() );
detail::output( os, p->polarization().theta() );
detail::output( os, p->polarization().phi() );
// since end_vertex is oftentimes null, this CREATES a null vertex
// in the map
detail::output( os, ( p->end_vertex() ? p->end_vertex()->barcode() : 0 ) );
os << ' ' << p->flow() << "\n";
return os;
}
} // HepMC
// Matt.Dobbs@Cern.CH, September 1999
// Updated: 7.1.2000 iterators complete and working!
// Updated: 10.1.2000 GenEvent::vertex, particle iterators are made
// constant WRT this event ... note that
// GenVertex::***_iterator is not const, since it must
// be able to return a mutable pointer to itself.
// Updated: 08.2.2000 the event now holds a set of all attached vertices
// rather than just the roots of the graph
// Event record for MC generators (for use at any stage of generation)
//////////////////////////////////////////////////////////////////////////
#include <iomanip>
#include "HepMC/GenEvent.h"
#include "HepMC/GenCrossSection.h"
#include "HepMC/Version.h"
#include "HepMC/StreamHelpers.h"
namespace HepMC {
GenEvent::GenEvent( int signal_process_id,
int event_number,
GenVertex* signal_vertex,
const WeightContainer& weights,
const std::vector<long>& random_states,
Units::MomentumUnit mom,
Units::LengthUnit len ) :
m_signal_process_id(signal_process_id),
m_event_number(event_number),
m_mpi(-1),
m_event_scale(-1),
m_alphaQCD(-1),
m_alphaQED(-1),
m_signal_process_vertex(signal_vertex),
m_beam_particle_1(0),
m_beam_particle_2(0),
m_weights(weights),
m_random_states(random_states),
m_vertex_barcodes(),
m_particle_barcodes(),
m_cross_section(0),
m_heavy_ion(0),
m_pdf_info(0),
m_momentum_unit(mom),
m_position_unit(len)
{
/// This constructor only allows null pointers to HeavyIon and PdfInfo
///
/// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
/// are as suggested in hep-ph/0109068, "Generic Interface..."
///
}
GenEvent::GenEvent( int signal_process_id, int event_number,
GenVertex* signal_vertex,
const WeightContainer& weights,
const std::vector<long>& random_states,
const HeavyIon& ion,
const PdfInfo& pdf,
Units::MomentumUnit mom,
Units::LengthUnit len ) :
m_signal_process_id(signal_process_id),
m_event_number(event_number),
m_mpi(-1),
m_event_scale(-1),
m_alphaQCD(-1),
m_alphaQED(-1),
m_signal_process_vertex(signal_vertex),
m_beam_particle_1(0),
m_beam_particle_2(0),
m_weights(weights),
m_random_states(random_states),
m_vertex_barcodes(),
m_particle_barcodes(),
m_cross_section(0),
m_heavy_ion( new HeavyIon(ion) ),
m_pdf_info( new PdfInfo(pdf) ),
m_momentum_unit(mom),
m_position_unit(len)
{
/// GenEvent makes its own copy of HeavyIon and PdfInfo
///
/// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
/// are as suggested in hep-ph/0109068, "Generic Interface..."
}
GenEvent::GenEvent( Units::MomentumUnit mom,
Units::LengthUnit len,
int signal_process_id,
int event_number,
GenVertex* signal_vertex,
const WeightContainer& weights,
const std::vector<long>& random_states ) :
m_signal_process_id(signal_process_id),
m_event_number(event_number),
m_mpi(-1),
m_event_scale(-1),
m_alphaQCD(-1),
m_alphaQED(-1),
m_signal_process_vertex(signal_vertex),
m_beam_particle_1(0),
m_beam_particle_2(0),
m_weights(weights),
m_random_states(random_states),
m_vertex_barcodes(),
m_particle_barcodes(),
m_cross_section(0),
m_heavy_ion(0),
m_pdf_info(0),
m_momentum_unit(mom),
m_position_unit(len)
{
/// constructor requiring units - all else is default
/// This constructor only allows null pointers to HeavyIon and PdfInfo
///
/// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
/// are as suggested in hep-ph/0109068, "Generic Interface..."
///
}
GenEvent::GenEvent( Units::MomentumUnit mom,
Units::LengthUnit len,
int signal_process_id, int event_number,
GenVertex* signal_vertex,
const WeightContainer& weights,
const std::vector<long>& random_states,
const HeavyIon& ion,
const PdfInfo& pdf ) :
m_signal_process_id(signal_process_id),
m_event_number(event_number),
m_mpi(-1),
m_event_scale(-1),
m_alphaQCD(-1),
m_alphaQED(-1),
m_signal_process_vertex(signal_vertex),
m_beam_particle_1(0),
m_beam_particle_2(0),
m_weights(weights),
m_random_states(random_states),
m_vertex_barcodes(),
m_particle_barcodes(),
m_cross_section(0),
m_heavy_ion( new HeavyIon(ion) ),
m_pdf_info( new PdfInfo(pdf) ),
m_momentum_unit(mom),
m_position_unit(len)
{
/// explicit constructor with units first that takes HeavyIon and PdfInfo
/// GenEvent makes its own copy of HeavyIon and PdfInfo
///
/// note: default values for m_event_scale, m_alphaQCD, m_alphaQED
/// are as suggested in hep-ph/0109068, "Generic Interface..."
}
GenEvent::GenEvent( const GenEvent& inevent )
: m_signal_process_id ( inevent.signal_process_id() ),
m_event_number ( inevent.event_number() ),
m_mpi ( inevent.mpi() ),
m_event_scale ( inevent.event_scale() ),
m_alphaQCD ( inevent.alphaQCD() ),
m_alphaQED ( inevent.alphaQED() ),
m_signal_process_vertex( /* inevent.m_signal_process_vertex */ ),
m_beam_particle_1 ( /* inevent.m_beam_particle_1 */ ),
m_beam_particle_2 ( /* inevent.m_beam_particle_2 */ ),
m_weights ( /* inevent.m_weights */ ),
m_random_states ( /* inevent.m_random_states */ ),
m_vertex_barcodes ( /* inevent.m_vertex_barcodes */ ),
m_particle_barcodes ( /* inevent.m_particle_barcodes */ ),
m_cross_section ( inevent.cross_section() ? new GenCrossSection(*inevent.cross_section()) : 0 ),
m_heavy_ion ( inevent.heavy_ion() ? new HeavyIon(*inevent.heavy_ion()) : 0 ),
m_pdf_info ( inevent.pdf_info() ? new PdfInfo(*inevent.pdf_info()) : 0 ),
m_momentum_unit ( inevent.momentum_unit() ),
m_position_unit ( inevent.length_unit() )
{
/// deep copy - makes a copy of all vertices!
//
// 1. create a NEW copy of all vertices from inevent
// taking care to map new vertices onto the vertices being copied
// and add these new vertices to this event.
// We do not use GenVertex::operator= because that would copy
// the attached particles as well.
std::map<const GenVertex*,GenVertex*> map_in_to_new;
for ( GenEvent::vertex_const_iterator v = inevent.vertices_begin();
v != inevent.vertices_end(); ++v ) {
GenVertex* newvertex = new GenVertex(
(*v)->position(), (*v)->id(), (*v)->weights() );
newvertex->suggest_barcode( (*v)->barcode() );
map_in_to_new[*v] = newvertex;
add_vertex( newvertex );
}
// 2. copy the signal process vertex info.
if ( inevent.signal_process_vertex() ) {
set_signal_process_vertex(
map_in_to_new[inevent.signal_process_vertex()] );
} else set_signal_process_vertex( 0 );
//
// 3. create a NEW copy of all particles from inevent
// taking care to attach them to the appropriate vertex
GenParticle* beam1(0);
GenParticle* beam2(0);
for ( GenEvent::particle_const_iterator p = inevent.particles_begin();
p != inevent.particles_end(); ++p )
{
GenParticle* oldparticle = *p;
GenParticle* newparticle = new GenParticle(*oldparticle);
if ( oldparticle->end_vertex() ) {
map_in_to_new[ oldparticle->end_vertex() ]->
add_particle_in(newparticle);
}
if ( oldparticle->production_vertex() ) {
map_in_to_new[ oldparticle->production_vertex() ]->
add_particle_out(newparticle);
}
if ( oldparticle == inevent.beam_particles().first ) beam1 = newparticle;
if ( oldparticle == inevent.beam_particles().second ) beam2 = newparticle;
}
set_beam_particles( beam1, beam2 );
//
// 4. now that vtx/particles are copied, copy weights and random states
set_random_states( inevent.random_states() );
weights() = inevent.weights();
}
void GenEvent::swap( GenEvent & other )
{
// if a container has a swap method, use that for improved performance
std::swap(m_signal_process_id , other.m_signal_process_id );
std::swap(m_event_number , other.m_event_number );
std::swap(m_mpi , other.m_mpi );
std::swap(m_event_scale , other.m_event_scale );
std::swap(m_alphaQCD , other.m_alphaQCD );
std::swap(m_alphaQED , other.m_alphaQED );
std::swap(m_signal_process_vertex, other.m_signal_process_vertex);
std::swap(m_beam_particle_1 , other.m_beam_particle_1 );
std::swap(m_beam_particle_2 , other.m_beam_particle_2 );
m_weights.swap( other.m_weights );
m_random_states.swap( other.m_random_states );
m_vertex_barcodes.swap( other.m_vertex_barcodes );
m_particle_barcodes.swap( other.m_particle_barcodes );
std::swap(m_cross_section , other.m_cross_section );
std::swap(m_heavy_ion , other.m_heavy_ion );
std::swap(m_pdf_info , other.m_pdf_info );
std::swap(m_momentum_unit , other.m_momentum_unit );
std::swap(m_position_unit , other.m_position_unit );
// must now adjust GenVertex back pointers
for ( GenEvent::vertex_const_iterator vthis = vertices_begin();
vthis != vertices_end(); ++vthis ) {
(*vthis)->change_parent_event_( this );
}
for ( GenEvent::vertex_const_iterator voth = other.vertices_begin();
voth != other.vertices_end(); ++voth ) {
(*voth)->change_parent_event_( &other );
}
}
GenEvent::~GenEvent()
{
/// Deep destructor.
/// deletes all vertices/particles in this GenEvent
/// deletes the associated HeavyIon and PdfInfo
delete_all_vertices();
delete m_cross_section;
delete m_heavy_ion;
delete m_pdf_info;
}
GenEvent& GenEvent::operator=( const GenEvent& inevent )
{
/// best practices implementation
GenEvent tmp( inevent );
swap( tmp );
return *this;
}
void GenEvent::print( std::ostream& ostr ) const {
/// dumps the content of this event to ostr
/// to dump to cout use: event.print();
/// if you want to write this event to file outfile.txt you could use:
/// std::ofstream outfile("outfile.txt"); event.print( outfile );
ostr << "________________________________________"
<< "________________________________________\n";
ostr << "GenEvent: #" << event_number()
<< " ID=" << signal_process_id()
<< " SignalProcessGenVertex Barcode: "
<< ( signal_process_vertex() ? signal_process_vertex()->barcode()
: 0 )
<< "\n";
write_units( ostr );
write_cross_section(ostr);
ostr << " Entries this event: " << vertices_size() << " vertices, "
<< particles_size() << " particles.\n";
if( m_beam_particle_1 && m_beam_particle_2 ) {
ostr << " Beam Particle barcodes: " << beam_particles().first->barcode() << " "
<< beam_particles().second->barcode() << " \n";
} else {
ostr << " Beam Particles are not defined.\n";
}
// Random State
ostr << " RndmState(" << m_random_states.size() << ")=";
for ( std::vector<long>::const_iterator rs
= m_random_states.begin();
rs != m_random_states.end(); ++rs ) { ostr << *rs << " "; }
ostr << "\n";
// Weights
ostr << " Wgts(" << weights().size() << ")=";
weights().print(ostr);
ostr << " EventScale " << event_scale()
<< " [energy] \t alphaQCD=" << alphaQCD()
<< "\t alphaQED=" << alphaQED() << std::endl;
// print a legend to describe the particle info
ostr << " GenParticle Legend\n";
ostr << " Barcode PDG ID "
<< "( Px, Py, Pz, E )"
<< " Stat DecayVtx\n";
ostr << "________________________________________"
<< "________________________________________\n";
// Print all Vertices
for ( GenEvent::vertex_const_iterator vtx = this->vertices_begin();
vtx != this->vertices_end(); ++vtx ) {
(*vtx)->print(ostr);
}
ostr << "________________________________________"
<< "________________________________________" << std::endl;
}
void GenEvent::print_version( std::ostream& ostr ) const {
ostr << "---------------------------------------------" << std::endl;
writeVersion( ostr );
ostr << "---------------------------------------------" << std::endl;
}
bool GenEvent::add_vertex( GenVertex* vtx ) {
/// returns true if successful - generally will only return false
/// if the inserted vertex is already included in the event.
if ( !vtx ) return false;
// if vtx previously pointed to another GenEvent, remove it from that
// GenEvent's list
if ( vtx->parent_event() && vtx->parent_event() != this ) {
bool remove_status = vtx->parent_event()->remove_vertex( vtx );
if ( !remove_status ) {
std::cerr << "GenEvent::add_vertex ERROR "
<< "GenVertex::parent_event points to \n"
<< "an event that does not point back to the "
<< "GenVertex. \n This probably indicates a deeper "
<< "problem. " << std::endl;
}
}
//
// setting the vertex parent also inserts the vertex into this
// event
vtx->set_parent_event_( this );
return ( m_vertex_barcodes.count(vtx->barcode()) ? true : false );
}
bool GenEvent::remove_vertex( GenVertex* vtx ) {
/// this removes vtx from the event but does NOT delete it.
/// returns True if an entry vtx existed in the table and was erased
if ( m_signal_process_vertex == vtx ) m_signal_process_vertex = 0;
if ( vtx->parent_event() == this ) vtx->set_parent_event_( 0 );
return ( m_vertex_barcodes.count(vtx->barcode()) ? false : true );
}
void GenEvent::clear()
{
/// remove all information from the event
/// deletes all vertices/particles in this evt
///
delete_all_vertices();
// remove existing objects and set pointers to null
delete m_cross_section;
m_cross_section = 0;
delete m_heavy_ion;
m_heavy_ion = 0;
delete m_pdf_info;
m_pdf_info = 0;
m_signal_process_id = 0;
m_beam_particle_1 = 0;
m_beam_particle_2 = 0;
m_event_number = 0;
m_mpi = -1;
m_event_scale = -1;
m_alphaQCD = -1;
m_alphaQED = -1;
m_weights = std::vector<double>();
m_random_states = std::vector<long>();
// resetting unit information
m_momentum_unit = Units::default_momentum_unit();
m_position_unit = Units::default_length_unit();
// error check just to be safe
if ( m_vertex_barcodes.size() != 0
|| m_particle_barcodes.size() != 0 ) {
std::cerr << "GenEvent::clear() strange result ... \n"
<< "either the particle and/or the vertex map isn't empty" << std::endl;
std::cerr << "Number vtx,particle the event after deleting = "
<< m_vertex_barcodes.size() << " "
<< m_particle_barcodes.size() << std::endl;
}
return;
}
void GenEvent::delete_all_vertices() {
/// deletes all vertices in the vertex container
/// (i.e. all vertices owned by this event)
/// The vertices are the "owners" of the particles, so as we delete
/// the vertices, the vertex desctructors are automatically
/// deleting their particles.
// delete each vertex individually (this deletes particles as well)
while ( !vertices_empty() ) {
GenVertex* vtx = ( m_vertex_barcodes.begin() )->second;
m_vertex_barcodes.erase( m_vertex_barcodes.begin() );
delete vtx;
}
//
// Error checking:
if ( !vertices_empty() || ! particles_empty() ) {
std::cerr << "GenEvent::delete_all_vertices strange result ... "
<< "after deleting all vertices, \nthe particle and "
<< "vertex maps aren't empty.\n This probably "
<< "indicates deeper problems or memory leak in the "
<< "code." << std::endl;
std::cerr << "Number vtx,particle the event after deleting = "
<< m_vertex_barcodes.size() << " "
<< m_particle_barcodes.size() << std::endl;
}
}
bool GenEvent::set_barcode( GenParticle* p, int suggested_barcode )
{
if ( p->parent_event() != this ) {
std::cerr << "GenEvent::set_barcode attempted, but the argument's"
<< "\n parent_event is not this ... request rejected."
<< std::endl;
return false;
}
// M.Dobbs Nov 4, 2002
// First we must check to see if the particle already has a
// barcode which is different from the suggestion. If yes, we
// remove it from the particle map.
if ( p->barcode() != 0 && p->barcode() != suggested_barcode ) {
if ( m_particle_barcodes.count(p->barcode()) &&
m_particle_barcodes[p->barcode()] == p ) {
m_particle_barcodes.erase( p->barcode() );
}
// At this point either the particle is NOT in
// m_particle_barcodes, or else it is in the map, but
// already with the suggested barcode.
}
//
// First case --- a valid barcode has been suggested
// (valid barcodes are numbers greater than zero)
bool insert_success = true;
if ( suggested_barcode > 0 ) {
if ( m_particle_barcodes.count(suggested_barcode) ) {
// the suggested_barcode is already used.
if ( m_particle_barcodes[suggested_barcode] == p ) {
// but it was used for this particle ... so everythings ok
p->set_barcode_( suggested_barcode );
return true;
}
insert_success = false;
suggested_barcode = 0;
} else { // suggested barcode is OK, proceed to insert
m_particle_barcodes[suggested_barcode] = p;
p->set_barcode_( suggested_barcode );
return true;
}
}
//
// Other possibility -- a valid barcode has not been suggested,
// so one is automatically generated
if ( suggested_barcode < 0 ) insert_success = false;
if ( suggested_barcode <= 0 ) {
if ( !m_particle_barcodes.empty() ) {
// in this case we find the highest barcode that was used,
// and increment it by 1
suggested_barcode = m_particle_barcodes.rbegin()->first;
++suggested_barcode;
}
// For the automatically assigned barcodes, the first one
// we use is 10001 ... this is just because when we read
// events from HEPEVT, we will suggest their index as the
// barcode, and that index has maximum value 10000.
// This way we will immediately be able to recognize the hepevt
// particles from the auto-assigned ones.
if ( suggested_barcode <= 10000 ) suggested_barcode = 10001;
}
// At this point we should have a valid barcode
if ( m_particle_barcodes.count(suggested_barcode) ) {
std::cerr << "GenEvent::set_barcode ERROR, this should never "
<< "happen \n report bug to matt.dobbs@cern.ch"
<< std::endl;
}
m_particle_barcodes[suggested_barcode] = p;
p->set_barcode_( suggested_barcode );
return insert_success;
}
bool GenEvent::set_barcode( GenVertex* v, int suggested_barcode )
{
if ( v->parent_event() != this ) {
std::cerr << "GenEvent::set_barcode attempted, but the argument's"
<< "\n parent_event is not this ... request rejected."
<< std::endl;
return false;
}
// M.Dobbs Nov 4, 2002
// First we must check to see if the vertex already has a
// barcode which is different from the suggestion. If yes, we
// remove it from the vertex map.
if ( v->barcode() != 0 && v->barcode() != suggested_barcode ) {
if ( m_vertex_barcodes.count(v->barcode()) &&
m_vertex_barcodes[v->barcode()] == v ) {
m_vertex_barcodes.erase( v->barcode() );
}
// At this point either the vertex is NOT in
// m_vertex_barcodes, or else it is in the map, but
// already with the suggested barcode.
}
//
// First case --- a valid barcode has been suggested
// (valid barcodes are numbers greater than zero)
bool insert_success = true;
if ( suggested_barcode < 0 ) {
if ( m_vertex_barcodes.count(suggested_barcode) ) {
// the suggested_barcode is already used.
if ( m_vertex_barcodes[suggested_barcode] == v ) {
// but it was used for this vertex ... so everythings ok
v->set_barcode_( suggested_barcode );
return true;
}
insert_success = false;
suggested_barcode = 0;
} else { // suggested barcode is OK, proceed to insert
m_vertex_barcodes[suggested_barcode] = v;
v->set_barcode_( suggested_barcode );
return true;
}
}
//
// Other possibility -- a valid barcode has not been suggested,
// so one is automatically generated
if ( suggested_barcode > 0 ) insert_success = false;
if ( suggested_barcode >= 0 ) {
if ( !m_vertex_barcodes.empty() ) {
// in this case we find the highest barcode that was used,
// and increment it by 1, (vertex barcodes are negative)
suggested_barcode = m_vertex_barcodes.rbegin()->first;
--suggested_barcode;
}
if ( suggested_barcode >= 0 ) suggested_barcode = -1;
}
// At this point we should have a valid barcode
if ( m_vertex_barcodes.count(suggested_barcode) ) {
std::cerr << "GenEvent::set_barcode ERROR, this should never "
<< "happen \n report bug to matt.dobbs@cern.ch"
<< std::endl;
}
m_vertex_barcodes[suggested_barcode] = v;
v->set_barcode_( suggested_barcode );
return insert_success;
}
/// test to see if we have two valid beam particles
bool GenEvent::valid_beam_particles() const {
bool have1 = false;
bool have2 = false;
// first check that both are defined
if(m_beam_particle_1 && m_beam_particle_2) {
// now look for a match with the particle "list"
for ( particle_const_iterator p = particles_begin();
p != particles_end(); ++p ) {
if( m_beam_particle_1 == *p ) have1 = true;
if( m_beam_particle_2 == *p ) have2 = true;
}
}
if( have1 && have2 ) return true;
return false;
}
/// construct the beam particle information using pointers to GenParticle
/// returns false if either GenParticle* is null
bool GenEvent::set_beam_particles(GenParticle* bp1, GenParticle* bp2) {
m_beam_particle_1 = bp1;
m_beam_particle_2 = bp2;
if( m_beam_particle_1 && m_beam_particle_2 ) return true;
return false;
}
/// construct the beam particle information using a std::pair of pointers to GenParticle
/// returns false if either GenParticle* is null
bool GenEvent::set_beam_particles(std::pair<HepMC::GenParticle*, HepMC::GenParticle*> const & bp) {
return set_beam_particles(bp.first,bp.second);
}
void GenEvent::write_units( std::ostream & os ) const {
os << " Momenutm units:" << std::setw(8) << name(momentum_unit());
os << " Position units:" << std::setw(8) << name(length_unit());
os << std::endl;
}
void GenEvent::write_cross_section( std::ostream& os ) const
{
// write the GenCrossSection information if the cross section was set
if( !cross_section() ) return;
if( cross_section()->is_set() ) {
os << " Cross Section: " << cross_section()->cross_section() ;
os << " +/- " << cross_section()->cross_section_error() ;
os << std::endl;
}
}
bool GenEvent::use_momentum_unit( Units::MomentumUnit newunit ) {
// currently not exception-safe.
// Easy to fix, though, if needed.
if ( m_momentum_unit != newunit ) {
const double factor = Units::conversion_factor( m_momentum_unit, newunit );
// multiply all momenta by 'factor',
// loop is entered only if particle list is not empty
for ( GenEvent::particle_iterator p = particles_begin();
p != particles_end(); ++p )
{
(*p)->convert_momentum(factor);
}
// ...
m_momentum_unit = newunit;
}
return true;
}
bool GenEvent::use_length_unit( Units::LengthUnit newunit ) {
// currently not exception-safe.
// Easy to fix, though, if needed.
if ( m_position_unit != newunit ) {
const double factor = Units::conversion_factor( m_position_unit, newunit );
// multiply all lengths by 'factor',
// loop is entered only if vertex list is not empty
for ( GenEvent::vertex_iterator vtx = vertices_begin();
vtx != vertices_end(); ++vtx ) {
(*vtx)->convert_position(factor);
}
// ...
m_position_unit = newunit;
}
return true;
}
bool GenEvent::use_momentum_unit( std::string& newunit ) {
if ( newunit == "MEV" ) return use_momentum_unit( Units::MEV );
else if( newunit == "GEV" ) return use_momentum_unit( Units::GEV );
else std::cerr << "GenEvent::use_momentum_unit ERROR: use either MEV or GEV\n";
return false;
}
bool GenEvent::use_length_unit( std::string& newunit ) {
if ( newunit == "MM" ) return use_length_unit( Units::MM );
else if( newunit == "CM" ) return use_length_unit( Units::CM );
else std::cerr << "GenEvent::use_length_unit ERROR: use either MEV or GEV\n";
return false;
}
bool GenEvent::is_valid() const {
/// A GenEvent is presumed valid if it has both associated
/// particles and vertices. No other information is checked.
if ( vertices_empty() ) return false;
if ( particles_empty() ) return false;
return true;
}
std::ostream & GenEvent::write_beam_particles(std::ostream & os,
std::pair<HepMC::GenParticle *,HepMC::GenParticle *> pr )
{
GenParticle* p = pr.first;
if(!p) {
detail::output( os, 0 );
} else {
detail::output( os, p->barcode() );
}
p = pr.second;
if(!p) {
detail::output( os, 0 );
} else {
detail::output( os, p->barcode() );
}
return os;
}
std::ostream & GenEvent::write_vertex(std::ostream & os, GenVertex const * v)
{
if ( !v || !os ) {
std::cerr << "GenEvent::write_vertex !v||!os, "
<< "v="<< v << " setting badbit" << std::endl;
os.clear(std::ios::badbit);
return os;
}
// First collect info we need
// count the number of orphan particles going into v
int num_orphans_in = 0;
for ( GenVertex::particles_in_const_iterator p1
= v->particles_in_const_begin();
p1 != v->particles_in_const_end(); ++p1 ) {
if ( !(*p1)->production_vertex() ) ++num_orphans_in;
}
//
os << 'V';
detail::output( os, v->barcode() ); // v's unique identifier
detail::output( os, v->id() );
detail::output( os, v->position().x() );
detail::output( os, v->position().y() );
detail::output( os, v->position().z() );
detail::output( os, v->position().t() );
detail::output( os, num_orphans_in );
detail::output( os, (int)v->particles_out_size() );
detail::output( os, (int)v->weights().size() );
for ( WeightContainer::const_iterator w = v->weights().begin();
w != v->weights().end(); ++w ) {
detail::output( os, *w );
}
detail::output( os,'\n');
// incoming particles
for ( GenVertex::particles_in_const_iterator p2
= v->particles_in_const_begin();
p2 != v->particles_in_const_end(); ++p2 ) {
if ( !(*p2)->production_vertex() ) {
write_particle( os, *p2 );
}
}
// outgoing particles
for ( GenVertex::particles_out_const_iterator p3
= v->particles_out_const_begin();
p3 != v->particles_out_const_end(); ++p3 ) {
write_particle( os, *p3 );
}
return os;
}
std::ostream & GenEvent::write_particle( std::ostream & os, GenParticle const * p )
{
if ( !p || !os ) {
std::cerr << "GenEvent::write_particle !p||!os, "
<< "p="<< p << " setting badbit" << std::endl;
os.clear(std::ios::badbit);
return os;
}
os << 'P';
detail::output( os, p->barcode() );
detail::output( os, p->pdg_id() );
detail::output( os, p->momentum().px() );
detail::output( os, p->momentum().py() );
detail::output( os, p->momentum().pz() );
detail::output( os, p->momentum().e() );
detail::output( os, p->generated_mass() );
detail::output( os, p->status() );
detail::output( os, p->polarization().theta() );
detail::output( os, p->polarization().phi() );
// since end_vertex is oftentimes null, this CREATES a null vertex
// in the map
detail::output( os, ( p->end_vertex() ? p->end_vertex()->barcode() : 0 ) );
os << ' ' << p->flow() << "\n";
return os;
}
} // HepMC