Index: trunk/Utilities/HepMC/src/CompareGenEvent.cc
===================================================================
--- trunk/Utilities/HepMC/src/CompareGenEvent.cc	(revision 349)
+++ trunk/Utilities/HepMC/src/CompareGenEvent.cc	(revision 571)
@@ -91,18 +91,7 @@
 
 bool compareWeights( GenEvent* e1, GenEvent* e2 ) {
-   WeightContainer w1 = e1->weights();
-   WeightContainer w2 = e2->weights();
-   if( w1.size() != w2.size() ) { 
-       std::cerr << "compareWeights: size of weight container differs " << std::endl;
-       return false; 
-   }
-   for( int i=0; i<w1.size(); ++i ) {
-       if( w1[i] != w2[i] ) {
-	   std::cerr << "compareWeights: weight container entry " 
-	             << i << " differs" << std::endl;
-	   return false; 
-       }
-   }
-   return true;
+   if( e1->weights() == e2->weights() ) return true;
+   std::cerr << "compareWeights: weight containers differ " << std::endl;
+   return false;
 }
 
Index: trunk/Utilities/HepMC/src/GenCrossSection.cc
===================================================================
--- trunk/Utilities/HepMC/src/GenCrossSection.cc	(revision 571)
+++ trunk/Utilities/HepMC/src/GenCrossSection.cc	(revision 571)
@@ -0,0 +1,106 @@
+//--------------------------------------------------------------------------
+//////////////////////////////////////////////////////////////////////////
+// garren@fnal.gov, January 2009
+//
+// The singleton GenCrossSection class holds run level information, 
+// such as the cross section.
+// 
+//////////////////////////////////////////////////////////////////////////
+//--------------------------------------------------------------------------
+
+#include <iostream>
+#include <sstream>
+
+#include "GenCrossSection.h"
+#include "IO_Exception.h"
+
+namespace HepMC {
+
+GenCrossSection::GenCrossSection( GenCrossSection const & orig )
+    : m_cross_section( orig.cross_section() ),
+      m_cross_section_error( orig.cross_section_error() ),
+      m_is_set( orig.is_set() )
+{}
+
+void GenCrossSection::swap( GenCrossSection & other)
+{
+    std::swap( m_cross_section, other.m_cross_section );
+    std::swap( m_cross_section_error, other.m_cross_section_error );
+    std::swap( m_is_set, other.m_is_set );
+}
+
+GenCrossSection &  GenCrossSection::operator = ( GenCrossSection const & rhs )
+{
+    GenCrossSection tmp( rhs );
+    swap( tmp );
+    return *this;
+}
+
+bool GenCrossSection::operator==( const GenCrossSection& rhs ) const
+{
+    if( rhs.cross_section() != this->cross_section() ) return false;
+    if( rhs.cross_section_error() != this->cross_section_error() ) return false;
+    return true;
+}
+
+bool GenCrossSection::operator!=( const GenCrossSection& rhs ) const
+{
+	return !( rhs == *this );
+}
+
+
+void GenCrossSection::clear() 
+{
+    m_cross_section       = 0.0;
+    m_cross_section_error = 0.0;
+    m_is_set              = false;
+}
+
+std::ostream & GenCrossSection::write( std::ostream & os ) const
+{
+    // make sure the stream is valid
+    if ( !os ) {
+	std::cerr << "GenCrossSection::print !os, setting badbit" << std::endl;
+	os.clear(std::ios::badbit); 
+	return os;
+    }
+    // write the GenCrossSection information if the cross section was set
+    if( is_set() ) {
+	os << "C " << m_cross_section 
+	   << " " << m_cross_section_error 
+	   << "\n";
+    }
+    return os;
+}
+
+std::istream & GenCrossSection::read( std::istream & is )
+{
+    // make sure the stream is valid
+    if ( !is ) { 
+      std::cerr << "GenCrossSection stream input setting badbit." << std::endl;
+      is.clear(std::ios::badbit); 
+      return is; 
+    }
+    // check to see if we have a GenCrossSection line
+    // This line is optional and may not exist
+    if ( is.peek()!='C' ) { 
+      return is; 
+    }
+    // get the GenCrossSection line
+    std::string line, firstc;
+    std::getline(is,line);
+    std::istringstream iline(line);
+    // Get first character and throw it away
+    iline >> firstc;
+    // Now get the numbers
+    double xs = 0., xserr = 0.;
+    iline >> xs ;
+    if(!iline) throw IO_Exception("GenCrossSection::read encounterd invalid data");
+    iline >> xserr ;
+    if(!iline) throw IO_Exception("GenCrossSection::read encounterd invalid data");
+    // set the data members
+    set_cross_section( xs, xserr );
+    return  is;
+}
+
+} // HepMC
Index: trunk/Utilities/HepMC/src/GenEvent.cc
===================================================================
--- trunk/Utilities/HepMC/src/GenEvent.cc	(revision 349)
+++ trunk/Utilities/HepMC/src/GenEvent.cc	(revision 571)
@@ -14,5 +14,7 @@
 
 #include "GenEvent.h"
+#include "GenCrossSection.h"
 #include "Version.h"
+#include "StreamHelpers.h"
 
 namespace HepMC {
@@ -38,4 +40,5 @@
 	m_vertex_barcodes(),
 	m_particle_barcodes(),
+	m_cross_section(0), 
 	m_heavy_ion(0), 
 	m_pdf_info(0),
@@ -71,4 +74,5 @@
 	m_vertex_barcodes(),
 	m_particle_barcodes(),
+	m_cross_section(0), 
 	m_heavy_ion( new HeavyIon(ion) ), 
 	m_pdf_info( new PdfInfo(pdf) ),
@@ -102,4 +106,5 @@
 	m_vertex_barcodes(),
 	m_particle_barcodes(),
+	m_cross_section(0), 
 	m_heavy_ion(0), 
 	m_pdf_info(0),
@@ -136,4 +141,5 @@
 	m_vertex_barcodes(),
 	m_particle_barcodes(),
+	m_cross_section(0), 
 	m_heavy_ion( new HeavyIon(ion) ), 
 	m_pdf_info( new PdfInfo(pdf) ),
@@ -162,4 +168,5 @@
 	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 ),
@@ -168,5 +175,4 @@
     {
 	/// deep copy - makes a copy of all vertices!
-	//++s_counter;
 	//
 
@@ -234,4 +240,5 @@
 	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             );
@@ -252,10 +259,10 @@
     {
 	/// Deep destructor.
-	/// deletes all vertices/particles in this evt
-	///
+	/// 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;
-	//--s_counter;
     }
 
@@ -281,9 +288,6 @@
 		  : 0 )
 	     << "\n";
-	//ostr << " Current Memory Usage: " 
-	//     << GenEvent::counter() << " events, "
-	//     << GenVertex::counter() << " vertices, "
-	//     << GenParticle::counter() << " particles.\n"; 
 	write_units( ostr );
+        write_cross_section(ostr);
 	ostr << " Entries this event: " << vertices_size() << " vertices, "
 	     << particles_size() << " particles.\n"; 
@@ -302,7 +306,5 @@
 	// Weights
 	ostr << " Wgts(" << weights().size() << ")=";
-	for ( WeightContainer::const_iterator wgt = weights().begin();
-	      wgt != weights().end(); ++wgt ) { ostr << *wgt << " "; }
-	ostr << "\n";
+	weights().print(ostr);
 	ostr << " EventScale " << event_scale() 
 	     << " [energy] \t alphaQCD=" << alphaQCD() 
@@ -367,10 +369,16 @@
 	///
 	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;
@@ -378,6 +386,7 @@
 	m_weights = std::vector<double>();
 	m_random_states = std::vector<long>();
-	//m_momentum_unit = ;
-	//m_position_unit = ;
+	// 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 
@@ -388,8 +397,4 @@
                       << m_vertex_barcodes.size() << "  " 
 		      << m_particle_barcodes.size() << std::endl;
-            //std::cerr << "Total Number vtx,particle in memory "
-             //         << "after method called = "
-             //         << GenVertex::counter() << "\t"
-		//      << GenParticle::counter() << std::endl;
         }
 	return;
@@ -420,8 +425,4 @@
                       << m_vertex_barcodes.size() << "  " 
 		      << m_particle_barcodes.size() << std::endl;
-            //std::cerr << "Total Number vtx,particle in memory "
-            //          << "after method called = "
-            //          << GenVertex::counter() << "\t"
-		//      << GenParticle::counter() << std::endl;
 	}
     }
@@ -602,4 +603,15 @@
     }
 
+    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. 
@@ -651,3 +663,106 @@
     }  
 
