JKM3Sim.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <fstream>
#include <vector>
#include <map>
 
#include "JLang/JException.hh"
#include "JMath/JConstants.hh"
#include "JPhysics/JConstants.hh"
#include "JPhysics/KM3NeT.hh"
#include "JPhysics/JDispersion.hh"
 
#include "Jeep/JStreamToolkit.hh"
#include "Jeep/JParser.hh"
 
namespace {
 
  using JLANG::JParseError;
 
  struct JKM3Sim {
 
    static constexpr const char* const detectorDepth_t     =  "detectorDepth";
    static constexpr const char* const PPCKOV_t            =  "PPCKOV";
    static constexpr const char* const RINDEX_WATER_t      =  "RINDEX_WATER";
    static constexpr const char* const ABSORPTION_WATER_t  =  "ABSORPTION_WATER"; 
    static constexpr const char* const RINDEX_GLASS_t      =  "RINDEX_GLASS";
    static constexpr const char* const ABSORPTION_GLASS_t  =  "ABSORPTION_GLASS";
    static constexpr const char* const RINDEX_GELL_t       =  "RINDEX_GELL";
    static constexpr const char* const RINDEX_AIR_t        =  "RINDEX_AIR";
    static constexpr const char* const ABSORPTION_AIR_t    =  "ABSORPTION_AIR";
    static constexpr const char* const RINDEX_CATH_t       =  "RINDEX_CATH";
    static constexpr const char* const ABSORPTION_CATH_t   =  "ABSORPTION_CATH";
    static constexpr const char* const Q_EFF_t             =  "Q_EFF";
    static constexpr const char* const MIE_WATER_t         =  "MIE_WATER";
    static constexpr const char* const MIE_PARAM_t         =  "MIE_PARAM";
    static constexpr const char* const COS_ANGLES_t        =  "COS_ANGLES";
    static constexpr const char* const ANG_ACCEPT_t        =  "ANG_ACCEPT";
 
