hepmc - Blame information for rev 128

Subversion Repositories:
Rev:
Rev Author Line No. Line
2 garren 1 //--------------------------------------------------------------------------
2  
3 //////////////////////////////////////////////////////////////////////////
4 // Matt.Dobbs@Cern.CH, January 2000
5 // event input/output in ascii format for machine reading
6 //////////////////////////////////////////////////////////////////////////
7  
8 #include "HepMC/IO_Ascii.h"
9 #include "HepMC/GenEvent.h"
10 #include "HepMC/ParticleDataTable.h"
117 garren 11 #include "HepMC/Version.h"
2 garren 12  
13 namespace HepMC {
14  
15     IO_Ascii::IO_Ascii( const char* filename, std::ios::openmode mode )
16         : m_mode(mode), m_file(filename, mode), m_finished_first_event_io(0)
17     {
128 garren 18         std::cout << "-------------------------------------------------------" << std::endl;
19         std::cout << "Use of HepMC/IO_Ascii is deprecated" << std::endl;
20         std::cout << "-------------------------------------------------------" << std::endl;
2 garren 21         if ( (m_mode&std::ios::out && m_mode&std::ios::in) ||
22              (m_mode&std::ios::app && m_mode&std::ios::in) ) {
23             std::cerr << "IO_Ascii::IO_Ascii Error, open of file requested "
24                       << "of input AND output type. Not allowed. Closing file."
25                       << std::endl;
26             m_file.close();
128 garren 27             return;
2 garren 28         }
29         // precision 16 (# digits following decimal point) is the minimum that
30         //  will capture the full information stored in a double
31         m_file.precision(16);
32         // we use decimal to store integers, because it is smaller than hex!
33         m_file.setf(std::ios::dec,std::ios::basefield);
34         m_file.setf(std::ios::scientific,std::ios::floatfield);
35     }
36  
37     IO_Ascii::~IO_Ascii() {
38         write_end_listing();
39         m_file.close();
40     }
41  
42     void IO_Ascii::print( std::ostream& ostr ) const {
43         ostr << "IO_Ascii: unformated ascii file IO for machine reading.\n"
44              << "\tFile openmode: " << m_mode
45              << " file state: " << m_file.rdstate()
46              << " bad:" << (m_file.rdstate()&std::ios::badbit)
47              << " eof:" << (m_file.rdstate()&std::ios::eofbit)
48              << " fail:" << (m_file.rdstate()&std::ios::failbit)
49              << " good:" << (m_file.rdstate()&std::ios::goodbit) << std::endl;
50     }
51  
52     void IO_Ascii::write_event( const GenEvent* evt ) {
65 garren 53         /// Writes evt to m_file. It does NOT delete the event after writing.
2 garren 54         //
55         // check the state of m_file is good, and that it is in output mode
56         if ( !evt || !m_file ) return;
57         if ( !m_mode&std::ios::out ) {
58             std::cerr << "HepMC::IO_Ascii::write_event "
59                       << " attempt to write to input file." << std::endl;
60             return;
61         }
62         //
63         // write event listing key before first event only.
64         if ( !m_finished_first_event_io ) {
65             m_finished_first_event_io = 1;
117 garren 66             m_file << "\n" << "HepMC::Version " << versionName();
2 garren 67             m_file << "\n" << "HepMC::IO_Ascii-START_EVENT_LISTING\n";
68         }
69         //
70         // output the event data including the number of primary vertices
71         //  and the total number of vertices
72         std::vector<long int> random_states = evt->random_states();
73         m_file << 'E';
74         output( evt->event_number() );
75         output( evt->event_scale() );
76         output( evt->alphaQCD() );
77         output( evt->alphaQED() );
78         output( evt->signal_process_id() );
79         output(   ( evt->signal_process_vertex() ?
80                     evt->signal_process_vertex()->barcode() : 0 )   );
81         output( evt->vertices_size() ); // total number of vertices.
82         output( (int)random_states.size() );
83         for ( std::vector<long int>::iterator rs = random_states.begin();
84               rs != random_states.end(); ++rs ) {
85             output( *rs );
86         }
87         output( (int)evt->weights().size() );
88         for ( WeightContainer::const_iterator w = evt->weights().begin();
89               w != evt->weights().end(); ++w ) {
90             output( *w );
91         }
92         output('\n');
93         //
94         // Output all of the vertices - note there is no real order.
95         for ( GenEvent::vertex_const_iterator v = evt->vertices_begin();
96               v != evt->vertices_end(); ++v ) {
97             write_vertex( *v );
98         }
99     }
100  
101     bool IO_Ascii::fill_next_event( GenEvent* evt ){
102         //
103         //
104         // test that evt pointer is not null
105         if ( !evt ) {
106             std::cerr
107                 << "IO_Ascii::fill_next_event error - passed null event."
108                 << std::endl;
109             return 0;
110         }
111         // check the state of m_file is good, and that it is in input mode
112         if ( !m_file ) return 0;
113         if ( !(m_mode&std::ios::in) ) {
114             std::cerr << "HepMC::IO_Ascii::fill_next_event "
115                       << " attempt to read from output file." << std::endl;
116             return 0;
117         }
118         //
119         // search for event listing key before first event only.
120         //
121         // skip through the file just after first occurence of the start_key
122         if ( !m_finished_first_event_io ) {
123             m_file.seekg( 0 ); // go to position zero in the file.
124             if (!search_for_key_end( m_file,
125                                      "HepMC::IO_Ascii-START_EVENT_LISTING\n")){
126                 std::cerr << "IO_Ascii::fill_next_event start key not found "
127                           << "setting badbit." << std::endl;
128                 m_file.clear(std::ios::badbit);
129                 return 0;
130             }
131             m_finished_first_event_io = 1;
132         }
133         //
134         // test to be sure the next entry is of type "E" then ignore it
135         if ( !m_file || m_file.peek()!='E' ) {
136             // if the E is not the next entry, then check to see if it is
137             // the end event listing key - if yes, search for another start key
138             if ( eat_key(m_file, "HepMC::IO_Ascii-END_EVENT_LISTING\n") ) {
139                 bool search_result = search_for_key_end(m_file,
140                                       "HepMC::IO_Ascii-START_EVENT_LISTING\n");
141                 if ( !search_result ) {
142                     // this is the only case where we set an EOF state
143                     m_file.clear(std::ios::eofbit);
144                     return 0;
145                 }
146             } else {
147                 std::cerr << "IO_Ascii::fill_next_event end key not found "
148                           << "setting badbit." << std::endl;
149                 m_file.clear(std::ios::badbit);
150                 return 0;
151             }
152         }
153         m_file.ignore();
154         // read values into temp variables, then create a new GenEvent
155         int event_number = 0, signal_process_id = 0, signal_process_vertex = 0,
156             num_vertices = 0, random_states_size = 0, weights_size = 0;
157         double eventScale = 0, alpha_qcd = 0, alpha_qed = 0;
158         m_file >> event_number >> eventScale >> alpha_qcd >> alpha_qed
159                >> signal_process_id >> signal_process_vertex
160                >> num_vertices >> random_states_size;
161         std::vector<long int> random_states(random_states_size);
162         for ( int i = 0; i < random_states_size; ++i ) {
163             m_file >> random_states[i];
164         }
165         m_file >> weights_size;
166         WeightContainer weights(weights_size);
167         for ( int ii = 0; ii < weights_size; ++ii ) m_file >> weights[ii];
168         m_file.ignore(2,'\n');
169         //
170         // fill signal_process_id, event_number, weights, random_states
171         evt->set_signal_process_id( signal_process_id );
172         evt->set_event_number( event_number );
173         evt->weights() = weights;
174         evt->set_random_states( random_states );
175         //
176         // the end vertices of the particles are not connected until
177         //  after the event is read --- we store the values in a map until then
178         std::map<GenParticle*,int> particle_to_end_vertex;
179         //
180         // read in the vertices
181         for ( int iii = 1; iii <= num_vertices; ++iii ) {
182             GenVertex* v = read_vertex(particle_to_end_vertex);
183             evt->add_vertex( v );
184         }
185         // set the signal process vertex
186         if ( signal_process_vertex ) {
187             evt->set_signal_process_vertex(
188                 evt->barcode_to_vertex(signal_process_vertex) );
189         }
190         //
191         // last connect particles to their end vertices
192         for ( std::map<GenParticle*,int>::iterator pmap
193                   = particle_to_end_vertex.begin();
194               pmap != particle_to_end_vertex.end(); ++pmap ) {
195             GenVertex* itsDecayVtx = evt->barcode_to_vertex(pmap->second);
196             if ( itsDecayVtx ) itsDecayVtx->add_particle_in( pmap->first );
197             else {
198                 std::cerr << "IO_Ascii::fill_next_event ERROR particle points"
199                           << "\n to null end vertex. " <<std::endl;
200             }
201         }
202         return 1;
203     }
204  
205     void IO_Ascii::write_comment( const std::string comment ) {
206         // check the state of m_file is good, and that it is in output mode
207         if ( !m_file ) return;
208         if ( !m_mode&std::ios::out ) {
209             std::cerr << "HepMC::IO_Ascii::write_particle_data_table "
210                       << " attempt to write to input file." << std::endl;
211             return;
212         }
213         // write end of event listing key if events have already been written
214         write_end_listing();
215         // insert the comment key before the comment
216         m_file << "\n" << "HepMC::IO_Ascii-COMMENT\n";
217         m_file << comment << std::endl;
218     }
219  
220     void IO_Ascii::write_particle_data_table( const ParticleDataTable* pdt) {
221         //
222         // check the state of m_file is good, and that it is in output mode
223         if ( !m_file ) return;
224         if ( !m_mode&std::ios::out ) {
225             std::cerr << "HepMC::IO_Ascii::write_particle_data_table "
226                       << " attempt to write to input file." << std::endl;
227             return;
228         }
229         // write end of event listing key if events have already been written
230         write_end_listing();
231         //
232         m_file << "\n" << "HepMC::IO_Ascii-START_PARTICLE_DATA\n";
233         for ( ParticleDataTable::const_iterator pd = pdt->begin();
234               pd != pdt->end(); pd++ ) {
235             write_particle_data( pd->second );
236         }
237         m_file << "HepMC::IO_Ascii-END_PARTICLE_DATA\n" << std::flush;
238     }
239  
240     bool IO_Ascii::fill_particle_data_table( ParticleDataTable* pdt ) {
241         //
242         // test that pdt pointer is not null
243         if ( !pdt ) {
244             std::cerr
245                 << "IO_Ascii::fill_particle_data_table - passed null table."
246                 << std::endl;
247             return 0;
248         }
249         //
250         // check the state of m_file is good, and that it is in input mode
251         if ( !m_file ) return 0;
252         if ( !m_mode&std::ios::in ) {
253             std::cerr << "HepMC::IO_Ascii::fill_particle_data_table "
254                       << " attempt to read from output file." << std::endl;
255             return 0;
256         }
257         // position to beginning of file
258         int initial_file_position = m_file.tellg();
259         std::ios::iostate initial_state = m_file.rdstate();
260         m_file.seekg( 0 );
261         // skip through the file just after first occurence of the start_key
262         if (!search_for_key_end( m_file,
263                                  "HepMC::IO_Ascii-START_PARTICLE_DATA\n")) {
264             m_file.seekg( initial_file_position );
265             std::cerr << "IO_Ascii::fill_particle_data_table start key not  "
266                       << "found setting badbit." << std::endl;
267             m_file.clear(std::ios::badbit);
268             return 0;
269         }
270         //
271         pdt->set_description("Read with IO_Ascii");
272         //
273         // read Individual GenParticle data entries
274         while ( read_particle_data( pdt ) );
275         //
276         // eat end_key
277         if ( !eat_key(m_file,"HepMC::IO_Ascii-END_PARTICLE_DATA\n") ){
278             std::cerr << "IO_Ascii::fill_particle_data_table end key not  "
279                       << "found setting badbit." << std::endl;
280             m_file.clear(std::ios::badbit);
281         }
282         // put the file back into its original state and position
283         m_file.clear( initial_state );
284         m_file.seekg( initial_file_position );
285         return 1;
286     }
287  
288     void IO_Ascii::write_vertex( GenVertex* v ) {
289         // assumes mode has already been checked
290         if ( !v || !m_file ) {
291             std::cerr << "IO_Ascii::write_vertex !v||!m_file, "
292                       << "v="<< v << " setting badbit" << std::endl;
293             m_file.clear(std::ios::badbit);
294             return;
295         }
296         // First collect info we need
297         // count the number of orphan particles going into v
298         int num_orphans_in = 0;
299         for ( GenVertex::particles_in_const_iterator p1
300                   = v->particles_in_const_begin();
301               p1 != v->particles_in_const_end(); ++p1 ) {
302             if ( !(*p1)->production_vertex() ) ++num_orphans_in;
303         }
304         //
305         m_file << 'V';
306         output( v->barcode() ); // v's unique identifier
307         output( v->id() );
308         output( v->position().x() );
309         output( v->position().y() );
310         output( v->position().z() );
311         output( v->position().t() );
312         output( num_orphans_in );
313         output( (int)v->particles_out_size() );
314         output( (int)v->weights().size() );
315         for ( WeightContainer::iterator w = v->weights().begin();
316               w != v->weights().end(); ++w ) {
317             output( *w );
318         }
319         output('\n');
320         for ( GenVertex::particles_in_const_iterator p2
321                   = v->particles_in_const_begin();
322               p2 != v->particles_in_const_end(); ++p2 ) {
323             if ( !(*p2)->production_vertex() ) {
324                 write_particle( *p2 );
325             }
326         }
327         for ( GenVertex::particles_out_const_iterator p3
328                   = v->particles_out_const_begin();
329               p3 != v->particles_out_const_end(); ++p3 ) {
330             write_particle( *p3 );
331         }
332     }
333  
334     void IO_Ascii::write_particle( GenParticle* p ) {
335         // assumes mode has already been checked
336         if ( !p || !m_file ) {
337             std::cerr << "IO_Ascii::write_vertex !p||!m_file, "
338                       << "v="<< p << " setting badbit" << std::endl;
339             m_file.clear(std::ios::badbit);
340             return;
341         }
342         m_file << 'P';
343         output( p->barcode() );
344         output( p->pdg_id() );
345         output( p->momentum().px() );
346         output( p->momentum().py() );
347         output( p->momentum().pz() );
348         output( p->momentum().e() );
349         output( p->status() );
350         output( p->polarization().theta() );
351         output( p->polarization().phi() );
352         // since end_vertex is oftentimes null, this CREATES a null vertex
353         // in the map
354         output(   ( p->end_vertex() ? p->end_vertex()->barcode() : 0 )  );
355         m_file << ' ' << p->flow() << "\n";
356     }
357  
358     void IO_Ascii::write_particle_data( const ParticleData* pdata ) {
359         // assumes mode has already been checked
360         if ( !pdata || !m_file ) {
361             std::cerr << "IO_Ascii::write_vertex !pdata||!m_file, "
362                       << "pdata="<< pdata << " setting badbit" << std::endl;
363             m_file.clear(std::ios::badbit);
364             return;
365         }
366         m_file << 'D';
367         output( pdata->pdg_id() );
368         output( pdata->charge() );
369         output( pdata->mass() );
370         output( pdata->clifetime() );
371         output( (int)(pdata->spin()*2.+.1) );
372         // writes the first 21 characters starting with 0
373         m_file << " " << pdata->name().substr(0,21) << "\n";
374     }
375  
376     GenVertex* IO_Ascii::read_vertex
377     ( std::map<GenParticle*,int>& particle_to_end_vertex )
378     {
379         // assumes mode has already been checked
380         //
381         // test to be sure the next entry is of type "V" then ignore it
382         if ( !m_file || m_file.peek()!='V' ) {
383             std::cerr << "IO_Ascii::read_vertex setting badbit." << std::endl;
384             m_file.clear(std::ios::badbit);
385             return 0;
386         }
387         m_file.ignore();
388         // read values into temp variables, then create a new GenVertex object
389         int identifier =0, id =0, num_orphans_in =0,
390             num_particles_out = 0, weights_size = 0;
391         double x = 0., y = 0., z = 0., t = 0.;
392         m_file >> identifier >> id >> x >> y >> z >> t
393                >> num_orphans_in >> num_particles_out >> weights_size;
394         WeightContainer weights(weights_size);
395         for ( int i1 = 0; i1 < weights_size; ++i1 ) m_file >> weights[i1];
396         m_file.ignore(2,'\n');
43 garren 397         GenVertex* v = new GenVertex( FourVector(x,y,z,t),
2 garren 398                                 id, weights);
399         v->suggest_barcode( identifier );
400         //
401         // read and create the associated particles. outgoing particles are
402         //  added to their production vertices immediately, while incoming
403         //  particles are added to a map and handles later.
404         for ( int i2 = 1; i2 <= num_orphans_in; ++i2 ) {
405             read_particle(particle_to_end_vertex);
406         }
407         for ( int i3 = 1; i3 <= num_particles_out; ++i3 ) {
408             v->add_particle_out( read_particle(particle_to_end_vertex) );
409         }
410         return v;
411     }
412  
413     GenParticle* IO_Ascii::read_particle(
414         std::map<GenParticle*,int>& particle_to_end_vertex ){
415         // assumes mode has already been checked
416         //
417         // test to be sure the next entry is of type "P" then ignore it
418         if ( !m_file || m_file.peek()!='P' ) {
419             std::cerr << "IO_Ascii::read_particle setting badbit."
420                       << std::endl;
421             m_file.clear(std::ios::badbit);
422             return 0;
423         }
424         m_file.ignore();
425         //
426         // declare variables to be read in to, and read everything except flow
427         double px = 0., py = 0., pz = 0., e = 0., theta = 0., phi = 0.;
428         int bar_code = 0, id = 0, status = 0, end_vtx_code = 0, flow_size = 0;
429         m_file >> bar_code >> id >> px >> py >> pz >> e >> status
430                >> theta >> phi >> end_vtx_code >> flow_size;
431         //
432         // read flow patterns if any exist
433         Flow flow;
434         int code_index, code;
435         for ( int i = 1; i <= flow_size; ++i ) {
436             m_file >> code_index >> code;
437             flow.set_icode( code_index,code);
438         }
439         m_file.ignore(2,'\n'); // '\n' at end of entry
43 garren 440         GenParticle* p = new GenParticle( FourVector(px,py,pz,e),
2 garren 441                                     id, status, flow,
442                                     Polarization(theta,phi) );
443         p->suggest_barcode( bar_code );
444         //
445         // all particles are connected to their end vertex separately
446         // after all particles and vertices have been created - so we keep
447         // a map of all particles that have end vertices
448         if ( end_vtx_code != 0 ) particle_to_end_vertex[p] = end_vtx_code;
449         return p;
450     }
451  
452     ParticleData* IO_Ascii::read_particle_data( ParticleDataTable* pdt ) {
453         // assumes mode has already been checked
454         //
455         // test to be sure the next entry is of type "D" then ignore it
456         if ( !m_file || m_file.peek()!='D' ) return 0;
457         m_file.ignore();
458         //
459         // read values into temp variables then create new ParticleData object
460         char its_name[22];
461         int its_id = 0, its_spin = 0;  
462         double its_charge = 0, its_mass = 0, its_clifetime = 0;
463         m_file >> its_id >> its_charge >> its_mass
464                >> its_clifetime >> its_spin;
465         m_file.ignore(1); // eat the " "
466         m_file.getline( its_name, 22, '\n' );
467         ParticleData* pdata = new ParticleData( its_name, its_id, its_charge,
468                                                 its_mass, its_clifetime,
469                                                 double(its_spin)/2.);
470         pdt->insert(pdata);
471         return pdata;
472     }
473  
474     bool IO_Ascii::write_end_listing() {
475         if ( m_finished_first_event_io && m_mode&std::ios::out ) {
476             m_file << "HepMC::IO_Ascii-END_EVENT_LISTING\n" << std::flush;
477             m_finished_first_event_io = 0;
478             return 1;
479         }
480         return 0;
481     }
482  
483     bool IO_Ascii::search_for_key_end( std::istream& in, const char* key ) {
65 garren 484         /// reads characters from in until the string of characters matching
485         /// key is found (success) or EOF is reached (failure).
486         /// It stops immediately thereafter. Returns T/F for success/fail
2 garren 487         //
488         char c[1];
489         unsigned int index = 0;
490         while ( in.get(c[0]) ) {
491             if ( c[0] == key[index] ) {
492                 ++index;
493             } else { index = 0; }
494             if ( index == strlen(key) ) return 1;
495         }
496         return 0;
497     }
498  
499     bool IO_Ascii::search_for_key_beginning( std::istream& in,
500                                              const char* key ) {
65 garren 501         /// not tested and NOT used anywhere!
2 garren 502         if ( search_for_key_end( in, key) ) {
503             int i = strlen(key);
504             while ( i>=0 ) in.putback(key[i--]);
505             return 1;
506         } else {
507             in.putback(EOF);
508             in.clear();
509             return 0;
510         }
511     }
512  
513     bool IO_Ascii::eat_key( std::iostream& in, const char* key ) {
65 garren 514         /// eats the character string key from istream in - only if the key
515         /// is the very next occurence in the stream
516         /// if the key is not the next occurence, it eats nothing ... i.e.
517         ///  it puts back whatever it would have eaten.
2 garren 518         int key_length = strlen(key);
519         // below is the only way I know of to get a variable length string
520         //  conforming to ansi standard.
521         char* c = new char[key_length +1];
522         int i=0;
523         // read the stream until get fails (because of EOF), a character
524         //  doesn't match a character in the string, or all characters in
525         //  the string have been checked and match.
526         while ( in.get(c[i]) && c[i]==key[i] && i<key_length ) {
527             ++i;
528         }
529         if ( i == key_length ) {
86 garren 530             delete [] c;
2 garren 531             return 1;
532         }
533         //
534         // if we get here, then we have eaten the wrong this and we must put it
535         // back
86 garren 536         //---------> non standard code --- probably does not work
2 garren 537         while ( i>=0 ) in.putback(c[i--]);
538         delete c;
539         return 0;
540     }
541  
542     int IO_Ascii::find_in_map( const std::map<GenVertex*,int>& m,
543                                GenVertex* v ) const {
544         std::map<GenVertex*,int>::const_iterator iter = m.find(v);
545         if ( iter == m.end() ) return 0;
546         return iter->second;
547     }
548  
549 } // HepMC
550  
551  
552  
553  
554  
555