+    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
Index: trunk/Utilities/HepMC/src/GenEventStreamIO.cc
===================================================================
--- trunk/Utilities/HepMC/src/GenEventStreamIO.cc	(revision 571)
+++ trunk/Utilities/HepMC/src/GenEventStreamIO.cc	(revision 571)
@@ -0,0 +1,819 @@
+//--------------------------------------------------------------------------
+//
+// GenEventStreamIO.cc
+// Author:  Lynn Garren
+//
+// Implement operator >> and operator <<
+//
+// ----------------------------------------------------------------------
+
+#include <iostream>
+#include <ostream>
+#include <istream>
+#include <sstream>
+
+#include "GenEvent.h"
+#include "GenCrossSection.h"
+#include "StreamInfo.h"
+#include "StreamHelpers.h"
+#include "Version.h"
+#include "IO_Exception.h"
+
+namespace HepMC {
+
+// ------------------------- local methods ----------------
+
+/// This method is called by the stream destructor.
+/// It does cleanup on stored user data (StreamInfo)
+/// and is registered by the first call to get_stream_info().
+void HepMCStreamCallback(std::ios_base::event e, std::ios_base& b, int i)
+{
+  // only clean up if the stream object is going away.
+  if(i!=0 && e!= std::ios_base::erase_event) return;
+
+  // retrieve the pointer to the object
+  StreamInfo* hd = (StreamInfo*)b.pword(i);
+  b.pword(i) = 0;
+  b.iword(i) = 0;
+#ifdef HEPMC_DEBUG
+  // the following line is just for sanity checking
+  if(hd) std::cerr << "deleted StreamInfo " << hd->stream_id() << "\n";
+#endif
+  delete hd;
+}
+
+// ------------------------- iomanip ----------------
+
+/// A custom iomanip that allows us to store and access user data (StreamInfo)
+/// associated with the stream.
+/// This method creates the StreamInfo object the first time it is called.
+template <class IO>
+StreamInfo& get_stream_info(IO& iost)
+{
+  if(iost.iword(0) == 0)
+    {
+      // make sure we add the callback if this is the first time through
+      iost.iword(0)=1;
+      iost.register_callback(&HepMCStreamCallback, 0);
+      // this is our special "context" record.
+      // there is one of these at the head of each IO block.
+      // allocate room for a StreamInfo in the userdata area
+      iost.pword(0) = new StreamInfo;
+#ifdef HEPMC_DEBUG
+      // the following line is just for sanity checking
+      std::cerr << "created StreamInfo " << ((StreamInfo*)iost.pword(0))->stream_id() << "\n";
+#endif
+    }
+  return *(StreamInfo*)iost.pword(0);
+}
+	
+// ------------------------- GenEvent member functions ----------------
+
+std::ostream& GenEvent::write( std::ostream& os )
+{
+    /// Writes evt to an output stream.
+
+    //
+    StreamInfo & info = get_stream_info(os);
+    //
+    // if this is the first event, set precision
+    if ( !info.finished_first_event() ) {
+	// precision 16 (# digits following decimal point) is the minimum that
+	//  will capture the full information stored in a double
+        //  However, we let the user set precision, since that is the expected functionality
+	// we use decimal to store integers, because it is smaller than hex!
+	os.setf(std::ios::dec,std::ios::basefield);
+	os.setf(std::ios::scientific,std::ios::floatfield);
+	//
+	info.set_finished_first_event(true);
+    }
+    //
+    // output the event data including the number of primary vertices
+    //  and the total number of vertices
+    //std::vector<long> random_states = random_states();
+    os << 'E';
+    detail::output( os, event_number() );
+    detail::output( os, mpi() );
+    detail::output( os, event_scale() );
+    detail::output( os, alphaQCD() );
+    detail::output( os, alphaQED() );
+    detail::output( os, signal_process_id() );
+    detail::output( os,   ( signal_process_vertex() ?
+		signal_process_vertex()->barcode() : 0 )   );
+    detail::output( os, vertices_size() ); // total number of vertices.
+    write_beam_particles( os, beam_particles() );
+    // random state
+    detail::output( os, (int)m_random_states.size() );
+    for ( std::vector<long>::iterator rs = m_random_states.begin(); 
+	  rs != m_random_states.end(); ++rs ) {
+	 detail::output( os, *rs );
+    }
+    // weights
+    // we need to iterate over the map so that the weights printed 
+    // here will be in the same order as the names printed next
+    os << ' ' << (int)weights().size() ;
+    for ( WeightContainer::const_map_iterator w = weights().map_begin(); 
+	  w != weights().map_end(); ++w ) {
+        detail::output( os, m_weights[w->second] );
+    }
+    detail::output( os,'\n');
+    // now add names for weights
+    // note that this prints a new line if and only if the weight container
+    // is not empty
+    if ( ! weights().empty() ) {
+	os << "N " << weights().size() << " " ;
+	for ( WeightContainer::const_map_iterator w = weights().map_begin(); 
+	      w != weights().map_end(); ++w ) {
+	    detail::output( os,'"');
+	    os << w->first;
+	    detail::output( os,'"');
+	    detail::output( os,' ');
+	}
+	detail::output( os,'\n');
+    }
+    //
+    // Units
+    os << "U " << name(momentum_unit());
+    os << " " << name(length_unit());
+    detail::output( os,'\n');
+    //
+    // write GenCrossSection if it has been set
+    if( m_cross_section ) m_cross_section->write(os);
+    //
+    // write HeavyIon and PdfInfo if they have been set
+    if( m_heavy_ion ) os << heavy_ion() ;
+    if( m_pdf_info ) os << pdf_info() ;
+    //
+    // Output all of the vertices - note there is no real order.
+    for ( GenEvent::vertex_const_iterator v = vertices_begin();
+	  v != vertices_end(); ++v ) {
+	write_vertex(os, *v);
+    }
+    return os;
+}
+
+std::istream& GenEvent::read( std::istream& is )
+{
+    /// read a GenEvent from streaming input
+    //
+    StreamInfo & info = get_stream_info(is);
+    clear();
+    //
+    // search for event listing key before first event only.
+    if ( !info.finished_first_event() ) {
+	//
+	find_file_type(is);
+	info.set_finished_first_event(true);
+    }
+    //
+    // make sure the stream is good
+    if ( !is ) {
+	std::cerr << "streaming input: end of stream found "
+		  << "setting badbit." << std::endl;
+	is.clear(std::ios::badbit); 
+        return is;
+    }
+
+    //
+    // test to be sure the next entry is of type "E" then ignore it
+    if ( is.peek()!='E' ) { 
+	// if the E is not the next entry, then check to see if it is
+	// the end event listing key - if yes, search for another start key
+	int ioendtype;
+	find_end_key(is,ioendtype);
+	if ( ioendtype == info.io_type() ) {
+	    find_file_type(is);
+	    // are we at the end of the file?
+	    if( !is ) return is;
+	} else if ( ioendtype > 0 ) {
+	    std::cerr << "streaming input: end key does not match start key "
+		      << "setting badbit." << std::endl;
+	    is.clear(std::ios::badbit); 
+	    return is;
+	} else if ( !info.has_key() ) {
+	    find_file_type(is);
+	    // are we at the end of the file?
+	    if( !is ) return is;
+	} else {
+	    std::cerr << "streaming input: end key not found "
+		      << "setting badbit." << std::endl;
+	    is.clear(std::ios::badbit); 
+	    return is;
+	}
+    } 
+
+    int signal_process_vertex = 0;
+    int num_vertices = 0, bp1 = 0, bp2 = 0;
+    // OK - now ready to start reading the event, so set the header flag
+    info.set_reading_event_header(true);
+    // The flag will be set to false when we reach the end of the header
+    while(info.reading_event_header()) {
+	switch(is.peek()) {
+	    case 'E':
+	    {	// deal with the event line
+		process_event_line( is, num_vertices, bp1, bp2, signal_process_vertex );
+	    } break;
+	    case 'N':
+	    {	// get weight names 
+	        read_weight_names( is );
+	    } break;
+	    case 'U':
+	    {	// get unit information if it exists
+		if( info.io_type() == gen ) {
+		    read_units( is );
+		}
+	    } break;
+	    case 'C':
+            {	// we have a GenCrossSection line
+		// create cross section
+		GenCrossSection xs;
+		// check for invalid data
+		try {
+		    // read the line
+		    xs.read(is);
+		}
+		catch (IO_Exception& e) {
+        	    detail::find_event_end( is );
+		}
+		if(xs.is_set()) { 
+		    set_cross_section( xs );
+		}
+            } break;
+	    case 'H':
+	    {	// we have a HeavyIon line OR an unexpected HepMC... line
+		if( info.io_type() == gen || info.io_type() == extascii ) {
+		    // get HeavyIon
+		    HeavyIon ion;
+		    // check for invalid data
+		    try {
+			is >> &ion;
+		    }
+		    catch (IO_Exception& e) {
+        		detail::find_event_end( is );
+		    }
+		    if(ion.is_valid()) { 
+			set_heavy_ion( ion );
+		    }
+		}
+	    } break;
+	    case 'F':
+	    {	// we have a PdfInfo line
+		if( info.io_type() == gen || info.io_type() == extascii ) {
+		    // get PdfInfo
+		    PdfInfo pdf;
+		    // check for invalid data
+		    try {
+			is >> &pdf;
+		    }
+		    catch (IO_Exception& e) {
+        		detail::find_event_end( is );
+		    }
+		    if(pdf.is_valid()) { 
+			set_pdf_info( pdf );
+		    }
+		}
+	    } break;
+	    case 'V':
+	    {	
+	        // this should be the first vertex line - exit this loop
+	        info.set_reading_event_header(false);
+	    } break;
+	    case 'P':
+	    {	// we should not find this line
+	        std::cerr << "streaming input: found unexpected line P" << std::endl;
+	        info.set_reading_event_header(false);
+	    } break;
+	    default:
+	        // ignore everything else
+	        break;
+	} // switch on line type
+    } // while reading_event_header
+    //
+    // the end vertices of the particles are not connected until
+    //  after the event is read --- we store the values in a map until then
+    TempParticleMap particle_to_end_vertex;
+    //
+    // read in the vertices
+    for ( int iii = 1; iii <= num_vertices; ++iii ) {
+	GenVertex* v = new GenVertex();
+	try {
+	    detail::read_vertex(is,particle_to_end_vertex,v);
+	}
+	catch (IO_Exception& e) {
+	    for( TempParticleMap::orderIterator it = particle_to_end_vertex.order_begin(); 
+	         it != particle_to_end_vertex.order_end(); ++it ) {
+		GenParticle* p = it->second;
+		// delete particles only if they are not already owned by a vertex
+        	if( p->production_vertex() ) {
+		} else if( p->end_vertex() ) {
+		} else {
+	             delete p;
+ 		}
+ 	    }
+	    delete v;
+            detail::find_event_end( is );
+	}
+	add_vertex( v );
+    }
+    // set the signal process vertex
+    if ( signal_process_vertex ) {
+	set_signal_process_vertex( 
+	    barcode_to_vertex(signal_process_vertex) );
+    }
+    //
+    // last connect particles to their end vertices
+    GenParticle* beam1(0);
+    GenParticle* beam2(0);
+    for ( TempParticleMap::orderIterator pmap 
+	      = particle_to_end_vertex.order_begin(); 
+	  pmap != particle_to_end_vertex.order_end(); ++pmap ) {
+	GenParticle* p =  pmap->second;
+	int vtx = particle_to_end_vertex.end_vertex( p );
+	GenVertex* itsDecayVtx = barcode_to_vertex(vtx);
+	if ( itsDecayVtx ) itsDecayVtx->add_particle_in( p );
+	else {
+	    std::cerr << "read_io_genevent: ERROR particle points"
+		      << " to null end vertex. " <<std::endl;
+	}
+	// also look for the beam particles
+	if( p->barcode() == bp1 ) beam1 = p;
+	if( p->barcode() == bp2 ) beam2 = p;
+    }
+    set_beam_particles(beam1,beam2);
+    return is;
+}
+
+// ------------------------- operator << and operator >> ----------------
+
+std::ostream & operator << (std::ostream & os, GenEvent & evt)
+{
+    /// Writes evt to an output stream.
+    evt.write(os);
+    return os;
+}
+
+std::istream & operator >> (std::istream & is, GenEvent & evt)
+{
+    evt.read(is);
+    return is;
+}
+
+// ------------------------- set units ----------------
+
+std::istream & set_input_units(std::istream & is, 
+                               Units::MomentumUnit mom,
+			       Units::LengthUnit len )
+{
+    //
+    StreamInfo & info = get_stream_info(is);
+    info.use_input_units( mom, len );
+    return is;
+}
+
+// ------------------------- begin and end block lines ----------------
+
+std::ostream & write_HepMC_IO_block_begin(std::ostream & os )
+{
+    //
+    StreamInfo & info = get_stream_info(os);
+
+    if( !info.finished_first_event() ) {
+    os << "\n" << "HepMC::Version " << versionName();
+    os << "\n";
+    os << info.IO_GenEvent_Key() << "\n";
+    }
+    return os;
+}
+
+std::ostream & write_HepMC_IO_block_end(std::ostream & os )
+{
+    //
+    StreamInfo & info = get_stream_info(os);
+
+    if( info.finished_first_event() ) {
+	os << info.IO_GenEvent_End() << "\n";
+	os << std::flush;
+    }
+    return os;
+}
+
+std::istream & GenEvent::process_event_line( std::istream & is, 
+                                             int & num_vertices,
+					     int & bp1, int & bp2,
+					     int & signal_process_vertex )
+{
+    //
+    if ( !is ) {
+	std::cerr << "GenEvent::process_event_line setting badbit." << std::endl;
+	is.clear(std::ios::badbit);
+	return is;
+    } 
+    //
+    StreamInfo & info = get_stream_info(is);
+    std::string line;
+    std::getline(is,line);
+    std::istringstream iline(line);
+    std::string firstc;
+    iline >> firstc;
+    //
+    // read values into temp variables, then fill GenEvent
+    int event_number = 0, signal_process_id = 0,
+	random_states_size = 0, nmpi = -1;
+    double eventScale = 0, alpha_qcd = 0, alpha_qed = 0;
+    iline >> event_number;
+    if(!iline) detail::find_event_end( is );
+    if( info.io_type() == gen || info.io_type() == extascii ) {
+        iline >> nmpi;
+        if(!iline) detail::find_event_end( is );
+        set_mpi( nmpi );
+    }
+    iline >> eventScale ;
+    if(!iline) detail::find_event_end( is );
+    iline >> alpha_qcd ;
+    if(!iline) detail::find_event_end( is );
+    iline >> alpha_qed;
+    if(!iline) detail::find_event_end( is );
+    iline >> signal_process_id ;
+    if(!iline) detail::find_event_end( is );
+    iline >> signal_process_vertex;
+    if(!iline) detail::find_event_end( is );
+    iline >> num_vertices;
+    if(!iline) detail::find_event_end( is );
+    if( info.io_type() == gen || info.io_type() == extascii ) {
+        iline >> bp1 ;
+        if(!iline) detail::find_event_end( is );
+	iline >> bp2;
+        if(!iline) detail::find_event_end( is );
+    }
+    iline >> random_states_size;
+    if(!iline) detail::find_event_end( is );
+    std::vector<long> random_states(random_states_size);
+    for ( int i = 0; i < random_states_size; ++i ) {
+	iline >> random_states[i];
+        if(!iline) detail::find_event_end( is );
+    }
+    WeightContainer::size_type weights_size = 0;
+    iline >> weights_size;
+    if(!iline) detail::find_event_end( is );
+    std::vector<double> wgt(weights_size);
+    for ( WeightContainer::size_type ii = 0; ii < weights_size; ++ii ) {
+        iline >> wgt[ii];
+        if(!iline) detail::find_event_end( is );
+    }
+    // weight names will be added later if they exist
+    if( weights_size > 0 ) m_weights = wgt;
+    // 
+    // fill signal_process_id, event_number, random_states, etc.
+    set_signal_process_id( signal_process_id );
+    set_event_number( event_number );
+    set_random_states( random_states );
+    set_event_scale( eventScale );
+    set_alphaQCD( alpha_qcd );
+    set_alphaQED( alpha_qed );
+    //
+    return is;
+}
+
+std::istream & GenEvent::read_weight_names( std::istream & is )
+{
+    // now check for a named weight line
+    if ( !is ) {
+	std::cerr << "GenEvent::read_weight_names setting badbit." << std::endl;
+	is.clear(std::ios::badbit);
+	return is;
+    } 
+    // Test to be sure the next entry is of type "N"
+    // If we have no named weight line, this is not an error
+    // releases prior to 2.06.00 do not have named weights
+    if ( is.peek() !='N') {
+	return is;
+    } 
+    // now get this line and process it
+    std::string line;
+    std::getline(is,line);
+    std::istringstream wline(line);
+    std::string firstc;
+    WeightContainer::size_type name_size = 0;
+    wline >> firstc >> name_size;
+    if(!wline) detail::find_event_end( is );
+    if( firstc != "N") { 
+        std::cout << "debug: first character of named weights is " << firstc << std::endl;
+        std::cout << "debug: We should never get here" << std::endl;
+	is.clear(std::ios::badbit);
+	return is;
+    }
+    if( m_weights.size() != name_size ) { 
+        std::cout << "debug: weight sizes do not match "<< std::endl;
+        std::cout << "debug: weight vector size is " << m_weights.size() << std::endl;
+        std::cout << "debug: weight name size is " << name_size << std::endl;
+	is.clear(std::ios::badbit);
+	return is;
+    }
+    std::string name;
+    std::string::size_type i1 = line.find("\"");
+    std::string::size_type i2;
+    std::string::size_type len = line.size();
+    WeightContainer namedWeight;
+    for ( WeightContainer::size_type ii = 0; ii < name_size; ++ii ) {
+        // weight names may contain blanks
+        if(i1 >= len) {
+            std::cout << "debug: attempting to read past the end of the named weight line " << std::endl;
+            std::cout << "debug: We should never get here" << std::endl;
+            std::cout << "debug: Looking for the end of this event" << std::endl;
+	    detail::find_event_end( is );
+	}
+	i2 = line.find("\"",i1+1);
+	name = line.substr(i1+1,i2-i1-1);
+	namedWeight[name] = m_weights[ii];
+	i1 = line.find("\"",i2+1);
+    }
+    m_weights = namedWeight;
+    return is;
+}
+
+std::istream & GenEvent::read_units( std::istream & is )
+{
+    //
+    if ( !is ) {
+	std::cerr << "GenEvent::read_units setting badbit." << std::endl;
+	is.clear(std::ios::badbit);
+	return is;
+    } 
+    //
+    StreamInfo & info = get_stream_info(is);
+    // test to be sure the next entry is of type "U" then ignore it
+    // if we have no units, this is not an error
+    // releases prior to 2.04.00 did not write unit information
+    if ( is.peek() !='U') {
+ 	use_units( info.io_momentum_unit(), 
+	               info.io_position_unit() );
+	return is;
+    } 
+    is.ignore();	// ignore the first character in the line
+    std::string mom, pos;
+    is >> mom >> pos;
+    is.ignore(1);      // eat the extra whitespace
+    use_units(mom,pos);
+    //
+    return is;
+}
+
+std::istream & GenEvent::find_file_type( std::istream & istr )
+{
+    //
+    // make sure the stream is good
+    if ( !istr ) return istr;
+
+    //
+    StreamInfo & info = get_stream_info(istr);
+
+    // if there is no input block line, then we assume this stream
+    // is in the IO_GenEvent format
+    if ( istr.peek()=='E' ) {
+	info.set_io_type( gen );
+	info.set_has_key(false);
+        return istr;
+    }
+    
+    std::string line;
+    while ( std::getline(istr,line) ) {
+	//
+	// search for event listing key before first event only.
+	//
+	if( line == info.IO_GenEvent_Key() ) {
+	    info.set_io_type( gen );
+	    info.set_has_key(true);
+	    return istr;
+	} else if( line == info.IO_Ascii_Key() ) {
+	    info.set_io_type( ascii );
+	    info.set_has_key(true);
+	    return istr;
+	} else if( line == info.IO_ExtendedAscii_Key() ) {
+	    info.set_io_type( extascii );
+	    info.set_has_key(true);
+	    return istr;
+	} else if( line == info.IO_Ascii_PDT_Key() ) {
+	    info.set_io_type( ascii_pdt );
+	    info.set_has_key(true);
+	    return istr;
+	} else if( line == info.IO_ExtendedAscii_PDT_Key() ) {
+	    info.set_io_type( extascii_pdt );
+	    info.set_has_key(true);
+	    return istr;
+	}
+    }
+    info.set_io_type( 0 );
+    info.set_has_key(false);
+    return istr;
+}
+
+std::istream & GenEvent::find_end_key( std::istream & istr, int & iotype )
+{
+    iotype = 0;
+    // peek at the first character before proceeding
+    if( istr.peek()!='H' ) return istr;
+    //
+    // we only check the next line
+    std::string line;
+    std::getline(istr,line);
+    //
+    StreamInfo & info = get_stream_info(istr);
+    //
+    // check to see if this is an end key
+    if( line == info.IO_GenEvent_End() ) {
+	iotype = gen;
+    } else if( line == info.IO_Ascii_End() ) {
+	iotype = ascii;
+    } else if( line == info.IO_ExtendedAscii_End() ) {
+	iotype = extascii;
+    } else if( line == info.IO_Ascii_PDT_End() ) {
+	iotype = ascii_pdt;
+    } else if( line == info.IO_ExtendedAscii_PDT_End() ) {
+	iotype = extascii_pdt;
+    }
+    if( iotype != 0 && info.io_type() != iotype ) {
+        std::cerr << "GenEvent::find_end_key: iotype keys have changed" << std::endl;
+    } else {
+        return istr;
+    }
+    //
+    // if we get here, then something has gotten badly confused
+    std::cerr << "GenEvent::find_end_key: MALFORMED INPUT" << std::endl;
+    istr.clear(std::ios::badbit); 
+    return istr;
+}
+
+std::ostream & establish_output_stream_info( std::ostream & os )
+{
+    StreamInfo & info = get_stream_info(os);
+    if ( !info.finished_first_event() ) {
+	// precision 16 (# digits following decimal point) is the minimum that
+	//  will capture the full information stored in a double
+	os.precision(16);
+	// we use decimal to store integers, because it is smaller than hex!
+	os.setf(std::ios::dec,std::ios::basefield);
+	os.setf(std::ios::scientific,std::ios::floatfield);
+    }
+    return os;
+}
+
+std::istream & establish_input_stream_info( std::istream & is )
+{
+    StreamInfo & info = get_stream_info(is);
+    if ( !info.finished_first_event() ) {
+	// precision 16 (# digits following decimal point) is the minimum that
+	//  will capture the full information stored in a double
+	is.precision(16);
+	// we use decimal to store integers, because it is smaller than hex!
+	is.setf(std::ios::dec,std::ios::basefield);
+	is.setf(std::ios::scientific,std::ios::floatfield);
+    }
+    return is;
+}
+
+
+// ------------------------- helper functions ----------------
+
+namespace detail {
+
+// The functions defined here need to use get_stream_info
+
+std::istream & read_particle( std::istream & is, 
+                              TempParticleMap & particle_to_end_vertex, 
+			      GenParticle * p )
+{
+    // get the next line
+    std::string line;
+    std::getline(is,line);
+    std::istringstream iline(line);
+    std::string firstc;
+    iline >> firstc;
+    if( firstc != "P" ) { 
+	std::cerr << "StreamHelpers::detail::read_particle invalid line type: " 
+	          << firstc << std::endl;
+	std::cerr << "StreamHelpers::detail::read_particle setting badbit." 
+		  << std::endl;
+	is.clear(std::ios::badbit); 
+	return is;
+    } 
+    //
+    StreamInfo & info = get_stream_info(is);
+    //testHepMC.cc
+    // declare variables to be read in to, and read everything except flow
+    double px = 0., py = 0., pz = 0., e = 0., m = 0., theta = 0., phi = 0.;
+    int bar_code = 0, id = 0, status = 0, end_vtx_code = 0, flow_size = 0;
+    // check that the input stream is still OK after reading item
+    iline >> bar_code ;
+    if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    iline >> id ;
+    if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    iline >> px ;
+    if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    iline >> py ;
+    if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    iline >> pz ;
+    if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    iline >> e ;
+    if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    if( info.io_type() != ascii ) {
+	iline >> m ;
+        if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    }
+    iline >> status ;
+    if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    iline >> theta ;
+    if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    iline >> phi ;
+    if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    iline >> end_vtx_code ;
+    if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    iline >> flow_size;
+    if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+    //
+    // read flow patterns if any exist
+    Flow flow;
+    int code_index, code;
+    for ( int i = 1; i <= flow_size; ++i ) {
+	iline >> code_index >> code;
+        if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
+	flow.set_icode( code_index,code);
+    }
+    p->set_momentum( FourVector(px,py,pz,e) );
+    p->set_pdg_id( id );
+    p->set_status( status );
+    p->set_flow( flow );
+    p->set_polarization( Polarization(theta,phi) );
+    if( info.io_type() == ascii ) {
+        p->set_generated_mass( p->momentum().m() );
+    } else {
+        p->set_generated_mass( m );
+    }
+    p->suggest_barcode( bar_code );
+    //
+    // all particles are connected to their end vertex separately 
+    // after all particles and vertices have been created - so we keep
+    // a map of all particles that have end vertices
+    if ( end_vtx_code != 0 ) {
+	particle_to_end_vertex.addEndParticle(p,end_vtx_code);
+    }
+    return is;
+}
+
+std::istream & read_units( std::istream & is, GenEvent & evt )
+{
+    //
+    if ( !is ) {
+	std::cerr << "StreamHelpers read_units setting badbit." << std::endl;
+	is.clear(std::ios::badbit);
+	return is;
+    } 
+    //
+    StreamInfo & info = get_stream_info(is);
+    // test to be sure the next entry is of type "U" then ignore it
+    // if we have no units, this is not an error
+    // releases prior to 2.04.00 did not write unit information
+    if ( is.peek() !='U') {
+ 	evt.use_units( info.io_momentum_unit(), 
+	               info.io_position_unit() );
+	return is;
+    } 
+    is.ignore();	// ignore the first character in the line
+    std::string mom, pos;
+    is >> mom >> pos;
+    is.ignore(1);      // eat the extra whitespace
+    evt.use_units(mom,pos);
+    //
+    return is;
+}
+
+std::ostream & establish_output_stream_info( std::ostream & os )
+{
+    StreamInfo & info = get_stream_info(os);
+    if ( !info.finished_first_event() ) {
+	// precision 16 (# digits following decimal point) is the minimum that
+	//  will capture the full information stored in a double
+	os.precision(16);
+	// we use decimal to store integers, because it is smaller than hex!
+	os.setf(std::ios::dec,std::ios::basefield);
+	os.setf(std::ios::scientific,std::ios::floatfield);
+    }
+    return os;
+}
+
+std::istream & establish_input_stream_info( std::istream & is )
+{
+    StreamInfo & info = get_stream_info(is);
+    if ( !info.finished_first_event() ) {
+	// precision 16 (# digits following decimal point) is the minimum that
+	//  will capture the full information stored in a double
+	is.precision(16);
+	// we use decimal to store integers, because it is smaller than hex!
+	is.setf(std::ios::dec,std::ios::basefield);
+	is.setf(std::ios::scientific,std::ios::floatfield);
+    }
+    return is;
+}
+
+} // detail
+
+} // HepMC
Index: trunk/Utilities/HepMC/src/GenParticle.cc
===================================================================
--- trunk/Utilities/HepMC/src/GenParticle.cc	(revision 349)
+++ trunk/Utilities/HepMC/src/GenParticle.cc	(revision 571)
@@ -39,5 +39,5 @@
 	m_pdg_id( inparticle.pdg_id() ), 
 	m_status( inparticle.status() ), 
