io_ascii.hh

Go to the documentation of this file.

#ifndef IO_ASCII_INCLUDED
#define IO_ASCII_INCLUDED
 
#include <string>
#include <istream>
#include <ostream>
#include <sstream>
 
#include "Evt.hh"
#include "Hit.hh"
#include "Trk.hh"
#include "Head.hh"
#include "Exception.hh"
 
#include "TDatabasePDG.h"
 
namespace mc_keys
{
  const char* const  auto_t          =  "auto";
  const char* const  start_run_t     =  "start_run:";
  const char* const  start_event_t   =  "start_event:";
  const char* const  hit_t           =  "hit:";
  const char* const  hit_raw_t       =  "hit_raw:";
  const char* const  track_in_t      =  "track_in:";
  const char* const  track_t         =  "track:";
  const char* const  track_fit_t     =  "track_fit:";
  const char* const  neutrino_t      =  "neutrino:";
  const char* const  track_primary_t =  "track_primary:";
  const char* const  primarylepton_t =  "primarylepton:";
  const char* const  weights_t       =  "weights:";
  const char* const  w2list_t        =  "w2list:";
  const char* const  w3list_t        =  "w3list:";
  const char* const  hourangle_t     =  "hourangle:";
  const char* const  center_on_can_t =  "center_on_can:";
  const char* const  muon_decay_t    =  "muon_decay:";
  const char* const  end_event_t     =  "end_event:";
}
 
namespace mc_usr_keys
{
  // track-level quantities
  const char* const energy_lost_in_can = "energy_lost_in_can";
  const char* const bx    = "bx";    // these follow naming from 
  const char* const by    = "by";    // 'updated tag list' document.
  const char* const ichan = "ichan";
  const char* const cc    = "cc";
 
  // event-level corsika variables
  const char* const muon_decay_x = "muon_decay_x";   
  const char* const muon_decay_y = "muon_decay_y";   
  const char* const muon_decay_z = "muon_decay_z";   
 
  const char* const center_on_can_x = "center_on_can_x";
  const char* const center_on_can_y = "center_on_can_y";
  const char* const center_on_can_z = "center_on_can_z";
  const char* const hourangle = "hourangle";
}
 
 
namespace io_stringutil
{
  /**
   * Check if string starts with given text.
   *
   * \param  a             input string
   * \param  b             text
   * \return               true if string start with text; else false
   */
  inline bool startswith( const std::string& a, const std::string& b ) 
  {
    if ( a.find( b ) == 0 ) return true;
    return false;
  }
 
  /**
   * Remove leading and trailing white spaces.
   *
   * \param  s             input string
   * \return               trimmed string
   */
  inline std::string trim(const std::string& s)
  {
    using namespace std;
 
    if ( s=="" ) return s;
 
    string::size_type i1;
    string::size_type i2;
 
    for (i1=0; i1< s.length(); i1++)
      {
        if ( !isspace (s[i1]) ) break;
      }
    for (i2 = s.length()-1 ; i2 >i1 ; i2--)
      {
        if ( !isspace (s[i2]) ) break;
      }
    return s.substr( i1, i2-i1+1 );
  }
 
}
 
 
/**
 * Convert Geant3 to PDG particle type.
 *
 * \param  geant3_code   Geant3 code
 * \return               PDG code
 */
inline int pdg_code( int geant3_code )
{  
  if ( geant3_code==-1 ) return -1; // used for k40 hits
  if ( geant3_code < 0 ) return pdg_code( -geant3_code ); // used for scattered
 
  return TDatabasePDG::Instance()->ConvertGeant3ToPdg( geant3_code );  
}
 
/**
 * Convert PDG to Geant3 particle type.
 *
 * \param  pdg_code      PDG code
 * \return               Geant3 code
 */
inline int geant3_code( int pdg_code )
{
  if (pdg_code == -1 )   return -1;  
  if (pdg_code == +311 ) return geant3_code( 130 ); // K0 -> K0long
  if (pdg_code == -311 ) return geant3_code( 130 ); // K0bar -> K0long
 
  return TDatabasePDG::Instance()->ConvertPdgToGeant3( pdg_code );  
}
 
 
/**
 * Read a Vec(tor) from a stream.
 *
 * \param  v             vector
 * \param  is            input stream
 * \return               true if correctly read; else false
 */
inline bool read ( Vec& v, std::istream& is)
{
  is >> v.x >> v.y >> v.z;
  return !is.fail();
}
 
/**
 * Write a Vec(tor) to a stream.
 *
 * \param  v             vector
 * \param  os            output stream
 * \return               true if correctly written; else false
 */
inline bool write( const Vec& v, std::ostream& os)
{
  os << v.x << ' ' << v.y << ' ' << v.z;
  return !os.fail();
}
 
