Fork me on GitHub

source: svn/trunk/Utilities/HepMC/src/IO_GenEvent.cc@ 364

Last change on this file since 364 was 349, checked in by severine ovyn, 16 years ago

first test
File size: 16.8 KB
RevLine 
[349]1//--------------------------------------------------------------------------
2
3//////////////////////////////////////////////////////////////////////////
4// garren@fnal.gov, July 2006
5// event input/output in ascii format for machine reading
6// IO_GenEvent format contains HeavyIon and PdfInfo classes
7//////////////////////////////////////////////////////////////////////////
8
9#include "IO_GenEvent.h"
10#include "GenEvent.h"
11#include "ParticleDataTable.h"
12#include "HeavyIon.h"
13#include "PdfInfo.h"
14#include "CommonIO.h"
15#include "Version.h"
16
17namespace HepMC {
18
19 IO_GenEvent::IO_GenEvent( const std::string& filename, std::ios::openmode mode )
20 : m_mode(mode),
21 m_file(filename.c_str(), mode),
22 m_ostr(0),
23 m_istr(0),
24 m_iostr(0),
25 m_finished_first_event_io(false),
26 m_have_file(false),
27 m_common_io()
28 {
29 if ( (m_mode&std::ios::out && m_mode&std::ios::in) ||
30 (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 << std::endl;
34 m_file.close();
35 return;
36 }
37 // precision 16 (# digits following decimal point) is the minimum that
38 // will capture the full information stored in a double
39 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 // now we set the streams
44 m_iostr = &m_file;
45 if ( m_mode&std::ios::in ) {
46 m_istr = &m_file;
47 m_ostr = NULL;
48 }
49 if ( m_mode&std::ios::out ) {
50 m_ostr = &m_file;
51 m_istr = NULL;
52 }
53 m_have_file = true;
54 }
55
56
57 IO_GenEvent::IO_GenEvent( std::istream & istr )
58 : m_ostr(0),
59 m_istr(&istr),
60 m_iostr(&istr),
61 m_finished_first_event_io(false),
62 m_have_file(false),
63 m_common_io()
64 { }
65
66 IO_GenEvent::IO_GenEvent( std::ostream & ostr )
67 : m_ostr(&ostr),
68 m_istr(0),
69 m_iostr(&ostr),
70 m_finished_first_event_io(false),
71 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 }
81
82 IO_GenEvent::~IO_GenEvent() {
83 write_end_listing();
84 if(m_have_file) m_file.close();
85 }
86
87 void IO_GenEvent::use_input_units( Units::MomentumUnit mom, Units::LengthUnit len ) {
88 m_common_io.use_input_units(mom,len);
89 }
90
91 void IO_GenEvent::print( std::ostream& ostr ) const {
92 ostr << "IO_GenEvent: unformated ascii file IO for machine reading.\n";
93 if(m_have_file) ostr << "\tFile openmode: " << m_mode ;
94 ostr << " stream state: " << m_ostr->rdstate()
95 << " bad:" << (m_ostr->rdstate()&std::ios::badbit)
96 << " eof:" << (m_ostr->rdstate()&std::ios::eofbit)
97 << " fail:" << (m_ostr->rdstate()&std::ios::failbit)
98 << " good:" << (m_ostr->rdstate()&std::ios::goodbit) << std::endl;
99 }
100
101 void IO_GenEvent::write_event( const GenEvent* evt ) {
102 /// Writes evt to output stream. It does NOT delete the event after writing.
103 //
104 // make sure the state is good, and that it is in output mode
105 if ( !evt ) return;
106 if ( m_ostr == NULL ) {
107 std::cerr << "HepMC::IO_GenEvent::write_event "
108 << " attempt to write to input file." << std::endl;
109 return;
110 }
111 //
112 // 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;
232 }
233
234 void IO_GenEvent::write_comment( const std::string comment ) {
235 // make sure the stream is good, and that it is in output mode
236 if ( !(*m_ostr) ) return;
237 if ( m_ostr == NULL ) {
238 std::cerr << "HepMC::IO_GenEvent::write_comment "
239 << " attempt to write to input file." << std::endl;
240 return;
241 }
242 // write end of event listing key if events have already been written
243 write_end_listing();
244 // insert the comment key before the comment
245 *m_ostr << "\n" << "HepMC::IO_GenEvent-COMMENT\n";
246 *m_ostr << comment << std::endl;
247 }
248
249 void IO_GenEvent::write_vertex( GenVertex* v ) {
250 // assumes mode has already been checked
251 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 need
258 // count the number of orphan particles going into v
259 int num_orphans_in = 0;
260 for ( GenVertex::particles_in_const_iterator p1
261 = 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 identifier
268 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 p2
282 = 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 p3
289 = 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 checked
314 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 default
322 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 checked
358 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 default
366 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 readable
400 *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 checked
407 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 vertex
425 // in the map
426 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 checked
432 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 0
445 *m_ostr << " " << pdata->name().substr(0,21) << "\n";
446 }
447
448 HeavyIon* IO_GenEvent::read_heavy_ion()
449 {
450 // assumes mode has already been checked
451 //
452 // test to be sure the next entry is of type "H" then ignore it
453 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 object
460 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 >> prot
464 >> 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 checked
477 //
478 // test to be sure the next entry is of type "P" then ignore it
479 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 flow
488 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 >> status
491 >> theta >> phi >> end_vtx_code >> flow_size;
492 //
493 // read flow patterns if any exist
494 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 entry
501 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 separately
508 // after all particles and vertices have been created - so we keep
509 // a map of all particles that have end vertices
510 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 checked
518 //
519 // test to be sure the next entry is of type "D" then ignore it
520 if ( !(*m_istr) || m_istr->peek()!='D' ) return false;
521 m_istr->ignore();
522 //
523 // read values into temp variables then create new ParticleData object
524 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_mass
528 >> 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
548} // HepMC
Note: See TracBrowser for help on using the repository browser.