JWifi.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
#include <chrono>
#include <map>
#include <sstream>
#include <cmath>
#include <vector>
 
#include <mutex>
#include <thread>
#include <condition_variable>
#include <vector>
#include <queue>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"
 
#include "JGizmo/JGizmoToolkit.hh"
 
#include "JPhysics/JPDF.hh"
#include "JPhysics/JPDFSupportkit.hh"
#include "JPhysics/JPhysicsSupportkit.hh"
#include "JPhysics/KM3NeT.hh"
 
#include "JMath/JConstants.hh"
 
#include "JTools/JRange.hh"
 
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JOmega3D.hh"
 
#include "JFit/JGradient.hh"
 
#include "JLang/JSTDObjectReader.hh"
 
#include "Jeep/JProperties.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  using JPHYSICS::JPDF_C;
  using JPHYSICS::JAbsorptionLength;
  using JPHYSICS::JScatteringLength;
  using JPHYSICS::JQE;
  using JFIT::JParameter_t;
 
 
  /**
   * PDF fixed parameters.
   */
  struct PDF {
    static constexpr int     numberOfPoints          =  25;
    static constexpr double  epsilon                 =  1.0e-10;
  };
 
 
  /**
   * Fit fixed parameters.
   */
  struct FIT {
    static constexpr int     number_of_iterations    = 10000;
    static constexpr int     number_of_extra_steps   =     1;
    static constexpr double  epsilon                 =   1.0e-4;
  };
 
 
  static const double PMIN     =   1.0e-6;             //!< minimal average probability 
  static const double TMIN_NS  =  -2.0;                //!< minimal arrival time [ns]
 
 
  /**
   * Auxiliary data structure for experiment.
   */
  struct experiment_type {
 
    JGIZMO::JRootObjectID Ha;         //!< scattered light
    JGIZMO::JRootObjectID Hb;         //!< direct    light
 
    /**
     * Read experiment from input stream.
     *
     * \param  in             input stream
     * \param  experiment     experiment
     * \return                input stream
     */
    friend inline std::istream& operator>>(std::istream& in, experiment_type& experiment)
    {
      return in >> experiment.Ha
                >> experiment.Hb;
    }
 
 
    /**
     * Write experiment to output stream.
     *
     * \param  out            output stream
     * \param  experiment     experiment
     * \return                output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const experiment_type& experiment)
    {
      return out << experiment.Ha << ' '
                 << experiment.Hb;
    }
  };
 
 
  typedef std::pair<double, double>  pair_type;     //!< bin edges
 
  /**
   * Auxiliary data structure for measurement.
   */
  struct bin_type {
 
    /**
     * Get chi2.
     *
     * \param  ps              average hit probability
     * \return                 chi2
     */
    double getChi2(const double ps) const
    {
      const double u  = (ratio -  ps/value) / error;
 
      return u*u;
    }
 
    pair_type R_m;                //!< distance [m]
    pair_type t_ns;               //!< time residual [ns]
    double    ratio;              //!< measured ratio
    double    error;              //!< measured error
    double    value;              //!< computed value
  };
 
 
  typedef std::vector<bin_type>                      data_type;
  typedef JLANG::JSTDObjectReader<const bin_type>    input_type;
 
  typedef JTOOLS::JRange<double>                     range_type;
 
 
  /**
   * Auxiliary class for editing value of fit parameter.
   */
  struct JEditor :
    public JParameter_t,
    public range_type
  {
    /**
     * Constructor.
     *
     * \param  ps                    parameter
     * \param  range                 allowed range
     */
    JEditor(double& ps, const range_type& range = range_type()) :
      ps(ps)
    {}
 
 
    /**
     * Apply step.
     *
     * \param  step                  step
     */
    virtual void apply(const double step) override
    {
      this->ps += step;
 
      this->ps  = constrain(this->ps);
    }
 
 
    /**
     * Get value.
     *
     * \return                       value
     */
    double getValue() const
    {
      return ps;
    }
 
  private:
    double& ps;                         //!< pointer to parameter
  };
 
 
  /**
   * Auxiliary data structure for correlated fit parameters.
   */
  struct JMixed {
 
    /**
     * Get factor of absorption length.
     *
     * \return               factor
     */
    static inline double getFactorAbsorptionLength()
    {
      return get_factor_a() * get_factor_b();
    }
 
 
    /**
     * Get offset of absorption length.
     *
     * \return               factor
     */
    static inline double getOffsetAbsorptionLength()
    {
      const double y = 1.0 / get_offset_a()  +  1.0 / get_offset_b();
 
      return 1.0 / y;
    }
 
    /**
     * Get absorption length.
     *
     * \param  lambda        wavelength of light [nm]
     * \return               absorption length   [m]
     */
    static inline double getAbsorptionLength(const double lambda)
    {
      const double y = 1.0 / (JAbsorptionLength::getAbsorptionLength(lambda) * getFactorAbsorptionLength())  +  1.0 / getOffsetAbsorptionLength();
 
      return 1.0 / y;
    }
 
 
    /**
     * Get scattering length.
     *
     * \param  lambda        wavelength of light [nm]
     * \return               scattering length   [m]
     */
    static inline double getScatteringLength(const double lambda)
    {
      return JScatteringLength::getScatteringLength(lambda);
    }
 
 
    /**
     * Get factor of quantum efficiency.
     *
     * \return               factor
     */
    static inline double getFactorQE()
    {
      return (get_factor_a() / get_factor_b())  *  ((ABSORPTION_LENGTH_M + get_offset_a()) / get_offset_a())  *  (get_offset_b() / (ABSORPTION_LENGTH_M + get_offset_b()));
    }
 
 
    /**
     * Get quantum efficiency.
     *
     * \param  lambda        wavelength of light [nm]
     * \return               quantum efficiency   [m]
     */
    static inline double getQE(const double lambda)
    {
      return JQE::getQE(lambda) * getFactorQE();
    }
 
 
    /**
     * Get reference to common factor.
     *
     * \return               factor
     */
    static double& get_factor_a()
    {
      static double a = 1.0;
 
      return a;
    }
 
 
    /**
     * Get reference to common factor.
     *
     * \return               factor
     */
    static double& get_factor_b()
    {
      static double b = 1.0;
 
      return b;
    }
 
 
    /**
     * Get reference to common offset.
     *
     * \return               offset
     */
    static double& get_offset_a()
    {
      static double a = 1.0e6;
 
      return a;
    }
 
 
    /**
     * Get reference to common offset.
     *
     * \return               offset
     */
    static double& get_offset_b()
    {
      static double b = 1.0e6;
 
      return b;
    }
 
    static constexpr double ABSORPTION_LENGTH_M  =  50.0;    //!< absorption length [m]
  };
 
 
  typedef JGEOMETRY3D::JAngle3D                      pmt_type;
  typedef std::vector<pmt_type>                      module_type;
 
 
  /**
   * Auxiliary data structure for chi2 function object.
   */
  struct JWifi {
 
    static const int ND  =  2;          //!< number of histograms
 
    /**
     * Check bin contents.
     *
     * \param  value                 value
     * \param  error                 error
     * \return                       true if bin contents are valid; else false
     */
    static bool check(const double value, const double error)
    {
      return (value > 0.0 &&
              error > 0.0);
    }
 
 
    /**
     * Auxiliary data structure for sampling of bins.
     */
    static struct bins_type {
 
      static const char SEPARATOR = 'x';
 
      /**
       * Read bins from input stream.
       *
       * \param  in             input stream
       * \param  bins           bins
       * \return                input stream
       */
      friend inline std::istream& operator>>(std::istream& in, bins_type& bins)
      {
        using namespace std;
 
        string buffer;
 
        getline(in,buffer);
 
        const size_t pos = buffer.find(SEPARATOR);
 
        if (pos != string::npos) {
 
          if (!(istringstream(buffer.substr(0,pos)) >> bins.R_m))  { in.setstate(ios::badbit); }
          if (!(istringstream(buffer.substr(pos+1)) >> bins.t_ns)) { in.setstate(ios::badbit); }
 
        } else {
 
          in.setstate(ios::badbit);
        }
 
        return in;
      }
 
 
      /**
       * Write bins to output stream.
       *
       * \param  out            output stream
       * \param  bins           bins
       * \return                output stream
       */
      friend inline std::ostream& operator<<(std::ostream& out, const bins_type& bins)
      {
        return out << bins.R_m  << SEPARATOR << bins.t_ns;
      }
 
      size_t R_m;
      size_t t_ns;
 
    } bins;
 
 
    /**
     * Get average hit probability at given time for optical module at given distance.
     *
     * \param  pdf                   PDF
     * \param  E_GeV                 muon energy [GeV]
     * \param  R_m                   distance [m]
     * \param  t_ns                  time residual [ns]
     * \param  module                module
     * \return                       probability
     */
    static double getPs(const JPDF_C&      pdf,
                        const double       E_GeV,
                        const pair_type&   R_m,
                        const pair_type&   t_ns,
                        const module_type& module)
    {
      // integrate time and average PMTs and distances
 
      const double dr = (R_m .second - R_m .first) / JWifi::bins.R_m;
      const double dt = (t_ns.second - t_ns.first) / JWifi::bins.t_ns;
 
      double ps = 0.0;
 
      for (const pmt_type& pmt : module) {
 
        for (double r1 = R_m.first + 0.5*dr; r1 < R_m.second; r1 += dr) {
 
          double npe = 0.0;
 
          for (double t1 = t_ns.first + 0.5*dt; t1 < t_ns.second; t1 += dt) {
            npe += pdf.getLightFromMuon(E_GeV, r1, pmt.getTheta(), pmt.getPhi(), t1) * dt;
          }
 
          ps += 1.0 - exp(-npe);
        }
      }
 
      return ps / (module.size() * JWifi::bins.R_m);
    }
 
 
    /**
     * Constructor.
     *
     * \param  E_GeV                 muon energy [GeV]
     * \param  ha                    histogram scattered light
     * \param  hb                    histogram direct    light
     * \param  angle_Deg             grid angle [deg]
     * \param  ns                    number of threads
     */
    JWifi(const double E_GeV,
          const TH2*   ha,
          const TH2*   hb,
          const double angle_Deg,
          const size_t ns) :
 
      pdf(KM3NET::getPhotocathodeArea(),
          JMixed::getQE,
          KM3NET::getAngularAcceptance,
          JMixed::getAbsorptionLength,
          JMixed::getScatteringLength,
          KM3NET::getScatteringProbability,
          KM3NET::getAmbientPressure(),
          JPHYSICS::getMinimalWavelength(),
          JPHYSICS::getMaximalWavelength(),
          PDF::numberOfPoints,
          PDF::epsilon),
      
      E_GeV(E_GeV),
      ns(ns)
    {
      using namespace JPP;
 
      // sets of PMTs conform 'stopping muon' definition of scattered (0) and direct (1) light
 
      const JOmega3D omega(angle_Deg * PI / 180.0);
      
      for (JOmega3D::const_iterator i = omega.begin(); i != omega.end(); ++i) {
 
        if      (i->getTheta() <= 0.5*PI && i->getPhi() <= 0.5*PI)
          module[0].push_back(*i);
        else if (i->getTheta() >= 0.5*PI && i->getPhi() >= 0.5*PI)
          module[1].push_back(*i);
      }
 
      // import data
 
      const TH2* h2[ND] = { ha, hb };
 
      for (int i = 0; i != ND; ++i) {
 
        for (Int_t ix = 1; ix <= h2[i]->GetXaxis()->GetNbins(); ++ix) {
          for (Int_t iy = 1; iy <= h2[i]->GetYaxis()->GetNbins(); ++iy) {
 
            pair_type R_m   =  { h2[i]->GetXaxis()->GetBinLowEdge(ix), h2[i]->GetXaxis()->GetBinUpEdge(ix) };
            pair_type t_ns  =  { h2[i]->GetYaxis()->GetBinLowEdge(iy), h2[i]->GetYaxis()->GetBinUpEdge(iy) };
 
            if (t_ns.first < TMIN_NS) {
              t_ns.first = TMIN_NS;
            }
 
            double value  =  h2[i]->GetBinContent(ix,iy);
            double error  =  h2[i]->GetBinError  (ix,iy);
 
            if (check(value, error)) {
 
              const double ps  =  getPs(pdf, E_GeV, R_m, t_ns, module[i]);
 
              if (ps >= PMIN) {
                buffer[i].push_back({R_m, t_ns, value, error, ps});
              }
            }
          }
        }
      }
    }
 
 
    /**
     * Get number of points.
     *
     * \return                       number of points
     */
    size_t getN() const
    {
      return buffer[0].size() + buffer[1].size();
    }
 
 
    /**
     * Get chi2.
     *
     * \param  option                option
     * \return                       chi2
     */
    double operator()(const int option) const
    {
      using namespace std;
 
      double chi2 = 0.0;
 
      for (int i = 0; i != ND; ++i) {
 
        if (ns == 0) {
 
          for (const auto& bin : buffer[i]) {
 
            const double ps =  getPs(pdf, E_GeV, bin.R_m, bin.t_ns, module[i]);
 
            chi2 += bin.getChi2(ps);
          }
 
        } else {
 
          {
            JPerth perth(pdf, E_GeV, buffer[i], module[i], ns);
          }
          
          chi2 += JPerth::chi2;
        }
      }
 
      return chi2;
    }
 
 
    const JPDF_C pdf;                //!< PDF
    const double E_GeV;              //!< muon energy
 
    const size_t ns;                 //!< number of threads
 
    data_type    buffer[ND];
    module_type  module[ND];
 
 
    /**
     * Thread pool for fits to data.
     */
    struct JPerth {
      /**
       * Constructor.
       *
       * \param  pdf                 PDF
       * \param  E_GeV               energy [GeV]
       * \param  data                data
       * \param  module              optical module
       * \param  ns                  number of threads
       */
      JPerth(const JPDF_C&       pdf,
             const double        E_GeV,
             const data_type&    data,
             const module_type&  module,
             const size_t        ns) :
        input(data),
        stop(false)
      {
        using namespace std;
        using namespace JPP;
 
        chi2 = 0.0;
 
        for (size_t i = 0; i < ns; ++i) {
 
          thread worker([this, pdf, E_GeV, module]() {
 
            for ( ; ; ) {
 
              bin_type bin;
 
              {
                unique_lock<mutex> lock(in);
 
                cv.wait(lock, [this]() { return stop || this->input.hasNext(); });
 
                if (stop && !this->input.hasNext()) {
                  return;
                }
 
                bin = *(this->input.next());
              }
 
              const double ps =  getPs(pdf, E_GeV, bin.R_m, bin.t_ns, module);
    
              {
                unique_lock<mutex> lock(out);
 
                chi2 += bin.getChi2(ps);
              }
            }
          });
 
          workers.emplace_back(std::move(worker));
        }
      }
 
 
      /**
       * Destructor.
       */   
      ~JPerth()
      {
        using namespace std;
 
        {
          unique_lock<mutex> lock(in);
 
          stop = true;
        }
      
        cv.notify_all();
 
        for (auto& worker : workers) {
          worker.join();
        }
      }
 
 
      static double             chi2;        //!< final result
 
    private:
      input_type                input;
      std::vector<std::thread>  workers;
      std::mutex                in;
      std::mutex                out;
      std::condition_variable   cv;
      bool                      stop;
    };
  };
 
  /**
   * Storage of bins.
   */
  JWifi::bins_type JWifi::bins = { 1, 1};
 
  /**
   * Storage of result of multi-threading.
   */
  double JWifi::JPerth::chi2  =  0.0;
 
 
  /**
   * Fit options.
   */
  const char* const absorption_factor_t  =  "absorption_factor";
  const char* const absorption_offset_t  =  "absorption_offset";
  const char* const scattering_factor_t  =  "scattering_factor";
  const char* const qe_factor_t          =  "qe_factor";
  const char* const mixed_factor_a_t     =  "mixed_factor_a";
  const char* const mixed_factor_b_t     =  "mixed_factor_b";
  const char* const mixed_offset_a_t     =  "mixed_offset_a";
  const char* const mixed_offset_b_t     =  "mixed_offset_b";
}
 
 
/**
 * Application to fit parameters to stopping muon data.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef std::pair<double, double>    parameter_type;
  typedef map<string, parameter_type>  option_type;
 
  experiment_type    setup;
  JAbsorptionLength  absorptionLength;
  JScatteringLength  scatteringLength;
  double             E_GeV;
  double             angle_Deg;
  size_t             threads;
  option_type        option;
  bool               normalise;
  vector<double>     xs;
  int                debug;
 
  try {
 
    JProperties properties;
 
    properties.insert(gmake_property(absorptionLength));
    properties.insert(gmake_property(scatteringLength));
 
    JParser<> zap;
 
    zap['@'] = make_field(properties, endl
                          << "possible options absorptionLength: "      << get_keys(absorptionLength)      << endl
                          << "possible options scatteringLength: "      << get_keys(scatteringLength)      << endl)  = JPARSER::initialised();
    zap['f'] = make_field(setup,
                          "douplet of histograms corresponding to scattered and direct light, "
                          << "each of which defined by <file name>:<histogram name>");
    zap['E'] = make_field(E_GeV,        "muon energy [GeV]")                             =  10.0;
    zap['G'] = make_field(angle_Deg,    "grid angle for PMT directions [deg]")           =  25.0;
    zap['M'] = make_field(JWifi::bins,  "number of points in distance and time bin")     =  JPARSER::initialised();
    zap['N'] = make_field(threads,      "number threads")                                =   0;
    zap['O'] = make_field(option,       "fit options: \"parameter <start value> <step size>\"" << endl
                          << "possible parameters: "
                          << absorption_factor_t  << ", "
                          << absorption_offset_t  << ", "
                          << scattering_factor_t  << ", "
                          << qe_factor_t          << ", "
                          << mixed_factor_a_t     << ", "
                          << mixed_factor_b_t     << ", "
                          << mixed_offset_a_t     << ", "
                          << mixed_offset_b_t);
    zap['n'] = make_field(normalise,  "normalise gradient");
    zap['x'] = make_field(xs,         "scan step sizes")                                 =  JPARSER::initialised();
    zap['d'] = make_field(debug)               = 2;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  const TH2* ha = dynamic_cast<const TH2*>(getObject(setup.Ha));
  const TH2* hb = dynamic_cast<const TH2*>(getObject(setup.Hb));
 
  if (ha == NULL) { FATAL("Missing histogram: " << setup.Ha << endl); }
  if (hb == NULL) { FATAL("Missing histogram: " << setup.Hb << endl); }
 
  const JWifi wifi(E_GeV, ha, hb, angle_Deg, threads);
 
  // fit
 
  JGradient fit(FIT::number_of_iterations, FIT::number_of_extra_steps, FIT::epsilon, 3);
 
  fit.normalise = normalise;
 
  auto fp = [&option, &fit](const string& key, double& value, const range_type& range) {
    if (option.count(key)) {
      fit.push_back(JModifier_t(key, new JEditor(value = option[key].first, range), option[key].second));
    }
  };
 
  fp(absorption_factor_t, JAbsorptionLength::get_factor(),   range_type(0.1, 1.0e1));
  fp(absorption_offset_t, JAbsorptionLength::get_offset(),   range_type(1.0, 1.0e6));
  fp(scattering_factor_t, JScatteringLength::get_factor(),   range_type(0.1, 1.0e1));
  fp(qe_factor_t,         JQE              ::get_factor(),   range_type(0.1, 1.0e1));
  fp(mixed_factor_a_t,    JMixed           ::get_factor_a(), range_type(0.1, 1.0e1));
  fp(mixed_factor_b_t,    JMixed           ::get_factor_b(), range_type(0.1, 1.0e1));
  fp(mixed_offset_a_t,    JMixed           ::get_offset_a(), range_type(1.0, 1.0e6));
  fp(mixed_offset_b_t,    JMixed           ::get_offset_b(), range_type(1.0, 1.0e6));
 
  const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
 
  const double chi2 = fit(wifi);
 
  const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
 
  cout << "chi2/NDF " << FIXED(9,5) << chi2 << "/" << setw(4) << (wifi.getN() - fit.size()) << endl;
 
  cout << "Elapsed time: " << setw(6) << (t1 - t0) / chrono::seconds(1) << " [s]" << endl;
 
  for (const auto& i : fit) {
    cout << left << setw(36) << i.name << ' ' << right << FIXED(13,3) << dynamic_cast<const JEditor*>(i.get())->getValue() << endl;
  }
 
  auto get_offset = [](const double x1, const double x2) { const double y = 1.0/x1 + 1.0/x2; return 1.0/y; };
  auto get_factor = [](const double x1, const double x2) { return x1 * x2; };
 
  cout << "Additional absorption length: " << FIXED(13,5) << get_offset(JAbsorptionLength::get_offset(), JMixed::getOffsetAbsorptionLength()) << " [m]" << endl;
  cout << "Scaling of absorption length: " << FIXED(13,5) << get_factor(JAbsorptionLength::get_factor(), JMixed::getFactorAbsorptionLength())           << endl;
  cout << "Scaling of scattering length: " << FIXED(13,5) << JScatteringLength::get_factor()                                                            << endl;
  cout << "Scaling of QE:                " << FIXED(13,5) << get_factor(JQE              ::get_factor(), JMixed::getFactorQE())                         << endl;
 
  if (!xs.empty()) {
    fit(cout, wifi, xs);
  }
}