/**
 * Read a hit from a stream.
 *
 * \param  h             hit
 * \param  is            input stream
 * \param  read_mc       option to read also type and origin
 * \return               true if correctly read; else false
 */
inline bool read ( Hit& h, std::istream& is, bool read_mc = false  )
{
  h.dom_id     = 0; // need a proper det file to
  h.channel_id = 0; // set these.
 
  is >> h.id >> h.pmt_id >> h.a >> h.t;
  if ( !read_mc ) 
    {
      return !is.fail();
    }
  else
    {
      is >> h.type >> h.origin;
    }
 
  // If we get here without failing, we will return true, but
  // will attempt to read also pure_a and pure_t which may
  // or may not be present, depending on the km3 version.
 
  bool r = !is.fail(); 
  is >> h.pure_a >> h.pure_t; 
  if (is.fail())
    {
      h.pure_a = h.pure_t =0 ;
      is.clear();
    }
  return r;
}
 
 
/**
 * Write a hit to a stream.
 *
 * \param  h             hit
 * \param  os            output stream
 * \param  tag           tag 
 * \return               true if correctly written; else false
 */
inline bool write( const Hit& h, std::ostream& os, const std::string& tag = mc_keys::hit_t)
{
  int om_id = h.pmt_id; // todo: deal with this better.
 
  os << tag << ' ' << h.id << ' ' << om_id << ' ' << h.a << ' ' << h.t;
  if ( tag != mc_keys::hit_raw_t ) {
    os << ' ' << h.type  << ' ' << h.origin << ' ' << h.pure_a << ' ' << h.pure_t; 
  }
  os << std::endl;
  return !os.fail();
}
 
 
/**
 * Read track from a stream.
 *
 * \param  t             track
 * \param  is            input stream
 * \param  tag           tag 
 * \return               true if correctly read; else false
 */
inline bool read ( Trk& t, std::istream& is , const std::string& tag = mc_keys::track_in_t)
{
  using namespace std;
 
  string line;
  getline( is, line );
  istringstream ii(line);
 
  ii >> t.id >> t.pos >> t.dir >> t.E >> t.t;
 
  if (!ii) return false;  
 
  t.len = 0.0;
  t.clearusr();
  t.comment = tag;
 
  if ( tag == mc_keys::track_in_t) 
    {
      ii >> t.type; if (!ii) return false;
      t.type = pdg_code( t.type );
      ii >> t.len;  
      if (!ii) return true; // missing length is not an error
 
      double eloss=0;
      ii >> eloss;
      if (!ii) return true; // missing eloss is not an error
 
      //dbg(eloss);
      t.setusr( mc_usr_keys::energy_lost_in_can, eloss); 
    }
 
  if (tag == mc_keys::primarylepton_t)
    {      
      ii >> t.type;  if (!ii) return false;
      t.type = pdg_code( t.type ); // treat the same as track_in
    }
 
 
  if ( tag == mc_keys::neutrino_t ) 
    {
      double bx,by;
      int cc, ichan;
      ii >> bx >> by >> ichan >> t.type >> cc;
 
      t.clearusr();
      t.setusr(mc_usr_keys::bx, bx);
      t.setusr(mc_usr_keys::by ,by);
      t.setusr(mc_usr_keys::ichan, ichan);
      t.setusr(mc_usr_keys::cc, cc);
    }
 
  if ( tag == mc_keys::track_primary_t )
    {
      ii >> t.type; // nucleus id (in pdg format or not?)
    }
 
  if ( ii.fail() )
    {
      ostream& out = Exception::getOstream();
 
      out << "Error reading trk ";
 
      t.print(out);
 
      throw Exception(static_cast<ostringstream&>(out).str());
    }
 
  return !ii.fail();
}
 
/**
 * Write track to a stream.
 *
 * \param  t             track
 * \param  os            output stream
 * \param  tag           tag 
 * \return               true if correctly read; else false
 */
inline bool write( const Trk& t, std::ostream& os, std::string tag = mc_keys::track_in_t)
{
  using namespace std;
 
  int type = t.type;
 
  if ( tag == mc_keys::auto_t )
    {
      tag = t.comment;
    }
 
  os << tag << ' ' << t.id <<' '<< t.pos <<' '<< t.dir <<' '<< t.E << ' ' << t.t;
 
  if ( tag == mc_keys::track_in_t)
    {
      os << ' ' << geant3_code(type) << ' ' << t.len;
 
      if ( t.haveusr( mc_usr_keys::energy_lost_in_can ) )
        {
          os << ' ' << t.getusr( mc_usr_keys::energy_lost_in_can );
        }
      os << endl;
    }
 
  else if (tag == mc_keys::track_primary_t)
    {
      os << ' ' << t.type << endl;
    }
 
  else if (tag == mc_keys::primarylepton_t)
    {
      os << ' ' << geant3_code(type) << endl;
    }
 
  else if ( tag == mc_keys::neutrino_t )
    {
      os << ' ' << t.getusr(mc_usr_keys::bx)
         << ' ' << t.getusr(mc_usr_keys::by)
         << ' ' << t.getusr(mc_usr_keys::ichan)
         << ' ' << t.type
         << ' ' << t.getusr(mc_usr_keys::cc)  
         << endl;
    }
 
  else
    os << endl;
 
  return !os.fail();
}
 
 
/**
 * Read data.
 *
 * \param  is            input stream
 * \return               data
 */
