hepmc - Blame information for rev 86

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