DDG4/plugins/Geant4EventReaderHepMC.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 /** \addtogroup Geant4EventReader
0015  *
0016  @{
0017   \package Geant4EventReaderHepMC
0018  * \brief Plugin to read HepMC2 ASCII files
0019  *
0020  *
0021 @}
0022  */
0023
0024 // Framework include files
0025 #include <DDG4/IoStreams.h>
0026 #include <DDG4/Geant4InputAction.h>
0027
0028 // C/C++ include files
0029
0030 /// Namespace for the AIDA detector description toolkit
0031 namespace dd4hep {
0032
0033   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0034   namespace sim {
0035
0036     /// HepMC namespace declaration
0037     namespace HepMC {
0038       /// HepMC EventStream class used internally by the Geant4EventReaderHepMC plugin
0039       class EventStream;
0040     }
0041
0042     /// Class to populate Geant4 primaries from HepMC(2) files.
0043     /**
0044      *  Class to populate Geant4 primary particles and vertices from a
0045      *  file in HepMC format (ASCII)
0046      *
0047      *  For details also see:
0048      *  http://hepmc.web.cern.ch/hepmc/ReaderAsciiHepMC2_8cc_source.html
0049      *
0050      *  \author  P.Kostka (main author)
0051      *  \author  M.Frank  (code reshuffeling into new DDG4 scheme)
0052      *  \version 1.0
0053      *  \ingroup DD4HEP_SIMULATION
0054      */
0055     class Geant4EventReaderHepMC : public Geant4EventReader  {
0056       typedef boost::iostreams::stream<dd4hep_file_source<int> > in_stream;
0057       //typedef boost::iostreams::stream<dd4hep_file_source<TFile*> > in_stream;
0058       typedef HepMC::EventStream EventStream;
0059     protected:
0060       in_stream    m_input;
0061       EventStream* m_events;
0062     public:
0063       /// Initializing constructor
0064       explicit Geant4EventReaderHepMC(const std::string& nam);
0065       /// Default destructor
0066       virtual ~Geant4EventReaderHepMC();
0067       /// Read an event and fill a vector of MCParticles.
0068       virtual EventReaderStatus readParticles(int event_number,
0069                                               Vertices& vertices,
0070                                               std::vector<Particle*>& particles)  override;
0071       virtual EventReaderStatus moveToEvent(int event_number)  override;
0072       virtual EventReaderStatus skipEvent() override { return EVENT_READER_OK; }
0073
0074     };
0075   }     /* End namespace sim   */
0076 }       /* End namespace dd4hep       */
0077
0078 //====================================================================
0079 //  AIDA Detector description implementation
0080 //--------------------------------------------------------------------
0081 //
0082 //====================================================================
0083 // #include <DDG4/Geant4EventReaderHepMC"
0084
0085 // Framework include files
0086 #include <DDG4/Factories.h>
0087 #include <DD4hep/Printout.h>
0088 #include <DDG4/Geant4Primary.h>
0089 #include <CLHEP/Units/SystemOfUnits.h>
0090 #include <CLHEP/Units/PhysicalConstants.h>
0091
0092 // C/C++ include files
0093 #include <cerrno>
0094 #include <climits>
0095 #include <algorithm>
0096
0097 using namespace dd4hep::sim;
0098 using PropertyMask = dd4hep::detail::ReferenceBitMask<int>;
0099
0100 // Factory entry
0101 DECLARE_GEANT4_EVENT_READER(Geant4EventReaderHepMC)
0102
0103 /// Namespace for the AIDA detector description toolkit
0104 namespace dd4hep {
0105
0106   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0107   namespace sim {
0108
0109     /// HepMC namespace declaration
0110     namespace HepMC {
0111
0112       /// HepMC EventHeader class used internally by the Geant4EventReaderHepMC plugin
0113       /*
0114        *  \author  P.Kostka (main author)
0115        *  \author  M.Frank  (code reshuffeling into new DDG4 scheme)
0116        *  \version 1.0
0117        *  \ingroup DD4HEP_SIMULATION
0118        */
0119       class EventHeader  {
0120       public:
0121         int   id;
0122         int   num_vertices;
0123         int   bp1, bp2;
0124         int   signal_process_id;
0125         int   signal_process_vertex;
0126         float scale;
0127         float alpha_qcd;
0128         float alpha_qed;
0129         std::vector<float>      weights;
0130         std::vector<long>       random;
0131         /// Default constructor
0132         EventHeader() : id(0), num_vertices(0), bp1(0), bp2(0),
0133                         signal_process_id(0), signal_process_vertex(0),
0134                         scale(0.0), alpha_qcd(0.0), alpha_qed(0.0), weights(), random() {}
0135       };
0136
0137       /// The known_io enum is used to track which type of input is being read
0138       enum known_io { gen=1, ascii, extascii, ascii_pdt, extascii_pdt };
0139
0140       /// HepMC EventStream class used internally by the Geant4EventReaderHepMC plugin
0141       /*
0142        *  \author  P.Kostka (main author)
0143        *  \author  M.Frank  (code reshuffeling into new DDG4 scheme)
0144        *  \version 1.0
0145        *  \ingroup DD4HEP_SIMULATION
0146        */
0147       class EventStream {
0148       public:
0149         typedef std::map<int,Geant4Vertex*> Vertices;
0150         typedef std::map<int,Geant4Particle*> Particles;
0151
0152         std::istream& instream;
0153
0154         // io information
0155         std::string key;
0156         double mom_unit, pos_unit;
0157         int    io_type;
0158
0159         float  xsection, xsection_err;
0160         EventHeader header;
0161         Vertices m_vertices;
0162         Particles m_particles;
0163
0164         /// Default constructor
0165         EventStream(std::istream& in) : instream(in), mom_unit(0.0), pos_unit(0.0),
0166                                    io_type(0), xsection(0.0), xsection_err(0.0)
0167         { use_default_units();                       }
0168         /// Check if data stream is in proper state and has data
0169         bool ok()  const;
0170         Geant4Vertex* vertex(int i);
0171         Particles& particles() { return m_particles; }
0172         Vertices&  vertices()  { return m_vertices;  }
0173         void set_io(int typ, const std::string& k)
0174         { io_type = typ;    key = k;                 }
0175         void use_default_units()
0176         { mom_unit = CLHEP::MeV;   pos_unit = CLHEP::mm;           }
0177         bool read();
0178         void clear();
0179       };
0180
0181       char get_input(std::istream& is, std::istringstream& iline);
0182       int read_until_event_end(std::istream & is);
0183       int read_weight_names(EventStream &, std::istringstream& iline);
0184       int read_particle(EventStream &info, std::istringstream& iline, Geant4Particle * p);
0185       int read_vertex(EventStream &info, std::istream& is, std::istringstream & iline);
0186       int read_event_header(EventStream &info, std::istringstream & input, EventHeader& header);
0187       int read_cross_section(EventStream &info, std::istringstream & input);
0188       int read_units(EventStream &info, std::istringstream & input);
0189       int read_heavy_ion(EventStream &, std::istringstream & input);
0190       int read_pdf(EventStream &, std::istringstream & input);
0191       Geant4Vertex* vertex(EventStream& info, int i);
0192       void fix_particles(EventStream &info);
0193     }
0194   }
0195 }
0196
0197 // For debugging specify here vertex numbers
0198 //#define DD4HEP_DEBUG_HEP_MC_VERTEX   -390
0199 //#define DD4HEP_DEBUG_HEP_MC_PARTICLE 418
0200
0201 /// Initializing constructor
0202 Geant4EventReaderHepMC::Geant4EventReaderHepMC(const std::string& nam)
0203   : Geant4EventReader(nam), m_input(), m_events(0)
0204 {
0205   // Now open the input file:
0206   m_input.open(nam.c_str(), BOOST_IOS::in|BOOST_IOS::binary);
0207   if ( not m_input.is_open() )   {
0208     except("+++ Failed to open input stream: %s Error:%s.", nam.c_str(), ::strerror(errno));
0209   }
0210   m_events = new HepMC::EventStream(m_input);
0211 }
0212
0213 /// Default destructor
0214 Geant4EventReaderHepMC::~Geant4EventReaderHepMC()    {
0215   delete m_events;
0216   m_events = 0;
0217   m_input.close();
0218 }
0219
0220 /// skipEvents if required
0221 Geant4EventReader::EventReaderStatus
0222 Geant4EventReaderHepMC::moveToEvent(int event_number) {
0223   if( m_currEvent < event_number && event_number != 0 ) {
0224     printout(INFO,"EventReaderHepMC::moveToEvent","Current event:%d Skipping the next %d events",
0225              m_currEvent, event_number);
0226     while ( m_currEvent < event_number ) {
0227       if ( not m_events->read() ) return EVENT_READER_ERROR;
0228       ++m_currEvent;
0229     }
0230   }
0231   printout(DEBUG,"EventReaderHepMC::moveToEvent","Current event number: %d",m_currEvent);
0232   return EVENT_READER_OK;
0233 }
0234
0235 /// Read an event and fill a vector of MCParticles.
0236 Geant4EventReaderHepMC::EventReaderStatus
0237 Geant4EventReaderHepMC::readParticles(int /* ev_id */,
0238                                       Vertices&  vertices,
0239                                       Particles& output) {
0240
0241   //fg: for now we create exactly one event vertex here ( as before )
0242   //    this needs revisiting as HepMC allows to have more than one vertex ...
0243   Geant4Vertex* primary_vertex = new Geant4Vertex ;
0244   vertices.emplace_back( primary_vertex );
0245   primary_vertex->x = 0;
0246   primary_vertex->y = 0;
0247   primary_vertex->z = 0;
0248
0249   if ( !m_events->ok() )  {
0250     vertices.clear();
0251     output.clear();
0252     return EVENT_READER_EOF;
0253   }
0254   else if ( m_events->read() )  {
0255     EventStream::Particles& parts = m_events->particles();
0256
0257     Position pos(primary_vertex->x,primary_vertex->y,primary_vertex->z);
0258
0259     output.reserve(parts.size());
0260     transform(parts.begin(),parts.end(),back_inserter(output),detail::reference2nd(parts));
0261     m_events->clear();
0262     if (pos.mag2() > std::numeric_limits<double>::epsilon() )  {
0263       for(Particles::iterator k=output.begin(); k != output.end(); ++k) {
0264         Geant4ParticleHandle p(*k);
0265         p->vsx += pos.x();
0266         p->vsy += pos.y();
0267         p->vsz += pos.z();
0268         p->vex += pos.x();
0269         p->vey += pos.y();
0270         p->vez += pos.z();
0271       }
0272     }
0273     for(Particles::const_iterator k=output.begin(); k != output.end(); ++k) {
0274       Geant4ParticleHandle p(*k);
0275       printout(VERBOSE,m_name,
0276                "+++ %s ID:%3d status:%08X typ:%9d Mom:(%+.2e,%+.2e,%+.2e)[MeV] "
0277                "time: %+.2e [ns] #Dau:%3d #Par:%1d",
0278                "",p->id,p->status,p->pdgID,
0279                p->psx/CLHEP::MeV,p->psy/CLHEP::MeV,p->psz/CLHEP::MeV,p->time/CLHEP::ns,
0280                p->daughters.size(),
0281                p->parents.size());
0282       //output.emplace_back(p);
0283
0284       //add particles to the 'primary vertex'
0285       if ( p->parents.size() == 0 )  {
0286         PropertyMask status(p->status);
0287         if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
0288           primary_vertex->in.insert(p->id);  // Beam particles and primary quarks etc.
0289         else
0290           primary_vertex->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
0291       }
0292     }
0293     ++m_currEvent;
0294     return EVENT_READER_OK;
0295   }
0296   vertices.clear();
0297   output.clear();
0298   return EVENT_READER_EOF;
0299 }
0300
0301 void HepMC::fix_particles(EventStream& info)  {
0302   EventStream::Particles& parts = info.particles();
0303   EventStream::Vertices&  verts = info.vertices();
0304   EventStream::Particles::iterator i;
0305   std::set<int>::const_iterator id, ip;
0306   for(i=parts.begin(); i != parts.end(); ++i)  {
0307     Geant4ParticleHandle p((*i).second);
0308     int end_vtx_id = p->secondaries;
0309     p->secondaries = 0;
0310     Geant4Vertex* v = vertex(info,end_vtx_id);
0311 #if defined(DD4HEP_DEBUG_HEP_MC_VERTEX)
0312     if ( end_vtx_id == DD4HEP_DEBUG_HEP_MC_VERTEX )   {
0313       std::cout << "End-vertex:" << end_vtx_id << std::endl;
0314     }
0315 #endif
0316     if ( v )   {
0317       p->vex = v->x;
0318       p->vey = v->y;
0319       p->vez = v->z;
0320       v->in.insert(p->id);
0321       for(id=v->out.begin(); id!=v->out.end();++id)    {
0322         EventStream::Particles::iterator ipp = parts.find(*id);
0323         Geant4Particle* dau = ipp != parts.end() ? (*ipp).second : 0;
0324         if ( !dau )
0325           std::cout << "ERROR: Invalid daughter particle: " << *id << std::endl;
0326         else
0327           dau->parents.insert(p->id);
0328         p->daughters.insert(*id);
0329       }
0330     }
0331   }
0332   for(const auto& iv : verts)   {
0333     Geant4Vertex* v = iv.second;
0334     for (int pout : v->out)   {
0335       EventStream::Particles::iterator ipp = parts.find(pout);
0336       if ( ipp != parts.end() )  {
0337         Geant4Particle* p = (*ipp).second;
0338         for (int d : v->in)   {
0339           p->parents.insert(d);
0340         }
0341       }
0342     }
0343   }
0344   /// Particles originating from the beam (=no parents) must be
0345   /// be stripped off their parents and the status set to G4PARTICLE_GEN_DECAYED!
0346   std::vector<Geant4Particle*> beam;
0347   for(const auto& ipart : parts)   {
0348     Geant4ParticleHandle p(ipart.second);
0349     if ( p->parents.size() == 0 )  {
0350       for(int d : p->daughters)   {
0351         Geant4Particle *pp = parts[d];
0352         beam.emplace_back(pp);
0353       }
0354     }
0355   }
0356   for(auto* ipp : beam)   {
0357     //std::cout << "Clear parents of " << (*ipp)->id << std::endl;
0358     ipp->parents.clear();
0359     ipp->status = G4PARTICLE_GEN_DECAYED;
0360   }
0361 }
0362
0363 Geant4Vertex* HepMC::vertex(EventStream& info, int i)   {
0364   EventStream::Vertices::iterator it=info.vertices().find(i);
0365   return (it==info.vertices().end()) ? 0 : (*it).second;
0366 }
0367
0368 char HepMC::get_input(std::istream& is, std::istringstream& iline)  {
0369   char value = is.peek();
0370   if ( !is ) {        // make sure the stream is valid
0371     std::cerr << "StreamHelpers: setting badbit." << std::endl;
0372     is.clear(std::ios::badbit);
0373     return -1;
0374   }
0375   std::string line, firstc;
0376   getline(is,line);
0377   if ( !is ) {        // make sure the stream is valid
0378     std::cerr << "StreamHelpers: setting badbit." << std::endl;
0379     is.clear(std::ios::badbit);
0380     return -1;
0381   }
0382   iline.clear();
0383   iline.str(line);
0384   if ( line.length()>0 && line[1] == ' ' ) iline >> firstc;
0385   return iline ? value : -1;
0386 }
0387
0388 int HepMC::read_until_event_end(std::istream & is) {
0389   std::string line;
0390   while ( is ) {
0391     char val = is.peek();
0392     if( val == 'E' ) {  // next event
0393       return 1;
0394     }
0395     getline(is,line);
0396     if ( !is.good() ) return 0;
0397   }
0398   return 0;
0399 }
0400
0401 int HepMC::read_weight_names(EventStream&, std::istringstream&)   {
0402 #if 0
0403   int HepMC::read_weight_names(EventStream& info, std::istringstream& iline)
0404     size_t name_size = 0;
0405   iline >> name_size;
0406   info.weights.names.clear();
0407   info.weights.weights.clear();
0408
0409   std::string name;
0410   WeightContainer namedWeight;
0411   std::string::size_type i1 = line.find("\""), i2, len = line.size();
0412   for(size_t ii = 0; ii < name_size; ++ii) {
0413     // weight names may contain blanks
0414     if(i1 >= len) {
0415       std::cout << "debug: attempting to read past the end of the named weight line " << std::endl;
0416       std::cout << "debug: We should never get here" << std::endl;
0417       std::cout << "debug: Looking for the end of this event" << std::endl;
0418       read_until_event_end(is);
0419     }
0420     i2 = line.find("\"",i1+1);
0421     name = line.substr(i1+1,i2-i1-1);
0422     if ( namedWeight.names.find(name) == namedWeight.names.end() )
0423       namedWeight.names[name] = namedWeight.names.size();
0424     namedWeight.weights[ii] = info.weights.weights[ii];
0425     i1 = line.find("\"",i2+1);
0426   }
0427
0428   info.weights.weights = namedWeight.weights;
0429   info.weights.names   = namedWeight.names;
0430 #endif
0431   return 1;
0432 }
0433
0434 int HepMC::read_particle(EventStream &info, std::istringstream& input, Geant4Particle * p)   {
0435   float ene = 0., theta = 0., phi = 0;
0436   int   size = 0, stat=0;
0437   PropertyMask status(p->status);
0438
0439   // check that the input stream is still OK after reading item
0440   input >> p->id >> p->pdgID >> p->psx >> p->psy >> p->psz >> ene;
0441   p->id = info.particles().size();
0442 #if defined(DD4HEP_DEBUG_HEP_MC_PARTICLE)
0443   if ( p->id == DD4HEP_DEBUG_HEP_MC_PARTICLE )   {
0444     std::cout << "Particle id: " << p->id << std::endl;
0445   }
0446 #endif
0447   p->charge = 0;
0448   p->psx *= info.mom_unit;
0449   p->psy *= info.mom_unit;
0450   p->psz *= info.mom_unit;
0451   ene *= info.mom_unit;
0452   if ( !input )
0453     return 0;
0454   else if ( info.io_type != ascii )  {
0455     input >> p->mass;
0456     p->mass *= info.mom_unit;
0457   }
0458   else   {
0459     p->mass = std::sqrt(fabs(ene*ene - (p->psx*p->psx + p->psy*p->psy + p->psz*p->psz)));
0460   }
0461   // Reuse here the secondaries to store the end-vertex ID
0462   input >> stat >> theta >> phi >> p->secondaries >> size;
0463   if( !input )   {
0464     return 0;
0465   }
0466   //
0467   //  Generator status
0468   //  Simulator status 0 until simulator acts on it
0469   status.clear();
0470   if ( stat == 0 )        status.set(G4PARTICLE_GEN_EMPTY);
0471   else if ( stat == 0x1 ) status.set(G4PARTICLE_GEN_STABLE);
0472   else if ( stat == 0x2 ) status.set(G4PARTICLE_GEN_DECAYED);
0473   else if ( stat == 0x3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
0474   else if ( stat == 0x4 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
0475   else if ( stat == 0xB ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
0476   else                    status.set(G4PARTICLE_GEN_OTHER);
0477   /// If there is an end vertex, the particle already decayed
0478   if ( p->secondaries != 0 )  {
0479     status.set(G4PARTICLE_GEN_DECAYED);
0480   }
0481   /// Keep a copy of the full generator status
0482   p->genStatus = stat&G4PARTICLE_GEN_STATUS_MASK;
0483
0484   // read flow patterns if any exist. Protect against tainted readings.
0485   size = std::min(size,100);
0486   for (int i = 0; i < size; ++i ) {
0487     input >> p->colorFlow[0] >> p->colorFlow[1];
0488     if(!input) return 0;
0489   }
0490   return 1;
0491 }
0492
0493 int HepMC::read_vertex(EventStream &info, std::istream& is, std::istringstream & input)    {
0494   int id=0, dummy = 0, num_orphans_in=0, num_particles_out=0, weights_size=0;
0495   std::vector<float> weights;
0496   Geant4Vertex* v = new Geant4Vertex();
0497   Geant4Particle* p;
0498
0499   if( !input ) {
0500     delete v;
0501     return 0;
0502   }
0503   input >> id >> dummy >> v->x >> v->y >> v->z >> v->time
0504         >> num_orphans_in >> num_particles_out >> weights_size;
0505   if( !input ) {
0506     delete v;
0507     return 0;
0508   }
0509   if ( weights_size < 0 || weights_size > USHRT_MAX )  {
0510     delete v;
0511     return 0;
0512   }
0513 #if defined(DD4HEP_DEBUG_HEP_MC_VERTEX)
0514   if ( id == DD4HEP_DEBUG_HEP_MC_VERTEX )   {
0515     printout(ALWAYS,"HepMC","++ Created Vertex ID=%d",id);
0516   }
0517 #endif
0518   v->x *= info.pos_unit;
0519   v->y *= info.pos_unit;
0520   v->z *= info.pos_unit;
0521   for (int i1 = 0; i1 < weights_size; ++i1) {
0522     float value = 0e0;
0523     input >> value;
0524     weights.emplace_back(value);
0525     if( !input ) {
0526       delete v;
0527       return 0;
0528     }
0529   }
0530   info.vertices().emplace(id,v);
0531   for(char value = is.peek(); value=='P'; value=is.peek())  {
0532     value = get_input(is,input);
0533     if( !input || value < 0 )
0534       return 0;
0535
0536     read_particle(info, input,p = new Geant4Particle());
0537     if( !input )   {
0538       printout(ERROR,"HepMC","++ Vertex %d Failed to daughter read particle!",id);
0539       delete p;
0540       return 0;
0541     }
0542     info.particles().emplace(p->id,p);
0543     p->pex = p->psx;
0544     p->pey = p->psy;
0545     p->pez = p->psz;
0546     if ( --num_orphans_in >= 0 )   {
0547       v->in.insert(p->id);
0548       p->vex = v->x;
0549       p->vey = v->y;
0550       p->vez = v->z;
0551 #if defined(DD4HEP_DEBUG_HEP_MC_VERTEX)
0552       if ( id == DD4HEP_DEBUG_HEP_MC_VERTEX )   {
0553         printout(ALWAYS,"HepMC","++ Vertex %d Add INGOING  Particle: %d",id,p->id);
0554       }
0555 #endif
0556     }
0557     else if ( num_particles_out >= 0 )   {
0558       v->out.insert(p->id);
0559       p->vsx = v->x;
0560       p->vsy = v->y;
0561       p->vsz = v->z;
0562 #if defined(DD4HEP_DEBUG_HEP_MC_VERTEX)
0563       if ( id == DD4HEP_DEBUG_HEP_MC_VERTEX )   {
0564         printout(ALWAYS,"HepMC","++ Vertex %d Add OUTGOING Particle: %d",id,p->id);
0565       }
0566 #endif
0567     }
0568     else  {
0569       delete p;
0570       except("HepMC", "Invalid number of particles....");
0571     }
0572   }
0573   return 1;
0574 }
0575
0576 int HepMC::read_event_header(EventStream &info, std::istringstream & input, EventHeader& header)   {
0577   // read values into temp variables, then fill GenEvent
0578   int size = 0;
0579   input >> header.id;
0580   if( info.io_type == gen || info.io_type == extascii ) {
0581     int nmpi = -1;
0582     input >> nmpi;
0583     if( input.fail() ) return 0;
0584     //MSF set_mpi( nmpi );
0585   }
0586   input >> header.scale;
0587   input >> header.alpha_qcd;
0588   input >> header.alpha_qed;
0589   input >> header.signal_process_id;
0590   input >> header.signal_process_vertex;
0591   input >> header.num_vertices;
0592   if( info.io_type == gen || info.io_type == extascii )
0593     input >> header.bp1 >> header.bp2;
0594
0595   printout(DEBUG,"HepMC","++ Event header: %s",input.str().c_str());
0596   input >> size;
0597   input.clear();
0598   if( input.fail() )
0599     return 0;
0600   if( size < 0 || size > USHRT_MAX )
0601     return 0;
0602
0603   for(int i = 0; i < size; ++i )  {
0604     long val = 0e0;
0605     input >> val;
0606     header.random.emplace_back(val);
0607     if( input.fail() ) return 0;
0608   }
0609
0610   input >> size;
0611   if( input.fail() )
0612     return 0;
0613   if( size < 0 || size > USHRT_MAX )
0614     return 0;
0615
0616   std::vector<float> wgt;
0617   for( int ii = 0; ii < size; ++ii )  {
0618     float val = 0e0;
0619     input >> val;
0620     wgt.emplace_back(val);
0621     if( input.fail() ) return 0;
0622   }
0623   // weight names will be added later if they exist
0624   if( !wgt.empty() ) header.weights = std::move(wgt);
0625   return 1;
0626 }
0627
0628 int HepMC::read_cross_section(EventStream &info, std::istringstream & input)   {
0629   input >> info.xsection >> info.xsection_err;
0630   return input.fail() ? 0 : 1;
0631 }
0632
0633 int HepMC::read_units(EventStream &info, std::istringstream & input)   {
0634   if( info.io_type == gen )  {
0635     std::string mom, pos;
0636     input >> mom >> pos;
0637     if ( !input.fail() )  {
0638       if ( mom == "KEV" ) info.mom_unit = CLHEP::keV;
0639       else if ( mom == "MEV" ) info.mom_unit = CLHEP::MeV;
0640       else if ( mom == "GEV" ) info.mom_unit = CLHEP::GeV;
0641       else if ( mom == "TEV" ) info.mom_unit = CLHEP::TeV;
0642
0643       if ( pos == "MM" ) info.pos_unit = CLHEP::mm;
0644       else if ( pos == "CM" ) info.pos_unit = CLHEP::cm;
0645       else if ( pos == "M"  ) info.pos_unit = CLHEP::m;
0646     }
0647   }
0648   return input.fail() ? 0 : 1;
0649 }
0650
0651 int HepMC::read_heavy_ion(EventStream &, std::istringstream & input)  {
0652   // read values into temp variables, then create a new HeavyIon object
0653   int nh =0, np =0, nt =0, nc =0,
0654     neut = 0, prot = 0, nw =0, nwn =0, nwnw =0;
0655   float impact = 0., plane = 0., xcen = 0., inel = 0.;
0656   input >> nh >> np >> nt >> nc >> neut >> prot >> nw >> nwn >> nwnw;
0657   input >> impact >> plane >> xcen >> inel;
0658   /*
0659     std::cerr << "Reading heavy ion, but igoring data!" << std::endl;
0660     ion->set_Ncoll_hard(nh);
0661     ion->set_Npart_proj(np);
0662     ion->set_Npart_targ(nt);
0663     ion->set_Ncoll(nc);
0664     ion->set_spectator_neutrons(neut);
0665     ion->set_spectator_protons(prot);
0666     ion->set_N_Nwounded_collisions(nw);
0667     ion->set_Nwounded_N_collisions(nwn);
0668     ion->set_Nwounded_Nwounded_collisions(nwnw);
0669     ion->set_impact_parameter(impact);
0670     ion->set_event_plane_angle(plane);
0671     ion->set_eccentricity(xcen);
0672     ion->set_sigma_inel_NN(inel);
0673   */
0674   return input.fail() ? 0 : 1;
0675 }
0676
0677 int HepMC::read_pdf(EventStream &, std::istringstream & input)  {
0678   // read values into temp variables, then create a new PdfInfo object
0679   int id1 =0, id2 =0;
0680   double  x1 = 0., x2 = 0., scale = 0., pdf1 = 0., pdf2 = 0.;
0681   input >> id1 ;
0682   if ( input.fail()  )
0683     return 0;
0684   // check now for empty PdfInfo line
0685   if( id1 == 0 )
0686     return 0;
0687   // continue reading
0688   input >> id2 >> x1 >> x2 >> scale >> pdf1 >> pdf2;
0689   if ( input.fail()  )
0690     return 0;
0691   /*
0692     std::cerr << "Reading pdf, but igoring data!" << std::endl;
0693     pdf->set_id1( id1 );
0694     pdf->set_id2( id2 );
0695     pdf->set_x1( x1 );
0696     pdf->set_x2( x2 );
0697     pdf->set_scalePDF( scale );
0698     pdf->set_pdf1( pdf1 );
0699     pdf->set_pdf2( pdf2 );
0700   */
0701   // check to see if we are at the end of the line
0702   if( !input.eof() )  {
0703     int pdf_id1=0, pdf_id2=0;
0704     input >> pdf_id1 >> pdf_id2;
0705     /*
0706       pdf->set_pdf_id1( pdf_id1 );
0707       pdf->set_pdf_id2( pdf_id2 );
0708     */
0709   }
0710   return input.fail() ? 0 : 1;
0711 }
0712
0713 /// Check if data stream is in proper state and has data
0714 bool HepMC::EventStream::ok()  const   {
0715   // make sure the stream is good
0716   if ( instream.eof() || instream.fail() )  {
0717     instream.clear(std::ios::badbit);
0718     return false;
0719   }
0720   return true;
0721 }
0722
0723 void HepMC::EventStream::clear()   {
0724   detail::releaseObjects(m_vertices);
0725   detail::releaseObjects(m_particles);
0726 }
0727
0728 bool HepMC::EventStream::read()   {
0729   EventStream& info = *this;
0730   bool event_read = false;
0731
0732   detail::releaseObjects(vertices());
0733   detail::releaseObjects(particles());
0734
0735   while( instream.good() ) {
0736     char value = instream.peek();
0737     std::istringstream input_line;
0738     if      ( value == 'E' && event_read )
0739       break;
0740     else if ( instream.eof() && event_read )
0741       goto Done;
0742     else if ( instream.eof() )
0743       return false;
0744     else if ( value=='#' || ::isspace(value) )  {
0745       get_input(instream,input_line);
0746       continue;
0747     }
0748     value = get_input(instream,input_line);
0749
0750     // On failure switch to end
0751     if( !input_line || value < 0 )
0752       goto Skip;
0753
0754     switch( value )   {
0755     case 'H':  {
0756       int iotype = 0;
0757       std::string key_value;
0758       input_line >> key_value;
0759       // search for event listing key before first event only.
0760       key_value = key_value.substr(0,key_value.find('\r'));
0761       if ( key_value == "H" && (this->io_type == gen || this->io_type == extascii) ) {
0762         read_heavy_ion(info, input_line);
0763         break;
0764       }
0765       else if( key_value == "HepMC::IO_GenEvent-START_EVENT_LISTING" )
0766         this->set_io(gen,key_value);
0767       else if( key_value == "HepMC::IO_Ascii-START_EVENT_LISTING" )
0768         this->set_io(ascii,key_value);
0769       else if( key_value == "HepMC::IO_ExtendedAscii-START_EVENT_LISTING" )
0770         this->set_io(extascii,key_value);
0771       else if( key_value == "HepMC::IO_Ascii-START_PARTICLE_DATA" )
0772         this->set_io(ascii_pdt,key_value);
0773       else if( key_value == "HepMC::IO_ExtendedAscii-START_PARTICLE_DATA" )
0774         this->set_io(extascii_pdt,key_value);
0775       else if( key_value == "HepMC::IO_GenEvent-END_EVENT_LISTING" )
0776         iotype = gen;
0777       else if( key_value == "HepMC::IO_Ascii-END_EVENT_LISTING" )
0778         iotype = ascii;
0779       else if( key_value == "HepMC::IO_ExtendedAscii-END_EVENT_LISTING" )
0780         iotype = extascii;
0781       else if( key_value == "HepMC::IO_Ascii-END_PARTICLE_DATA" )
0782         iotype = ascii_pdt;
0783       else if( key_value == "HepMC::IO_ExtendedAscii-END_PARTICLE_DATA" )
0784         iotype = extascii_pdt;
0785
0786       if( iotype != 0 && this->io_type != iotype )  {
0787         std::cerr << "GenEvent::find_end_key: iotype keys have changed. "
0788                   << "MALFORMED INPUT" << std::endl;
0789         instream.clear(std::ios::badbit);
0790         return false;
0791       }
0792       else if ( iotype != 0 )  {
0793         break;
0794       }
0795       continue;
0796     }
0797     case 'E':           // deal with the event line
0798       if ( !read_event_header(info, input_line, this->header) )
0799         goto Skip;
0800       event_read = true;
0801       continue;
0802
0803     case 'N':           // get weight names
0804       if ( !read_weight_names(info, input_line) )
0805         goto Skip;
0806       continue;
0807
0808     case 'U':           // get unit information if it exists
0809       if ( !read_units(info, input_line) )
0810         goto Skip;
0811       continue;
0812
0813     case 'C':           // we have a GenCrossSection line
0814       if ( !read_cross_section(info, input_line) )
0815         goto Skip;
0816       continue;
0817
0818     case 'V':           // Read vertex with particles
0819       if ( !read_vertex(info, instream, input_line) )
0820         goto Skip;
0821       continue;
0822
0823     case 'F':           // Read PDF
0824       if ( !read_pdf(info, input_line) )
0825         goto Skip;
0826       continue;
0827
0828     case 'P':           // we should not find this line
0829       std::cerr << "streaming input: found unexpected Particle line." << std::endl;
0830       continue;
0831
0832     default:            // ignore everything else
0833       continue;
0834     }
0835     continue;
0836   Skip:
0837     printout(WARNING,"HepMC::EventStream","+++ Skip event with ID: %d",this->header.id);
0838     detail::releaseObjects(vertices());
0839     detail::releaseObjects(particles());
0840     read_until_event_end(instream);
0841     event_read = false;
0842     if ( instream.eof() ) return false;
0843   }
0844
0845   if( not instream.good() ) return false;
0846  Done:
0847   fix_particles(info);
0848   detail::releaseObjects(vertices());
0849   return true;
0850 }
0851