inline std::vector<double> read_line_to_vector( std::istream& is )
{
  using namespace std;
 
  vector<double> r;
 
  string ss;
  getline(is,ss);
  istringstream il(ss);
  for ( double x; il >> x ; ) r.push_back( x );
 
  return r;
}
 
 
/**
 * Put value in front of data.
 *
 * \param  vec           data
 * \param  value         value
 */
template<typename T>
inline void push_front( std::vector<T>& vec, T& value )
{
  vec.insert( vec.begin(), value );
}
 
 
/**
 * Read event from a stream.
 *
 * \param  evt           event
 * \param  is            input stream
 * \param  skip_hits     option to skip reading of hits
 * \return               true if correctly read; else false
 */
inline bool read ( Evt& evt, std::istream& is, bool skip_hits = false )
{
  using namespace std;
 
  string w;
 
  // find next start_event
 
  while ( w != mc_keys::start_event_t && is.good() ) is >> w;
 
  int mc_event_type;   // dummy - this is always 1 in all files.
  is >> evt.mc_id >> mc_event_type;
 
  Trk trk_nu, trk_primary, trk_primary_lepton; 
  bool have_trk_nu(false), have_trk_primary(false), have_trk_primary_lepton(false); 
 
 
  evt.mc_trks.clear();
  evt.hits.clear();
  evt.mc_hits.clear();
 
  string w_old;
 
  while ( is )
    {
      static Hit h;
      static Trk t;
 
      is >> w;
 
      if (skip_hits && ( w == mc_keys::hit_t || w == mc_keys::hit_raw_t) )
        {
          is.ignore( 1000, '\n' );
          continue;
        }
 
      if ( w == mc_keys::hit_t ) 
        {         
          read( h, is, true ); 
          evt.mc_hits.push_back( h );
        }
      else if ( w == mc_keys::hit_raw_t )
        {
          read( h, is, false); 
          evt.hits.push_back( h );
        }
      else if ( w == mc_keys::track_in_t )
        {
          read( t, is , w);
          evt.mc_trks.push_back( t );
        }
      else if ( w == mc_keys::track_t )
        {
          read( t, is );
          evt.trks.push_back( t );
        }
      else if ( w == mc_keys::neutrino_t)
        {
          read( trk_nu ,is,  w);
          have_trk_nu = true;
        } 
      else if ( w == mc_keys::track_primary_t)
        {         
          read( trk_primary, is,  w);
          have_trk_primary = true;
        }
      else if ( w == mc_keys::primarylepton_t)
        {         
          read( trk_primary_lepton, is,  w);
          have_trk_primary_lepton = true;
        }      
      else if ( w == mc_keys::weights_t)
        {
          evt.w = read_line_to_vector( is );    
        }
      else if ( w == mc_keys::w2list_t) 
        {
          evt.w2list = read_line_to_vector( is );       
        }
      else if ( w == mc_keys::w3list_t) 
        {
          evt.w3list = read_line_to_vector( is );       
        }
      else if ( w == mc_keys::hourangle_t )
        {
          double ha;
          is >> ha;
          evt.setusr(mc_usr_keys::hourangle, ha );
        }
 
 
      else if ( w == mc_keys::center_on_can_t) 
        {
          // in corsika files, there is the (undocumented?) center_on_can tag,
          // which denoets the projection of the primary on the can. The direction
          // of the center_on_can 'track' is by defintion the direction of the
          // primary. We record the position in the usr data.
 
          vector<double> v = read_line_to_vector( is ); 
          if ( v.size() > 3 )
            {
              evt.setusr(mc_usr_keys::center_on_can_x, v[1] );
              evt.setusr(mc_usr_keys::center_on_can_y, v[2] );
              evt.setusr(mc_usr_keys::center_on_can_z, v[3] );
            }
        }
 
      else if ( w == mc_keys::muon_decay_t) 
        {
          // in km3sim files, there are additional tags, including this one
          vector<double> v = read_line_to_vector( is ); 
          if ( v.size() > 3 )
            {
              evt.setusr(mc_usr_keys::muon_decay_x, v[2] );
              evt.setusr(mc_usr_keys::muon_decay_y, v[3] );
              evt.setusr(mc_usr_keys::muon_decay_z, v[4] );
            }
        }
 
 
      else if ( w == mc_keys::end_event_t)
        {
          // finalize
 
          Trk none;
 
          if ( have_trk_primary_lepton )
            {
              Trk* p_primlep = 0 ;
 
              double min_edif = 1e100;
              for( auto& t: evt.mc_trks )
                {
                  if ( t.type != trk_primary_lepton.type ) continue;
 
                  // since we are reading ascii, we cannot demand an exact match
                  // in the floating point... so take the one that best matches
                  double edif = fabs( (t.E / trk_primary_lepton.E) -1 );
                  if (edif < min_edif )
                    {
                      min_edif = edif;
                      p_primlep = &t;
                    }
                }
 
              if ( p_primlep ) p_primlep -> comment =mc_keys:: primarylepton_t;
            }
 
 
          // If there is both a primary, and a neutrino, then the primary
          // will go second (mc_trks[1]) and the neutrino will be mc_trks[0].
          // The primarylepton tag is not stored as a seperate Trk.
 
          if ( have_trk_primary  )
            {
              trk_primary.id = 0;
              push_front( evt.mc_trks, trk_primary);
            }
 
          if ( have_trk_nu )
            {
              trk_nu.id = 0;
              push_front( evt.mc_trks, trk_nu);
            }
 
 
          return true;
        }
      else
        {
          //cout << "unknown tag " << w << endl;
          //cout << "previous = " << w_old << endl;
          is.ignore( 1000, '\n' );
        }
      w_old = w;
    }
 
  if (!is.eof()) 
    {
      THROW(Exception, "Error while reading ascii event" << w << ' ' << evt.id);
    }
 
  //else cout << "eof"<< endl;
  return false; 
}
 
 
/**
 * Write event to a stream.
 *
 * \param  evt           event
 * \param  os            output stream
 * \return               true if correctly read; else false
 */