-	m_flow(this),
+	m_flow(inparticle.flow()),
 	m_polarization( inparticle.polarization() ),
 	m_production_vertex(0), 
@@ -188,4 +188,7 @@
     /// Dump this particle's full info to ostr
     std::ostream& operator<<( std::ostream& ostr, const GenParticle& part ) {
+        // find the current stream state
+	std::ios_base::fmtflags orig = ostr.flags();
+ 	std::streamsize prec = ostr.precision();
 	ostr << " ";
 	ostr.width(9);
@@ -199,11 +202,8 @@
 	ostr << part.momentum().px() << ",";
 	ostr.width(9);
-        ostr.precision(2);
 	ostr << part.momentum().py() << ",";
 	ostr.width(9);
-        ostr.precision(2);
 	ostr << part.momentum().pz() << ",";
 	ostr.width(9);
-        ostr.precision(2);
 	ostr << part.momentum().e() << " ";
         ostr.setf(std::ios::fmtflags(0), std::ios::floatfield);
@@ -227,4 +227,7 @@
 	    ostr << (void*)part.end_vertex();
 	}
+        // restore the stream state
+        ostr.flags(orig);
+        ostr.precision(prec);
 	return ostr;
     }
Index: trunk/Utilities/HepMC/src/GenRanges.cc
===================================================================
--- trunk/Utilities/HepMC/src/GenRanges.cc	(revision 571)
+++ trunk/Utilities/HepMC/src/GenRanges.cc	(revision 571)
@@ -0,0 +1,86 @@
+//--------------------------------------------------------------------------
+//////////////////////////////////////////////////////////////////////////
+// garren@fnal.gov, March 2010
+// 
+//////////////////////////////////////////////////////////////////////////
+//--------------------------------------------------------------------------
+
+#include <iostream>
+
+#include "GenRanges.h"
+#include "GenEvent.h"
+#include "GenVertex.h"
+
+namespace HepMC {
+
+GenEventVertexRange GenEvent::vertex_range()
+{
+    return GenEventVertexRange(*this);
+}
+
+ConstGenEventVertexRange GenEvent::vertex_range() const
+{
+    return ConstGenEventVertexRange(*this);
+}
+
+GenEventParticleRange GenEvent::particle_range()
+{
+    return GenEventParticleRange(*this);
+}
+
+ConstGenEventParticleRange GenEvent::particle_range() const
+{
+    return ConstGenEventParticleRange(*this);
+}
+
+GenVertexParticleRange GenVertex::particles( IteratorRange range )
+{
+    return GenVertexParticleRange(*this,range);
+}
+
+GenParticleProductionRange GenVertex::particles_in( GenParticle& p, IteratorRange range )
+{
+    return GenParticleProductionRange(p,range);
+}
+
+ConstGenParticleProductionRange GenVertex:: particles_in( GenParticle const & p, IteratorRange range ) const
+{
+    return ConstGenParticleProductionRange(p,range);
+}
+
+GenParticleEndRange GenVertex::particles_out( GenParticle& p, IteratorRange range )
+{
+    return GenParticleEndRange(p,range);
+}
+
+ConstGenParticleEndRange GenVertex::particles_out( GenParticle const & p, IteratorRange range ) const
+{
+    return ConstGenParticleEndRange(p,range);
+}
+
+GenParticleProductionRange GenParticle::particles_in( IteratorRange range )
+{
+    return GenParticleProductionRange(*this,range);
+}
+
+
+ConstGenParticleProductionRange GenParticle::particles_in( IteratorRange range ) const
+{
+    return ConstGenParticleProductionRange(*this,range);
+}
+
+
+GenParticleEndRange GenParticle::particles_out( IteratorRange range )
+{
+    return GenParticleEndRange(*this,range);
+}
+
+
+ConstGenParticleEndRange GenParticle::particles_out( IteratorRange range ) const
+{
+    return ConstGenParticleEndRange(*this,range);
+}
+
+
+
+} // HepMC
Index: trunk/Utilities/HepMC/src/GenVertex.cc
===================================================================
--- trunk/Utilities/HepMC/src/GenVertex.cc	(revision 349)
+++ trunk/Utilities/HepMC/src/GenVertex.cc	(revision 571)
@@ -144,4 +144,7 @@
 
     void GenVertex::print( std::ostream& ostr ) const {
+        // find the current stream state
+	std::ios_base::fmtflags orig = ostr.flags();
+ 	std::streamsize prec = ostr.precision();
 	if ( barcode()!=0 ) {
 	    if ( position() != FourVector(0,0,0,0) ) {
@@ -159,11 +162,8 @@
 		ostr << position().x() << ",";
 		ostr.width(9);
-                ostr.precision(2);
 		ostr << position().y() << ",";
 		ostr.width(9);
-                ostr.precision(2);
 		ostr << position().z() << ",";
 		ostr.width(9);
-                ostr.precision(2);
 		ostr << position().t();
                 ostr.setf(std::ios::fmtflags(0), std::ios::floatfield);
@@ -198,11 +198,8 @@
 		ostr << position().x();
 		ostr.width(9);
-                ostr.precision(2);
 		ostr << position().y();
 		ostr.width(9);
-                ostr.precision(2);
 		ostr << position().z();
 		ostr.width(9);
-                ostr.precision(2);
 		ostr << position().t();
                 ostr.setf(std::ios::fmtflags(0), std::ios::floatfield);
@@ -249,4 +246,7 @@
 	    ostr << **part2 << std::endl;
 	}
+        // restore the stream state
+        ostr.flags(orig);
+        ostr.precision(prec);
     }
 
Index: trunk/Utilities/HepMC/src/HeavyIon.cc
===================================================================
--- trunk/Utilities/HepMC/src/HeavyIon.cc	(revision 571)
+++ trunk/Utilities/HepMC/src/HeavyIon.cc	(revision 571)
@@ -0,0 +1,144 @@
+//--------------------------------------------------------------------------
+//
+// HeavyIon.cc
+// Author:  Lynn Garren
+//
+// Implement operator >> and operator <<
+//
+// ----------------------------------------------------------------------
+
+#include <iostream>
+#include <ostream>
+#include <istream>
+#include <sstream>
+
+#include "HeavyIon.h"
+#include "StreamHelpers.h"
+#include "IO_Exception.h"
+
+namespace HepMC {
+
+/// Write the contents of HeavyIon to an output stream.
+/// GenEvent stores a pointer to a HeavyIon.
+std::ostream & operator << (std::ostream & os, HeavyIon const * ion)
+{
+    if ( !os ) {
+	std::cerr << "HeavyIon output stream !os, "
+		  << " setting badbit" << std::endl;
+	os.clear(std::ios::badbit); 
+	return os;
+    }
+    os << 'H';
+    // HeavyIon* is set to 0 by default
+    if ( !ion  ) {
+	detail::output( os, 0 );
+	detail::output( os, 0 );
+	detail::output( os, 0 );
+	detail::output( os, 0 );
+	detail::output( os, 0 );
+	detail::output( os, 0 );
+	detail::output( os, 0 );
+	detail::output( os, 0 );
+	detail::output( os, 0 );
+	detail::output( os, 0. );
+	detail::output( os, 0. );
+	detail::output( os, 0. );
+	detail::output( os, 0. );
+	detail::output( os,'\n');
+	return os;
+    }
+    //
+    detail::output( os, ion->Ncoll_hard() );
+    detail::output( os, ion->Npart_proj() );
+    detail::output( os, ion->Npart_targ() );
+    detail::output( os, ion->Ncoll() );
+    detail::output( os, ion->spectator_neutrons() );
+    detail::output( os, ion->spectator_protons() );
+    detail::output( os, ion->N_Nwounded_collisions() );
+    detail::output( os, ion->Nwounded_N_collisions() );
+    detail::output( os, ion->Nwounded_Nwounded_collisions() );
+    detail::output( os, ion->impact_parameter() );
+    detail::output( os, ion->event_plane_angle() );
+    detail::output( os, ion->eccentricity() );
+    detail::output( os, ion->sigma_inel_NN() );
+    detail::output( os,'\n');
+
+    return os;
+}
+
+/// Read the contents of HeavyIon from an input stream.
+/// GenEvent stores a pointer to a HeavyIon.
+std::istream & operator >> (std::istream & is, HeavyIon * ion)
+{
+    // make sure the stream is valid
+    if ( !is ) { 
+      std::cerr << "HeavyIon input stream setting badbit." << std::endl;
+      is.clear(std::ios::badbit); 
+      return is; 
+    }
+    // get the HeavyIon line
+    std::string line;
+    std::getline(is,line);
+    std::istringstream iline(line);
+    std::string firstc;
+    iline >> firstc;
+    // test to be sure the next entry is of type "H" 
+    if( firstc != "H" ) { 
+	std::cerr << "HeavyIon input stream invalid line type: " 
+	          << firstc << std::endl;
+	// The most likely problem is that we have found a HepMC block line
+	throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    } 
+    // read values into temp variables, then create a new HeavyIon object
+    int nh =0, np =0, nt =0, nc =0, 
+        neut = 0, prot = 0, nw =0, nwn =0, nwnw =0;
+    float impact = 0., plane = 0., xcen = 0., inel = 0.; 
+    iline >> nh ;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> np ;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> nt ;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> nc ;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> neut ;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> prot;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> nw ;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> nwn ;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> nwnw ;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> impact ;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> plane ;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> xcen ;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    iline >> inel;
+    if(!iline) throw IO_Exception("HeavyIon input stream encounterd invalid data");
+    if( nh == 0 ) {
+        return is;
+    }
+    
+    ion->set_Ncoll_hard(nh);
+    ion->set_Npart_proj(np);
+    ion->set_Npart_targ(nt);
+    ion->set_Ncoll(nc);
+    ion->set_spectator_neutrons(neut);
+    ion->set_spectator_protons(prot);
+    ion->set_N_Nwounded_collisions(nw);
+    ion->set_Nwounded_N_collisions(nwn);
+    ion->set_Nwounded_Nwounded_collisions(nwnw);
+    ion->set_impact_parameter(impact);
+    ion->set_event_plane_angle(plane);
+    ion->set_eccentricity(xcen);
+    ion->set_sigma_inel_NN(inel);
+
+    return is;
+}
+
+
+} // HepMC
Index: trunk/Utilities/HepMC/src/IO_AsciiParticles.cc
===================================================================
--- trunk/Utilities/HepMC/src/IO_AsciiParticles.cc	(revision 349)
+++ trunk/Utilities/HepMC/src/IO_AsciiParticles.cc	(revision 571)
@@ -11,5 +11,4 @@
 #include "IO_AsciiParticles.h"
 #include "GenEvent.h"
-#include "ParticleDataTable.h"
 #include "Version.h"
 
Index: trunk/Utilities/HepMC/src/IO_GenEvent.cc
===================================================================
--- trunk/Utilities/HepMC/src/IO_GenEvent.cc	(revision 349)
+++ trunk/Utilities/HepMC/src/IO_GenEvent.cc	(revision 571)
@@ -8,10 +8,7 @@
 
 #include "IO_GenEvent.h"
+#include "IO_Exception.h"
 #include "GenEvent.h"
-#include "ParticleDataTable.h"
-#include "HeavyIon.h"
-#include "PdfInfo.h"
-#include "CommonIO.h"
-#include "Version.h"
+#include "StreamHelpers.h"
 
 namespace HepMC {
@@ -23,22 +20,16 @@
       m_istr(0),
       m_iostr(0),
-      m_finished_first_event_io(false),
       m_have_file(false),
-      m_common_io()
+      m_error_type(IO_Exception::OK),
+      m_error_message()
     {
 	if ( (m_mode&std::ios::out && m_mode&std::ios::in) ||
 	     (m_mode&std::ios::app && m_mode&std::ios::in) ) {
-	    std::cerr << "IO_GenEvent::IO_GenEvent Error, open of file requested "
-		      << "of input AND output type. Not allowed. Closing file."
-		      << std::endl;
+            m_error_type = IO_Exception::InputAndOutput;
+	    m_error_message ="IO_GenEvent::IO_GenEvent Error, open of file requested of input AND output type. Not allowed. Closing file.";
+	    std::cerr << m_error_message << std::endl;
 	    m_file.close();
 	    return;
 	}
-	// precision 16 (# digits following decimal point) is the minimum that
-	//  will capture the full information stored in a double
-	m_file.precision(16);
-	// we use decimal to store integers, because it is smaller than hex!
-	m_file.setf(std::ios::dec,std::ios::basefield);
-	m_file.setf(std::ios::scientific,std::ios::floatfield);
 	// now we set the streams
 	m_iostr = &m_file;
@@ -46,8 +37,10 @@
 	    m_istr = &m_file;
 	    m_ostr = NULL;
+	    detail::establish_input_stream_info(m_file);
 	}
 	if ( m_mode&std::ios::out ) {
 	    m_ostr = &m_file;
 	    m_istr = NULL;
+	    detail::establish_output_stream_info(m_file);
 	}
 	m_have_file = true;
@@ -59,8 +52,10 @@
       m_istr(&istr),
       m_iostr(&istr),
-      m_finished_first_event_io(false),
       m_have_file(false),
-      m_common_io()
-    { }
+      m_error_type(IO_Exception::OK),
+      m_error_message()
+    { 
+        detail::establish_input_stream_info( istr );
+    }
 
     IO_GenEvent::IO_GenEvent( std::ostream & ostr )
@@ -68,23 +63,23 @@
       m_istr(0),
       m_iostr(&ostr),
-      m_finished_first_event_io(false),
       m_have_file(false),
-      m_common_io()
-    {
-	// precision 16 (# digits following decimal point) is the minimum that
-	//  will capture the full information stored in a double
-	m_ostr->precision(16);
-	// we use decimal to store integers, because it is smaller than hex!
-	m_ostr->setf(std::ios::dec,std::ios::basefield);
-	m_ostr->setf(std::ios::scientific,std::ios::floatfield);
-    }
+      m_error_type(IO_Exception::OK),
+      m_error_message()
+   {
+        detail::establish_output_stream_info( ostr );
+   }
 
     IO_GenEvent::~IO_GenEvent() {
-	write_end_listing();
+    	if ( m_ostr != NULL ) {
+	    write_HepMC_IO_block_end(*m_ostr);
+	}
 	if(m_have_file) m_file.close();
     }
 
-    void IO_GenEvent::use_input_units( Units::MomentumUnit mom, Units::LengthUnit len ) {
-	m_common_io.use_input_units(mom,len);
+    void IO_GenEvent::use_input_units( Units::MomentumUnit mom, 
+                                       Units::LengthUnit len ) {
+        if( m_istr != NULL ) {
+            set_input_units( *m_istr, mom, len );
+	}
     }
 
@@ -99,4 +94,51 @@
     }
 
+    void IO_GenEvent::precision( int size )  { 
+        if( size > 16 ) { 
+	    std::cerr << "IO_GenEvent::precision Error, "
+		      << "precision is greater than 16. "
+		      << "Not allowed. Using default precision of 16."
+		      << std::endl;
+            size = 16;
+	}
+        if(m_ostr) {
+	    m_ostr->precision(size);
+	}
+    }
+	
+    bool IO_GenEvent::fill_next_event( GenEvent* evt ){
+	//
+	// reset error type
+        m_error_type = IO_Exception::OK;
+	//
+	// test that evt pointer is not null
+	if ( !evt ) {
+            m_error_type = IO_Exception::NullEvent;
+	    m_error_message = "IO_GenEvent::fill_next_event error - passed null event.";
+	    std::cerr << m_error_message << std::endl;
+	    return false;
+	}
+	// make sure the stream is good, and that it is in input mode
+	if ( !(*m_istr) ) return false;
+	if ( !m_istr ) {
+            m_error_type = IO_Exception::WrongFileType;
+	    m_error_message = "HepMC::IO_GenEvent::fill_next_event attempt to read from output file.";
+	    std::cerr << m_error_message << std::endl;
+	    return false;
+	}
+	// use streaming input
+        try {
+	    *m_istr >> *evt;
+	}
+        catch (IO_Exception& e) {
+            m_error_type = IO_Exception::InvalidData;
+	    m_error_message = e.what();
+	    evt->clear();
+ 	    return false;
+        }
+	if( evt->is_valid() ) return true;
+	return false;
+    }
+
     void IO_GenEvent::write_event( const GenEvent* evt ) {
 	/// Writes evt to output stream. It does NOT delete the event after writing.
@@ -105,129 +147,15 @@
 	if ( !evt  ) return;
 	if ( m_ostr == NULL ) {
-	    std::cerr << "HepMC::IO_GenEvent::write_event "
-		      << " attempt to write to input file." << std::endl;
+            m_error_type = IO_Exception::WrongFileType;
+	    m_error_message = "HepMC::IO_GenEvent::write_event attempt to write to input file.";
+	    std::cerr << m_error_message << std::endl;
 	    return;
 	}
 	//
 	// write event listing key before first event only.
-	if ( !m_finished_first_event_io ) {
-	    m_finished_first_event_io = 1;
-	    *m_ostr << "\n" << "HepMC::Version " << versionName();
-	    *m_ostr << "\n";
-	    m_common_io.write_IO_GenEvent_Key(*m_ostr);
-	}
-	//
-	// output the event data including the number of primary vertices
-	//  and the total number of vertices
-	std::vector<long> random_states = evt->random_states();
-	*m_ostr << 'E';
-	output( evt->event_number() );
-	output( evt->mpi() );
-	output( evt->event_scale() );
-	output( evt->alphaQCD() );
-	output( evt->alphaQED() );
-	output( evt->signal_process_id() );
-	output(   ( evt->signal_process_vertex() ?
-		    evt->signal_process_vertex()->barcode() : 0 )   );
-	output( evt->vertices_size() ); // total number of vertices.
-	write_beam_particles( evt->beam_particles() );
-	output( (int)random_states.size() );
-	for ( std::vector<long>::iterator rs = random_states.begin(); 
-	      rs != random_states.end(); ++rs ) {
-	    output( *rs );
-	}
-	output( (int)evt->weights().size() );
-	for ( WeightContainer::const_iterator w = evt->weights().begin(); 
-	      w != evt->weights().end(); ++w ) {
-	    output( *w );
-	}
-	output('\n');
-	write_unit_info( evt );
-	write_heavy_ion( evt->heavy_ion() );
-	write_pdf_info( evt->pdf_info() );
-	//
-	// Output all of the vertices - note there is no real order.
-	for ( GenEvent::vertex_const_iterator v = evt->vertices_begin();
-	      v != evt->vertices_end(); ++v ) {
-	    write_vertex( *v );
-	}
-    }
-
-    bool IO_GenEvent::fill_next_event( GenEvent* evt ){
-	//
-	//
-	// test that evt pointer is not null
-	if ( !evt ) {
-	    std::cerr 
-		<< "IO_GenEvent::fill_next_event error - passed null event." 
-		<< std::endl;
-	    return false;
-	}
-	// make sure the stream is good, and that it is in input mode
-	if ( !(*m_istr) ) return false;
-	if ( !m_istr ) {
-	    std::cerr << "HepMC::IO_GenEvent::fill_next_event "
-		      << " attempt to read from output file." << std::endl;
-	    return false;
-	}
-	//
-	// search for event listing key before first event only.
-	//
-	// skip through the file just after first occurence of the start_key
-	int iotype = 0;
-	if ( !m_finished_first_event_io ) {
-	    iotype = m_common_io.find_file_type(*m_istr);
-	    if( iotype <= 0 ) {
-		std::cerr << "IO_GenEvent::fill_next_event start key not found "
-			  << "setting badbit." << std::endl;
-		m_istr->clear(std::ios::badbit); 
-		return false;
-	    }
-	    m_finished_first_event_io = 1;
-	}
-	//
-	// test to be sure the next entry is of type "E" then ignore it
-	if ( !(*m_istr) ) { 
-		std::cerr << "IO_GenEvent::fill_next_event end of stream found "
-			  << "setting badbit." << std::endl;
-		m_istr->clear(std::ios::badbit); 
-		return false;
-	}
-	if ( !(*m_istr) || m_istr->peek()!='E' ) { 
-	    // if the E is not the next entry, then check to see if it is
-	    // the end event listing key - if yes, search for another start key
-	    int ioendtype = m_common_io.find_end_key(*m_istr);
-	    if ( ioendtype == iotype ) {
-		iotype = m_common_io.find_file_type(*m_istr);
-		if( iotype <= 0 ) {
-		    // this is the only case where we set an EOF state
-		    m_istr->clear(std::ios::eofbit);
-		    return false;
-		}
-	    } else if ( ioendtype > 0 ) {
-		std::cerr << "IO_GenEvent::fill_next_event end key does not match start key "
-			  << "setting badbit." << std::endl;
-		m_istr->clear(std::ios::badbit); 
-		return false;
-	    } else {
-		std::cerr << "IO_GenEvent::fill_next_event end key not found "
-			  << "setting badbit." << std::endl;
-		m_istr->clear(std::ios::badbit); 
-		return false;
-	    }
-	}
-	m_istr->ignore();
-	// call the appropriate read method
-	if( m_common_io.io_type() == gen ) {
-	    return m_common_io.read_io_genevent(m_istr, evt);
-	} else if( m_common_io.io_type() == ascii ) { 
-	    return m_common_io.read_io_ascii(m_istr, evt);
-	} else if( m_common_io.io_type() == extascii ) { 
-	    return m_common_io.read_io_extendedascii(m_istr, evt);
-	} else if( m_common_io.io_type() == ascii_pdt ) { 
-	} else if( m_common_io.io_type() == extascii_pdt ) { 
-	}
-	// should not get to this statement
-	return false;
+	write_HepMC_IO_block_begin(*m_ostr);
+	// explicit cast is necessary
+	GenEvent e = *evt;
+	*m_ostr << e ;
     }
 
@@ -236,313 +164,15 @@
 	if ( !(*m_ostr) ) return;
 	if ( m_ostr == NULL ) {
-	    std::cerr << "HepMC::IO_GenEvent::write_comment "
-		      << " attempt to write to input file." << std::endl;
+            m_error_type = IO_Exception::WrongFileType;
+	    m_error_message = "HepMC::IO_GenEvent::write_event attempt to write to input file.";
+	    std::cerr << m_error_message << std::endl;
 	    return;
 	}
 	// write end of event listing key if events have already been written
-	write_end_listing();
+	write_HepMC_IO_block_end(*m_ostr);
 	// insert the comment key before the comment
 	*m_ostr << "\n" << "HepMC::IO_GenEvent-COMMENT\n";
 	*m_ostr << comment << std::endl;
     }
-
-    void IO_GenEvent::write_vertex( GenVertex* v ) {
-	// assumes mode has already been checked
-	if ( !v || !(*m_ostr) ) {
-	    std::cerr << "IO_GenEvent::write_vertex !v||!(*m_ostr), "
-		      << "v="<< v << " setting badbit" << std::endl;
-	    m_ostr->clear(std::ios::badbit); 
-	    return;
-	}
-	// 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;
-	}
-	//
-	*m_ostr << 'V';
-	output( v->barcode() ); // v's unique identifier
-	output( v->id() );
-	output( v->position().x() );
-	output( v->position().y() );
-	output( v->position().z() );
-	output( v->position().t() );
-	output( num_orphans_in );
-	output( (int)v->particles_out_size() );
-	output( (int)v->weights().size() );
-	for ( WeightContainer::iterator w = v->weights().begin(); 
-	      w != v->weights().end(); ++w ) {
-	    output( *w );
-	}
-	output('\n');
-	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( *p2 );
-	    }
-	}
-	for ( GenVertex::particles_out_const_iterator p3 
-		  = v->particles_out_const_begin();
-	      p3 != v->particles_out_const_end(); ++p3 ) {
-	    write_particle( *p3 );
-	}
-    }
-
-    void IO_GenEvent::write_beam_particles( 
-        std::pair<HepMC::GenParticle *,HepMC::GenParticle *> pr ) {
-        GenParticle* p = pr.first;
-	//m_file << 'B';
-	if(!p) {
-	   output( 0 );
-	} else {
-	   output( p->barcode() );
-	}
-        p = pr.second;
-	if(!p) {
-	   output( 0 );
-	} else {
-	   output( p->barcode() );
-	}
-    }
-
-    void IO_GenEvent::write_heavy_ion( HeavyIon const * ion ) {
-	// assumes mode has already been checked
-	if ( !(*m_ostr) ) {
-	    std::cerr << "IO_GenEvent::write_heavy_ion !(*m_ostr), "
-		      << " setting badbit" << std::endl;
-	    m_ostr->clear(std::ios::badbit); 
-	    return;
-	}
-	*m_ostr << 'H';
-	// HeavyIon* is set to 0 by default
-	if ( !ion  ) {
-	    output( 0 );
-	    output( 0 );
-	    output( 0 );
-	    output( 0 );
-	    output( 0 );
-	    output( 0 );
-	    output( 0 );
-	    output( 0 );
-	    output( 0 );
-	    output( 0. );
-	    output( 0. );
-	    output( 0. );
-	    output( 0. );
-	    output('\n');
-	    return;
-	}
-	//
-	output( ion->Ncoll_hard() );
-	output( ion->Npart_proj() );
-	output( ion->Npart_targ() );
-	output( ion->Ncoll() );
-	output( ion->spectator_neutrons() );
-	output( ion->spectator_protons() );
-	output( ion->N_Nwounded_collisions() );
-	output( ion->Nwounded_N_collisions() );
-	output( ion->Nwounded_Nwounded_collisions() );
-	output( ion->impact_parameter() );
-	output( ion->event_plane_angle() );
-	output( ion->eccentricity() );
-	output( ion->sigma_inel_NN() );
-	output('\n');
-    }
-
-    void IO_GenEvent::write_pdf_info( PdfInfo const * pdf ) {
-	// assumes mode has already been checked
-	if ( !(*m_ostr) ) {
-	    std::cerr << "IO_GenEvent::write_pdf_info !(*m_ostr), "
-		      << " setting badbit" << std::endl;
-	    m_ostr->clear(std::ios::badbit); 
-	    return;
-	}
-	*m_ostr << 'F';
-	// PdfInfo* is set to 0 by default
-	if ( !pdf ) {
-	    output( 0 );
-	    output( 0 );
-	    output( 0. );
-	    output( 0. );
-	    output( 0. );
-	    output( 0. );
-	    output( 0. );
-	    output( 0 );
-	    output( 0 );
-	    output('\n');
-	    return;
-	}
-	//
-	output( pdf->id1() );
-	output( pdf->id2() );
-	output( pdf->x1() );
-	output( pdf->x2() );
-	output( pdf->scalePDF() );
-	output( pdf->pdf1() );
-	output( pdf->pdf2() );
-	output( pdf->pdf_id1() );
-	output( pdf->pdf_id2() );
-	output('\n');
-    }
-
-    void IO_GenEvent::write_unit_info( const GenEvent* evt ) {
-	if ( !(*m_ostr) ) {
-	    std::cerr << "IO_GenEvent::write_unit_info !(*m_ostr), "
-		      << " setting badbit" << std::endl;
-	    m_ostr->clear(std::ios::badbit); 
-	    return;
-	}
-	// could write enums here, but strings are more readable
-	*m_ostr << "U " << name(evt->momentum_unit());
-	*m_ostr << " " << name(evt->length_unit());
-	output('\n');
-    }
-
-    void IO_GenEvent::write_particle( GenParticle* p ) {
-	// assumes mode has already been checked
-	if ( !p || !(*m_ostr) ) {
-	    std::cerr << "IO_GenEvent::write_particle !p||!(*m_ostr), "
-		      << "p="<< p << " setting badbit" << std::endl;
-	    m_ostr->clear(std::ios::badbit); 
-	    return;
-	}
-	*m_ostr << 'P';
-	output( p->barcode() );
-	output( p->pdg_id() );
-	output( p->momentum().px() );
-	output( p->momentum().py() );
-	output( p->momentum().pz() );
-	output( p->momentum().e() );
-	output( p->generated_mass() );
-	output( p->status() );
-	output( p->polarization().theta() );
-	output( p->polarization().phi() );
-	// since end_vertex is oftentimes null, this CREATES a null vertex
-	// in the map
-	output(   ( p->end_vertex() ? p->end_vertex()->barcode() : 0 )  );
-	*m_ostr << ' ' << p->flow() << "\n";
-    }
-
-    void IO_GenEvent::write_particle_data( const ParticleData* pdata ) {
-	// assumes mode has already been checked
-	if ( !pdata || !(*m_ostr) ) {
-	    std::cerr << "IO_GenEvent::write_particle_data !pdata||!(*m_ostr), "
-		      << "pdata="<< pdata << " setting badbit" << std::endl;
-	    m_ostr->clear(std::ios::badbit); 
-	    return;
-	}
-	*m_ostr << 'D';
-	output( pdata->pdg_id() );
-	output( pdata->charge() );
-	output( pdata->mass() );
-	output( pdata->clifetime() );
-	output( (int)(pdata->spin()*2.+.1) );
-	// writes the first 21 characters starting with 0
-	*m_ostr << " " << pdata->name().substr(0,21) << "\n";
-    }
-
-    HeavyIon* IO_GenEvent::read_heavy_ion()
-    {
-	// assumes mode has already been checked
-	//
-	// test to be sure the next entry is of type "H" then ignore it
-	if ( !(*m_istr) || m_istr->peek()!='H' ) {
-	    std::cerr << "IO_GenEvent::read_heavy_ion setting badbit." << std::endl;
-	    m_istr->clear(std::ios::badbit); 
-	    return false;
-	} 
-	m_istr->ignore();
-	// read values into temp variables, then create a new HeavyIon object
-	int nh =0, np =0, nt =0, nc =0, 
-            neut = 0, prot = 0, nw =0, nwn =0, nwnw =0;
-	float impact = 0., plane = 0., xcen = 0., inel = 0.; 
-	*m_istr >> nh >> np >> nt >> nc >> neut >> prot
-	       >> nw >> nwn >> nwnw >> impact >> plane >> xcen >> inel;
-	m_istr->ignore(2,'\n');
-	if( nh == 0 ) return false;
-	HeavyIon* ion = new HeavyIon(nh, np, nt, nc, neut, prot,
-	                             nw, nwn, nwnw, 
-				     impact, plane, xcen, inel );
-	//
-	return ion;
-    }
-
-    GenParticle* IO_GenEvent::read_particle(
-	TempParticleMap& particle_to_end_vertex ){
-	// assumes mode has already been checked
-	//
-	// test to be sure the next entry is of type "P" then ignore it
-	if ( !(*m_istr) || m_istr->peek()!='P' ) { 
-	    std::cerr << "IO_GenEvent::read_particle setting badbit." 
-		      << std::endl;
-	    m_istr->clear(std::ios::badbit); 
-	    return false;
-	} 
-	m_istr->ignore();
-	//
-	// declare variables to be read in to, and read everything except flow
-	double px = 0., py = 0., pz = 0., e = 0., m = 0., theta = 0., phi = 0.;
-	int bar_code = 0, id = 0, status = 0, end_vtx_code = 0, flow_size = 0;
-	*m_istr >> bar_code >> id >> px >> py >> pz >> e >> m >> status 
-	       >> theta >> phi >> end_vtx_code >> flow_size;
-	//
-	// read flow patterns if any exist
-	Flow flow;
-	int code_index, code;
-	for ( int i = 1; i <= flow_size; ++i ) {
-	    *m_istr >> code_index >> code;
-	    flow.set_icode( code_index,code);
-	}
-	m_istr->ignore(2,'\n'); // '\n' at end of entry
-	GenParticle* p = new GenParticle( FourVector(px,py,pz,e), 
-				    id, status, flow, 
-				    Polarization(theta,phi) );
-        p->set_generated_mass( m );
-	p->suggest_barcode( bar_code );
-	//
-	// all particles are connected to their end vertex separately 
-	// after all particles and vertices have been created - so we keep
-	// a map of all particles that have end vertices
-	if ( end_vtx_code != 0 ) {
-	    particle_to_end_vertex.addEndParticle(p,end_vtx_code);
-	}
-	return p;
-    }
-
-    ParticleData* IO_GenEvent::read_particle_data( ParticleDataTable* pdt ) {
-	// assumes mode has already been checked
-	//
-	// test to be sure the next entry is of type "D" then ignore it
-	if ( !(*m_istr) || m_istr->peek()!='D' ) return false;
-	m_istr->ignore();
-	//
-	// read values into temp variables then create new ParticleData object
-	char its_name[22];
-	int its_id = 0, its_spin = 0;  
-	double its_charge = 0, its_mass = 0, its_clifetime = 0;
-	*m_istr >> its_id >> its_charge >> its_mass 
-	       >> its_clifetime >> its_spin;
-	m_istr->ignore(1); // eat the " "
-	m_istr->getline( its_name, 22, '\n' );
-	ParticleData* pdata = new ParticleData( its_name, its_id, its_charge, 
-						its_mass, its_clifetime, 
-						double(its_spin)/2.);
-	pdt->insert(pdata);
-	return pdata;
-    }
-
-    bool IO_GenEvent::write_end_listing() {
-    	if ( m_finished_first_event_io && (m_ostr != NULL) ) {
-	    m_common_io.write_IO_GenEvent_End(*m_ostr);
-	    *m_ostr << std::flush;
-	    m_finished_first_event_io = 0;
-	    return true;
-	}
-	return false;
-    }
 	
 } // HepMC
