hepmc - Blame information for rev 537

Subversion Repositories:

Rev 521 - Details - Compare with Previous - Last modification - View Log -

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