HepMC Reference Documentation

HepMC

GenEventStreamIO.cc

Go to the documentation of this file.
00001 //--------------------------------------------------------------------------
00002 //
00003 // GenEventStreamIO.cc
00004 // Author:  Lynn Garren
00005 //
00006 // Implement operator >> and operator <<
00007 //
00008 // ----------------------------------------------------------------------
00009 
00010 #include <iostream>
00011 #include <ostream>
00012 #include <istream>
00013 #include <sstream>
00014 
00015 #include "HepMC/GenEvent.h"
00016 #include "HepMC/GenCrossSection.h"
00017 #include "HepMC/StreamInfo.h"
00018 #include "HepMC/StreamHelpers.h"
00019 #include "HepMC/Version.h"
00020 #include "HepMC/IO_Exception.h"
00021 
00022 namespace HepMC {
00023 
00024 // ------------------------- local methods ----------------
00025 
00029 void HepMCStreamCallback(std::ios_base::event e, std::ios_base& b, int i)
00030 {
00031   // only clean up if the stream object is going away.
00032   if(i!=0 && e!= std::ios_base::erase_event) return;
00033 
00034   // retrieve the pointer to the object
00035   StreamInfo* hd = (StreamInfo*)b.pword(i);
00036   b.pword(i) = 0;
00037   b.iword(i) = 0;
00038 #ifdef HEPMC_DEBUG
00039   // the following line is just for sanity checking
00040   if(hd) std::cerr << "deleted StreamInfo " << hd->stream_id() << "\n";
00041 #endif
00042   delete hd;
00043 }
00044 
00045 // ------------------------- iomanip ----------------
00046 
00050 template <class IO>
00051 StreamInfo& get_stream_info(IO& iost)
00052 {
00053   if(iost.iword(0) == 0)
00054     {
00055       // make sure we add the callback if this is the first time through
00056       iost.iword(0)=1;
00057       iost.register_callback(&HepMCStreamCallback, 0);
00058       // this is our special "context" record.
00059       // there is one of these at the head of each IO block.
00060       // allocate room for a StreamInfo in the userdata area
00061       iost.pword(0) = new StreamInfo;
00062 #ifdef HEPMC_DEBUG
00063       // the following line is just for sanity checking
00064       std::cerr << "created StreamInfo " << ((StreamInfo*)iost.pword(0))->stream_id() << "\n";
00065 #endif
00066     }
00067   return *(StreamInfo*)iost.pword(0);
00068 }
00069         
00070 // ------------------------- GenEvent member functions ----------------
00071 
00072 std::ostream& GenEvent::write( std::ostream& os )
00073 {
00075 
00076     //
00077     StreamInfo & info = get_stream_info(os);
00078     //
00079     // if this is the first event, set precision
00080     if ( !info.finished_first_event() ) {
00081         // precision 16 (# digits following decimal point) is the minimum that
00082         //  will capture the full information stored in a double
00083         //  However, we let the user set precision, since that is the expected functionality
00084         // we use decimal to store integers, because it is smaller than hex!
00085         os.setf(std::ios::dec,std::ios::basefield);
00086         os.setf(std::ios::scientific,std::ios::floatfield);
00087         //
00088         info.set_finished_first_event(true);
00089     }
00090     //
00091     // output the event data including the number of primary vertices
00092     //  and the total number of vertices
00093     //std::vector<long> random_states = random_states();
00094     os << 'E';
00095     detail::output( os, event_number() );
00096     detail::output( os, mpi() );
00097     detail::output( os, event_scale() );
00098     detail::output( os, alphaQCD() );
00099     detail::output( os, alphaQED() );
00100     detail::output( os, signal_process_id() );
00101     detail::output( os,   ( signal_process_vertex() ?
00102                 signal_process_vertex()->barcode() : 0 )   );
00103     detail::output( os, vertices_size() ); // total number of vertices.
00104     write_beam_particles( os, beam_particles() );
00105     // random state
00106     detail::output( os, (int)m_random_states.size() );
00107     for ( std::vector<long>::iterator rs = m_random_states.begin(); 
00108           rs != m_random_states.end(); ++rs ) {
00109          detail::output( os, *rs );
00110     }
00111     // weights
00112     // we need to iterate over the map so that the weights printed 
00113     // here will be in the same order as the names printed next
00114     os << ' ' << (int)weights().size() ;
00115     for ( WeightContainer::const_map_iterator w = weights().map_begin(); 
00116           w != weights().map_end(); ++w ) {
00117         detail::output( os, m_weights[w->second] );
00118     }
00119     detail::output( os,'\n');
00120     // now add names for weights
00121     // note that this prints a new line if and only if the weight container
00122     // is not empty
00123     if ( ! weights().empty() ) {
00124         os << "N " << weights().size() << " " ;
00125         for ( WeightContainer::const_map_iterator w = weights().map_begin(); 
00126               w != weights().map_end(); ++w ) {
00127             detail::output( os,'"');
00128             os << w->first;
00129             detail::output( os,'"');
00130             detail::output( os,' ');
00131         }
00132         detail::output( os,'\n');
00133     }
00134     //
00135     // Units
00136     os << "U " << name(momentum_unit());
00137     os << " " << name(length_unit());
00138     detail::output( os,'\n');
00139     //
00140     // write GenCrossSection if it has been set
00141     if( m_cross_section ) m_cross_section->write(os);
00142     //
00143     // write HeavyIon and PdfInfo if they have been set
00144     if( m_heavy_ion ) os << heavy_ion() ;
00145     if( m_pdf_info ) os << pdf_info() ;
00146     //
00147     // Output all of the vertices - note there is no real order.
00148     for ( GenEvent::vertex_const_iterator v = vertices_begin();
00149           v != vertices_end(); ++v ) {
00150         write_vertex(os, *v);
00151     }
00152     return os;
00153 }
00154 
00155 std::istream& GenEvent::read( std::istream& is )
00156 {
00158     //
00159     StreamInfo & info = get_stream_info(is);
00160     clear();
00161     //
00162     // search for event listing key before first event only.
00163     if ( !info.finished_first_event() ) {
00164         //
00165         find_file_type(is);
00166         info.set_finished_first_event(true);
00167     }
00168     //
00169     // make sure the stream is good
00170     if ( !is ) {
00171         std::cerr << "streaming input: end of stream found "
00172                   << "setting badbit." << std::endl;
00173         is.clear(std::ios::badbit); 
00174         return is;
00175     }
00176 
00177     //
00178     // test to be sure the next entry is of type "E" then ignore it
00179     if ( is.peek()!='E' ) { 
00180         // if the E is not the next entry, then check to see if it is
00181         // the end event listing key - if yes, search for another start key
00182         int ioendtype;
00183         find_end_key(is,ioendtype);
00184         if ( ioendtype == info.io_type() ) {
00185             find_file_type(is);
00186             // are we at the end of the file?
00187             if( !is ) return is;
00188         } else if ( ioendtype > 0 ) {
00189             std::cerr << "streaming input: end key does not match start key "
00190                       << "setting badbit." << std::endl;
00191             is.clear(std::ios::badbit); 
00192             return is;
00193         } else if ( !info.has_key() ) {
00194             find_file_type(is);
00195             // are we at the end of the file?
00196             if( !is ) return is;
00197         } else {
00198             std::cerr << "streaming input: end key not found "
00199                       << "setting badbit." << std::endl;
00200             is.clear(std::ios::badbit); 
00201             return is;
00202         }
00203     } 
00204 
00205     int signal_process_vertex = 0;
00206     int num_vertices = 0, bp1 = 0, bp2 = 0;
00207     bool units_line = false;
00208     // OK - now ready to start reading the event, so set the header flag
00209     info.set_reading_event_header(true);
00210     // The flag will be set to false when we reach the end of the header
00211     while(info.reading_event_header()) {
00212         switch(is.peek()) {
00213             case 'E':
00214             {   // deal with the event line
00215                 process_event_line( is, num_vertices, bp1, bp2, signal_process_vertex );
00216             } break;
00217             case 'N':
00218             {   // get weight names 
00219                 read_weight_names( is );
00220             } break;
00221             case 'U':
00222             {   // get unit information if it exists
00223                 units_line = true;
00224                 if( info.io_type() == gen ) {
00225                     read_units( is );
00226                 }
00227             } break;
00228             case 'C':
00229             {   // we have a GenCrossSection line
00230                 // create cross section
00231                 GenCrossSection xs;
00232                 // check for invalid data
00233                 try {
00234                     // read the line
00235                     xs.read(is);
00236                 }
00237                 catch (IO_Exception& e) {
00238                     detail::find_event_end( is );
00239                 }
00240                 if(xs.is_set()) { 
00241                     set_cross_section( xs );
00242                 }
00243             } break;
00244             case 'H':
00245             {   // we have a HeavyIon line OR an unexpected HepMC... line
00246                 if( info.io_type() == gen || info.io_type() == extascii ) {
00247                     // get HeavyIon
00248                     HeavyIon ion;
00249                     // check for invalid data
00250                     try {
00251                         is >> &ion;
00252                     }
00253                     catch (IO_Exception& e) {
00254                         detail::find_event_end( is );
00255                     }
00256                     if(ion.is_valid()) { 
00257                         set_heavy_ion( ion );
00258                     }
00259                 }
00260             } break;
00261             case 'F':
00262             {   // we have a PdfInfo line
00263                 if( info.io_type() == gen || info.io_type() == extascii ) {
00264                     // get PdfInfo
00265                     PdfInfo pdf;
00266                     // check for invalid data
00267                     try {
00268                         is >> &pdf;
00269                     }
00270                     catch (IO_Exception& e) {
00271                         detail::find_event_end( is );
00272                     }
00273                     if(pdf.is_valid()) { 
00274                         set_pdf_info( pdf );
00275                     }
00276                 }
00277             } break;
00278             case 'V':
00279             {   
00280                 // this should be the first vertex line - exit this loop
00281                 info.set_reading_event_header(false);
00282             } break;
00283             case 'P':
00284             {   // we should not find this line
00285                 std::cerr << "streaming input: found unexpected line P" << std::endl;
00286                 info.set_reading_event_header(false);
00287             } break;
00288             default:
00289                 // ignore everything else
00290                 break;
00291         } // switch on line type
00292     } // while reading_event_header
00293     // before proceeding - did we find a units line?
00294     if( !units_line ) {
00295         use_units( info.io_momentum_unit(), 
00296                        info.io_position_unit() );
00297     }
00298     //
00299     // the end vertices of the particles are not connected until
00300     //  after the event is read --- we store the values in a map until then
00301     TempParticleMap particle_to_end_vertex;
00302     //
00303     // read in the vertices
00304     for ( int iii = 1; iii <= num_vertices; ++iii ) {
00305         GenVertex* v = new GenVertex();
00306         try {
00307             detail::read_vertex(is,particle_to_end_vertex,v);
00308         }
00309         catch (IO_Exception& e) {
00310             for( TempParticleMap::orderIterator it = particle_to_end_vertex.order_begin(); 
00311                  it != particle_to_end_vertex.order_end(); ++it ) {
00312                 GenParticle* p = it->second;
00313                 // delete particles only if they are not already owned by a vertex
00314                 if( p->production_vertex() ) {
00315                 } else if( p->end_vertex() ) {
00316                 } else {
00317                      delete p;
00318                 }
00319             }
00320             delete v;
00321             detail::find_event_end( is );
00322         }
00323         add_vertex( v );
00324     }
00325     // set the signal process vertex
00326     if ( signal_process_vertex ) {
00327         set_signal_process_vertex( 
00328             barcode_to_vertex(signal_process_vertex) );
00329     }
00330     //
00331     // last connect particles to their end vertices
00332     GenParticle* beam1(0);
00333     GenParticle* beam2(0);
00334     for ( TempParticleMap::orderIterator pmap 
00335               = particle_to_end_vertex.order_begin(); 
00336           pmap != particle_to_end_vertex.order_end(); ++pmap ) {
00337         GenParticle* p =  pmap->second;
00338         int vtx = particle_to_end_vertex.end_vertex( p );
00339         GenVertex* itsDecayVtx = barcode_to_vertex(vtx);
00340         if ( itsDecayVtx ) itsDecayVtx->add_particle_in( p );
00341         else {
00342             std::cerr << "read_io_genevent: ERROR particle points"
00343                       << " to null end vertex. " <<std::endl;
00344         }
00345         // also look for the beam particles
00346         if( p->barcode() == bp1 ) beam1 = p;
00347         if( p->barcode() == bp2 ) beam2 = p;
00348     }
00349     set_beam_particles(beam1,beam2);
00350     return is;
00351 }
00352 
00353 // ------------------------- operator << and operator >> ----------------
00354 
00355 std::ostream & operator << (std::ostream & os, GenEvent & evt)
00356 {
00358     evt.write(os);
00359     return os;
00360 }
00361 
00362 std::istream & operator >> (std::istream & is, GenEvent & evt)
00363 {
00364     evt.read(is);
00365     return is;
00366 }
00367 
00368 // ------------------------- set units ----------------
00369 
00370 std::istream & set_input_units(std::istream & is, 
00371                                Units::MomentumUnit mom,
00372                                Units::LengthUnit len )
00373 {
00374     //
00375     StreamInfo & info = get_stream_info(is);
00376     info.use_input_units( mom, len );
00377     return is;
00378 }
00379 
00380 // ------------------------- begin and end block lines ----------------
00381 
00382 std::ostream & write_HepMC_IO_block_begin(std::ostream & os )
00383 {
00384     //
00385     StreamInfo & info = get_stream_info(os);
00386 
00387     if( !info.finished_first_event() ) {
00388     os << "\n" << "HepMC::Version " << versionName();
00389     os << "\n";
00390     os << info.IO_GenEvent_Key() << "\n";
00391     }
00392     return os;
00393 }
00394 
00395 std::ostream & write_HepMC_IO_block_end(std::ostream & os )
00396 {
00397     //
00398     StreamInfo & info = get_stream_info(os);
00399 
00400     if( info.finished_first_event() ) {
00401         os << info.IO_GenEvent_End() << "\n";
00402         os << std::flush;
00403     }
00404     return os;
00405 }
00406 
00407 std::istream & GenEvent::process_event_line( std::istream & is, 
00408                                              int & num_vertices,
00409                                              int & bp1, int & bp2,
00410                                              int & signal_process_vertex )
00411 {
00412     //
00413     if ( !is ) {
00414         std::cerr << "GenEvent::process_event_line setting badbit." << std::endl;
00415         is.clear(std::ios::badbit);
00416         return is;
00417     } 
00418     //
00419     StreamInfo & info = get_stream_info(is);
00420     std::string line;
00421     std::getline(is,line);
00422     std::istringstream iline(line);
00423     std::string firstc;
00424     iline >> firstc;
00425     //
00426     // read values into temp variables, then fill GenEvent
00427     int event_number = 0, signal_process_id = 0,
00428         random_states_size = 0, nmpi = -1;
00429     double eventScale = 0, alpha_qcd = 0, alpha_qed = 0;
00430     iline >> event_number;
00431     if(!iline) detail::find_event_end( is );
00432     if( info.io_type() == gen || info.io_type() == extascii ) {
00433         iline >> nmpi;
00434         if(!iline) detail::find_event_end( is );
00435         set_mpi( nmpi );
00436     }
00437     iline >> eventScale ;
00438     if(!iline) detail::find_event_end( is );
00439     iline >> alpha_qcd ;
00440     if(!iline) detail::find_event_end( is );
00441     iline >> alpha_qed;
00442     if(!iline) detail::find_event_end( is );
00443     iline >> signal_process_id ;
00444     if(!iline) detail::find_event_end( is );
00445     iline >> signal_process_vertex;
00446     if(!iline) detail::find_event_end( is );
00447     iline >> num_vertices;
00448     if(!iline) detail::find_event_end( is );
00449     if( info.io_type() == gen || info.io_type() == extascii ) {
00450         iline >> bp1 ;
00451         if(!iline) detail::find_event_end( is );
00452         iline >> bp2;
00453         if(!iline) detail::find_event_end( is );
00454     }
00455     iline >> random_states_size;
00456     if(!iline) detail::find_event_end( is );
00457     std::vector<long> random_states(random_states_size);
00458     for ( int i = 0; i < random_states_size; ++i ) {
00459         iline >> random_states[i];
00460         if(!iline) detail::find_event_end( is );
00461     }
00462     WeightContainer::size_type weights_size = 0;
00463     iline >> weights_size;
00464     if(!iline) detail::find_event_end( is );
00465     std::vector<double> wgt(weights_size);
00466     for ( WeightContainer::size_type ii = 0; ii < weights_size; ++ii ) {
00467         iline >> wgt[ii];
00468         if(!iline) detail::find_event_end( is );
00469     }
00470     // weight names will be added later if they exist
00471     if( weights_size > 0 ) m_weights = wgt;
00472     // 
00473     // fill signal_process_id, event_number, random_states, etc.
00474     set_signal_process_id( signal_process_id );
00475     set_event_number( event_number );
00476     set_random_states( random_states );
00477     set_event_scale( eventScale );
00478     set_alphaQCD( alpha_qcd );
00479     set_alphaQED( alpha_qed );
00480     //
00481     return is;
00482 }
00483 
00484 std::istream & GenEvent::read_weight_names( std::istream & is )
00485 {
00486     // now check for a named weight line
00487     if ( !is ) {
00488         std::cerr << "GenEvent::read_weight_names setting badbit." << std::endl;
00489         is.clear(std::ios::badbit);
00490         return is;
00491     } 
00492     // Test to be sure the next entry is of type "N"
00493     // If we have no named weight line, this is not an error
00494     // releases prior to 2.06.00 do not have named weights
00495     if ( is.peek() !='N') {
00496         return is;
00497     } 
00498     // now get this line and process it
00499     std::string line;
00500     std::getline(is,line);
00501     std::istringstream wline(line);
00502     std::string firstc;
00503     WeightContainer::size_type name_size = 0;
00504     wline >> firstc >> name_size;
00505     if(!wline) detail::find_event_end( is );
00506     if( firstc != "N") { 
00507         std::cout << "debug: first character of named weights is " << firstc << std::endl;
00508         std::cout << "debug: We should never get here" << std::endl;
00509         is.clear(std::ios::badbit);
00510         return is;
00511     }
00512     if( m_weights.size() != name_size ) { 
00513         std::cout << "debug: weight sizes do not match "<< std::endl;
00514         std::cout << "debug: weight vector size is " << m_weights.size() << std::endl;
00515         std::cout << "debug: weight name size is " << name_size << std::endl;
00516         is.clear(std::ios::badbit);
00517         return is;
00518     }
00519     std::string name;
00520     std::string::size_type i1 = line.find("\"");
00521     std::string::size_type i2;
00522     std::string::size_type len = line.size();
00523     WeightContainer namedWeight;
00524     for ( WeightContainer::size_type ii = 0; ii < name_size; ++ii ) {
00525         // weight names may contain blanks
00526         if(i1 >= len) {
00527             std::cout << "debug: attempting to read past the end of the named weight line " << std::endl;
00528             std::cout << "debug: We should never get here" << std::endl;
00529             std::cout << "debug: Looking for the end of this event" << std::endl;
00530             detail::find_event_end( is );
00531         }
00532         i2 = line.find("\"",i1+1);
00533         name = line.substr(i1+1,i2-i1-1);
00534         namedWeight[name] = m_weights[ii];
00535         i1 = line.find("\"",i2+1);
00536     }
00537     m_weights = namedWeight;
00538     return is;
00539 }
00540 
00541 std::istream & GenEvent::read_units( std::istream & is )
00542 {
00543     //
00544     if ( !is ) {
00545         std::cerr << "GenEvent::read_units setting badbit." << std::endl;
00546         is.clear(std::ios::badbit);
00547         return is;
00548     } 
00549     //
00550     StreamInfo & info = get_stream_info(is);
00551     // test to be sure the next entry is of type "U" then ignore it
00552     // if we have no units, this is not an error
00553     // releases prior to 2.04.00 did not write unit information
00554     if ( is.peek() !='U') {
00555         use_units( info.io_momentum_unit(), 
00556                        info.io_position_unit() );
00557         return is;
00558     } 
00559     is.ignore();        // ignore the first character in the line
00560     std::string mom, pos;
00561     is >> mom >> pos;
00562     is.ignore(1);      // eat the extra whitespace
00563     use_units(mom,pos);
00564     //
00565     return is;
00566 }
00567 
00568 std::istream & GenEvent::find_file_type( std::istream & istr )
00569 {
00570     //
00571     // make sure the stream is good
00572     if ( !istr ) return istr;
00573 
00574     //
00575     StreamInfo & info = get_stream_info(istr);
00576 
00577     // if there is no input block line, then we assume this stream
00578     // is in the IO_GenEvent format
00579     if ( istr.peek()=='E' ) {
00580         info.set_io_type( gen );
00581         info.set_has_key(false);
00582         return istr;
00583     }
00584     
00585     std::string line;
00586     while ( std::getline(istr,line) ) {
00587         //
00588         // search for event listing key before first event only.
00589         //
00590         if( line == info.IO_GenEvent_Key() ) {
00591             info.set_io_type( gen );
00592             info.set_has_key(true);
00593             return istr;
00594         } else if( line == info.IO_Ascii_Key() ) {
00595             info.set_io_type( ascii );
00596             info.set_has_key(true);
00597             return istr;
00598         } else if( line == info.IO_ExtendedAscii_Key() ) {
00599             info.set_io_type( extascii );
00600             info.set_has_key(true);
00601             return istr;
00602         } else if( line == info.IO_Ascii_PDT_Key() ) {
00603             info.set_io_type( ascii_pdt );
00604             info.set_has_key(true);
00605             return istr;
00606         } else if( line == info.IO_ExtendedAscii_PDT_Key() ) {
00607             info.set_io_type( extascii_pdt );
00608             info.set_has_key(true);
00609             return istr;
00610         }
00611     }
00612     info.set_io_type( 0 );
00613     info.set_has_key(false);
00614     return istr;
00615 }
00616 
00617 std::istream & GenEvent::find_end_key( std::istream & istr, int & iotype )
00618 {
00619     iotype = 0;
00620     // peek at the first character before proceeding
00621     if( istr.peek()!='H' ) return istr;
00622     //
00623     // we only check the next line
00624     std::string line;
00625     std::getline(istr,line);
00626     //
00627     StreamInfo & info = get_stream_info(istr);
00628     //
00629     // check to see if this is an end key
00630     if( line == info.IO_GenEvent_End() ) {
00631         iotype = gen;
00632     } else if( line == info.IO_Ascii_End() ) {
00633         iotype = ascii;
00634     } else if( line == info.IO_ExtendedAscii_End() ) {
00635         iotype = extascii;
00636     } else if( line == info.IO_Ascii_PDT_End() ) {
00637         iotype = ascii_pdt;
00638     } else if( line == info.IO_ExtendedAscii_PDT_End() ) {
00639         iotype = extascii_pdt;
00640     }
00641     if( iotype != 0 && info.io_type() != iotype ) {
00642         std::cerr << "GenEvent::find_end_key: iotype keys have changed" << std::endl;
00643     } else {
00644         return istr;
00645     }
00646     //
00647     // if we get here, then something has gotten badly confused
00648     std::cerr << "GenEvent::find_end_key: MALFORMED INPUT" << std::endl;
00649     istr.clear(std::ios::badbit); 
00650     return istr;
00651 }
00652 
00653 std::ostream & establish_output_stream_info( std::ostream & os )
00654 {
00655     StreamInfo & info = get_stream_info(os);
00656     if ( !info.finished_first_event() ) {
00657         // precision 16 (# digits following decimal point) is the minimum that
00658         //  will capture the full information stored in a double
00659         os.precision(16);
00660         // we use decimal to store integers, because it is smaller than hex!
00661         os.setf(std::ios::dec,std::ios::basefield);
00662         os.setf(std::ios::scientific,std::ios::floatfield);
00663     }
00664     return os;
00665 }
00666 
00667 std::istream & establish_input_stream_info( std::istream & is )
00668 {
00669     StreamInfo & info = get_stream_info(is);
00670     if ( !info.finished_first_event() ) {
00671         // precision 16 (# digits following decimal point) is the minimum that
00672         //  will capture the full information stored in a double
00673         is.precision(16);
00674         // we use decimal to store integers, because it is smaller than hex!
00675         is.setf(std::ios::dec,std::ios::basefield);
00676         is.setf(std::ios::scientific,std::ios::floatfield);
00677     }
00678     return is;
00679 }
00680 
00681 
00682 // ------------------------- helper functions ----------------
00683 
00684 namespace detail {
00685 
00686 // The functions defined here need to use get_stream_info
00687 
00688 std::istream & read_particle( std::istream & is, 
00689                               TempParticleMap & particle_to_end_vertex, 
00690                               GenParticle * p )
00691 {
00692     // get the next line
00693     std::string line;
00694     std::getline(is,line);
00695     std::istringstream iline(line);
00696     std::string firstc;
00697     iline >> firstc;
00698     if( firstc != "P" ) { 
00699         std::cerr << "StreamHelpers::detail::read_particle invalid line type: " 
00700                   << firstc << std::endl;
00701         std::cerr << "StreamHelpers::detail::read_particle setting badbit." 
00702                   << std::endl;
00703         is.clear(std::ios::badbit); 
00704         return is;
00705     } 
00706     //
00707     StreamInfo & info = get_stream_info(is);
00708     //testHepMC.cc
00709     // declare variables to be read in to, and read everything except flow
00710     double px = 0., py = 0., pz = 0., e = 0., m = 0., theta = 0., phi = 0.;
00711     int bar_code = 0, id = 0, status = 0, end_vtx_code = 0, flow_size = 0;
00712     // check that the input stream is still OK after reading item
00713     iline >> bar_code ;
00714     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00715     iline >> id ;
00716     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00717     iline >> px ;
00718     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00719     iline >> py ;
00720     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00721     iline >> pz ;
00722     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00723     iline >> e ;
00724     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00725     if( info.io_type() != ascii ) {
00726         iline >> m ;
00727         if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00728     }
00729     iline >> status ;
00730     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00731     iline >> theta ;
00732     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00733     iline >> phi ;
00734     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00735     iline >> end_vtx_code ;
00736     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00737     iline >> flow_size;
00738     if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00739     //
00740     // read flow patterns if any exist
00741     Flow flow;
00742     int code_index, code;
00743     for ( int i = 1; i <= flow_size; ++i ) {
00744         iline >> code_index >> code;
00745         if(!iline) {  delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
00746         flow.set_icode( code_index,code);
00747     }
00748     p->set_momentum( FourVector(px,py,pz,e) );
00749     p->set_pdg_id( id );
00750     p->set_status( status );
00751     p->set_flow( flow );
00752     p->set_polarization( Polarization(theta,phi) );
00753     if( info.io_type() == ascii ) {
00754         p->set_generated_mass( p->momentum().m() );
00755     } else {
00756         p->set_generated_mass( m );
00757     }
00758     p->suggest_barcode( bar_code );
00759     //
00760     // all particles are connected to their end vertex separately 
00761     // after all particles and vertices have been created - so we keep
00762     // a map of all particles that have end vertices
00763     if ( end_vtx_code != 0 ) {
00764         particle_to_end_vertex.addEndParticle(p,end_vtx_code);
00765     }
00766     return is;
00767 }
00768 
00769 std::ostream & establish_output_stream_info( std::ostream & os )
00770 {
00771     StreamInfo & info = get_stream_info(os);
00772     if ( !info.finished_first_event() ) {
00773         // precision 16 (# digits following decimal point) is the minimum that
00774         //  will capture the full information stored in a double
00775         os.precision(16);
00776         // we use decimal to store integers, because it is smaller than hex!
00777         os.setf(std::ios::dec,std::ios::basefield);
00778         os.setf(std::ios::scientific,std::ios::floatfield);
00779     }
00780     return os;
00781 }
00782 
00783 std::istream & establish_input_stream_info( std::istream & is )
00784 {
00785     StreamInfo & info = get_stream_info(is);
00786     if ( !info.finished_first_event() ) {
00787         // precision 16 (# digits following decimal point) is the minimum that
00788         //  will capture the full information stored in a double
00789         is.precision(16);
00790         // we use decimal to store integers, because it is smaller than hex!
00791         is.setf(std::ios::dec,std::ios::basefield);
00792         is.setf(std::ios::scientific,std::ios::floatfield);
00793     }
00794     return is;
00795 }
00796 
00797 } // detail
00798 
00799 } // HepMC

Generated on Fri Feb 17 00:31:25 2012 for HepMC by  doxygen 1.4.7