    static constexpr const char* const END_t               =  "END";
    static constexpr char              SEPARATOR           =  '*';
 
 
    /**
     * Default constructor.
     */
    JKM3Sim()
    {
      using namespace JPP;
 
      detectorDepth =  0.0;
      mie_param     =  0.0;
      R             =  4.7462;                      // PMT radius [cm]
      A             =  PI*R*R;                      // PMT area   [cm^2]
    }
 
 
    /**
     * Type definition of value and unit.
     */
    typedef std::pair<double, std::string> pair_type;
 
 
    /**
     * Get value and unit.
     *
     * \param  buffer       input
     * \return              value and input
     */
    static inline pair_type parse(const std::string& buffer)
    {
      using namespace std;
 
      const string::size_type pos = buffer.find(SEPARATOR);
 
      if (pos != string::npos) {
 
        double value;
 
        istringstream(buffer.substr(0,pos)) >> value;
 
        return make_pair(value, buffer.substr(pos+1));
 
      } else {
 
        THROW(JParseError, "Error parsing <" << buffer << ">.");
      }
    }
 
 
    /**
     * Read KM3Sim inputs from input strea.
     *
     * \param  in             input stream
     * \param  object         KM3Sim inputs
     * \return                input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JKM3Sim& object)
    {
      using namespace std;
 
      for (string key; in >> key; ) {
 
        string  buffer[2];
        size_t  N;
        double  value;
 
        if        (key == detectorDepth_t) {
 
          if (in >> buffer[0]) {
            object.detectorDepth = parse(buffer[0]).first;
          }
 
        } else if (key == PPCKOV_t) {
 
          if (in >> N) {
            for (size_t i = 0; i != N; ++i) {
              if (in >> buffer[0]) {
                object.ppckov.push_back(parse(buffer[0]).first);
              }
            }
          }
 
        } else if (key == RINDEX_WATER_t) {
 
          for (size_t i = 0; i != object.ppckov.size(); ++i) {
            if (in >> value) {
              object.rindex_water.push_back(value);
            }
          }
 
        } else if (key == ABSORPTION_WATER_t) {
 
          for (in >> buffer[0]; in >> buffer[0] >> buffer[1] && buffer[0] != END_t; ) {
            object.absorption_water[parse(buffer[0]).first] = parse(buffer[1]).first;
          }
 
        } else if (key == MIE_WATER_t) {
 
          for (size_t i = 0; i != object.ppckov.size(); ++i) {
            if (in >> buffer[0]) {
              object.mie_water.push_back(parse(buffer[0]).first);
            }
          }
 
        } else if (key == Q_EFF_t) {
 
          double QE;
 
          in >> QE;
      
          for (size_t i = 0; i != object.ppckov.size(); ++i) {
            if (in >> value) {
              object.q_eff.push_back(value*QE);
            }
          }
 
        } else if (key == MIE_PARAM_t) {
 
          in >> object.mie_param;
 
        } else if (key == COS_ANGLES_t) {
 
          if (in >> N) {
            for (size_t i = 0; i != N; ++i) {
              if (in >> value) {
                object.cos_angles.push_back(value);
              }
            }
          }
 
        } else if (key == ANG_ACCEPT_t) {
 
          for (size_t i = 0; i != object.cos_angles.size(); ++i) {
            if (in >> value) {
              object.ang_accept.push_back(value);
            }
          }
        }
      }
 
      return in;
    }
 
 
    /**
     * Write KM3Sim inputs to output strea.
     *
     * \param  out            output stream
     * \param  object         KM3Sim inputs
     * \return                output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JKM3Sim& object)
    {
      using namespace std;
      using namespace JPP;
 
      out << detectorDepth_t     << ' ' << object.detectorDepth                << endl;
      out << PPCKOV_t            << ' ' << JEEPZ() << object.ppckov            << endl;
      out << RINDEX_WATER_t      << ' ' << JEEPZ() << object.rindex_water      << endl;
      out << ABSORPTION_WATER_t  << ' ' << JEEPZ() << object.absorption_water  << endl;
      out << MIE_WATER_t         << ' ' << JEEPZ() << object.mie_water         << endl;
      out << Q_EFF_t             << ' ' << JEEPZ() << object.q_eff             << endl;
      out << MIE_PARAM_t         << ' ' << object.mie_param                    << endl;
      out << COS_ANGLES_t        << ' ' << JEEPZ() << object.cos_angles        << endl;
      out << ANG_ACCEPT_t        << ' ' << JEEPZ() << object.ang_accept        << endl;
 
      return out;
    }
 
 
    double R;                                 // PMT radius [cm]
    double A;                                 // PMT area   [cm^2]
  
    double                    detectorDepth;
    std::vector<double>       ppckov;
    std::vector<double>       rindex_water;
    std::map<double, double>  absorption_water;
    std::vector<double>       q_eff;
    std::vector<double>       mie_water;
    double                    mie_param;
    std::vector<double>       cos_angles;
    std::vector<double>       ang_accept;
  };
 
 
  /**
   * Auxiliary data structure for map with default value.
   */
  struct map_type :
    public std::map<std::string, double>
  {
    /**
     * Constructor.
     *
     * \param  value        default value
     */
    map_type(const double value) :
      value(value)
    {}
 
 
    /**
     * Get value for given key.
     *
     * \param  key          key
     * \return              corresponding value if key is present; else default value
     */
    double get(const std::string key) const
    {
      const_iterator p = this->find(key);
      
      if (p != this->end())
        return p->second;
      else
        return value;
    }
 
    friend inline std::istream& operator>>(std::istream& in, map_type& object) { return JPP::readObject(in, static_cast<std::map<std::string, double>&>(object)); }
    friend inline std::ostream& operator<<(std::ostream& out, const map_type& object) { return JPP::writeObject(out, static_cast<const std::map<std::string, double>&>(object)); }
 
    const double value;
  };
}
 
 
/**
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  string    inputFile;
  map_type  precision (0.0);
  map_type  correction(1.0);
  int       debug;
 
  precision[JKM3Sim::RINDEX_WATER_t]      =   1.0e-4;
  precision[JKM3Sim::ABSORPTION_WATER_t]  =   1.0e-4;
  precision[JKM3Sim::MIE_WATER_t]         =   0.2;
  precision[JKM3Sim::ANG_ACCEPT_t]        =   0.5;
  precision[JKM3Sim::Q_EFF_t]             =   0.01;
  
  try {
 
    JParser<> zap;
 
    zap['f'] = make_field(inputFile);
    zap['p'] = make_field(precision)   = JPARSER::initialised();
    zap['c'] = make_field(correction)  = JPARSER::initialised();
    zap['d'] = make_field(debug)       = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  JKM3Sim km3sim;
 
  ifstream in(inputFile.c_str());
 
  in >> km3sim;
 
  in.close();
 
  DEBUG(km3sim);
  
  using namespace KM3NET;
 
 
  {
    /*
      const double g = 9.80665;                                                   // [m/s^2]
      const double P = detectorDepth * DENSITY_SEA_WATER * 1.0e-3 * g * 1.0e1;    // [Atm]
    */
    const double P = 315.0;                                                     // [Atm]
 
