Subversion Repositories hepmc

Rev

Blame | Last modification | View Log | Download | RSS feed

  1. //--------------------------------------------------------------------------
  2. //
  3. // CommonIO.cc
  4. // Author:  Lynn Garren
  5. //
  6. //  Allowed keys used at the beginning and end of HepMC data dumps
  7. //
  8. // ----------------------------------------------------------------------
  9.  
  10. #include "HepMC/CommonIO.h"
  11. #include "HepMC/GenEvent.h"
  12. #include "HepMC/HeavyIon.h"
  13. #include "HepMC/PdfInfo.h"
  14. #include "HepMC/TempParticleMap.h"
  15.  
  16. namespace HepMC {
  17.  
  18. int CommonIO::find_file_type( std::istream& istr )
  19. {
  20.     std::string line;
  21.     while ( std::getline(istr,line) ) {
  22.         //
  23.         // search for event listing key before first event only.
  24.         //
  25.         if( line == m_io_genevent_start ) {
  26.             std::cout << "begin IO_GenEvent" << std::endl;
  27.             return gen;
  28.         } else if( line == m_io_ascii_start ) {
  29.             std::cout << "begin IO_Ascii" << std::endl;
  30.             return ascii;
  31.         } else if( line == m_io_extendedascii_start ) {
  32.             std::cout << "begin IO_ExtendedAscii" << std::endl;
  33.             return extascii;
  34.         }
  35.     }
  36.     return -1;
  37. }
  38.  
  39. int CommonIO::find_end_key( std::istream& istr )
  40. {
  41.     // peek at the first character before proceeding
  42.     if( istr.peek()!='H' ) return false;
  43.     //
  44.     // we only check the next line
  45.     std::string line;
  46.     std::getline(istr,line);
  47.     //
  48.     // check to see if this is an end key
  49.     if( line == m_io_genevent_end ) {
  50.         std::cout << "end IO_GenEvent" << std::endl;
  51.         return gen;
  52.     } else if( line == m_io_ascii_end ) {
  53.         std::cout << "end IO_Ascii" << std::endl;
  54.         return ascii;
  55.     } else if( line == m_io_extendedascii_end ) {
  56.         std::cout << "end IO_ExtendedAscii" << std::endl;
  57.         return extascii;
  58.     }
  59.     //
  60.     // if we get here, then something has gotten badly confused
  61.     std::cerr << "CommonIO::find_end_key: MALFORMED INPUT" << std::endl;
  62.     return -1;
  63. }
  64.  
  65. bool CommonIO::read_io_ascii_event( std::istream* is, GenEvent* evt )
  66. {
  67.     return false;
  68. }
  69.  
  70. bool CommonIO::read_io_extendedascii_event( std::istream* is, GenEvent* evt )
  71. {
  72.     return false;
  73. }
  74.  
  75. bool CommonIO::read_io_genevent_event( std::istream* is, GenEvent* evt )
  76. {
  77.     /// this method ONLY works if called from fill_next_event
  78.     //
  79.     // assume that we've been handed a valid stream and event pointer
  80.     //
  81.     // read values into temp variables, then fill GenEvent
  82.     int event_number = 0, signal_process_id = 0, signal_process_vertex = 0,
  83.         num_vertices = 0, random_states_size = 0, weights_size = 0,
  84.         nmpi = 0, bp1 = 0, bp2 = 0;
  85.     double eventScale = 0, alpha_qcd = 0, alpha_qed = 0;
  86.     *is >> event_number >> nmpi >> eventScale >> alpha_qcd >> alpha_qed
  87.            >> signal_process_id >> signal_process_vertex
  88.            >> num_vertices >> bp1 >> bp2 >> random_states_size;
  89.     std::vector<long int> random_states(random_states_size);
  90.     for ( int i = 0; i < random_states_size; ++i ) {
  91.         *is >> random_states[i];
  92.     }
  93.     *is >> weights_size;
  94.     WeightContainer weights(weights_size);
  95.     for ( int ii = 0; ii < weights_size; ++ii ) *is >> weights[ii];
  96.     is->ignore(2,'\n');
  97.     //
  98.     // fill signal_process_id, event_number, weights, random_states, etc.
  99.     evt->set_signal_process_id( signal_process_id );
  100.     evt->set_event_number( event_number );
  101.     evt->set_mpi( nmpi );
  102.     evt->weights() = weights;
  103.     evt->set_random_states( random_states );
  104.     evt->set_event_scale( eventScale );
  105.     evt->set_alphaQCD( alpha_qcd );
  106.     evt->set_alphaQED( alpha_qed );
  107.     // get HeavyIon and PdfInfo
  108.     HeavyIon* ion = read_heavy_ion(is);
  109.     if(ion) evt->set_heavy_ion( *ion );
  110.     PdfInfo* pdf = read_pdf_info(is);
  111.     if(pdf) evt->set_pdf_info( *pdf );
  112.     //
  113.     // the end vertices of the particles are not connected until
  114.     //  after the event is read --- we store the values in a map until then
  115.     TempParticleMap particle_to_end_vertex;
  116.     //
  117.     // read in the vertices
  118.     for ( int iii = 1; iii <= num_vertices; ++iii ) {
  119.         GenVertex* v = read_vertex(is,particle_to_end_vertex);
  120.         evt->add_vertex( v );
  121.     }
  122.     // set the signal process vertex
  123.     if ( signal_process_vertex ) {
  124.         evt->set_signal_process_vertex(
  125.             evt->barcode_to_vertex(signal_process_vertex) );
  126.     }
  127.     //
  128.     // last connect particles to their end vertices
  129.     GenParticle* beam1(0);
  130.     GenParticle* beam2(0);
  131.     for ( std::map<int,GenParticle*>::iterator pmap
  132.               = particle_to_end_vertex.order_begin();
  133.           pmap != particle_to_end_vertex.order_end(); ++pmap ) {
  134.         GenParticle* p =  pmap->second;
  135.         int vtx = particle_to_end_vertex.end_vertex( p );
  136.         GenVertex* itsDecayVtx = evt->barcode_to_vertex(vtx);
  137.         if ( itsDecayVtx ) itsDecayVtx->add_particle_in( p );
  138.         else {
  139.             std::cerr << "read_io_genevent_event: ERROR particle points"
  140.                       << "\n to null end vertex. " <<std::endl;
  141.         }
  142.         // also look for the beam particles
  143.         if( p->barcode() == bp1 ) beam1 = p;
  144.         if( p->barcode() == bp2 ) beam2 = p;
  145.     }
  146.     evt->set_beam_particles(beam1,beam2);
  147.     return true;
  148. }
  149.  
  150. HeavyIon* CommonIO::read_heavy_ion(std::istream* is)
  151. {
  152.     // assumes mode has already been checked
  153.     //
  154.     // test to be sure the next entry is of type "H" then ignore it
  155.     if ( !(*is) || is->peek()!='H' ) {
  156.         std::cerr << "IO_GenEvent::read_heavy_ion setting badbit." << std::endl;
  157.         is->clear(std::ios::badbit);
  158.         return false;
  159.     }
  160.     is->ignore();
  161.     // read values into temp variables, then create a new HeavyIon object
  162.     int nh =0, np =0, nt =0, nc =0,
  163.         neut = 0, prot = 0, nw =0, nwn =0, nwnw =0;
  164.     float impact = 0., plane = 0., xcen = 0., inel = 0.;
  165.     *is >> nh >> np >> nt >> nc >> neut >> prot
  166.            >> nw >> nwn >> nwnw >> impact >> plane >> xcen >> inel;
  167.     is->ignore(2,'\n');
  168.     if( nh == 0 ) return false;
  169.     HeavyIon* ion = new HeavyIon(nh, np, nt, nc, neut, prot,
  170.                                  nw, nwn, nwnw,
  171.                                  impact, plane, xcen, inel );
  172.     //
  173.     return ion;
  174. }
  175.  
  176. PdfInfo* CommonIO::read_pdf_info(std::istream* is)
  177. {
  178.     // assumes mode has already been checked
  179.     //
  180.     // test to be sure the next entry is of type "F" then ignore it
  181.     if ( !(*is) || is->peek() !='F') {
  182.         std::cerr << "IO_GenEvent::read_pdf_info setting badbit." << std::endl;
  183.         is->clear(std::ios::badbit);
  184.         return false;
  185.     }
  186.     is->ignore();
  187.     // read values into temp variables, then create a new PdfInfo object
  188.     int id1 =0, id2 =0;
  189.     double  x1 = 0., x2 = 0., scale = 0., pdf1 = 0., pdf2 = 0.;
  190.     *is >> id1 >> id2 >> x1 >> x2 >> scale >> pdf1 >> pdf2;
  191.     is->ignore(2,'\n');
  192.     if( id1 == 0 ) return false;
  193.     PdfInfo* pdf = new PdfInfo( id1, id2, x1, x2, scale, pdf1, pdf2);
  194.     //
  195.     return pdf;
  196. }
  197.  
  198. GenVertex* CommonIO::read_vertex( std::istream* is, TempParticleMap& particle_to_end_vertex )
  199. {
  200.     // assumes mode has already been checked
  201.     //
  202.     // test to be sure the next entry is of type "V" then ignore it
  203.     if ( !(*is) || is->peek()!='V' ) {
  204.         std::cerr << "IO_GenEvent::read_vertex setting badbit." << std::endl;
  205.         is->clear(std::ios::badbit);
  206.         return false;
  207.     }
  208.     is->ignore();
  209.     // read values into temp variables, then create a new GenVertex object
  210.     int identifier =0, id =0, num_orphans_in =0,
  211.         num_particles_out = 0, weights_size = 0;
  212.     double x = 0., y = 0., z = 0., t = 0.;
  213.     *is >> identifier >> id >> x >> y >> z >> t
  214.            >> num_orphans_in >> num_particles_out >> weights_size;
  215.     WeightContainer weights(weights_size);
  216.     for ( int i1 = 0; i1 < weights_size; ++i1 ) *is >> weights[i1];
  217.     is->ignore(2,'\n');
  218.     GenVertex* v = new GenVertex( FourVector(x,y,z,t),
  219.                             id, weights);
  220.     v->suggest_barcode( identifier );
  221.     //
  222.     // read and create the associated particles. outgoing particles are
  223.     //  added to their production vertices immediately, while incoming
  224.     //  particles are added to a map and handles later.
  225.     for ( int i2 = 1; i2 <= num_orphans_in; ++i2 ) {
  226.         read_particle(is,particle_to_end_vertex);
  227.     }
  228.     for ( int i3 = 1; i3 <= num_particles_out; ++i3 ) {
  229.         v->add_particle_out( read_particle(is,particle_to_end_vertex) );
  230.     }
  231.     return v;
  232. }
  233.  
  234. GenParticle* CommonIO::read_particle(std::istream* is,
  235.     TempParticleMap& particle_to_end_vertex ){
  236.     // assumes mode has already been checked
  237.     //
  238.     // test to be sure the next entry is of type "P" then ignore it
  239.     if ( !(*is) || is->peek()!='P' ) {
  240.         std::cerr << "IO_GenEvent::read_particle setting badbit."
  241.                   << std::endl;
  242.         is->clear(std::ios::badbit);
  243.         return false;
  244.     }
  245.     is->ignore();
  246.     //
  247.     // declare variables to be read in to, and read everything except flow
  248.     double px = 0., py = 0., pz = 0., e = 0., m = 0., theta = 0., phi = 0.;
  249.     int bar_code = 0, id = 0, status = 0, end_vtx_code = 0, flow_size = 0;
  250.     *is >> bar_code >> id >> px >> py >> pz >> e >> m >> status
  251.            >> theta >> phi >> end_vtx_code >> flow_size;
  252.     //
  253.     // read flow patterns if any exist
  254.     Flow flow;
  255.     int code_index, code;
  256.     for ( int i = 1; i <= flow_size; ++i ) {
  257.         *is >> code_index >> code;
  258.         flow.set_icode( code_index,code);
  259.     }
  260.     is->ignore(2,'\n'); // '\n' at end of entry
  261.     GenParticle* p = new GenParticle( FourVector(px,py,pz,e),
  262.                                 id, status, flow,
  263.                                 Polarization(theta,phi) );
  264.     p->set_generated_mass( m );
  265.     p->suggest_barcode( bar_code );
  266.     //
  267.     // all particles are connected to their end vertex separately
  268.     // after all particles and vertices have been created - so we keep
  269.     // a map of all particles that have end vertices
  270.     if ( end_vtx_code != 0 ) {
  271.         particle_to_end_vertex.addEndParticle(p,end_vtx_code);
  272.     }
  273.     return p;
  274. }
  275.  
  276. }       // end namespace HepMC
  277.  
  278.