Fork me on GitHub

source: svn/trunk/Utilities/HepMC/src/GenEventStreamIO.cc@ 609

Last change on this file since 609 was 571, checked in by cp3-support, 13 years ago

upgrade HepMC to version 2.06.05

File size: 25.6 KB
Line 
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 "GenEvent.h"
16#include "GenCrossSection.h"
17#include "StreamInfo.h"
18#include "StreamHelpers.h"
19#include "Version.h"
20#include "IO_Exception.h"
21
22namespace 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().
29void 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.
50template <class IO>
51StreamInfo& 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
72std::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
155std::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 // OK - now ready to start reading the event, so set the header flag
208 info.set_reading_event_header(true);
209 // The flag will be set to false when we reach the end of the header
210 while(info.reading_event_header()) {
211 switch(is.peek()) {
212 case 'E':
213 { // deal with the event line
214 process_event_line( is, num_vertices, bp1, bp2, signal_process_vertex );
215 } break;
216 case 'N':
217 { // get weight names
218 read_weight_names( is );
219 } break;
220 case 'U':
221 { // get unit information if it exists
222 if( info.io_type() == gen ) {
223 read_units( is );
224 }
225 } break;
226 case 'C':
227 { // we have a GenCrossSection line
228 // create cross section
229 GenCrossSection xs;
230 // check for invalid data
231 try {
232 // read the line
233 xs.read(is);
234 }
235 catch (IO_Exception& e) {
236 detail::find_event_end( is );
237 }
238 if(xs.is_set()) {
239 set_cross_section( xs );
240 }
241 } break;
242 case 'H':
243 { // we have a HeavyIon line OR an unexpected HepMC... line
244 if( info.io_type() == gen || info.io_type() == extascii ) {
245 // get HeavyIon
246 HeavyIon ion;
247 // check for invalid data
248 try {
249 is >> &ion;
250 }
251 catch (IO_Exception& e) {
252 detail::find_event_end( is );
253 }
254 if(ion.is_valid()) {
255 set_heavy_ion( ion );
256 }
257 }
258 } break;
259 case 'F':
260 { // we have a PdfInfo line
261 if( info.io_type() == gen || info.io_type() == extascii ) {
262 // get PdfInfo
263 PdfInfo pdf;
264 // check for invalid data
265 try {
266 is >> &pdf;
267 }
268 catch (IO_Exception& e) {
269 detail::find_event_end( is );
270 }
271 if(pdf.is_valid()) {
272 set_pdf_info( pdf );
273 }
274 }
275 } break;
276 case 'V':
277 {
278 // this should be the first vertex line - exit this loop
279 info.set_reading_event_header(false);
280 } break;
281 case 'P':
282 { // we should not find this line
283 std::cerr << "streaming input: found unexpected line P" << std::endl;
284 info.set_reading_event_header(false);
285 } break;
286 default:
287 // ignore everything else
288 break;
289 } // switch on line type
290 } // while reading_event_header
291 //
292 // the end vertices of the particles are not connected until
293 // after the event is read --- we store the values in a map until then
294 TempParticleMap particle_to_end_vertex;
295 //
296 // read in the vertices
297 for ( int iii = 1; iii <= num_vertices; ++iii ) {
298 GenVertex* v = new GenVertex();
299 try {
300 detail::read_vertex(is,particle_to_end_vertex,v);
301 }
302 catch (IO_Exception& e) {
303 for( TempParticleMap::orderIterator it = particle_to_end_vertex.order_begin();
304 it != particle_to_end_vertex.order_end(); ++it ) {
305 GenParticle* p = it->second;
306 // delete particles only if they are not already owned by a vertex
307 if( p->production_vertex() ) {
308 } else if( p->end_vertex() ) {
309 } else {
310 delete p;
311 }
312 }
313 delete v;
314 detail::find_event_end( is );
315 }
316 add_vertex( v );
317 }
318 // set the signal process vertex
319 if ( signal_process_vertex ) {
320 set_signal_process_vertex(
321 barcode_to_vertex(signal_process_vertex) );
322 }
323 //
324 // last connect particles to their end vertices
325 GenParticle* beam1(0);
326 GenParticle* beam2(0);
327 for ( TempParticleMap::orderIterator pmap
328 = particle_to_end_vertex.order_begin();
329 pmap != particle_to_end_vertex.order_end(); ++pmap ) {
330 GenParticle* p = pmap->second;
331 int vtx = particle_to_end_vertex.end_vertex( p );
332 GenVertex* itsDecayVtx = barcode_to_vertex(vtx);
333 if ( itsDecayVtx ) itsDecayVtx->add_particle_in( p );
334 else {
335 std::cerr << "read_io_genevent: ERROR particle points"
336 << " to null end vertex. " <<std::endl;
337 }
338 // also look for the beam particles
339 if( p->barcode() == bp1 ) beam1 = p;
340 if( p->barcode() == bp2 ) beam2 = p;
341 }
342 set_beam_particles(beam1,beam2);
343 return is;
344}
345
346// ------------------------- operator << and operator >> ----------------
347
348std::ostream & operator << (std::ostream & os, GenEvent & evt)
349{
350 /// Writes evt to an output stream.
351 evt.write(os);
352 return os;
353}
354
355std::istream & operator >> (std::istream & is, GenEvent & evt)
356{
357 evt.read(is);
358 return is;
359}
360
361// ------------------------- set units ----------------
362
363std::istream & set_input_units(std::istream & is,
364 Units::MomentumUnit mom,
365 Units::LengthUnit len )
366{
367 //
368 StreamInfo & info = get_stream_info(is);
369 info.use_input_units( mom, len );
370 return is;
371}
372
373// ------------------------- begin and end block lines ----------------
374
375std::ostream & write_HepMC_IO_block_begin(std::ostream & os )
376{
377 //
378 StreamInfo & info = get_stream_info(os);
379
380 if( !info.finished_first_event() ) {
381 os << "\n" << "HepMC::Version " << versionName();
382 os << "\n";
383 os << info.IO_GenEvent_Key() << "\n";
384 }
385 return os;
386}
387
388std::ostream & write_HepMC_IO_block_end(std::ostream & os )
389{
390 //
391 StreamInfo & info = get_stream_info(os);
392
393 if( info.finished_first_event() ) {
394 os << info.IO_GenEvent_End() << "\n";
395 os << std::flush;
396 }
397 return os;
398}
399
400std::istream & GenEvent::process_event_line( std::istream & is,
401 int & num_vertices,
402 int & bp1, int & bp2,
403 int & signal_process_vertex )
404{
405 //
406 if ( !is ) {
407 std::cerr << "GenEvent::process_event_line setting badbit." << std::endl;
408 is.clear(std::ios::badbit);
409 return is;
410 }
411 //
412 StreamInfo & info = get_stream_info(is);
413 std::string line;
414 std::getline(is,line);
415 std::istringstream iline(line);
416 std::string firstc;
417 iline >> firstc;
418 //
419 // read values into temp variables, then fill GenEvent
420 int event_number = 0, signal_process_id = 0,
421 random_states_size = 0, nmpi = -1;
422 double eventScale = 0, alpha_qcd = 0, alpha_qed = 0;
423 iline >> event_number;
424 if(!iline) detail::find_event_end( is );
425 if( info.io_type() == gen || info.io_type() == extascii ) {
426 iline >> nmpi;
427 if(!iline) detail::find_event_end( is );
428 set_mpi( nmpi );
429 }
430 iline >> eventScale ;
431 if(!iline) detail::find_event_end( is );
432 iline >> alpha_qcd ;
433 if(!iline) detail::find_event_end( is );
434 iline >> alpha_qed;
435 if(!iline) detail::find_event_end( is );
436 iline >> signal_process_id ;
437 if(!iline) detail::find_event_end( is );
438 iline >> signal_process_vertex;
439 if(!iline) detail::find_event_end( is );
440 iline >> num_vertices;
441 if(!iline) detail::find_event_end( is );
442 if( info.io_type() == gen || info.io_type() == extascii ) {
443 iline >> bp1 ;
444 if(!iline) detail::find_event_end( is );
445 iline >> bp2;
446 if(!iline) detail::find_event_end( is );
447 }
448 iline >> random_states_size;
449 if(!iline) detail::find_event_end( is );
450 std::vector<long> random_states(random_states_size);
451 for ( int i = 0; i < random_states_size; ++i ) {
452 iline >> random_states[i];
453 if(!iline) detail::find_event_end( is );
454 }
455 WeightContainer::size_type weights_size = 0;
456 iline >> weights_size;
457 if(!iline) detail::find_event_end( is );
458 std::vector<double> wgt(weights_size);
459 for ( WeightContainer::size_type ii = 0; ii < weights_size; ++ii ) {
460 iline >> wgt[ii];
461 if(!iline) detail::find_event_end( is );
462 }
463 // weight names will be added later if they exist
464 if( weights_size > 0 ) m_weights = wgt;
465 //
466 // fill signal_process_id, event_number, random_states, etc.
467 set_signal_process_id( signal_process_id );
468 set_event_number( event_number );
469 set_random_states( random_states );
470 set_event_scale( eventScale );
471 set_alphaQCD( alpha_qcd );
472 set_alphaQED( alpha_qed );
473 //
474 return is;
475}
476
477std::istream & GenEvent::read_weight_names( std::istream & is )
478{
479 // now check for a named weight line
480 if ( !is ) {
481 std::cerr << "GenEvent::read_weight_names setting badbit." << std::endl;
482 is.clear(std::ios::badbit);
483 return is;
484 }
485 // Test to be sure the next entry is of type "N"
486 // If we have no named weight line, this is not an error
487 // releases prior to 2.06.00 do not have named weights
488 if ( is.peek() !='N') {
489 return is;
490 }
491 // now get this line and process it
492 std::string line;
493 std::getline(is,line);
494 std::istringstream wline(line);
495 std::string firstc;
496 WeightContainer::size_type name_size = 0;
497 wline >> firstc >> name_size;
498 if(!wline) detail::find_event_end( is );
499 if( firstc != "N") {
500 std::cout << "debug: first character of named weights is " << firstc << std::endl;
501 std::cout << "debug: We should never get here" << std::endl;
502 is.clear(std::ios::badbit);
503 return is;
504 }
505 if( m_weights.size() != name_size ) {
506 std::cout << "debug: weight sizes do not match "<< std::endl;
507 std::cout << "debug: weight vector size is " << m_weights.size() << std::endl;
508 std::cout << "debug: weight name size is " << name_size << std::endl;
509 is.clear(std::ios::badbit);
510 return is;
511 }
512 std::string name;
513 std::string::size_type i1 = line.find("\"");
514 std::string::size_type i2;
515 std::string::size_type len = line.size();
516 WeightContainer namedWeight;
517 for ( WeightContainer::size_type ii = 0; ii < name_size; ++ii ) {
518 // weight names may contain blanks
519 if(i1 >= len) {
520 std::cout << "debug: attempting to read past the end of the named weight line " << std::endl;
521 std::cout << "debug: We should never get here" << std::endl;
522 std::cout << "debug: Looking for the end of this event" << std::endl;
523 detail::find_event_end( is );
524 }
525 i2 = line.find("\"",i1+1);
526 name = line.substr(i1+1,i2-i1-1);
527 namedWeight[name] = m_weights[ii];
528 i1 = line.find("\"",i2+1);
529 }
530 m_weights = namedWeight;
531 return is;
532}
533
534std::istream & GenEvent::read_units( std::istream & is )
535{
536 //
537 if ( !is ) {
538 std::cerr << "GenEvent::read_units setting badbit." << std::endl;
539 is.clear(std::ios::badbit);
540 return is;
541 }
542 //
543 StreamInfo & info = get_stream_info(is);
544 // test to be sure the next entry is of type "U" then ignore it
545 // if we have no units, this is not an error
546 // releases prior to 2.04.00 did not write unit information
547 if ( is.peek() !='U') {
548 use_units( info.io_momentum_unit(),
549 info.io_position_unit() );
550 return is;
551 }
552 is.ignore(); // ignore the first character in the line
553 std::string mom, pos;
554 is >> mom >> pos;
555 is.ignore(1); // eat the extra whitespace
556 use_units(mom,pos);
557 //
558 return is;
559}
560
561std::istream & GenEvent::find_file_type( std::istream & istr )
562{
563 //
564 // make sure the stream is good
565 if ( !istr ) return istr;
566
567 //
568 StreamInfo & info = get_stream_info(istr);
569
570 // if there is no input block line, then we assume this stream
571 // is in the IO_GenEvent format
572 if ( istr.peek()=='E' ) {
573 info.set_io_type( gen );
574 info.set_has_key(false);
575 return istr;
576 }
577
578 std::string line;
579 while ( std::getline(istr,line) ) {
580 //
581 // search for event listing key before first event only.
582 //
583 if( line == info.IO_GenEvent_Key() ) {
584 info.set_io_type( gen );
585 info.set_has_key(true);
586 return istr;
587 } else if( line == info.IO_Ascii_Key() ) {
588 info.set_io_type( ascii );
589 info.set_has_key(true);
590 return istr;
591 } else if( line == info.IO_ExtendedAscii_Key() ) {
592 info.set_io_type( extascii );
593 info.set_has_key(true);
594 return istr;
595 } else if( line == info.IO_Ascii_PDT_Key() ) {
596 info.set_io_type( ascii_pdt );
597 info.set_has_key(true);
598 return istr;
599 } else if( line == info.IO_ExtendedAscii_PDT_Key() ) {
600 info.set_io_type( extascii_pdt );
601 info.set_has_key(true);
602 return istr;
603 }
604 }
605 info.set_io_type( 0 );
606 info.set_has_key(false);
607 return istr;
608}
609
610std::istream & GenEvent::find_end_key( std::istream & istr, int & iotype )
611{
612 iotype = 0;
613 // peek at the first character before proceeding
614 if( istr.peek()!='H' ) return istr;
615 //
616 // we only check the next line
617 std::string line;
618 std::getline(istr,line);
619 //
620 StreamInfo & info = get_stream_info(istr);
621 //
622 // check to see if this is an end key
623 if( line == info.IO_GenEvent_End() ) {
624 iotype = gen;
625 } else if( line == info.IO_Ascii_End() ) {
626 iotype = ascii;
627 } else if( line == info.IO_ExtendedAscii_End() ) {
628 iotype = extascii;
629 } else if( line == info.IO_Ascii_PDT_End() ) {
630 iotype = ascii_pdt;
631 } else if( line == info.IO_ExtendedAscii_PDT_End() ) {
632 iotype = extascii_pdt;
633 }
634 if( iotype != 0 && info.io_type() != iotype ) {
635 std::cerr << "GenEvent::find_end_key: iotype keys have changed" << std::endl;
636 } else {
637 return istr;
638 }
639 //
640 // if we get here, then something has gotten badly confused
641 std::cerr << "GenEvent::find_end_key: MALFORMED INPUT" << std::endl;
642 istr.clear(std::ios::badbit);
643 return istr;
644}
645
646std::ostream & establish_output_stream_info( std::ostream & os )
647{
648 StreamInfo & info = get_stream_info(os);
649 if ( !info.finished_first_event() ) {
650 // precision 16 (# digits following decimal point) is the minimum that
651 // will capture the full information stored in a double
652 os.precision(16);
653 // we use decimal to store integers, because it is smaller than hex!
654 os.setf(std::ios::dec,std::ios::basefield);
655 os.setf(std::ios::scientific,std::ios::floatfield);
656 }
657 return os;
658}
659
660std::istream & establish_input_stream_info( std::istream & is )
661{
662 StreamInfo & info = get_stream_info(is);
663 if ( !info.finished_first_event() ) {
664 // precision 16 (# digits following decimal point) is the minimum that
665 // will capture the full information stored in a double
666 is.precision(16);
667 // we use decimal to store integers, because it is smaller than hex!
668 is.setf(std::ios::dec,std::ios::basefield);
669 is.setf(std::ios::scientific,std::ios::floatfield);
670 }
671 return is;
672}
673
674
675// ------------------------- helper functions ----------------
676
677namespace detail {
678
679// The functions defined here need to use get_stream_info
680
681std::istream & read_particle( std::istream & is,
682 TempParticleMap & particle_to_end_vertex,
683 GenParticle * p )
684{
685 // get the next line
686 std::string line;
687 std::getline(is,line);
688 std::istringstream iline(line);
689 std::string firstc;
690 iline >> firstc;
691 if( firstc != "P" ) {
692 std::cerr << "StreamHelpers::detail::read_particle invalid line type: "
693 << firstc << std::endl;
694 std::cerr << "StreamHelpers::detail::read_particle setting badbit."
695 << std::endl;
696 is.clear(std::ios::badbit);
697 return is;
698 }
699 //
700 StreamInfo & info = get_stream_info(is);
701 //testHepMC.cc
702 // declare variables to be read in to, and read everything except flow
703 double px = 0., py = 0., pz = 0., e = 0., m = 0., theta = 0., phi = 0.;
704 int bar_code = 0, id = 0, status = 0, end_vtx_code = 0, flow_size = 0;
705 // check that the input stream is still OK after reading item
706 iline >> bar_code ;
707 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
708 iline >> id ;
709 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
710 iline >> px ;
711 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
712 iline >> py ;
713 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
714 iline >> pz ;
715 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
716 iline >> e ;
717 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
718 if( info.io_type() != ascii ) {
719 iline >> m ;
720 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
721 }
722 iline >> status ;
723 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
724 iline >> theta ;
725 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
726 iline >> phi ;
727 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
728 iline >> end_vtx_code ;
729 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
730 iline >> flow_size;
731 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
732 //
733 // read flow patterns if any exist
734 Flow flow;
735 int code_index, code;
736 for ( int i = 1; i <= flow_size; ++i ) {
737 iline >> code_index >> code;
738 if(!iline) { delete p; throw IO_Exception("read_particle input stream encounterd invalid data"); }
739 flow.set_icode( code_index,code);
740 }
741 p->set_momentum( FourVector(px,py,pz,e) );
742 p->set_pdg_id( id );
743 p->set_status( status );
744 p->set_flow( flow );
745 p->set_polarization( Polarization(theta,phi) );
746 if( info.io_type() == ascii ) {
747 p->set_generated_mass( p->momentum().m() );
748 } else {
749 p->set_generated_mass( m );
750 }
751 p->suggest_barcode( bar_code );
752 //
753 // all particles are connected to their end vertex separately
754 // after all particles and vertices have been created - so we keep
755 // a map of all particles that have end vertices
756 if ( end_vtx_code != 0 ) {
757 particle_to_end_vertex.addEndParticle(p,end_vtx_code);
758 }
759 return is;
760}
761
762std::istream & read_units( std::istream & is, GenEvent & evt )
763{
764 //
765 if ( !is ) {
766 std::cerr << "StreamHelpers read_units setting badbit." << std::endl;
767 is.clear(std::ios::badbit);
768 return is;
769 }
770 //
771 StreamInfo & info = get_stream_info(is);
772 // test to be sure the next entry is of type "U" then ignore it
773 // if we have no units, this is not an error
774 // releases prior to 2.04.00 did not write unit information
775 if ( is.peek() !='U') {
776 evt.use_units( info.io_momentum_unit(),
777 info.io_position_unit() );
778 return is;
779 }
780 is.ignore(); // ignore the first character in the line
781 std::string mom, pos;
782 is >> mom >> pos;
783 is.ignore(1); // eat the extra whitespace
784 evt.use_units(mom,pos);
785 //
786 return is;
787}
788
789std::ostream & establish_output_stream_info( std::ostream & os )
790{
791 StreamInfo & info = get_stream_info(os);
792 if ( !info.finished_first_event() ) {
793 // precision 16 (# digits following decimal point) is the minimum that
794 // will capture the full information stored in a double
795 os.precision(16);
796 // we use decimal to store integers, because it is smaller than hex!
797 os.setf(std::ios::dec,std::ios::basefield);
798 os.setf(std::ios::scientific,std::ios::floatfield);
799 }
800 return os;
801}
802
803std::istream & establish_input_stream_info( std::istream & is )
804{
805 StreamInfo & info = get_stream_info(is);
806 if ( !info.finished_first_event() ) {
807 // precision 16 (# digits following decimal point) is the minimum that
808 // will capture the full information stored in a double
809 is.precision(16);
810 // we use decimal to store integers, because it is smaller than hex!
811 is.setf(std::ios::dec,std::ios::basefield);
812 is.setf(std::ios::scientific,std::ios::floatfield);
813 }
814 return is;
815}
816
817} // detail
818
819} // HepMC
Note: See TracBrowser for help on using the repository browser.