00001
00002
00003
00004
00005
00006
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
00025
00029 void HepMCStreamCallback(std::ios_base::event e, std::ios_base& b, int i)
00030 {
00031
00032 if(i!=0 && e!= std::ios_base::erase_event) return;
00033
00034
00035 StreamInfo* hd = (StreamInfo*)b.pword(i);
00036 b.pword(i) = 0;
00037 b.iword(i) = 0;
00038 #ifdef HEPMC_DEBUG
00039
00040 if(hd) std::cerr << "deleted StreamInfo " << hd->stream_id() << "\n";
00041 #endif
00042 delete hd;
00043 }
00044
00045
00046
00050 template <class IO>
00051 StreamInfo& get_stream_info(IO& iost)
00052 {
00053 if(iost.iword(0) == 0)
00054 {
00055
00056 iost.iword(0)=1;
00057 iost.register_callback(&HepMCStreamCallback, 0);
00058
00059
00060
00061 iost.pword(0) = new StreamInfo;
00062 #ifdef HEPMC_DEBUG
00063
00064 std::cerr << "created StreamInfo " << ((StreamInfo*)iost.pword(0))->stream_id() << "\n";
00065 #endif
00066 }
00067 return *(StreamInfo*)iost.pword(0);
00068 }
00069
00070
00071
00072 std::ostream& GenEvent::write( std::ostream& os )
00073 {
00075
00076
00077 StreamInfo & info = get_stream_info(os);
00078
00079
00080 if ( !info.finished_first_event() ) {
00081
00082
00083
00084
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
00092
00093
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() );
00104 write_beam_particles( os, beam_particles() );
00105
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
00112
00113
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
00121
00122
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
00136 os << "U " << name(momentum_unit());
00137 os << " " << name(length_unit());
00138 detail::output( os,'\n');
00139
00140
00141 if( m_cross_section ) m_cross_section->write(os);
00142
00143
00144 if( m_heavy_ion ) os << heavy_ion() ;
00145 if( m_pdf_info ) os << pdf_info() ;
00146
00147
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
00163 if ( !info.finished_first_event() ) {
00164
00165 find_file_type(is);
00166 info.set_finished_first_event(true);
00167 }
00168
00169
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
00179 if ( is.peek()!='E' ) {
00180
00181
00182 int ioendtype;
00183 find_end_key(is,ioendtype);
00184 if ( ioendtype == info.io_type() ) {
00185 find_file_type(is);
00186
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
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
00209 info.set_reading_event_header(true);
00210
00211 while(info.reading_event_header()) {
00212 switch(is.peek()) {
00213 case 'E':
00214 {
00215 process_event_line( is, num_vertices, bp1, bp2, signal_process_vertex );
00216 } break;
00217 case 'N':
00218 {
00219 read_weight_names( is );
00220 } break;
00221 case 'U':
00222 {
00223 units_line = true;
00224 if( info.io_type() == gen ) {
00225 read_units( is );
00226 }
00227 } break;
00228 case 'C':
00229 {
00230
00231 GenCrossSection xs;
00232
00233 try {
00234
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 {
00246 if( info.io_type() == gen || info.io_type() == extascii ) {
00247
00248 HeavyIon ion;
00249
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 {
00263 if( info.io_type() == gen || info.io_type() == extascii ) {
00264
00265 PdfInfo pdf;
00266
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
00281 info.set_reading_event_header(false);
00282 } break;
00283 case 'P':
00284 {
00285 std::cerr << "streaming input: found unexpected line P" << std::endl;
00286 info.set_reading_event_header(false);
00287 } break;
00288 default:
00289
00290 break;
00291 }
00292 }
00293
00294 if( !units_line ) {
00295 use_units( info.io_momentum_unit(),
00296 info.io_position_unit() );
00297 }
00298
00299
00300
00301 TempParticleMap particle_to_end_vertex;
00302
00303
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
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
00326 if ( signal_process_vertex ) {
00327 set_signal_process_vertex(
00328 barcode_to_vertex(signal_process_vertex) );
00329 }
00330
00331
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
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
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
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
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
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
00471 if( weights_size > 0 ) m_weights = wgt;
00472
00473
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
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
00493
00494
00495 if ( is.peek() !='N') {
00496 return is;
00497 }
00498
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
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
00552
00553
00554 if ( is.peek() !='U') {
00555 use_units( info.io_momentum_unit(),
00556 info.io_position_unit() );
00557 return is;
00558 }
00559 is.ignore();
00560 std::string mom, pos;
00561 is >> mom >> pos;
00562 is.ignore(1);
00563 use_units(mom,pos);
00564
00565 return is;
00566 }
00567
00568 std::istream & GenEvent::find_file_type( std::istream & istr )
00569 {
00570
00571
00572 if ( !istr ) return istr;
00573
00574
00575 StreamInfo & info = get_stream_info(istr);
00576
00577
00578
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
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
00621 if( istr.peek()!='H' ) return istr;
00622
00623
00624 std::string line;
00625 std::getline(istr,line);
00626
00627 StreamInfo & info = get_stream_info(istr);
00628
00629
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
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
00658
00659 os.precision(16);
00660
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
00672
00673 is.precision(16);
00674
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
00683
00684 namespace detail {
00685
00686
00687
00688 std::istream & read_particle( std::istream & is,
00689 TempParticleMap & particle_to_end_vertex,
00690 GenParticle * p )
00691 {
00692
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
00709
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
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
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
00761
00762
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
00774
00775 os.precision(16);
00776
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
00788
00789 is.precision(16);
00790
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 }
00798
00799 }