Index: trunk/Utilities/HepMC/src/PdfInfo.cc
===================================================================
--- trunk/Utilities/HepMC/src/PdfInfo.cc	(revision 571)
+++ trunk/Utilities/HepMC/src/PdfInfo.cc	(revision 571)
@@ -0,0 +1,121 @@
+//--------------------------------------------------------------------------
+//
+// PdfInfo.cc
+// Author:  Lynn Garren
+//
+// Implement operator >> and operator <<
+//
+// ----------------------------------------------------------------------
+
+#include <iostream>
+#include <ostream>
+#include <istream>
+#include <sstream>
+
+#include "PdfInfo.h"
+#include "StreamHelpers.h"
+#include "IO_Exception.h"
+
+namespace HepMC {
+
+std::ostream & operator << ( std::ostream & os, PdfInfo const * pdf)
+{
+    if ( !os ) {
+	std::cerr << "operator << for PdfInfo: !os, "
+		  << " setting badbit" << std::endl;
+	os.clear(std::ios::badbit); 
+	return os;
+    }
+    os << 'F';
+    // PdfInfo* is set to 0 by default
+    if ( !pdf ) {
+	detail::output( os, 0 );
+	detail::output( os, 0 );
+	detail::output( os, 0. );
+	detail::output( os, 0. );
+	detail::output( os, 0. );
+	detail::output( os, 0. );
+	detail::output( os, 0. );
+	detail::output( os, 0 );
+	detail::output( os, 0 );
+	detail::output( os,'\n');
+	return os;
+    }
+    //
+    detail::output( os, pdf->id1() );
+    detail::output( os, pdf->id2() );
+    detail::output( os, pdf->x1() );
+    detail::output( os, pdf->x2() );
+    detail::output( os, pdf->scalePDF() );
+    detail::output( os, pdf->pdf1() );
+    detail::output( os, pdf->pdf2() );
+    detail::output( os, pdf->pdf_id1() );
+    detail::output( os, pdf->pdf_id2() );
+    detail::output( os,'\n');
+
+    return os;
+}
+
+std::istream & operator >> (std::istream & is, PdfInfo * pdf)
+{
+    // make sure the stream is valid
+    if ( !is ) {
+	std::cerr << "PdfInfo input stream setting badbit." << std::endl;
+	is.clear(std::ios::badbit); 
+	return is;
+    } 
+    //
+    // get the PdfInfo line
+    std::string line;
+    std::getline(is,line);
+    std::istringstream iline(line);
+    std::string firstc;
+    iline >> firstc;
+    // test to be sure the next entry is of type "F" then ignore it
+    if ( firstc != "F" ) {
+	std::cerr << "PdfInfo input stream invalid line type: " 
+	          << firstc << std::endl;
+	// this is non-recoverable, so throw here 
+	throw IO_Exception("PdfInfo input stream encounterd invalid data");
+    } 
+    // read values into temp variables, then create a new PdfInfo object
+    int id1 =0, id2 =0, pdf_id1=0, pdf_id2=0;
+    double  x1 = 0., x2 = 0., scale = 0., pdf1 = 0., pdf2 = 0.; 
+    iline >> id1 ;
+    if(!iline) throw IO_Exception("PdfInfo input stream encounterd invalid data");
+    // check now for empty PdfInfo line
+    if( id1 == 0 ) return is;
+    // continue reading
+    iline >> id2 ;
+    if(!iline) throw IO_Exception("PdfInfo input stream encounterd invalid data");
+    iline >> x1 ;
+    if(!iline) throw IO_Exception("PdfInfo input stream encounterd invalid data");
+    iline >> x2 ;
+    if(!iline) throw IO_Exception("PdfInfo input stream encounterd invalid data");
+    iline >> scale ;
+    if(!iline) throw IO_Exception("PdfInfo input stream encounterd invalid data");
+    iline >> pdf1 ;
+    if(!iline) throw IO_Exception("PdfInfo input stream encounterd invalid data");
+    iline >> pdf2;
+    if(!iline) throw IO_Exception("PdfInfo input stream encounterd invalid data");
+    // check to see if we are at the end of the line
+    if( !iline.eof() ) {
+        iline >> pdf_id1 ;
+        if(!iline) throw IO_Exception("PdfInfo input stream encounterd invalid data");
+	iline >> pdf_id2;
+        if(!iline) throw IO_Exception("PdfInfo input stream encounterd invalid data");
+    }
+    pdf->set_id1( id1 );
+    pdf->set_id2( id2 );
+    pdf->set_pdf_id1( pdf_id1 );
+    pdf->set_pdf_id2( pdf_id2 );
+    pdf->set_x1( x1 );
+    pdf->set_x2( x2 );
+    pdf->set_scalePDF( scale );
+    pdf->set_pdf1( pdf1 );
+    pdf->set_pdf2( pdf2 );
+
+    return is;
+}
+
+} // HepMC
Index: trunk/Utilities/HepMC/src/Polarization.cc
===================================================================
--- trunk/Utilities/HepMC/src/Polarization.cc	(revision 349)
+++ trunk/Utilities/HepMC/src/Polarization.cc	(revision 571)
@@ -9,17 +9,26 @@
 namespace HepMC {
 
+    Polarization::Polarization( )
+    : m_theta( 0. ),
+      m_phi( 0. ),
+      m_defined( false )
+    { }
+
     Polarization::Polarization( double theta, double phi )
     : m_theta( valid_theta(theta) ),
-      m_phi  ( valid_phi(phi) )
+      m_phi  ( valid_phi(phi) ),
+      m_defined( true )
     { }
 
     Polarization::Polarization( const Polarization& inpolar )
     : m_theta( valid_theta( inpolar.theta() ) ),
-      m_phi  ( valid_phi(   inpolar.phi()   ) )
+      m_phi  ( valid_phi(   inpolar.phi()   ) ),
+      m_defined( inpolar.is_defined() )
     { }
 
     Polarization::Polarization( const ThreeVector& vec3in ) 
     : m_theta( valid_theta( vec3in.theta() ) ),
-      m_phi  ( valid_phi(   vec3in.phi()   ) )
+      m_phi  ( valid_phi(   vec3in.phi()   ) ),
+      m_defined( true )
     { }
 
@@ -28,4 +37,5 @@
  	std::swap( m_theta, other.m_theta );
  	std::swap( m_phi,   other.m_phi   );
+	std::swap( m_defined, other.m_defined );
     }
 
