Changeset 571 in svn for trunk/Utilities/HepMC/src/IO_GenEvent.cc
- Timestamp:
- Nov 2, 2011, 5:06:22 PM (13 years ago)
- File:
-
- 1 edited
-
trunk/Utilities/HepMC/src/IO_GenEvent.cc (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/Utilities/HepMC/src/IO_GenEvent.cc
r349 r571 8 8 9 9 #include "IO_GenEvent.h" 10 #include "IO_Exception.h" 10 11 #include "GenEvent.h" 11 #include "ParticleDataTable.h" 12 #include "HeavyIon.h" 13 #include "PdfInfo.h" 14 #include "CommonIO.h" 15 #include "Version.h" 12 #include "StreamHelpers.h" 16 13 17 14 namespace HepMC { … … 23 20 m_istr(0), 24 21 m_iostr(0), 25 m_finished_first_event_io(false),26 22 m_have_file(false), 27 m_common_io() 23 m_error_type(IO_Exception::OK), 24 m_error_message() 28 25 { 29 26 if ( (m_mode&std::ios::out && m_mode&std::ios::in) || 30 27 (m_mode&std::ios::app && m_mode&std::ios::in) ) { 31 std::cerr << "IO_GenEvent::IO_GenEvent Error, open of file requested " 32 << "of input AND output type. Not allowed. Closing file."33 28 m_error_type = IO_Exception::InputAndOutput; 29 m_error_message ="IO_GenEvent::IO_GenEvent Error, open of file requested of input AND output type. Not allowed. Closing file."; 30 std::cerr << m_error_message << std::endl; 34 31 m_file.close(); 35 32 return; 36 33 } 37 // precision 16 (# digits following decimal point) is the minimum that38 // will capture the full information stored in a double39 m_file.precision(16);40 // we use decimal to store integers, because it is smaller than hex!41 m_file.setf(std::ios::dec,std::ios::basefield);42 m_file.setf(std::ios::scientific,std::ios::floatfield);43 34 // now we set the streams 44 35 m_iostr = &m_file; … … 46 37 m_istr = &m_file; 47 38 m_ostr = NULL; 39 detail::establish_input_stream_info(m_file); 48 40 } 49 41 if ( m_mode&std::ios::out ) { 50 42 m_ostr = &m_file; 51 43 m_istr = NULL; 44 detail::establish_output_stream_info(m_file); 52 45 } 53 46 m_have_file = true; … … 59 52 m_istr(&istr), 60 53 m_iostr(&istr), 61 m_finished_first_event_io(false),62 54 m_have_file(false), 63 m_common_io() 64 { } 55 m_error_type(IO_Exception::OK), 56 m_error_message() 57 { 58 detail::establish_input_stream_info( istr ); 59 } 65 60 66 61 IO_GenEvent::IO_GenEvent( std::ostream & ostr ) … … 68 63 m_istr(0), 69 64 m_iostr(&ostr), 70 m_finished_first_event_io(false),71 65 m_have_file(false), 72 m_common_io() 73 { 74 // precision 16 (# digits following decimal point) is the minimum that 75 // will capture the full information stored in a double 76 m_ostr->precision(16); 77 // we use decimal to store integers, because it is smaller than hex! 78 m_ostr->setf(std::ios::dec,std::ios::basefield); 79 m_ostr->setf(std::ios::scientific,std::ios::floatfield); 80 } 66 m_error_type(IO_Exception::OK), 67 m_error_message() 68 { 69 detail::establish_output_stream_info( ostr ); 70 } 81 71 82 72 IO_GenEvent::~IO_GenEvent() { 83 write_end_listing(); 73 if ( m_ostr != NULL ) { 74 write_HepMC_IO_block_end(*m_ostr); 75 } 84 76 if(m_have_file) m_file.close(); 85 77 } 86 78 87 void IO_GenEvent::use_input_units( Units::MomentumUnit mom, Units::LengthUnit len ) { 88 m_common_io.use_input_units(mom,len); 79 void IO_GenEvent::use_input_units( Units::MomentumUnit mom, 80 Units::LengthUnit len ) { 81 if( m_istr != NULL ) { 82 set_input_units( *m_istr, mom, len ); 83 } 89 84 } 90 85 … … 99 94 } 100 95 96 void IO_GenEvent::precision( int size ) { 97 if( size > 16 ) { 98 std::cerr << "IO_GenEvent::precision Error, " 99 << "precision is greater than 16. " 100 << "Not allowed. Using default precision of 16." 101 << std::endl; 102 size = 16; 103 } 104 if(m_ostr) { 105 m_ostr->precision(size); 106 } 107 } 108 109 bool IO_GenEvent::fill_next_event( GenEvent* evt ){ 110 // 111 // reset error type 112 m_error_type = IO_Exception::OK; 113 // 114 // test that evt pointer is not null 115 if ( !evt ) { 116 m_error_type = IO_Exception::NullEvent; 117 m_error_message = "IO_GenEvent::fill_next_event error - passed null event."; 118 std::cerr << m_error_message << std::endl; 119 return false; 120 } 121 // make sure the stream is good, and that it is in input mode 122 if ( !(*m_istr) ) return false; 123 if ( !m_istr ) { 124 m_error_type = IO_Exception::WrongFileType; 125 m_error_message = "HepMC::IO_GenEvent::fill_next_event attempt to read from output file."; 126 std::cerr << m_error_message << std::endl; 127 return false; 128 } 129 // use streaming input 130 try { 131 *m_istr >> *evt; 132 } 133 catch (IO_Exception& e) { 134 m_error_type = IO_Exception::InvalidData; 135 m_error_message = e.what(); 136 evt->clear(); 137 return false; 138 } 139 if( evt->is_valid() ) return true; 140 return false; 141 } 142 101 143 void IO_GenEvent::write_event( const GenEvent* evt ) { 102 144 /// Writes evt to output stream. It does NOT delete the event after writing. … … 105 147 if ( !evt ) return; 106 148 if ( m_ostr == NULL ) { 107 std::cerr << "HepMC::IO_GenEvent::write_event " 108 << " attempt to write to input file." << std::endl; 149 m_error_type = IO_Exception::WrongFileType; 150 m_error_message = "HepMC::IO_GenEvent::write_event attempt to write to input file."; 151 std::cerr << m_error_message << std::endl; 109 152 return; 110 153 } 111 154 // 112 155 // write event listing key before first event only. 113 if ( !m_finished_first_event_io ) { 114 m_finished_first_event_io = 1; 115 *m_ostr << "\n" << "HepMC::Version " << versionName(); 116 *m_ostr << "\n"; 117 m_common_io.write_IO_GenEvent_Key(*m_ostr); 118 } 119 // 120 // output the event data including the number of primary vertices 121 // and the total number of vertices 122 std::vector<long> random_states = evt->random_states(); 123 *m_ostr << 'E'; 124 output( evt->event_number() ); 125 output( evt->mpi() ); 126 output( evt->event_scale() ); 127 output( evt->alphaQCD() ); 128 output( evt->alphaQED() ); 129 output( evt->signal_process_id() ); 130 output( ( evt->signal_process_vertex() ? 131 evt->signal_process_vertex()->barcode() : 0 ) ); 132 output( evt->vertices_size() ); // total number of vertices. 133 write_beam_particles( evt->beam_particles() ); 134 output( (int)random_states.size() ); 135 for ( std::vector<long>::iterator rs = random_states.begin(); 136 rs != random_states.end(); ++rs ) { 137 output( *rs ); 138 } 139 output( (int)evt->weights().size() ); 140 for ( WeightContainer::const_iterator w = evt->weights().begin(); 141 w != evt->weights().end(); ++w ) { 142 output( *w ); 143 } 144 output('\n'); 145 write_unit_info( evt ); 146 write_heavy_ion( evt->heavy_ion() ); 147 write_pdf_info( evt->pdf_info() ); 148 // 149 // Output all of the vertices - note there is no real order. 150 for ( GenEvent::vertex_const_iterator v = evt->vertices_begin(); 151 v != evt->vertices_end(); ++v ) { 152 write_vertex( *v ); 153 } 154 } 155 156 bool IO_GenEvent::fill_next_event( GenEvent* evt ){ 157 // 158 // 159 // test that evt pointer is not null 160 if ( !evt ) { 161 std::cerr 162 << "IO_GenEvent::fill_next_event error - passed null event." 163 << std::endl; 164 return false; 165 } 166 // make sure the stream is good, and that it is in input mode 167 if ( !(*m_istr) ) return false; 168 if ( !m_istr ) { 169 std::cerr << "HepMC::IO_GenEvent::fill_next_event " 170 << " attempt to read from output file." << std::endl; 171 return false; 172 } 173 // 174 // search for event listing key before first event only. 175 // 176 // skip through the file just after first occurence of the start_key 177 int iotype = 0; 178 if ( !m_finished_first_event_io ) { 179 iotype = m_common_io.find_file_type(*m_istr); 180 if( iotype <= 0 ) { 181 std::cerr << "IO_GenEvent::fill_next_event start key not found " 182 << "setting badbit." << std::endl; 183 m_istr->clear(std::ios::badbit); 184 return false; 185 } 186 m_finished_first_event_io = 1; 187 } 188 // 189 // test to be sure the next entry is of type "E" then ignore it 190 if ( !(*m_istr) ) { 191 std::cerr << "IO_GenEvent::fill_next_event end of stream found " 192 << "setting badbit." << std::endl; 193 m_istr->clear(std::ios::badbit); 194 return false; 195 } 196 if ( !(*m_istr) || m_istr->peek()!='E' ) { 197 // if the E is not the next entry, then check to see if it is 198 // the end event listing key - if yes, search for another start key 199 int ioendtype = m_common_io.find_end_key(*m_istr); 200 if ( ioendtype == iotype ) { 201 iotype = m_common_io.find_file_type(*m_istr); 202 if( iotype <= 0 ) { 203 // this is the only case where we set an EOF state 204 m_istr->clear(std::ios::eofbit); 205 return false; 206 } 207 } else if ( ioendtype > 0 ) { 208 std::cerr << "IO_GenEvent::fill_next_event end key does not match start key " 209 << "setting badbit." << std::endl; 210 m_istr->clear(std::ios::badbit); 211 return false; 212 } else { 213 std::cerr << "IO_GenEvent::fill_next_event end key not found " 214 << "setting badbit." << std::endl; 215 m_istr->clear(std::ios::badbit); 216 return false; 217 } 218 } 219 m_istr->ignore(); 220 // call the appropriate read method 221 if( m_common_io.io_type() == gen ) { 222 return m_common_io.read_io_genevent(m_istr, evt); 223 } else if( m_common_io.io_type() == ascii ) { 224 return m_common_io.read_io_ascii(m_istr, evt); 225 } else if( m_common_io.io_type() == extascii ) { 226 return m_common_io.read_io_extendedascii(m_istr, evt); 227 } else if( m_common_io.io_type() == ascii_pdt ) { 228 } else if( m_common_io.io_type() == extascii_pdt ) { 229 } 230 // should not get to this statement 231 return false; 156 write_HepMC_IO_block_begin(*m_ostr); 157 // explicit cast is necessary 158 GenEvent e = *evt; 159 *m_ostr << e ; 232 160 } 233 161 … … 236 164 if ( !(*m_ostr) ) return; 237 165 if ( m_ostr == NULL ) { 238 std::cerr << "HepMC::IO_GenEvent::write_comment " 239 << " attempt to write to input file." << std::endl; 166 m_error_type = IO_Exception::WrongFileType; 167 m_error_message = "HepMC::IO_GenEvent::write_event attempt to write to input file."; 168 std::cerr << m_error_message << std::endl; 240 169 return; 241 170 } 242 171 // write end of event listing key if events have already been written 243 write_ end_listing();172 write_HepMC_IO_block_end(*m_ostr); 244 173 // insert the comment key before the comment 245 174 *m_ostr << "\n" << "HepMC::IO_GenEvent-COMMENT\n"; 246 175 *m_ostr << comment << std::endl; 247 176 } 248 249 void IO_GenEvent::write_vertex( GenVertex* v ) {250 // assumes mode has already been checked251 if ( !v || !(*m_ostr) ) {252 std::cerr << "IO_GenEvent::write_vertex !v||!(*m_ostr), "253 << "v="<< v << " setting badbit" << std::endl;254 m_ostr->clear(std::ios::badbit);255 return;256 }257 // First collect info we need258 // count the number of orphan particles going into v259 int num_orphans_in = 0;260 for ( GenVertex::particles_in_const_iterator p1261 = v->particles_in_const_begin();262 p1 != v->particles_in_const_end(); ++p1 ) {263 if ( !(*p1)->production_vertex() ) ++num_orphans_in;264 }265 //266 *m_ostr << 'V';267 output( v->barcode() ); // v's unique identifier268 output( v->id() );269 output( v->position().x() );270 output( v->position().y() );271 output( v->position().z() );272 output( v->position().t() );273 output( num_orphans_in );274 output( (int)v->particles_out_size() );275 output( (int)v->weights().size() );276 for ( WeightContainer::iterator w = v->weights().begin();277 w != v->weights().end(); ++w ) {278 output( *w );279 }280 output('\n');281 for ( GenVertex::particles_in_const_iterator p2282 = v->particles_in_const_begin();283 p2 != v->particles_in_const_end(); ++p2 ) {284 if ( !(*p2)->production_vertex() ) {285 write_particle( *p2 );286 }287 }288 for ( GenVertex::particles_out_const_iterator p3289 = v->particles_out_const_begin();290 p3 != v->particles_out_const_end(); ++p3 ) {291 write_particle( *p3 );292 }293 }294 295 void IO_GenEvent::write_beam_particles(296 std::pair<HepMC::GenParticle *,HepMC::GenParticle *> pr ) {297 GenParticle* p = pr.first;298 //m_file << 'B';299 if(!p) {300 output( 0 );301 } else {302 output( p->barcode() );303 }304 p = pr.second;305 if(!p) {306 output( 0 );307 } else {308 output( p->barcode() );309 }310 }311 312 void IO_GenEvent::write_heavy_ion( HeavyIon const * ion ) {313 // assumes mode has already been checked314 if ( !(*m_ostr) ) {315 std::cerr << "IO_GenEvent::write_heavy_ion !(*m_ostr), "316 << " setting badbit" << std::endl;317 m_ostr->clear(std::ios::badbit);318 return;319 }320 *m_ostr << 'H';321 // HeavyIon* is set to 0 by default322 if ( !ion ) {323 output( 0 );324 output( 0 );325 output( 0 );326 output( 0 );327 output( 0 );328 output( 0 );329 output( 0 );330 output( 0 );331 output( 0 );332 output( 0. );333 output( 0. );334 output( 0. );335 output( 0. );336 output('\n');337 return;338 }339 //340 output( ion->Ncoll_hard() );341 output( ion->Npart_proj() );342 output( ion->Npart_targ() );343 output( ion->Ncoll() );344 output( ion->spectator_neutrons() );345 output( ion->spectator_protons() );346 output( ion->N_Nwounded_collisions() );347 output( ion->Nwounded_N_collisions() );348 output( ion->Nwounded_Nwounded_collisions() );349 output( ion->impact_parameter() );350 output( ion->event_plane_angle() );351 output( ion->eccentricity() );352 output( ion->sigma_inel_NN() );353 output('\n');354 }355 356 void IO_GenEvent::write_pdf_info( PdfInfo const * pdf ) {357 // assumes mode has already been checked358 if ( !(*m_ostr) ) {359 std::cerr << "IO_GenEvent::write_pdf_info !(*m_ostr), "360 << " setting badbit" << std::endl;361 m_ostr->clear(std::ios::badbit);362 return;363 }364 *m_ostr << 'F';365 // PdfInfo* is set to 0 by default366 if ( !pdf ) {367 output( 0 );368 output( 0 );369 output( 0. );370 output( 0. );371 output( 0. );372 output( 0. );373 output( 0. );374 output( 0 );375 output( 0 );376 output('\n');377 return;378 }379 //380 output( pdf->id1() );381 output( pdf->id2() );382 output( pdf->x1() );383 output( pdf->x2() );384 output( pdf->scalePDF() );385 output( pdf->pdf1() );386 output( pdf->pdf2() );387 output( pdf->pdf_id1() );388 output( pdf->pdf_id2() );389 output('\n');390 }391 392 void IO_GenEvent::write_unit_info( const GenEvent* evt ) {393 if ( !(*m_ostr) ) {394 std::cerr << "IO_GenEvent::write_unit_info !(*m_ostr), "395 << " setting badbit" << std::endl;396 m_ostr->clear(std::ios::badbit);397 return;398 }399 // could write enums here, but strings are more readable400 *m_ostr << "U " << name(evt->momentum_unit());401 *m_ostr << " " << name(evt->length_unit());402 output('\n');403 }404 405 void IO_GenEvent::write_particle( GenParticle* p ) {406 // assumes mode has already been checked407 if ( !p || !(*m_ostr) ) {408 std::cerr << "IO_GenEvent::write_particle !p||!(*m_ostr), "409 << "p="<< p << " setting badbit" << std::endl;410 m_ostr->clear(std::ios::badbit);411 return;412 }413 *m_ostr << 'P';414 output( p->barcode() );415 output( p->pdg_id() );416 output( p->momentum().px() );417 output( p->momentum().py() );418 output( p->momentum().pz() );419 output( p->momentum().e() );420 output( p->generated_mass() );421 output( p->status() );422 output( p->polarization().theta() );423 output( p->polarization().phi() );424 // since end_vertex is oftentimes null, this CREATES a null vertex425 // in the map426 output( ( p->end_vertex() ? p->end_vertex()->barcode() : 0 ) );427 *m_ostr << ' ' << p->flow() << "\n";428 }429 430 void IO_GenEvent::write_particle_data( const ParticleData* pdata ) {431 // assumes mode has already been checked432 if ( !pdata || !(*m_ostr) ) {433 std::cerr << "IO_GenEvent::write_particle_data !pdata||!(*m_ostr), "434 << "pdata="<< pdata << " setting badbit" << std::endl;435 m_ostr->clear(std::ios::badbit);436 return;437 }438 *m_ostr << 'D';439 output( pdata->pdg_id() );440 output( pdata->charge() );441 output( pdata->mass() );442 output( pdata->clifetime() );443 output( (int)(pdata->spin()*2.+.1) );444 // writes the first 21 characters starting with 0445 *m_ostr << " " << pdata->name().substr(0,21) << "\n";446 }447 448 HeavyIon* IO_GenEvent::read_heavy_ion()449 {450 // assumes mode has already been checked451 //452 // test to be sure the next entry is of type "H" then ignore it453 if ( !(*m_istr) || m_istr->peek()!='H' ) {454 std::cerr << "IO_GenEvent::read_heavy_ion setting badbit." << std::endl;455 m_istr->clear(std::ios::badbit);456 return false;457 }458 m_istr->ignore();459 // read values into temp variables, then create a new HeavyIon object460 int nh =0, np =0, nt =0, nc =0,461 neut = 0, prot = 0, nw =0, nwn =0, nwnw =0;462 float impact = 0., plane = 0., xcen = 0., inel = 0.;463 *m_istr >> nh >> np >> nt >> nc >> neut >> prot464 >> nw >> nwn >> nwnw >> impact >> plane >> xcen >> inel;465 m_istr->ignore(2,'\n');466 if( nh == 0 ) return false;467 HeavyIon* ion = new HeavyIon(nh, np, nt, nc, neut, prot,468 nw, nwn, nwnw,469 impact, plane, xcen, inel );470 //471 return ion;472 }473 474 GenParticle* IO_GenEvent::read_particle(475 TempParticleMap& particle_to_end_vertex ){476 // assumes mode has already been checked477 //478 // test to be sure the next entry is of type "P" then ignore it479 if ( !(*m_istr) || m_istr->peek()!='P' ) {480 std::cerr << "IO_GenEvent::read_particle setting badbit."481 << std::endl;482 m_istr->clear(std::ios::badbit);483 return false;484 }485 m_istr->ignore();486 //487 // declare variables to be read in to, and read everything except flow488 double px = 0., py = 0., pz = 0., e = 0., m = 0., theta = 0., phi = 0.;489 int bar_code = 0, id = 0, status = 0, end_vtx_code = 0, flow_size = 0;490 *m_istr >> bar_code >> id >> px >> py >> pz >> e >> m >> status491 >> theta >> phi >> end_vtx_code >> flow_size;492 //493 // read flow patterns if any exist494 Flow flow;495 int code_index, code;496 for ( int i = 1; i <= flow_size; ++i ) {497 *m_istr >> code_index >> code;498 flow.set_icode( code_index,code);499 }500 m_istr->ignore(2,'\n'); // '\n' at end of entry501 GenParticle* p = new GenParticle( FourVector(px,py,pz,e),502 id, status, flow,503 Polarization(theta,phi) );504 p->set_generated_mass( m );505 p->suggest_barcode( bar_code );506 //507 // all particles are connected to their end vertex separately508 // after all particles and vertices have been created - so we keep509 // a map of all particles that have end vertices510 if ( end_vtx_code != 0 ) {511 particle_to_end_vertex.addEndParticle(p,end_vtx_code);512 }513 return p;514 }515 516 ParticleData* IO_GenEvent::read_particle_data( ParticleDataTable* pdt ) {517 // assumes mode has already been checked518 //519 // test to be sure the next entry is of type "D" then ignore it520 if ( !(*m_istr) || m_istr->peek()!='D' ) return false;521 m_istr->ignore();522 //523 // read values into temp variables then create new ParticleData object524 char its_name[22];525 int its_id = 0, its_spin = 0;526 double its_charge = 0, its_mass = 0, its_clifetime = 0;527 *m_istr >> its_id >> its_charge >> its_mass528 >> its_clifetime >> its_spin;529 m_istr->ignore(1); // eat the " "530 m_istr->getline( its_name, 22, '\n' );531 ParticleData* pdata = new ParticleData( its_name, its_id, its_charge,532 its_mass, its_clifetime,533 double(its_spin)/2.);534 pdt->insert(pdata);535 return pdata;536 }537 538 bool IO_GenEvent::write_end_listing() {539 if ( m_finished_first_event_io && (m_ostr != NULL) ) {540 m_common_io.write_IO_GenEvent_End(*m_ostr);541 *m_ostr << std::flush;542 m_finished_first_event_io = 0;543 return true;544 }545 return false;546 }547 177 548 178 } // HepMC
Note:
See TracChangeset
for help on using the changeset viewer.
