TEvtGen/HepMC/GenEventStreamIO.cc

   1 //--------------------------------------------------------------------------
   2 //
   3 // GenEventStreamIO.cc
   4 // Author:  Lynn Garren
   5 //
   6 // Implement operator >> and operator <<
   7 //
   8 // ----------------------------------------------------------------------
   9
  10 #include <iostream>
  11 #include <ostream>
  12 #include <istream>
  13 #include <sstream>
  14
  15 #include "HepMC/GenEvent.h"
  16 #include "HepMC/GenCrossSection.h"
  17 #include "HepMC/StreamInfo.h"
  18 #include "HepMC/StreamHelpers.h"
  19 #include "HepMC/Version.h"
  20 #include "HepMC/IO_Exception.h"
  21
  22 namespace HepMC {
  23
  24 // ------------------------- local methods ----------------
  25
  26 /// This method is called by the stream destructor.
  27 /// It does cleanup on stored user data (StreamInfo)
  28 /// and is registered by the first call to get_stream_info().
  29 void HepMCStreamCallback(std::ios_base::event e, std::ios_base& b, int i)
  30 {
  31   // only clean up if the stream object is going away.
  32   if(i!=0 && e!= std::ios_base::erase_event) return;
  33
  34   // retrieve the pointer to the object
  35   StreamInfo* hd = (StreamInfo*)b.pword(i);
  36   b.pword(i) = 0;
  37   b.iword(i) = 0;
  38 #ifdef HEPMC_DEBUG
  39   // the following line is just for sanity checking
  40   if(hd) std::cerr << "deleted StreamInfo " << hd->stream_id() << "\n";
  41 #endif
  42   delete hd;
  43 }
  44
  45 // ------------------------- iomanip ----------------
  46
  47 /// A custom iomanip that allows us to store and access user data (StreamInfo)
  48 /// associated with the stream.
  49 /// This method creates the StreamInfo object the first time it is called.
  50 template <class IO>
  51 StreamInfo& get_stream_info(IO& iost)
  52 {
  53   if(iost.iword(0) == 0)
  54     {
  55       // make sure we add the callback if this is the first time through
  56       iost.iword(0)=1;
  57       iost.register_callback(&HepMCStreamCallback, 0);
  58       // this is our special "context" record.
  59       // there is one of these at the head of each IO block.
  60       // allocate room for a StreamInfo in the userdata area
  61       iost.pword(0) = new StreamInfo;
  62 #ifdef HEPMC_DEBUG
  63       // the following line is just for sanity checking
  64       std::cerr << "created StreamInfo " << ((StreamInfo*)iost.pword(0))->stream_id() << "\n";
  65 #endif
  66     }
  67   return *(StreamInfo*)iost.pword(0);
  68 }
  69
  70 // ------------------------- GenEvent member functions ----------------
  71
  72 std::ostream& GenEvent::write( std::ostream& os )
  73 {
  74     /// Writes evt to an output stream.
  75
  76     //
  77     StreamInfo & info = get_stream_info(os);
  78     //
  79     // if this is the first event, set precision
  80     if ( !info.finished_first_event() ) {
  81         // precision 16 (# digits following decimal point) is the minimum that
  82         //  will capture the full information stored in a double
  83         //  However, we let the user set precision, since that is the expected functionality
  84         // we use decimal to store integers, because it is smaller than hex!
  85         os.setf(std::ios::dec,std::ios::basefield);
  86         os.setf(std::ios::scientific,std::ios::floatfield);
  87         //
  88         info.set_finished_first_event(true);
  89     }
  90     //
  91     // output the event data including the number of primary vertices
  92     //  and the total number of vertices
  93     //std::vector<long> random_states = random_states();
  94     os << 'E';
  95     detail::output( os, event_number() );
  96     detail::output( os, mpi() );
  97     detail::output( os, event_scale() );
  98     detail::output( os, alphaQCD() );
  99     detail::output( os, alphaQED() );
 100     detail::output( os, signal_process_id() );
 101     detail::output( os,   ( signal_process_vertex() ?
 102                 signal_process_vertex()->barcode() : 0 )   );
 103     detail::output( os, vertices_size() ); // total number of vertices.
 104     write_beam_particles( os, beam_particles() );
 105     // random state
 106     detail::output( os, (int)m_random_states.size() );
 107     for ( std::vector<long>::iterator rs = m_random_states.begin();
 108           rs != m_random_states.end(); ++rs ) {
 109          detail::output( os, *rs );
 110     }
 111     // weights
 112     // we need to iterate over the map so that the weights printed
 113     // here will be in the same order as the names printed next
 114     os << ' ' << (int)weights().size() ;
 115     for ( WeightContainer::const_map_iterator w = weights().map_begin();
 116           w != weights().map_end(); ++w ) {
 117         detail::output( os, m_weights[w->second] );
 118     }
 119     detail::output( os,'\n');
 120     // now add names for weights
 121     // note that this prints a new line if and only if the weight container
 122     // is not empty
 123     if ( ! weights().empty() ) {
 124         os << "N " << weights().size() << " " ;
 125         for ( WeightContainer::const_map_iterator w = weights().map_begin();
 126               w != weights().map_end(); ++w ) {
 127             detail::output( os,'"');
 128             os << w->first;
 129             detail::output( os,'"');
 130             detail::output( os,' ');
 131         }
 132         detail::output( os,'\n');
 133     }
 134     //
 135     // Units
 136     os << "U " << name(momentum_unit());
 137     os << " " << name(length_unit());
 138     detail::output( os,'\n');
 139     //
 140     // write GenCrossSection if it has been set
 141     if( m_cross_section ) m_cross_section->write(os);
 142     //
 143     // write HeavyIon and PdfInfo if they have been set
 144     if( m_heavy_ion ) os << heavy_ion() ;
 145     if( m_pdf_info ) os << pdf_info() ;
 146     //
 147     // Output all of the vertices - note there is no real order.
 148     for ( GenEvent::vertex_const_iterator v = vertices_begin();
 149           v != vertices_end(); ++v ) {
 150         write_vertex(os, *v);
 151     }
 152     return os;
 153 }
 154
 155 std::istream& GenEvent::read( std::istream& is )
 156 {
 157     /// read a GenEvent from streaming input
 158     //
 159     StreamInfo & info = get_stream_info(is);
 160     clear();
 161     //
 162     // search for event listing key before first event only.
 163     if ( !info.finished_first_event() ) {
 164         //
 165         find_file_type(is);
 166         info.set_finished_first_event(true);
 167     }
 168     //
 169     // make sure the stream is good
 170     if ( !is ) {
 171         std::cerr << "streaming input: end of stream found "
 172                   << "setting badbit." << std::endl;
 173         is.clear(std::ios::badbit);
 174         return is;
 175     }
 176
 177     //
 178     // test to be sure the next entry is of type "E" then ignore it
 179     if ( is.peek()!='E' ) {
 180         // if the E is not the next entry, then check to see if it is
 181         // the end event listing key - if yes, search for another start key
 182         int ioendtype;
 183         find_end_key(is,ioendtype);
 184         if ( ioendtype == info.io_type() ) {
 185             find_file_type(is);
 186             // are we at the end of the file?
 187             if( !is ) return is;
 188         } else if ( ioendtype > 0 ) {
 189             std::cerr << "streaming input: end key does not match start key "
 190                       << "setting badbit." << std::endl;
 191             is.clear(std::ios::badbit);
 192             return is;
 193         } else if ( !info.has_key() ) {
 194             find_file_type(is);
 195             // are we at the end of the file?
 196             if( !is ) return is;
 197         } else {
 198             std::cerr << "streaming input: end key not found "
 199                       << "setting badbit." << std::endl;
 200             is.clear(std::ios::badbit);
 201             return is;
 202         }
 203     }
 204
 205     int signal_process_vertex = 0;
 206     int num_vertices = 0, bp1 = 0, bp2 = 0;
 207     bool units_line = false;
 208     // OK - now ready to start reading the event, so set the header flag
 209     info.set_reading_event_header(true);
 210     // The flag will be set to false when we reach the end of the header
 211     while(info.reading_event_header()) {
 212         switch(is.peek()) {
 213             case 'E':
 214             {   // deal with the event line
 215                 process_event_line( is, num_vertices, bp1, bp2, signal_process_vertex );
 216             } break;
 217             case 'N':
 218             {   // get weight names
 219                 read_weight_names( is );
 220             } break;
 221             case 'U':
 222             {   // get unit information if it exists
 223                 units_line = true;
 224                 if( info.io_type() == gen ) {
 225                     read_units( is );
 226                 }
 227             } break;
 228             case 'C':
 229             {   // we have a GenCrossSection line
 230                 // create cross section
 231                 GenCrossSection xs;
 232                 // check for invalid data
 233                 try {
 234                     // read the line
 235                     xs.read(is);
 236                 }
 237                 catch (IO_Exception& e) {
 238                     detail::find_event_end( is );
 239                 }
 240                 if(xs.is_set()) {
 241                     set_cross_section( xs );
 242                 }
 243             } break;
 244             case 'H':
 245             {   // we have a HeavyIon line OR an unexpected HepMC... line
 246                 if( info.io_type() == gen || info.io_type() == extascii ) {
 247                     // get HeavyIon
 248                     HeavyIon ion;
 249                     // check for invalid data
 250                     try {
 251                         is >> &ion;
 252                     }
 253                     catch (IO_Exception& e) {
 254                         detail::find_event_end( is );
 255                     }
 256                     if(ion.is_valid()) {
 257                         set_heavy_ion( ion );
 258                     }
 259                 }
 260             } break;
 261             case 'F':
 262             {   // we have a PdfInfo line
 263                 if( info.io_type() == gen || info.io_type() == extascii ) {
 264                     // get PdfInfo
 265                     PdfInfo pdf;
 266                     // check for invalid data
 267                     try {
 268                         is >> &pdf;
 269                     }
 270                     catch (IO_Exception& e) {
 271                         detail::find_event_end( is );
 272                     }
 273                     if(pdf.is_valid()) {
 274                         set_pdf_info( pdf );
 275                     }
 276                 }
 277             } break;
 278             case 'V':
 279             {
 280                 // this should be the first vertex line - exit this loop
 281                 info.set_reading_event_header(false);
 282             } break;
 283             case 'P':
 284             {   // we should not find this line
 285                 std::cerr << "streaming input: found unexpected line P" << std::endl;
 286                 info.set_reading_event_header(false);
 287             } break;
 288             default:
 289                 // ignore everything else
 290                 break;
 291         } // switch on line type
 292     } // while reading_event_header
 293     // before proceeding - did we find a units line?
 294     if( !units_line ) {
 295         use_units( info.io_momentum_unit(),
 296                        info.io_position_unit() );
 297     }
 298     //
 299     // the end vertices of the particles are not connected until
 300     //  after the event is read --- we store the values in a map until then
 301     TempParticleMap particle_to_end_vertex;
 302     //
 303     // read in the vertices
 304     for ( int iii = 1; iii <= num_vertices; ++iii ) {
 305         GenVertex* v = new GenVertex();
 306         try {
 307             detail::read_vertex(is,particle_to_end_vertex,v);
 308         }
 309         catch (IO_Exception& e) {
 310             for( TempParticleMap::orderIterator it = particle_to_end_vertex.order_begin();
 311                  it != particle_to_end_vertex.order_end(); ++it ) {
 312                 GenParticle* p = it->second;
 313                 // delete particles only if they are not already owned by a vertex
 314                 if( p->production_vertex() ) {
 315                 } else if( p->end_vertex() ) {
 316                 } else {
 317                      delete p;
 318                 }
 319             }
 320             delete v;
 321             detail::find_event_end( is );
 322         }
 323         add_vertex( v );
 324     }
 325     // set the signal process vertex
 326     if ( signal_process_vertex ) {
 327         set_signal_process_vertex(
 328             barcode_to_vertex(signal_process_vertex) );
 329     }
 330     //
 331     // last connect particles to their end vertices
 332     GenParticle* beam1(0);
 333     GenParticle* beam2(0);
 334     for ( TempParticleMap::orderIterator pmap
 335               = particle_to_end_vertex.order_begin();
 336           pmap != particle_to_end_vertex.order_end(); ++pmap ) {
 337         GenParticle* p =  pmap->second;
 338         int vtx = particle_to_end_vertex.end_vertex( p );
 339         GenVertex* itsDecayVtx = barcode_to_vertex(vtx);
 340         if ( itsDecayVtx ) itsDecayVtx->add_particle_in( p );
 341         else {
 342             std::cerr << "read_io_genevent: ERROR particle points"
 343                       << " to null end vertex. " <<std::endl;
 344         }
 345         // also look for the beam particles
 346         if( p->barcode() == bp1 ) beam1 = p;
 347         if( p->barcode() == bp2 ) beam2 = p;
 348     }
 349     set_beam_particles(beam1,beam2);
 350     return is;
 351 }
 352
 353 // ------------------------- operator << and operator >> ----------------
 354
 355 std::ostream & operator << (std::ostream & os, GenEvent & evt)
 356 {
 357     /// Writes evt to an output stream.
 358     evt.write(os);
 359     return os;
 360 }
 361
 362 std::istream & operator >> (std::istream & is, GenEvent & evt)
 363 {
 364     evt.read(is);
 365     return is;
 366 }
 367
 368 // ------------------------- set units ----------------
 369
 370 std::istream & set_input_units(std::istream & is,
 371                                Units::MomentumUnit mom,
 372                                Units::LengthUnit len )
 373 {
 374     //
 375     StreamInfo & info = get_stream_info(is);
 376     info.use_input_units( mom, len );
 377     return is;
 378 }
 379
 380 // ------------------------- begin and end block lines ----------------
 381
 382 std::ostream & write_HepMC_IO_block_begin(std::ostream & os )
 383 {
 384     //
 385     StreamInfo & info = get_stream_info(os);
 386
 387     if( !info.finished_first_event() ) {
 388     os << "\n" << "HepMC::Version " << versionName();
 389     os << "\n";
 390     os << info.IO_GenEvent_Key() << "\n";
 391     }
 392     return os;
 393 }
 394
 395 std::ostream & write_HepMC_IO_block_end(std::ostream & os )
 396 {
 397     //
 398     StreamInfo & info = get_stream_info(os);
 399
 400     if( info.finished_first_event() ) {
 401         os << info.IO_GenEvent_End() << "\n";
 402         os << std::flush;
 403     }
 404     return os;
 405 }
 406
 407 std::istream & GenEvent::process_event_line( std::istream & is,
 408                                              int & num_vertices,
 409                                              int & bp1, int & bp2,
 410                                              int & signal_process_vertex )
 411 {
 412     //
 413     if ( !is ) {
 414         std::cerr << "GenEvent::process_event_line setting badbit." << std::endl;
 415         is.clear(std::ios::badbit);
 416         return is;
 417     }
 418     //
 419     StreamInfo & info = get_stream_info(is);
 420     std::string line;
 421     std::getline(is,line);
 422     std::istringstream iline(line);
 423     std::string firstc;
 424     iline >> firstc;
 425     //
 426     // read values into temp variables, then fill GenEvent
 427     int event_number = 0, signal_process_id = 0,
 428         random_states_size = 0, nmpi = -1;
 429     double eventScale = 0, alpha_qcd = 0, alpha_qed = 0;
 430     iline >> event_number;
 431     if(!iline) detail::find_event_end( is );
 432     if( info.io_type() == gen || info.io_type() == extascii ) {
 433         iline >> nmpi;
 434         if(!iline) detail::find_event_end( is );
 435         set_mpi( nmpi );
 436     }
 437     iline >> eventScale ;
 438     if(!iline) detail::find_event_end( is );
 439     iline >> alpha_qcd ;
 440     if(!iline) detail::find_event_end( is );
 441     iline >> alpha_qed;
 442     if(!iline) detail::find_event_end( is );
 443     iline >> signal_process_id ;
 444     if(!iline) detail::find_event_end( is );
 445     iline >> signal_process_vertex;
 446     if(!iline) detail::find_event_end( is );
 447     iline >> num_vertices;
 448     if(!iline) detail::find_event_end( is );
 449     if( info.io_type() == gen || info.io_type() == extascii ) {
 450         iline >> bp1 ;
 451         if(!iline) detail::find_event_end( is );
 452         iline >> bp2;
 453         if(!iline) detail::find_event_end( is );
 454     }
 455     iline >> random_states_size;
 456     if(!iline) detail::find_event_end( is );
 457     std::vector<long> random_states(random_states_size);
 458     for ( int i = 0; i < random_states_size; ++i ) {
 459         iline >> random_states[i];
 460         if(!iline) detail::find_event_end( is );
 461     }
 462     WeightContainer::size_type weights_size = 0;
 463     iline >> weights_size;
 464     if(!iline) detail::find_event_end( is );
 465     std::vector<double> wgt(weights_size);
 466     for ( WeightContainer::size_type ii = 0; ii < weights_size; ++ii ) {
 467         iline >> wgt[ii];
 468         if(!iline) detail::find_event_end( is );
 469     }
 470     // weight names will be added later if they exist
 471     if( weights_size > 0 ) m_weights = wgt;
 472     //
 473     // fill signal_process_id, event_number, random_states, etc.
 474     set_signal_process_id( signal_process_id );
 475     set_event_number( event_number );
 476     set_random_states( random_states );
 477     set_event_scale( eventScale );
 478     set_alphaQCD( alpha_qcd );
 479     set_alphaQED( alpha_qed );
 480     //
 481     return is;
 482 }
 483
 484 std::istream & GenEvent::read_weight_names( std::istream & is )
 485 {
 486     // now check for a named weight line
 487     if ( !is ) {
 488         std::cerr << "GenEvent::read_weight_names setting badbit." << std::endl;
 489         is.clear(std::ios::badbit);
 490         return is;
 491     }
 492     // Test to be sure the next entry is of type "N"
 493     // If we have no named weight line, this is not an error
 494     // releases prior to 2.06.00 do not have named weights
 495     if ( is.peek() !='N') {
 496         return is;
 497     }
 498     // now get this line and process it
 499     std::string line;
 500     std::getline(is,line);
 501     std::istringstream wline(line);
 502     std::string firstc;
 503     WeightContainer::size_type name_size = 0;
 504     wline >> firstc >> name_size;
 505     if(!wline) detail::find_event_end( is );
 506     if( firstc != "N") {
 507         std::cout << "debug: first character of named weights is " << firstc << std::endl;
 508         std::cout << "debug: We should never get here" << std::endl;
 509         is.clear(std::ios::badbit);
 510         return is;
 511     }
 512     if( m_weights.size() != name_size ) {
 513         std::cout << "debug: weight sizes do not match "<< std::endl;
 514         std::cout << "debug: weight vector size is " << m_weights.size() << std::endl;
 515         std::cout << "debug: weight name size is " << name_size << std::endl;
 516         is.clear(std::ios::badbit);
 517         return is;
 518     }
 519     std::string name;
 520     std::string::size_type i1 = line.find("\"");
 521     std::string::size_type i2;
 522     std::string::size_type len = line.size();
 523     WeightContainer namedWeight;
 524     for ( WeightContainer::size_type ii = 0; ii < name_size; ++ii ) {
 525         // weight names may contain blanks
 526         if(i1 >= len) {
 527             std::cout << "debug: attempting to read past the end of the named weight line " << std::endl;
 528             std::cout << "debug: We should never get here" << std::endl;
 529             std::cout << "debug: Looking for the end of this event" << std::endl;
 530             detail::find_event_end( is );
 531         }
 532         i2 = line.find("\"",i1+1);
 533         name = line.substr(i1+1,i2-i1-1);
 534         namedWeight[name] = m_weights[ii];
 535         i1 = line.find("\"",i2+1);
 536     }
 537     m_weights = namedWeight;
 538     return is;
 539 }
 540
 541 std::istream & GenEvent::read_units( std::istream & is )
 542 {
 543     //
 544     if ( !is ) {
 545         std::cerr << "GenEvent::read_units setting badbit." << std::endl;
 546         is.clear(std::ios::badbit);
 547         return is;
 548     }
 549     //
 550     StreamInfo & info = get_stream_info(is);
 551     // test to be sure the next entry is of type "U" then ignore it
 552     // if we have no units, this is not an error
 553     // releases prior to 2.04.00 did not write unit information
 554     if ( is.peek() !='U') {
 555         use_units( info.io_momentum_unit(),
 556                        info.io_position_unit() );
 557         return is;
 558     }
 559     is.ignore();        // ignore the first character in the line
 560     std::string mom, pos;
 561     is >> mom >> pos;
 562     is.ignore(1);      // eat the extra whitespace
 563     use_units(mom,pos);
 564     //
 565     return is;
 566 }
 567
 568 std::istream & GenEvent::find_file_type( std::istream & istr )
 569 {
 570     //
 571     // make sure the stream is good
 572     if ( !istr ) return istr;
 573
 574     //
 575     StreamInfo & info = get_stream_info(istr);
 576
 577     // if there is no input block line, then we assume this stream
 578     // is in the IO_GenEvent format
 579     if ( istr.peek()=='E' ) {
 580         info.set_io_type( gen );
 581         info.set_has_key(false);
 582         return istr;
 583     }
 584
 585     std::string line;
 586     while ( std::getline(istr,line) ) {
 587         //
 588         // search for event listing key before first event only.
 589         //
 590         if( line == info.IO_GenEvent_Key() ) {
 591             info.set_io_type( gen );
 592             info.set_has_key(true);
 593             return istr;
 594         } else if( line == info.IO_Ascii_Key() ) {
 595             info.set_io_type( ascii );
 596             info.set_has_key(true);
 597             return istr;
 598         } else if( line == info.IO_ExtendedAscii_Key() ) {
 599             info.set_io_type( extascii );
 600             info.set_has_key(true);
 601             return istr;
 602         } else if( line == info.IO_Ascii_PDT_Key() ) {
 603             info.set_io_type( ascii_pdt );
 604             info.set_has_key(true);
 605             return istr;
 606         } else if( line == info.IO_ExtendedAscii_PDT_Key() ) {
 607             info.set_io_type( extascii_pdt );
 608             info.set_has_key(true);
 609             return istr;
 610         }
 611     }
 612     info.set_io_type( 0 );
 613     info.set_has_key(false);
 614     return istr;
 615 }
 616
 617 std::istream & GenEvent::find_end_key( std::istream & istr, int & iotype )
 618 {
 619     iotype = 0;
 620     // peek at the first character before proceeding
 621     if( istr.peek()!='H' ) return istr;
 622     //
 623     // we only check the next line
 624     std::string line;
 625     std::getline(istr,line);
 626     //
 627     StreamInfo & info = get_stream_info(istr);
 628     //
 629     // check to see if this is an end key
 630     if( line == info.IO_GenEvent_End() ) {
 631         iotype = gen;
 632     } else if( line == info.IO_Ascii_End() ) {
 633         iotype = ascii;
 634     } else if( line == info.IO_ExtendedAscii_End() ) {
 635         iotype = extascii;
 636     } else if( line == info.IO_Ascii_PDT_End() ) {
 637         iotype = ascii_pdt;
 638     } else if( line == info.IO_ExtendedAscii_PDT_End() ) {
 639         iotype = extascii_pdt;
 640     }
 641     if( iotype != 0 && info.io_type() != iotype ) {
 642         std::cerr << "GenEvent::find_end_key: iotype keys have changed" << std::endl;
 643     } else {
 644         return istr;
 645     }
 646     //
 647     // if we get here, then something has gotten badly confused
 648     std::cerr << "GenEvent::find_end_key: MALFORMED INPUT" << std::endl;
 649     istr.clear(std::ios::badbit);
 650     return istr;
 651 }
 652
 653 std::ostream & establish_output_stream_info( std::ostream & os )
 654 {
 655     StreamInfo & info = get_stream_info(os);
 656     if ( !info.finished_first_event() ) {
 657         // precision 16 (# digits following decimal point) is the minimum that
 658         //  will capture the full information stored in a double
 659         os.precision(16);
 660         // we use decimal to store integers, because it is smaller than hex!
 661         os.setf(std::ios::dec,std::ios::basefield);
 662         os.setf(std::ios::scientific,std::ios::floatfield);
 663     }
 664     return os;
 665 }
 666
 667 std::istream & establish_input_stream_info( std::istream & is )
 668 {
 669     StreamInfo & info = get_stream_info(is);
 670     if ( !info.finished_first_event() ) {
 671         // precision 16 (# digits following decimal point) is the minimum that
 672         //  will capture the full information stored in a double
 673         is.precision(16);
 674         // we use decimal to store integers, because it is smaller than hex!
 675         is.setf(std::ios::dec,std::ios::basefield);
 676         is.setf(std::ios::scientific,std::ios::floatfield);
 677     }
 678     return is;
 679 }
 680
 681
 682 // ------------------------- helper functions ----------------
 683
 684 namespace detail {
 685
 686 // The functions defined here need to use get_stream_info
 687
 688 std::istream & read_particle( std::istream & is,
 689                               TempParticleMap & particle_to_end_vertex,
 690                               GenParticle * p )
 691 {
 692     // get the next line
 693     std::string line;
 694     std::getline(is,line);
 695     std::istringstream iline(line);
 696     std::string firstc;
 697     iline >> firstc;
 698     if( firstc != "P" ) {
 699         std::cerr << "StreamHelpers::detail::read_particle invalid line type: "
 700                   << firstc << std::endl;
 701         std::cerr << "StreamHelpers::detail::read_particle setting badbit."
 702                   << std::endl;
 703         is.clear(std::ios::badbit);
 704         return is;
 705     }
 706     //
 707     StreamInfo & info = get_stream_info(is);
 708     //testHepMC.cc
 709     // declare variables to be read in to, and read everything except flow
 710     double px = 0., py = 0., pz = 0., e = 0., m = 0., theta = 0., phi = 0.;
 711     int bar_code = 0, id = 0, status = 0, end_vtx_code = 0, flow_size = 0;
 712     // check that the input stream is still OK after reading item
 713     iline >> bar_code ;
 714     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 715     iline >> id ;
 716     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 717     iline >> px ;
 718     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 719     iline >> py ;
 720     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 721     iline >> pz ;
 722     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 723     iline >> e ;
 724     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 725     if( info.io_type() != ascii ) {
 726         iline >> m ;
 727         if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 728     }
 729     iline >> status ;
 730     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 731     iline >> theta ;
 732     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 733     iline >> phi ;
 734     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 735     iline >> end_vtx_code ;
 736     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 737     iline >> flow_size;
 738     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 739     //
 740     // read flow patterns if any exist
 741     Flow flow;
 742     int code_index, code;
 743     for ( int i = 1; i <= flow_size; ++i ) {
 744         iline >> code_index >> code;
 745         if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
 746         flow.set_icode( code_index,code);
 747     }
 748     p->set_momentum( FourVector(px,py,pz,e) );
 749     p->set_pdg_id( id );
 750     p->set_status( status );
 751     p->set_flow( flow );
 752     p->set_polarization( Polarization(theta,phi) );
 753     if( info.io_type() == ascii ) {
 754         p->set_generated_mass( p->momentum().m() );
 755     } else {
 756         p->set_generated_mass( m );
 757     }
 758     p->suggest_barcode( bar_code );
 759     //
 760     // all particles are connected to their end vertex separately
 761     // after all particles and vertices have been created - so we keep
 762     // a map of all particles that have end vertices
 763     if ( end_vtx_code != 0 ) {
 764         particle_to_end_vertex.addEndParticle(p,end_vtx_code);
 765     }
 766     return is;
 767 }
 768
 769 std::ostream & establish_output_stream_info( std::ostream & os )
 770 {
 771     StreamInfo & info = get_stream_info(os);
 772     if ( !info.finished_first_event() ) {
 773         // precision 16 (# digits following decimal point) is the minimum that
 774         //  will capture the full information stored in a double
 775         os.precision(16);
 776         // we use decimal to store integers, because it is smaller than hex!
 777         os.setf(std::ios::dec,std::ios::basefield);
 778         os.setf(std::ios::scientific,std::ios::floatfield);
 779     }
 780     return os;
 781 }
 782
 783 std::istream & establish_input_stream_info( std::istream & is )
 784 {
 785     StreamInfo & info = get_stream_info(is);
 786     if ( !info.finished_first_event() ) {
 787         // precision 16 (# digits following decimal point) is the minimum that
 788         //  will capture the full information stored in a double
 789         is.precision(16);
 790         // we use decimal to store integers, because it is smaller than hex!
 791         is.setf(std::ios::dec,std::ios::basefield);
 792         is.setf(std::ios::scientific,std::ios::floatfield);
 793     }
 794     return is;
 795 }
 796
 797 } // detail
 798
 799 } // HepMC