inline bool write( const Evt& evt, std::ostream& os)
{
  using namespace std;
 
  os << mc_keys::start_event_t << ' ' << evt.mc_id << ' ' << 1 << endl;
 
  for( auto& trk: evt.mc_trks )  write ( trk, os, mc_keys::auto_t );
  for( auto& trk: evt.trks    )  write ( trk, os, mc_keys::track_fit_t );
  for( auto& hit: evt.mc_hits )  write ( hit, os, mc_keys::hit_t);
  for( auto& hit: evt.hits    )  write ( hit, os, mc_keys::hit_raw_t);
 
  os << mc_keys::weights_t; for (auto&  w: evt.w      ) os << ' ' << w; os << endl;
  os << mc_keys::w2list_t;  for (auto&  w: evt.w2list ) os << ' ' << w; os << endl;
  os << mc_keys::w3list_t;  for (auto&  w: evt.w3list ) os << ' ' << w; os << endl;
  os << mc_keys::end_event_t << endl;
 
  return true;
}
 
 
 
/**
 * Read header from a stream.
 *
 * The stream may be positioned at any point before the tag mc_keys::start_run_t,
 * which marks the beginning of the header.
 * Information before the header (if any) is disgarded.
 *
 * \param  hdr           header
 * \param  is            input stream
 * \return               true if correctly read; else false
 */
inline bool read( Head& hdr, std::istream& is )
{
  using namespace std;
 
  string line;
 
  bool start = false;
 
  while (getline( is, line ))
    {
      if ( io_stringutil::startswith(line, mc_keys::end_event_t))
        {
          break;
        }
 
      if ( io_stringutil::startswith(line, mc_keys::start_run_t))
        {
          start = true;
        }
 
      if ( io_stringutil::startswith(line, mc_keys::start_event_t))
        {
          // this should not happen:
          THROW(Exception, "Unexpected tag " << mc_keys::start_event_t << " found while reading header at " << line << " (could mean the evt file has no header)");
        }
 
      if (!start) continue; 
 
      vector<string> v = splitstring( line, ':' );
      if (v.size() < 2 )
        {
          THROW(Exception, "Invalid line found when reading header at " << line);
          continue;
        }
 
      // make sure key is unique
      string key = v[0];
      for (int i=1; hdr.find(key) != hdr.end() ; i++)
        {
          key = v[0]+"_"+to_string(i);
        }
 
      hdr[key] = io_stringutil::trim(v[1]);
    }
 
  if (!start)
    {
      THROW(Exception, "Reading of MC header terminated before finding a start_run: tag. Please check your file");
    }
 
  return start;
}
 
 
/**
 * Write header to a stream.
 *
 * \param  hdr           header
 * \param  os            output stream
 * \return               true if correctly written; else false
 */
inline bool write( const Head& hdr, std::ostream& os )
{
  hdr.print(os);
  return true;
}
 
#endif