@@ -64,8 +74,19 @@
 	return m_phi = valid_phi( phi );
     }
+    
+    bool Polarization::is_defined( ) const {
+        return m_defined;
+    }
+    
+    void Polarization::set_undefined() {
+        m_defined = false;
+	m_theta = 0.;
+	m_phi = 0.;
+    }
 
     void Polarization::set_theta_phi( double theta, double phi ) {
 	set_theta( theta );
 	set_phi( phi ) ;
+	m_defined = true;
     }
 
@@ -73,4 +94,5 @@
 	set_theta( vec3in.theta() );
 	set_phi( vec3in.phi() );
+	m_defined = true;
 	return vec3in;
     }
Index: trunk/Utilities/HepMC/src/StreamHelpers.cc
===================================================================
--- trunk/Utilities/HepMC/src/StreamHelpers.cc	(revision 571)
+++ trunk/Utilities/HepMC/src/StreamHelpers.cc	(revision 571)
@@ -0,0 +1,123 @@
+//--------------------------------------------------------------------------
+//
+// GenEventStreamIO.cc
+// Author:  Lynn Garren
+//
+// helper functions used by streaming IO
+//
+// ----------------------------------------------------------------------
+
+#include <ostream>
+#include <istream>
+#include <sstream>
+
+#include "GenVertex.h"
+#include "GenParticle.h"
+#include "StreamHelpers.h"
+#include "IO_Exception.h"
+
+namespace HepMC {
+
+namespace detail {
+
+std::istream & read_vertex( std::istream & is, 
+                            TempParticleMap & particle_to_end_vertex, 
+			    GenVertex * v )
+{
+    //
+    // make sure the stream is valid
+    if ( !is ) {
+	std::cerr << "StreamHelpers::detail::read_vertex setting badbit." << std::endl;
+	is.clear(std::ios::badbit); 
+	return is;
+    } 
+    //
+    // get the vertex line
+    std::string line;
+    std::getline(is,line);
+    std::istringstream iline(line);
+    std::string firstc;
+    iline >> firstc;
+    //
+    // test to be sure the next entry is of type "V" 
+    if ( firstc != "V"  ) {
+	std::cerr << "StreamHelpers::detail::read_vertex invalid line type: " 
+	          << firstc << std::endl;
+	std::cerr << "StreamHelpers::detail::read_vertex setting badbit." << std::endl;
+	is.clear(std::ios::badbit); 
+	return is;
+    } 
+    // read values into temp variables, then create a new GenVertex object
+    int identifier =0, id =0, num_orphans_in =0, 
+        num_particles_out = 0, weights_size = 0;
+    double x = 0., y = 0., z = 0., t = 0.; 
+    iline >> identifier ;
+    if(!iline) { throw IO_Exception("read_vertex input stream encounterd invalid data"); }
+    iline >> id ;
+    if(!iline) { throw IO_Exception("read_vertex input stream encounterd invalid data"); }
+    iline >> x ;
+    if(!iline) { throw IO_Exception("read_vertex input stream encounterd invalid data"); }
+    iline >> y ;
+    if(!iline) { throw IO_Exception("read_vertex input stream encounterd invalid data"); }
+    iline >> z ;
+    if(!iline) { throw IO_Exception("read_vertex input stream encounterd invalid data"); }
+    iline >> t;
+    if(!iline) { throw IO_Exception("read_vertex input stream encounterd invalid data"); }
+    iline >> num_orphans_in ;
+    if(!iline) { throw IO_Exception("read_vertex input stream encounterd invalid data"); }
+    iline >> num_particles_out ;
+    if(!iline) { throw IO_Exception("read_vertex input stream encounterd invalid data"); }
+    iline >> weights_size;
+    if(!iline) { throw IO_Exception("read_vertex input stream encounterd invalid data"); }
+    WeightContainer weights(weights_size);
+    for ( int i1 = 0; i1 < weights_size; ++i1 ) {
+        iline >> weights[i1];
+        if(!iline) { throw IO_Exception("read_vertex input stream encounterd invalid data"); }
+    }
+    v->set_position( FourVector(x,y,z,t) );
+    v->set_id( id );
+    v->weights() = weights;
+    v->suggest_barcode( identifier );
+    //
+    // read and create the associated particles. outgoing particles are
+    //  added to their production vertices immediately, while incoming
+    //  particles are added to a map and handled later.
+    for ( int i2 = 1; i2 <= num_orphans_in; ++i2 ) {
+        GenParticle* p1 = new GenParticle( ); 
+	detail::read_particle(is,particle_to_end_vertex,p1);
+    }
+    for ( int i3 = 1; i3 <= num_particles_out; ++i3 ) {
+        GenParticle* p2 = new GenParticle( ); 
+	detail::read_particle(is,particle_to_end_vertex,p2);
+	v->add_particle_out( p2 );
+    }
+
+    return is;
+}
+
+std::istream & find_event_end( std::istream & is ) {
+    // since there is no end of event flag, 
+    // read one line at time until we find the next event 
+    // or the end of event block
+    // don't throw until we find the end of the event
+    std::string line, firstc;
+    while ( is ) { 
+	is >> firstc;
+        if( firstc=="E" ) {	// next event
+	    is.unget();
+            throw IO_Exception("input stream encountered invalid data");
+	    return is;
+	} else if( firstc.size() > 1 ) { // no more events in this block
+            throw IO_Exception("input stream encountered invalid data, now at end of event block");
+	    return is;
+	}
+        std::getline(is,line);
+    }
+    // the stream is bad 
+    throw IO_Exception("input stream encountered invalid data, stream is now corrupt");
+    return is;
+}
+
+} // detail
+
+} // HepMC
Index: trunk/Utilities/HepMC/src/StreamInfo.cc
===================================================================
--- trunk/Utilities/HepMC/src/StreamInfo.cc	(revision 571)
+++ trunk/Utilities/HepMC/src/StreamInfo.cc	(revision 571)
@@ -0,0 +1,59 @@
+//--------------------------------------------------------------------------
+//
+// StreamInfo.cc
+// Author:  Lynn Garren
+//
+// ----------------------------------------------------------------------
+
+#include <string>
+#include "StreamInfo.h"
+
+namespace HepMC {
+
+StreamInfo::StreamInfo( )
+: m_finished_first_event_io(false),
+  m_io_genevent_start("HepMC::IO_GenEvent-START_EVENT_LISTING"),
+  m_io_ascii_start("HepMC::IO_Ascii-START_EVENT_LISTING"),
+  m_io_extendedascii_start("HepMC::IO_ExtendedAscii-START_EVENT_LISTING"),
+  m_io_genevent_end("HepMC::IO_GenEvent-END_EVENT_LISTING"),
+  m_io_ascii_end("HepMC::IO_Ascii-END_EVENT_LISTING"),
+  m_io_extendedascii_end("HepMC::IO_ExtendedAscii-END_EVENT_LISTING"),
+  m_io_ascii_pdt_start("HepMC::IO_Ascii-START_PARTICLE_DATA"),
+  m_io_extendedascii_pdt_start("HepMC::IO_ExtendedAscii-START_PARTICLE_DATA"),
+  m_io_ascii_pdt_end("HepMC::IO_Ascii-END_PARTICLE_DATA"),
+  m_io_extendedascii_pdt_end("HepMC::IO_ExtendedAscii-END_PARTICLE_DATA"),
+  m_io_type(0),
+  m_has_key(true),
+  m_io_momentum_unit(Units::default_momentum_unit()),
+  m_io_position_unit(Units::default_length_unit()),
+  m_stream_id(m_stream_counter),
+  m_reading_event_header(false)
+{
+    ++m_stream_counter;
+}
+
+/// static counter 
+unsigned int StreamInfo::m_stream_counter = 0; 
+
+void StreamInfo::use_input_units( Units::MomentumUnit mom, Units::LengthUnit len ) {
+    m_io_momentum_unit = mom;
+    m_io_position_unit = len;
+}
+
+void StreamInfo::set_io_type( int io ) {
+    m_io_type = io;
+}
+
+void StreamInfo::set_has_key( bool io ) {
+    m_has_key = io;
+}
+
+bool StreamInfo::reading_event_header() {
+    return m_reading_event_header;
+}
+
+void StreamInfo::set_reading_event_header(bool tf) {
+    m_reading_event_header = tf;
+}
+
+} // HepMC
Index: trunk/Utilities/HepMC/src/Units.cc
===================================================================
--- trunk/Utilities/HepMC/src/Units.cc	(revision 349)
+++ trunk/Utilities/HepMC/src/Units.cc	(revision 571)
@@ -52,5 +52,4 @@
     // if this function fails to compile, rerun configure using --with-length_units
     LengthUnit default_length_unit() {
-      //return @HEPMC_DEFAULT_LEN_UNIT@ ;
       return MM ;
     }
