hepmc - Rev 129

Subversion Repositories:
Rev:
  1. //--------------------------------------------------------------------------
  2.  
  3. //////////////////////////////////////////////////////////////////////////
  4. // garren@fnal.gov, July 2006
  5. // event input/output in ascii format for machine reading
  6. // extended format contains HeavyIon and PdfInfo classes
  7. //////////////////////////////////////////////////////////////////////////
  8.  
  9. #include "HepMC/IO_GenEvent.h"
  10. #include "HepMC/GenEvent.h"
  11. #include "HepMC/ParticleDataTable.h"
  12. #include "HepMC/HeavyIon.h"
  13. #include "HepMC/PdfInfo.h"
  14. #include "HepMC/Version.h"
  15.  
  16. namespace HepMC {
  17.  
  18.     IO_GenEvent::IO_GenEvent( const char* filename, std::ios::openmode mode )
  19.     : m_mode(mode), m_file(filename, mode), m_finished_first_event_io(0)
  20.     {
  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_GenEvent::IO_GenEvent Error, open of file requested "
  24.                       << "of input AND output type. Not allowed. Closing file."
  25.                       << std::endl;
  26.             m_file.close();
  27.             return;
  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.         // now we set the streams
  36.         m_iostr = &m_file;
  37.         if ( m_mode&std::ios::in ) {
  38.             m_istr = &m_file;
  39.             m_ostr = NULL;
  40.         }
  41.         if ( m_mode&std::ios::out ) {
  42.             m_ostr = &m_file;
  43.             m_istr = NULL;
  44.         }
  45.         m_have_file = true;
  46.     }
  47.  
  48.     IO_GenEvent::IO_GenEvent( std::istream * istr )
  49.     : m_istr(istr),
  50.       m_iostr(istr),
  51.       m_finished_first_event_io(0)
  52.     {
  53.         m_ostr = NULL;
  54.         m_have_file = false;
  55.     }
  56.  
  57.     IO_GenEvent::IO_GenEvent( std::ostream * ostr )
  58.     : m_ostr(ostr),
  59.       m_iostr(ostr),
  60.       m_finished_first_event_io(0)
  61.     {
  62.         m_istr = NULL;
  63.         m_have_file = false;
  64.         // precision 16 (# digits following decimal point) is the minimum that
  65.         //  will capture the full information stored in a double
  66.         m_ostr->precision(16);
  67.         // we use decimal to store integers, because it is smaller than hex!
  68.         m_ostr->setf(std::ios::dec,std::ios::basefield);
  69.         m_ostr->setf(std::ios::scientific,std::ios::floatfield);
  70.     }
  71.  
  72.     IO_GenEvent::~IO_GenEvent() {
  73.         write_end_listing();
  74.         if(m_have_file) m_file.close();
  75.     }
  76.  
  77.     void IO_GenEvent::print( std::ostream& ostr ) const {
  78.         ostr << "IO_GenEvent: unformated ascii file IO for machine reading.\n";
  79.         if(m_have_file)    ostr  << "\tFile openmode: " << m_mode ;
  80.         ostr << " stream state: " << m_ostr->rdstate()
  81.              << " bad:" << (m_ostr->rdstate()&std::ios::badbit)
  82.              << " eof:" << (m_ostr->rdstate()&std::ios::eofbit)
  83.              << " fail:" << (m_ostr->rdstate()&std::ios::failbit)
  84.              << " good:" << (m_ostr->rdstate()&std::ios::goodbit) << std::endl;
  85.     }
  86.  
  87.     void IO_GenEvent::write_event( const GenEvent* evt ) {
  88.         /// Writes evt to m_file. It does NOT delete the event after writing.
  89.         //
  90.         // make sure the state is good, and that it is in output mode
  91.         if ( !evt  ) return;
  92.         if ( m_ostr == NULL ) {
  93.             std::cerr << "HepMC::IO_GenEvent::write_event "
  94.                       << " attempt to write to input file." << std::endl;
  95.             return;
  96.         }
  97.         //
  98.         // write event listing key before first event only.
  99.         if ( !m_finished_first_event_io ) {
  100.             m_finished_first_event_io = 1;
  101.             *m_ostr << "\n" << "HepMC::Version " << versionName();
  102.             *m_ostr << "\n" << "HepMC::IO_GenEvent-START_EVENT_LISTING\n";
  103.         }
  104.         //
  105.         // output the event data including the number of primary vertices
  106.         //  and the total number of vertices
  107.         std::vector<long int> random_states = evt->random_states();
  108.         *m_ostr << 'E';
  109.         output( evt->event_number() );
  110.         output( evt->mpi() );
  111.         output( evt->event_scale() );
  112.         output( evt->alphaQCD() );
  113.         output( evt->alphaQED() );
  114.         output( evt->signal_process_id() );
  115.         output(   ( evt->signal_process_vertex() ?
  116.                     evt->signal_process_vertex()->barcode() : 0 )   );
  117.         output( evt->vertices_size() ); // total number of vertices.
  118.         write_beam_particles( evt->beam_particles() );
  119.         output( (int)random_states.size() );
  120.         for ( std::vector<long int>::iterator rs = random_states.begin();
  121.               rs != random_states.end(); ++rs ) {
  122.             output( *rs );
  123.         }
  124.         output( (int)evt->weights().size() );
  125.         for ( WeightContainer::const_iterator w = evt->weights().begin();
  126.               w != evt->weights().end(); ++w ) {
  127.             output( *w );
  128.         }
  129.         output('\n');
  130.         write_heavy_ion( evt->heavy_ion() );
  131.         write_pdf_info( evt->pdf_info() );
  132.         //
  133.         // Output all of the vertices - note there is no real order.
  134.         for ( GenEvent::vertex_const_iterator v = evt->vertices_begin();
  135.               v != evt->vertices_end(); ++v ) {
  136.             write_vertex( *v );
  137.         }
  138.     }
  139.  
  140.     bool IO_GenEvent::fill_next_event( GenEvent* evt ){
  141.         //
  142.         //
  143.         // test that evt pointer is not null
  144.         if ( !evt ) {
  145.             std::cerr
  146.                 << "IO_GenEvent::fill_next_event error - passed null event."
  147.                 << std::endl;
  148.             return false;
  149.         }
  150.         // make sure the stream is good, and that it is in input mode
  151.         if ( !(*m_istr) ) return false;
  152.         if ( !m_istr ) {
  153.             std::cerr << "HepMC::IO_GenEvent::fill_next_event "
  154.                       << " attempt to read from output file." << std::endl;
  155.             return false;
  156.         }
  157.         //
  158.         // search for event listing key before first event only.
  159.         //
  160.         // skip through the file just after first occurence of the start_key
  161.         if ( !m_finished_first_event_io ) {
  162.             //m_file.seekg( 0 ); // go to position zero in the file.
  163.             if (!search_for_key_end( m_file,
  164.                                      "HepMC::IO_GenEvent-START_EVENT_LISTING\n")){
  165.                 std::cerr << "IO_GenEvent::fill_next_event start key not found "
  166.                           << "setting badbit." << std::endl;
  167.                 m_file.clear(std::ios::badbit);
  168.                 return false;
  169.             }
  170.             m_finished_first_event_io = 1;
  171.         }
  172.         //
  173.         // test to be sure the next entry is of type "E" then ignore it
  174.         if ( !(*m_istr) || m_istr->peek()!='E' ) {
  175.             // if the E is not the next entry, then check to see if it is
  176.             // the end event listing key - if yes, search for another start key
  177.             if ( eat_key(*m_istr, "HepMC::IO_GenEvent-END_EVENT_LISTING\n") ) {
  178.                 bool search_result = search_for_key_end(*m_istr,
  179.                                       "HepMC::IO_GenEvent-START_EVENT_LISTING\n");
  180.                 if ( !search_result ) {
  181.                     // this is the only case where we set an EOF state
  182.                     m_istr->clear(std::ios::eofbit);
  183.                     return false;
  184.                 }
  185.             } else {
  186.                 std::cerr << "IO_GenEvent::fill_next_event end key not found "
  187.                           << "setting badbit." << std::endl;
  188.                 m_istr->clear(std::ios::badbit);
  189.                 return false;
  190.             }
  191.         }
  192.         m_istr->ignore();
  193.         // read values into temp variables, then create a new GenEvent
  194.         int event_number = 0, signal_process_id = 0, signal_process_vertex = 0,
  195.             num_vertices = 0, random_states_size = 0, weights_size = 0,
  196.             nmpi = 0, bp1 = 0, bp2 = 0;
  197.         double eventScale = 0, alpha_qcd = 0, alpha_qed = 0;
  198.         *m_istr >> event_number >> nmpi >> eventScale >> alpha_qcd >> alpha_qed
  199.                >> signal_process_id >> signal_process_vertex
  200.                >> num_vertices >> bp1 >> bp2 >> random_states_size;
  201.         std::vector<long int> random_states(random_states_size);
  202.         for ( int i = 0; i < random_states_size; ++i ) {
  203.             *m_istr >> random_states[i];
  204.         }
  205.         *m_istr >> weights_size;
  206.         WeightContainer weights(weights_size);
  207.         for ( int ii = 0; ii < weights_size; ++ii ) *m_istr >> weights[ii];
  208.         m_istr->ignore(2,'\n');
  209.         //
  210.         // fill signal_process_id, event_number, weights, random_states
  211.         evt->set_signal_process_id( signal_process_id );
  212.         evt->set_event_number( event_number );
  213.         evt->set_mpi( nmpi );
  214.         evt->weights() = weights;
  215.         evt->set_random_states( random_states );
  216.         // get HeavyIon and PdfInfo
  217.         HeavyIon* ion = read_heavy_ion();
  218.         if(ion) evt->set_heavy_ion( *ion );
  219.         PdfInfo* pdf = read_pdf_info();
  220.         if(pdf) evt->set_pdf_info( *pdf );
  221.         //
  222.         // the end vertices of the particles are not connected until
  223.         //  after the event is read --- we store the values in a map until then
  224.         std::map<GenParticle*,int> particle_to_end_vertex;
  225.         //
  226.         // read in the vertices
  227.         for ( int iii = 1; iii <= num_vertices; ++iii ) {
  228.             GenVertex* v = read_vertex(particle_to_end_vertex);
  229.             evt->add_vertex( v );
  230.         }
  231.         // set the signal process vertex
  232.         if ( signal_process_vertex ) {
  233.             evt->set_signal_process_vertex(
  234.                 evt->barcode_to_vertex(signal_process_vertex) );
  235.         }
  236.         //
  237.         // last connect particles to their end vertices
  238.         GenParticle* beam1(0);
  239.         GenParticle* beam2(0);
  240.         for ( std::map<GenParticle*,int>::iterator pmap
  241.                   = particle_to_end_vertex.begin();
  242.               pmap != particle_to_end_vertex.end(); ++pmap ) {
  243.             GenVertex* itsDecayVtx = evt->barcode_to_vertex(pmap->second);
  244.             if ( itsDecayVtx ) itsDecayVtx->add_particle_in( pmap->first );
  245.             else {
  246.                 std::cerr << "IO_GenEvent::fill_next_event ERROR particle points"
  247.                           << "\n to null end vertex. " <<std::endl;
  248.             }
  249.             // also look for the beam particles
  250.             GenParticle* p =  pmap->first;
  251.             if( p->barcode() == bp1 ) beam1 = p;
  252.             if( p->barcode() == bp2 ) beam2 = p;
  253.         }
  254.         evt->set_beam_particles(beam1,beam2);
  255.         return true;
  256.     }
  257.  
  258.     void IO_GenEvent::write_comment( const std::string comment ) {
  259.         // make sure the stream is good, and that it is in output mode
  260.         if ( !(*m_ostr) ) return;
  261.         if ( m_ostr == NULL ) {
  262.             std::cerr << "HepMC::IO_GenEvent::write_comment "
  263.                       << " attempt to write to input file." << std::endl;
  264.             return;
  265.         }
  266.         // write end of event listing key if events have already been written
  267.         write_end_listing();
  268.         // insert the comment key before the comment
  269.         *m_ostr << "\n" << "HepMC::IO_GenEvent-COMMENT\n";
  270.         *m_ostr << comment << std::endl;
  271.     }
  272.  
  273.     void IO_GenEvent::write_particle_data_table( const ParticleDataTable* pdt) {
  274.         //
  275.         // make sure the stream is good, and that it is in output mode
  276.         if ( !(*m_ostr) ) return;
  277.         if ( m_ostr == NULL ) {
  278.             std::cerr << "HepMC::IO_GenEvent::write_particle_data_table "
  279.                       << " attempt to write to input file." << std::endl;
  280.             return;
  281.         }
  282.         // write end of event listing key if events have already been written
  283.         write_end_listing();
  284.         //
  285.         *m_ostr << "\n" << "HepMC::IO_GenEvent-START_PARTICLE_DATA\n";
  286.         for ( ParticleDataTable::const_iterator pd = pdt->begin();
  287.               pd != pdt->end(); pd++ ) {
  288.             write_particle_data( pd->second );
  289.         }
  290.         *m_ostr << "HepMC::IO_GenEvent-END_PARTICLE_DATA\n" << std::flush;
  291.     }
  292.  
  293.     bool IO_GenEvent::fill_particle_data_table( ParticleDataTable* pdt ) {
  294.         //
  295.         // test that pdt pointer is not null
  296.         if ( !pdt ) {
  297.             std::cerr
  298.                 << "IO_GenEvent::fill_particle_data_table - passed null table."
  299.                 << std::endl;
  300.             return false;
  301.         }
  302.         //
  303.         // check the state of m_file is good, and that it is in input mode
  304.         if ( !(*m_istr) ) return false;
  305.         if ( m_istr == NULL ) {
  306.             std::cerr << "HepMC::IO_GenEvent::fill_particle_data_table "
  307.                       << " attempt to read from output file." << std::endl;
  308.             return false;
  309.         }
  310.         // position to beginning of file
  311.         int initial_file_position = m_istr->tellg();
  312.         std::ios::iostate initial_state = m_istr->rdstate();
  313.         m_istr->seekg( 0 );
  314.         // skip through the file just after first occurence of the start_key
  315.         if (!search_for_key_end( m_file,
  316.                                  "HepMC::IO_GenEvent-START_PARTICLE_DATA\n")) {
  317.             m_istr->seekg( initial_file_position );
  318.             std::cerr << "IO_GenEvent::fill_particle_data_table start key not  "
  319.                       << "found setting badbit." << std::endl;
  320.             m_istr->clear(std::ios::badbit);
  321.             return false;
  322.         }
  323.         //
  324.         pdt->set_description("Read with IO_GenEvent");
  325.         //
  326.         // read Individual GenParticle data entries
  327.         while ( read_particle_data( pdt ) );
  328.         //
  329.         // eat end_key
  330.         if ( !eat_key(m_file,"HepMC::IO_GenEvent-END_PARTICLE_DATA\n") ){
  331.             std::cerr << "IO_GenEvent::fill_particle_data_table end key not  "
  332.                       << "found setting badbit." << std::endl;
  333.             m_istr->clear(std::ios::badbit);
  334.         }
  335.         // put the file back into its original state and position
  336.         m_istr->clear( initial_state );
  337.         m_istr->seekg( initial_file_position );
  338.         return true;
  339.     }
  340.  
  341.     void IO_GenEvent::write_vertex( GenVertex* v ) {
  342.         // assumes mode has already been checked
  343.         if ( !v || !(*m_ostr) ) {
  344.             std::cerr << "IO_GenEvent::write_vertex !v||!(*m_ostr), "
  345.                       << "v="<< v << " setting badbit" << std::endl;
  346.             m_ostr->clear(std::ios::badbit);
  347.             return;
  348.         }
  349.         // First collect info we need
  350.         // count the number of orphan particles going into v
  351.         int num_orphans_in = 0;
  352.         for ( GenVertex::particles_in_const_iterator p1
  353.                   = v->particles_in_const_begin();
  354.               p1 != v->particles_in_const_end(); ++p1 ) {
  355.             if ( !(*p1)->production_vertex() ) ++num_orphans_in;
  356.         }
  357.         //
  358.         *m_ostr << 'V';
  359.         output( v->barcode() ); // v's unique identifier
  360.         output( v->id() );
  361.         output( v->position().x() );
  362.         output( v->position().y() );
  363.         output( v->position().z() );
  364.         output( v->position().t() );
  365.         output( num_orphans_in );
  366.         output( (int)v->particles_out_size() );
  367.         output( (int)v->weights().size() );
  368.         for ( WeightContainer::iterator w = v->weights().begin();
  369.               w != v->weights().end(); ++w ) {
  370.             output( *w );
  371.         }
  372.         output('\n');
  373.         for ( GenVertex::particles_in_const_iterator p2
  374.                   = v->particles_in_const_begin();
  375.               p2 != v->particles_in_const_end(); ++p2 ) {
  376.             if ( !(*p2)->production_vertex() ) {
  377.                 write_particle( *p2 );
  378.             }
  379.         }
  380.         for ( GenVertex::particles_out_const_iterator p3
  381.                   = v->particles_out_const_begin();
  382.               p3 != v->particles_out_const_end(); ++p3 ) {
  383.             write_particle( *p3 );
  384.         }
  385.     }
  386.  
  387.     void IO_GenEvent::write_beam_particles(
  388.         std::pair<GenParticle *,GenParticle *> pr ) {
  389.         GenParticle* p = pr.first;
  390.         //m_file << 'B';
  391.         if(!p) {
  392.            output( 0 );
  393.         } else {
  394.            output( p->barcode() );
  395.         }
  396.         p = pr.second;
  397.         if(!p) {
  398.            output( 0 );
  399.         } else {
  400.            output( p->barcode() );
  401.         }
  402.     }
  403.  
  404.     void IO_GenEvent::write_heavy_ion( HeavyIon* ion ) {
  405.         // assumes mode has already been checked
  406.         if ( !(*m_ostr) ) {
  407.             std::cerr << "IO_GenEvent::write_heavy_ion !(*m_ostr), "
  408.                       << " setting badbit" << std::endl;
  409.             m_ostr->clear(std::ios::badbit);
  410.             return;
  411.         }
  412.         *m_ostr << 'H';
  413.         // HeavyIon* is set to 0 by default
  414.         if ( !ion  ) {
  415.             output( 0 );
  416.             output( 0 );
  417.             output( 0 );
  418.             output( 0 );
  419.             output( 0 );
  420.             output( 0 );
  421.             output( 0 );
  422.             output( 0 );
  423.             output( 0 );
  424.             output( 0. );
  425.             output( 0. );
  426.             output( 0. );
  427.             output( 0. );
  428.             output('\n');
  429.             return;
  430.         }
  431.         //
  432.         output( ion->Ncoll_hard() );
  433.         output( ion->Npart_proj() );
  434.         output( ion->Npart_targ() );
  435.         output( ion->Ncoll() );
  436.         output( ion->spectator_neutrons() );
  437.         output( ion->spectator_protons() );
  438.         output( ion->N_Nwounded_collisions() );
  439.         output( ion->Nwounded_N_collisions() );
  440.         output( ion->Nwounded_Nwounded_collisions() );
  441.         output( ion->impact_parameter() );
  442.         output( ion->event_plane_angle() );
  443.         output( ion->eccentricity() );
  444.         output( ion->sigma_inel_NN() );
  445.         output('\n');
  446.     }
  447.  
  448.     void IO_GenEvent::write_pdf_info( PdfInfo* pdf ) {
  449.         // assumes mode has already been checked
  450.         if ( !(*m_ostr) ) {
  451.             std::cerr << "IO_GenEvent::write_pdf_info !(*m_ostr), "
  452.                       << " setting badbit" << std::endl;
  453.             m_ostr->clear(std::ios::badbit);
  454.             return;
  455.         }
  456.         *m_ostr << 'F';
  457.         // PdfInfo* is set to 0 by default
  458.         if ( !pdf ) {
  459.             output( 0 );
  460.             output( 0 );
  461.             output( 0. );
  462.             output( 0. );
  463.             output( 0. );
  464.             output( 0. );
  465.             output( 0. );
  466.             output('\n');
  467.             return;
  468.         }
  469.         //
  470.         output( pdf->id1() );
  471.         output( pdf->id2() );
  472.         output( pdf->x1() );
  473.         output( pdf->x2() );
  474.         output( pdf->scalePDF() );
  475.         output( pdf->pdf1() );
  476.         output( pdf->pdf2() );
  477.         output('\n');
  478.     }
  479.  
  480.     void IO_GenEvent::write_particle( GenParticle* p ) {
  481.         // assumes mode has already been checked
  482.         if ( !p || !(*m_ostr) ) {
  483.             std::cerr << "IO_GenEvent::write_particle !p||!(*m_ostr), "
  484.                       << "v="<< p << " setting badbit" << std::endl;
  485.             m_ostr->clear(std::ios::badbit);
  486.             return;
  487.         }
  488.         *m_ostr << 'P';
  489.         output( p->barcode() );
  490.         output( p->pdg_id() );
  491.         output( p->momentum().px() );
  492.         output( p->momentum().py() );
  493.         output( p->momentum().pz() );
  494.         output( p->momentum().e() );
  495.         output( p->generated_mass() );
  496.         output( p->status() );
  497.         output( p->polarization().theta() );
  498.         output( p->polarization().phi() );
  499.         // since end_vertex is oftentimes null, this CREATES a null vertex
  500.         // in the map
  501.         output(   ( p->end_vertex() ? p->end_vertex()->barcode() : 0 )  );
  502.         *m_ostr << ' ' << p->flow() << "\n";
  503.     }
  504.  
  505.     void IO_GenEvent::write_particle_data( const ParticleData* pdata ) {
  506.         // assumes mode has already been checked
  507.         if ( !pdata || !(*m_ostr) ) {
  508.             std::cerr << "IO_GenEvent::write_particle_data !pdata||!(*m_ostr), "
  509.                       << "pdata="<< pdata << " setting badbit" << std::endl;
  510.             m_ostr->clear(std::ios::badbit);
  511.             return;
  512.         }
  513.         *m_ostr << 'D';
  514.         output( pdata->pdg_id() );
  515.         output( pdata->charge() );
  516.         output( pdata->mass() );
  517.         output( pdata->clifetime() );
  518.         output( (int)(pdata->spin()*2.+.1) );
  519.         // writes the first 21 characters starting with 0
  520.         *m_ostr << " " << pdata->name().substr(0,21) << "\n";
  521.     }
  522.  
  523.     GenVertex* IO_GenEvent::read_vertex
  524.     ( std::map<GenParticle*,int>& particle_to_end_vertex )
  525.     {
  526.         // assumes mode has already been checked
  527.         //
  528.         // test to be sure the next entry is of type "V" then ignore it
  529.         if ( !(*m_istr) || m_istr->peek()!='V' ) {
  530.             std::cerr << "IO_GenEvent::read_vertex setting badbit." << std::endl;
  531.             m_istr->clear(std::ios::badbit);
  532.             return false;
  533.         }
  534.         m_istr->ignore();
  535.         // read values into temp variables, then create a new GenVertex object
  536.         int identifier =0, id =0, num_orphans_in =0,
  537.             num_particles_out = 0, weights_size = 0;
  538.         double x = 0., y = 0., z = 0., t = 0.;
  539.         *m_istr >> identifier >> id >> x >> y >> z >> t
  540.                >> num_orphans_in >> num_particles_out >> weights_size;
  541.         WeightContainer weights(weights_size);
  542.         for ( int i1 = 0; i1 < weights_size; ++i1 ) *m_istr >> weights[i1];
  543.         m_istr->ignore(2,'\n');
  544.         GenVertex* v = new GenVertex( FourVector(x,y,z,t),
  545.                                 id, weights);
  546.         v->suggest_barcode( identifier );
  547.         //
  548.         // read and create the associated particles. outgoing particles are
  549.         //  added to their production vertices immediately, while incoming
  550.         //  particles are added to a map and handles later.
  551.         for ( int i2 = 1; i2 <= num_orphans_in; ++i2 ) {
  552.             read_particle(particle_to_end_vertex);
  553.         }
  554.         for ( int i3 = 1; i3 <= num_particles_out; ++i3 ) {
  555.             v->add_particle_out( read_particle(particle_to_end_vertex) );
  556.         }
  557.         return v;
  558.     }
  559.  
  560.     HeavyIon* IO_GenEvent::read_heavy_ion()
  561.     {
  562.         // assumes mode has already been checked
  563.         //
  564.         // test to be sure the next entry is of type "H" then ignore it
  565.         if ( !(*m_istr) || m_istr->peek()!='H' ) {
  566.             std::cerr << "IO_GenEvent::read_heavy_ion setting badbit." << std::endl;
  567.             m_istr->clear(std::ios::badbit);
  568.             return false;
  569.         }
  570.         m_istr->ignore();
  571.         // read values into temp variables, then create a new HeavyIon object
  572.         int nh =0, np =0, nt =0, nc =0,
  573.             neut = 0, prot = 0, nw =0, nwn =0, nwnw =0;
  574.         float impact = 0., plane = 0., xcen = 0., inel = 0.;
  575.         *m_istr >> nh >> np >> nt >> nc >> neut >> prot
  576.                >> nw >> nwn >> nwnw >> impact >> plane >> xcen >> inel;
  577.         m_istr->ignore(2,'\n');
  578.         if( nh == 0 ) return false;
  579.         HeavyIon* ion = new HeavyIon(nh, np, nt, nc, neut, prot,
  580.                                      nw, nwn, nwnw,
  581.                                      impact, plane, xcen, inel );
  582.         //
  583.         return ion;
  584.     }
  585.  
  586.     PdfInfo* IO_GenEvent::read_pdf_info()
  587.     {
  588.         // assumes mode has already been checked
  589.         //
  590.         // test to be sure the next entry is of type "F" then ignore it
  591.         if ( !(*m_istr) || m_istr->peek() !='F') {
  592.             std::cerr << "IO_GenEvent::read_pdf_info setting badbit." << std::endl;
  593.             m_istr->clear(std::ios::badbit);
  594.             return false;
  595.         }
  596.         m_istr->ignore();
  597.         // read values into temp variables, then create a new PdfInfo object
  598.         int id1 =0, id2 =0;
  599.         double  x1 = 0., x2 = 0., scale = 0., pdf1 = 0., pdf2 = 0.;
  600.         *m_istr >> id1 >> id2 >> x1 >> x2 >> scale >> pdf1 >> pdf2;
  601.         m_istr->ignore(2,'\n');
  602.         if( id1 == 0 ) return false;
  603.         PdfInfo* pdf = new PdfInfo( id1, id2, x1, x2, scale, pdf1, pdf2);
  604.         //
  605.         return pdf;
  606.     }
  607.  
  608.     GenParticle* IO_GenEvent::read_particle(
  609.         std::map<GenParticle*,int>& particle_to_end_vertex ){
  610.         // assumes mode has already been checked
  611.         //
  612.         // test to be sure the next entry is of type "P" then ignore it
  613.         if ( !(*m_istr) || m_istr->peek()!='P' ) {
  614.             std::cerr << "IO_GenEvent::read_particle setting badbit."
  615.                       << std::endl;
  616.             m_istr->clear(std::ios::badbit);
  617.             return false;
  618.         }
  619.         m_istr->ignore();
  620.         //
  621.         // declare variables to be read in to, and read everything except flow
  622.         double px = 0., py = 0., pz = 0., e = 0., m = 0., theta = 0., phi = 0.;
  623.         int bar_code = 0, id = 0, status = 0, end_vtx_code = 0, flow_size = 0;
  624.         *m_istr >> bar_code >> id >> px >> py >> pz >> e >> m >> status
  625.                >> theta >> phi >> end_vtx_code >> flow_size;
  626.         //
  627.         // read flow patterns if any exist
  628.         Flow flow;
  629.         int code_index, code;
  630.         for ( int i = 1; i <= flow_size; ++i ) {
  631.             *m_istr >> code_index >> code;
  632.             flow.set_icode( code_index,code);
  633.         }
  634.         m_istr->ignore(2,'\n'); // '\n' at end of entry
  635.         GenParticle* p = new GenParticle( FourVector(px,py,pz,e),
  636.                                     id, status, flow,
  637.                                     Polarization(theta,phi) );
  638.         p->set_generated_mass( m );
  639.         p->suggest_barcode( bar_code );
  640.         //
  641.         // all particles are connected to their end vertex separately
  642.         // after all particles and vertices have been created - so we keep
  643.         // a map of all particles that have end vertices
  644.         if ( end_vtx_code != 0 ) particle_to_end_vertex[p] = end_vtx_code;
  645.         return p;
  646.     }
  647.  
  648.     ParticleData* IO_GenEvent::read_particle_data( ParticleDataTable* pdt ) {
  649.         // assumes mode has already been checked
  650.         //
  651.         // test to be sure the next entry is of type "D" then ignore it
  652.         if ( !(*m_istr) || m_istr->peek()!='D' ) return false;
  653.         m_istr->ignore();
  654.         //
  655.         // read values into temp variables then create new ParticleData object
  656.         char its_name[22];
  657.         int its_id = 0, its_spin = 0;  
  658.         double its_charge = 0, its_mass = 0, its_clifetime = 0;
  659.         *m_istr >> its_id >> its_charge >> its_mass
  660.                >> its_clifetime >> its_spin;
  661.         m_istr->ignore(1); // eat the " "
  662.         m_istr->getline( its_name, 22, '\n' );
  663.         ParticleData* pdata = new ParticleData( its_name, its_id, its_charge,
  664.                                                 its_mass, its_clifetime,
  665.                                                 double(its_spin)/2.);
  666.         pdt->insert(pdata);
  667.         return pdata;
  668.     }
  669.  
  670.     bool IO_GenEvent::write_end_listing() {
  671.         if ( m_finished_first_event_io && (m_ostr != NULL) ) {
  672.             *m_ostr << "HepMC::IO_GenEvent-END_EVENT_LISTING\n" << std::flush;
  673.             m_finished_first_event_io = 0;
  674.             return true;
  675.         }
  676.         return false;
  677.     }
  678.  
  679.     bool IO_GenEvent::search_for_key_end( std::istream& in, const char* key ) {
  680.         /// reads characters from in until the string of characters matching
  681.         /// key is found (success) or EOF is reached (failure).
  682.         /// It stops immediately thereafter. Returns T/F for success/fail
  683.         //
  684.         char c[1];
  685.         unsigned int index = 0;
  686.         while ( in.get(c[0]) ) {
  687.             if ( c[0] == key[index] ) {
  688.                 ++index;
  689.             } else { index = 0; }
  690.             if ( index == strlen(key) ) return true;
  691.         }
  692.         return false;
  693.     }
  694.  
  695.     bool IO_GenEvent::search_for_key_beginning( std::istream& in,
  696.                                              const char* key ) {
  697.         /// not tested and NOT used anywhere!
  698.         if ( search_for_key_end( in, key) ) {
  699.             int i = strlen(key);
  700.             while ( i>=0 ) in.putback(key[i--]);
  701.             return true;
  702.         } else {
  703.             in.putback(EOF);
  704.             in.clear();
  705.             return false;
  706.         }
  707.     }
  708.  
  709.     bool IO_GenEvent::eat_key( std::istream& in, const char* key ) {
  710.         /// eats the character string key from istream in - only if the key
  711.         /// is the very next occurence in the stream
  712.         /// if the key is not the next occurence, it eats nothing ... i.e.
  713.         ///  it puts back whatever it would have eaten.
  714.         int key_length = strlen(key);
  715.         // below is the only way I know of to get a variable length string
  716.         //  conforming to ansi standard.
  717.         char* c = new char[key_length +1];
  718.         int i=0;
  719.         // read the stream until get fails (because of EOF), a character
  720.         //  doesn't match a character in the string, or all characters in
  721.         //  the string have been checked and match.
  722.         while ( in.get(c[i]) && c[i]==key[i] && i<key_length ) {
  723.             ++i;
  724.         }
  725.         if ( i == key_length ) {
  726.             delete [] c;
  727.             return true;
  728.         }
  729.         //
  730.         // if we get here, then we have eaten the wrong this and we must put it
  731.         // back
  732.         while ( i>=0 ) in.putback(c[i--]);
  733.         delete c;
  734.         return false;
  735.     }
  736.  
  737.     int IO_GenEvent::find_in_map( const std::map<GenVertex*,int>& m,
  738.                                GenVertex* v ) const {
  739.         std::map<GenVertex*,int>::const_iterator iter = m.find(v);
  740.         if ( iter == m.end() ) return false;
  741.         return iter->second;
  742.     }
  743.        
  744. } // HepMC
  745.