    DEBUG("Pressure " << P << " [Atm]" << endl);
 
    const JDispersion dispersion(P);
 
    for (size_t i = 0; i != km3sim.ppckov.size(); ++i) {
 
      const double x = km3sim.ppckov[i];
 
      ASSERT(fabs(km3sim.rindex_water[i] - dispersion.getIndexOfRefractionPhase(x) * correction.get(JKM3Sim::RINDEX_WATER_t)) <= precision.get(JKM3Sim::RINDEX_WATER_t),
             "Test of index of refraction at " << FIXED(6,1) << x << ' ' << FIXED(7,4) << km3sim.rindex_water[i] << ' ' << FIXED(7,4) << dispersion.getIndexOfRefractionPhase(x) * correction.get(JKM3Sim::RINDEX_WATER_t));
    }
  }
  {
    for (const auto& i : km3sim.absorption_water) {
      ASSERT(fabs(i.second - getAbsorptionLength(i.first) * correction.get(JKM3Sim::ABSORPTION_WATER_t)) <= precision.get(JKM3Sim::ABSORPTION_WATER_t),
             "Test of absorption length at " << FIXED(6,1) << i.first << ' ' << FIXED(7,3) << i.second << ' ' << FIXED(7,3) << getAbsorptionLength(i.first) * correction.get(JKM3Sim::ABSORPTION_WATER_t));
    }
  }
  {
    for (size_t i = 0; i != km3sim.ppckov.size(); ++i) {
 
      const double x = km3sim.ppckov[i];
 
      ASSERT(fabs(km3sim.mie_water[i] - getScatteringLength(x) * correction.get(JKM3Sim::MIE_WATER_t)) <= precision.get(JKM3Sim::MIE_WATER_t),
             "Test of scattering length at " << FIXED(6,1) << x << ' ' << FIXED(7,3) << km3sim.mie_water[i] << ' ' << FIXED(7,3) << getScatteringLength(x) * correction.get(JKM3Sim::MIE_WATER_t));
    }
  }
  {
    for (size_t i = 0; i != km3sim.ppckov.size(); ++i) {
 
      const double x = km3sim.ppckov[i];
 
      ASSERT(fabs(km3sim.q_eff[i] - getQE(x) * correction.get(JKM3Sim::Q_EFF_t)) <= precision.get(JKM3Sim::Q_EFF_t),
             "Test of QE at " << FIXED(6,1) << x << ' ' << FIXED(7,4) << km3sim.q_eff[i] << ' ' << FIXED(7,4) << getQE(x) * correction.get(JKM3Sim::Q_EFF_t));
    }
  }
  {
    const double U = 1.0e4;                                                     // m^2 -> cm^2
 
    for (size_t i = 0; i != km3sim.cos_angles.size(); ++i) {
 
      const double x = km3sim.cos_angles[i];
 
      ASSERT(fabs(km3sim.ang_accept[i]*km3sim.A - getPhotocathodeArea() * getAngularAcceptance(x) * U * correction.get(km3sim.ANG_ACCEPT_t)) <= precision.get(km3sim.ANG_ACCEPT_t),
             "Test of PMT acceptance at " << FIXED(6,2) << x << ' ' << FIXED(9,5) << km3sim.ang_accept[i]*km3sim.A << ' ' << FIXED(9,5) << getPhotocathodeArea() * getAngularAcceptance(x) * U * correction.get(km3sim.ANG_ACCEPT_t));
    }
  }
}