@@ -58,6 +57,5 @@
     // if this function fails to compile, rerun configure using --with-momentum_units
     MomentumUnit default_momentum_unit() {
-      //return @HEPMC_DEFAULT_MOM_UNIT@ ;
-      return GEV;
+      return GEV ;
     }
 
Index: trunk/Utilities/HepMC/src/WeightContainer.cc
===================================================================
--- trunk/Utilities/HepMC/src/WeightContainer.cc	(revision 571)
+++ trunk/Utilities/HepMC/src/WeightContainer.cc	(revision 571)
@@ -0,0 +1,139 @@
+//////////////////////////////////////////////////////////////////////////
+// Matt.Dobbs@Cern.CH, November 2000, refer to:
+// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
+// High Energy Physics", Computer Physics Communications (to be published).
+//
+// Container for the Weights associated with an event or vertex.
+// Basically just an interface to STL vector with extra map-like attributes
+//////////////////////////////////////////////////////////////////////////
+
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+#include <vector>
+#include <string>
+#include <map>
+#include <stdexcept>
+
+#include "WeightContainer.h"
+
+namespace HepMC {
+
+WeightContainer::WeightContainer( size_type n, double value ) 
+    : m_weights(n,value), m_names()
+{ set_default_names(n); }
+
+WeightContainer::WeightContainer( const std::vector<double>& wgts )
+    : m_weights(wgts), m_names()
+{ set_default_names(size()); }
+
+void WeightContainer::set_default_names( size_type n )
+{
+    // internal program used by the constructors
+    std::ostringstream name;
+    for ( size_type count = 0; count<n; ++count ) 
+    { 
+	name.str(std::string());
+	name << count;
+	m_names[name.str()] = count;
+    }
+}
+
+void WeightContainer::push_back( const double& value) 
+{ 
+    size_type count = m_weights.size();
+    m_weights.push_back(value); 
+    std::ostringstream name;
+    name << count;
+    m_names[name.str()] = count;
+}
+
+void WeightContainer::pop_back() 
+{
+    // this needs to remove the last entry in the vector 
+    // and ALSO the associated map entry
+    size_type vit = size() - 1;
+    for ( map_iterator m = m_names.begin(); m != m_names.end(); ++m ) 
+    { 
+        if( m->second == vit ) { 
+	    m_names.erase(m->first); 
+	    continue;
+	}
+    }
+    m_weights.pop_back(); 
+}
+
+double& WeightContainer::operator[]( const std::string& s ) 
+{ 
+    const_map_iterator m = m_names.find(s);
+    if( m != m_names.end() ) {
+        return m_weights[m->second]; 
+    }
+    // doesn't exist - have to create it
+    size_type count = m_weights.size();
+    m_weights.push_back(0); 
+    m_names[s] = count;
+    return m_weights.back(); 
+}
+
+
+const double& WeightContainer::operator[]( const std::string& s ) const
+{ 
+    const_map_iterator m = m_names.find(s);
+    if( m != m_names.end() ) {
+        return m_weights[m->second]; 
+    }
+    // doesn't exist and we cannot create it
+    // note that std::map does not support this (const) operator
+    // throw an appropriate error, we choose the error thrown by std::vector
+    throw std::out_of_range("const WeightContainer::operator[] ERROR: string "+s+" not found in  WeightContainer" );
+}
+
+bool WeightContainer::operator==( const WeightContainer & other ) const
+{
+   if( size() != other.size() ) { return false; }
+   if( m_names != other.m_names ) { return false; }
+   if( m_weights != other.m_weights ) { return false; }
+   return true;
+}
+
+bool WeightContainer::operator!=( const WeightContainer & other ) const
+{
+   return !(*this == other );
+}
+
+bool WeightContainer::has_key( const std::string& s ) const
+{
+    // look up the name in the map
+    return m_names.find(s) != m_names.end();
+}
+
+void WeightContainer::print( std::ostream& ostr ) const 
+{ 
+    // print a name, weight pair
+    for ( const_map_iterator m = map_begin(); m != map_end(); ++m )
+    {
+	ostr << "(" << m->first << "," << m_weights[m->second] << ") ";
+    }
+    ostr << std::endl; 
+}
+
+void WeightContainer::write( std::ostream& ostr ) const 
+{ 
+    size_type count = 0;
+    for ( const_iterator w = begin(); w != end(); ++w ) 
+    { 
+	std::string name;
+	for ( const_map_iterator m = map_begin(); m != map_end(); ++m )
+	{
+	    if( m->second == count ) name = m->first;
+	}
+	ostr << "Weight " << std::setw(4) << count 
+	     << " with name " << std::setw(10) <<  name
+	     << " is " << *w << std::endl;
+	++count;
+    }
+}
+
+} // HepMC
+
Index: trunk/Utilities/HepMC/src/filterEvent.cc
===================================================================
--- trunk/Utilities/HepMC/src/filterEvent.cc	(revision 571)
+++ trunk/Utilities/HepMC/src/filterEvent.cc	(revision 571)
@@ -0,0 +1,82 @@
+// from Andy Buckley
+
+#include "GenEvent.h"
+
+  void filterEvent(HepMC::GenEvent* ge) {
+    // We treat beam particles a bit specially
+    const std::pair<HepMC::GenParticle*, HepMC::GenParticle*> beams = ge->beam_particles();
+
+    // Make list of non-physical particle entries
+    std::vector<HepMC::GenParticle*> unphys_particles;
+    for (HepMC::GenEvent::particle_const_iterator pi = ge->particles_begin(); 
+         pi != ge->particles_end(); ++pi) {
+      // Beam particles might not have status = 4, but we want them anyway
+      if (beams.first == *pi || beams.second == *pi) continue;
+      // Filter by status
+      const int status = (*pi)->status();
+      if (status != 1 && status != 2 && status != 4) {
+        unphys_particles.push_back(*pi);
+      }      
+    }
+
+    // Remove each unphysical particle from the list
+    while (unphys_particles.size()) {
+      HepMC::GenParticle* gp = unphys_particles.back();
+
+      // Get start and end vertices
+      HepMC::GenVertex* vstart = gp->production_vertex();
+      HepMC::GenVertex* vend = gp->end_vertex();
+
+      if (vend == vstart) {
+        // Deal with loops
+        vstart->remove_particle(gp);
+      } else {
+
+        // Connect decay particles from end vertex to start vertex
+        /// @todo Have to build a list, since the GV::add_particle_out method modifies the end vertex!
+        if (vend && vend->particles_out_size()) {
+          std::vector<HepMC::GenParticle*> end_particles;
+          for (HepMC::GenVertex::particles_out_const_iterator gpe = vend->particles_out_const_begin(); 
+               gpe != vend->particles_out_const_end(); ++gpe) {
+            end_particles.push_back(*gpe);
+          }
+          // Reset production vertices of child particles to bypass unphysical particle
+          for (std::vector<HepMC::GenParticle*>::const_iterator gpe = end_particles.begin();
+               gpe != end_particles.end(); ++gpe) {
+            //std::cout << vstart << ", " << vend << std::endl;
+            if (vstart) vstart->add_particle_out(*gpe);
+          }
+        } else {
+          // If null end_vertex... stable unphysical particle?
+        }
+
+        // Delete unphysical particle and its orphaned end vertex
+        delete vend;
+        if (vstart) {
+          delete vstart->remove_particle(gp);
+        }// else {
+        /// @todo Why does this cause an error?
+        //  delete gp;
+        //}
+      }
+
+      // Remove deleted particle from list
+      unphys_particles.pop_back();
+      //std::cout << unphys_particles.size() << std::endl;
+    }
+
+    // Delete any orphaned vertices
+    std::vector<HepMC::GenVertex*> orphaned_vtxs;
+    for (HepMC::GenEvent::vertex_const_iterator vi = ge->vertices_begin(); 
+         vi != ge->vertices_end(); ++vi) {
+      if ((*vi)->particles_in_size() == 0 && (*vi)->particles_out_size() == 0) {
+        orphaned_vtxs.push_back(*vi);
+      }
+    }
+    while (orphaned_vtxs.size()) {
+      delete orphaned_vtxs.back();
+      orphaned_vtxs.pop_back();
+    }
+  }
+
+
