JFitK40.hh

Go to the documentation of this file.

#ifndef __JCALIBRATE_JFITK40__
#define __JCALIBRATE_JFITK40__
 
#include <vector>
#include <array>
#include <map>
#include <memory>
#include <limits>
#include <cmath>
 
#include "km3net-dataformat/online/JDAQ.hh"
 
#include "JLang/JException.hh"
#include "JLang/JManip.hh"
 
#include "Jeep/JMessage.hh"
 
#include "JTools/JRange.hh"
 
#include "JDetector/JModule.hh"
 
#include "JMath/JVectorND.hh"
#include "JMath/JMatrixNS.hh"
#include "JMath/JMath.hh"
#include "JMath/JBell.hh"
#include "JMath/JMathToolkit.hh"
 
#include "JFit/JMEstimator.hh"
 
#include "JCalibrate/JCalibrateK40.hh"
#include "JCalibrate/JTDC_t.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JCALIBRATE {}
namespace JPP { using namespace JCALIBRATE; }
 
namespace JCALIBRATE {
 
  using KM3NETDAQ::NUMBER_OF_PMTS;
  using JLANG::JValueOutOfRange;
  using JDETECTOR::JModule;
  using JMATH::JMath;
  using JFIT::JMEstimator;
 
 
  /**
   * Fit options.
   */
  enum JOption_t {
    FIT_PMTS_t                         =   1,     //!< fit parameters of PMTs
    FIT_PMTS_AND_ANGULAR_DEPENDENCE_t  =   2,     //!< fit parameters of PMTs and angular dependence of K40 rate
    FIT_PMTS_AND_BACKGROUND_t          =   3,     //!< fit parameters of PMTs and background
    FIT_PMTS_QE_FIXED_t                =   4,     //!< fit parameters of PMTs with QE fixed
    FIT_MODEL_t                        =   5      //!< fit parameters of K40 rate and TTSs of PMTs
  };
 
  static const int INVALID_INDEX       =  -1;     //!< invalid index
 
  static const int NUMBER_OF_RINGS                   =  6;      //!< number of rings in optical module.
  
  static double  TEROSTAT_DZ  =  0.40;    //!< maximal PMT inclination
  static double  TEROSTAT_R1  =  1.00;    //!< scaling factor
  static double  BELL_SHAPE   =  1.55;    //!< Bell shape
 
 
  /**
   * Data structure for measured coincidence rate of pair of PMTs.
   */
  struct rate_type {
    /**
     * Default constructor.
     */
    rate_type() :
      dt_ns(0.0),
      value(0.0),
      error(0.0)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  dt_ns       time difference [ns]
     * \param  value       value of rate [Hz/ns]
     * \param  error       error of rate [Hz/ns]
     */
    rate_type(double  dt_ns,
              double  value,
              double  error) :
      dt_ns(dt_ns),
      value(value),
      error(error)
    {}
 
    double  dt_ns;    //!< time difference [ns]
    double  value;    //!< value of rate [Hz/ns]
    double  error;    //!< error of rate [Hz/ns]
  };
 
 
  /**
   * Data structure for measured coincidence rates of all pairs of PMTs in optical module.
   */
  struct data_type :
    public std::map<pair_type, std::vector<rate_type> >
  {};
 
 
  /**
   * Auxiliary class for fit parameter with optional limits.
   */
  class JParameter_t :
    public JMath<JParameter_t>
  {
  public:
    /**
     * Fit options.
     */
    enum FIT_t {
      FREE_t  =  0,                    //!< free
      FIXED_t                          //!< fixed
    };
 
 
    /**
     * Type definition for range of parameter values.
     */
    typedef JTOOLS::JRange<double>  range_type;
 
 
    /**
     * Default constructor.
     */
    JParameter_t()
    {
      set(0.0);
    }
 
 
    /**
     * Constructor.
     *
     * \param  value           value
     * \param  range           range
     */
    JParameter_t(const double      value,
                 const range_type& range = range_type::DEFAULT_RANGE()) :
      range(range)
    {
      set(value);
    }
 
 
    /** 
     * Negate parameter.
     *
     * \return                 this parameter
     */
    JParameter_t& negate()
    {
      set(-get());
 
      return *this;
    }
 
 
    /** 
     * Add parameter.
     *
     * \param  parameter       parameter
     * \return                 this parameter
     */
    JParameter_t& add(const JParameter_t& parameter)
    {
      set(get() + parameter.get());
 
      return *this;
    }
 
 
    /** 
     * Subtract parameter.
     *
     * \param  parameter       parameter
     * \return                 this parameter
     */
    JParameter_t& sub(const JParameter_t& parameter)
    {
      set(get() - parameter.get());
 
      return *this;
    }
 
 
    /** 
     * Scale parameter.
     *
     * \param  factor          multiplication factor
     * \return                 this parameter
     */
    JParameter_t& mul(const double factor)
    {
      set(get() * factor);
 
      return *this;
    }
 
 
    /** 
     * Scale parameter.
     *
     * \param  factor          division factor
     * \return                 this parameter
     */
    JParameter_t& div(const double factor)
    {
      set(get() / factor);
 
      return *this;
    }
 
 
    /** 
     * Scale parameter.
     *
     * \param  first           first  parameter
     * \param  second          second parameter
     * \return                 this parameter
     */
    JParameter_t& mul(const JParameter_t& first, const JParameter_t& second)
    {
      set(first.get() * second.get());
 
      return *this;
    }
 
 
    /**
     * Check if parameter is free.
     *
     * \return                 true if free; else false
     */
    bool isFree() const
    {
      return option == FREE_t;
    }
 
 
    /**
     * Check if parameter is fixed.
     *
     * \return                 true if fixed; else false
     */
    bool isFixed() const
    {
      return option == FIXED_t;
    }
 
 
    /**
     * Check if parameter is bound.
     *
     * \return                 true if bound; else false
     */
    bool isBound() const
    {
      return range.is_valid();
    }
 
 
    /**
     * Set current value.
     */
    void set()
    {
      option = FREE_t;
    }
 
 
    /**
     * Fix current value.
     */
    void fix()
    {
      option = FIXED_t;
    }
 
 
    /**
     * Get value.
     *
     * \return                 value
     */
    double get() const
    {
      if (isBound()) 
        return range.getLowerLimit()  +  0.5 * range.getLength() * (sin(value) + 1.0);
      else
        return value;
    }
 
 
    /**
     * Set value.
     *
     * \param  value           value
     */
    void set(const double value)
    {
      if (isBound()) 
        this->value = asin(2.0 * (range.constrain(value) - range.getLowerLimit()) / range.getLength() - 1.0);
      else
        this->value = value;
 
      set();
    }
 
 
    /**
     * Set limits.
     *
     * \param xmin             minimal value
     * \param xmax             maximal value
     */
    void setLimits(const double xmin, const double xmax)
    {
      const double x = get();
 
      range = range_type(xmin, xmax);
 
      set(x);
    }
 
 
    /**
     * Relax limits.
     */
    void relax()
    {
      const double x = get();
 
      range = range_type::DEFAULT_RANGE();
 
      set(x);
    }
 
 
    /**
     * Check if parameter is at limit.
     *
     * \param precision        precision
     * \return                 true if at limit; else false
     */
    bool atLimit(const double precision) const
    {
      if (isBound())
        return (get() - range.getLowerLimit() <= precision ||
                range.getUpperLimit() - get() <= precision);
      else
        return false;
    }
 
 
    /**
     * Fix value.
     *
     * \param  value           value
     */
    void fix(const double value)
    {
      set(value);
 
      fix();
    }
 
 
    /**
     * Get derivative of value.
     *
     * \return                 derivative of value
     */
    double getDerivative() const
    {
      if (isBound())
        return 1.0 / (0.5 * range.getLength() * cos(value));
      else
        return 1.0;
    }
 
 
    /**
     * Type conversion operator.
     *
     * \return                 value
     */
    double operator()() const
    {
      return get();
    }
 
 
    /**
     * Type conversion operator.
     *
     * \return                 value
     */
    operator double() const
    {
      return get();
    }
 
 
    /**
     * Assignment operator.
     *
     * \param  value           value
     * \return                 this parameter
     */
    JParameter_t& operator=(double value)
    {
      set(value);
 
      return *this;
    }
 
 
    /**
     * Read parameter from input stream.
     *
     * \param  in             input stream
     * \param  object          parameter
     * \return                 input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JParameter_t& object)
    {
      return in >> object.value;
    }
 
 
    /**
     * Write parameter to output stream.
     *
     * \param  out             output stream
     * \param  object          parameter
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JParameter_t& object)
    {
      using namespace std;
 
      out << FIXED(12,6) <<  object.get()                       << ' '
          << setw(5)     << (object.isFixed() ? "fixed" : " ")  << ' ';
 
      if (object.isBound()) {
        out << "[" << FIXED(12,6) << object.range.getLowerLimit() << "," << FIXED(12,6) << object.range.getUpperLimit() << "]";
      }
 
      return out;
    }
 
 
    double     value   =  0.0;
    FIT_t      option  =  FREE_t;
    range_type range   =  range_type::DEFAULT_RANGE();
  };
 
 
  /**
   * Fit parameters for single PMT.
   */
  struct JPMTParameters_t {
 
    static constexpr double  QE_MIN  =  0.0;   //!< minimal QE
    static constexpr double  QE_MAX  =  2.0;   //!< maximal QE
    static constexpr double  TTS_NS  =  2.0;   //!< start value transition-time spread [ns]
 
    /**
     * Default constructor.
     */
    JPMTParameters_t()
    {
      reset();
    }
 
 
    /**
     * Get default values.
     *
     * \return              parameters
     */
    static const JPMTParameters_t& getInstance()
    {
      static JPMTParameters_t parameters;
 
      parameters.QE.range = JParameter_t::range_type(QE_MIN, QE_MAX);
 
      parameters.status   = true;      
 
      parameters.QE .set(1.0);
      parameters.TTS.set(TTS_NS);
      parameters.t0 .set(0.0);
      parameters.bg .set(0.0);
 
      return parameters;
    }
 
 
    /**
     * Reset.
     */
    void reset()
    {
      status = true;
 
      QE .set(0.0);
      TTS.set(0.0);
      t0 .set(0.0);
      bg .set(0.0);
    }
 
 
    /**
     * Set parameters that are free to given values.
     *
     * \param  parameters   parameters
     */
    void set(const JPMTParameters_t& parameters)
    {
      if (QE .isFree()) { QE .set(parameters.QE); }
      if (TTS.isFree()) { TTS.set(parameters.TTS); }
      if (t0 .isFree()) { t0 .set(parameters.t0); }
      if (bg .isFree()) { bg .set(parameters.bg); }
    }
 
    
    /**
     * Get number of fit parameters.
     *
     * \return               number of parameters
     */
    inline size_t getN() const
    {
      return ((QE. isFree() ? 1 : 0) +
              (TTS.isFree() ? 1 : 0) +
              (t0 .isFree() ? 1 : 0) +
              (bg .isFree() ? 1 : 0));
    }
 
 
    /**
     * Disable PMT.
     */
    void disable()
    {
      status = false;
 
      QE .fix(0.0);
      TTS.fix(TTS_NS);
      t0 .fix(0.0);
      bg .fix(0.0);
    }
 
 
    /**
     * Enable PMT.
     */
    void enable()
    {
      status = true;
 
      QE .set();
      TTS.set();
      t0 .set();
      bg .set();
    }
 
 
    /**
     * Write PMT parameters to output stream.
     *
     * \param  out             output stream
     * \param  object          PMT parameters
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JPMTParameters_t& object)
    {
      using namespace std;
 
      out << "QE  " << FIXED(7,3) << object.QE  << endl;
      out << "TTS " << FIXED(7,3) << object.TTS << endl;
      out << "t0  " << FIXED(7,3) << object.t0  << endl;
      out << "bg  " << FIXED(7,3) << object.bg  << endl;
 
      return out;
    }
 
 
    bool          status;         //!< status
    JParameter_t  QE;             //!< relative quantum efficiency [unit]
    JParameter_t  TTS;            //!< transition-time spread [ns]
    JParameter_t  t0;             //!< time offset [ns]
    JParameter_t  bg;             //!< background  [Hz/ns]
  };
 
 
  /**
   * Fit parameters for two-fold coincidence rate due to K40.
   */
  struct JK40Parameters_t {
    /**
     * Default constructor.
     */
    JK40Parameters_t()
    {
      reset();
    }
 
 
    /**
     * Get K40 parameters.
     *
     * \return              K40 parameters
     */
    const JK40Parameters_t& getK40Parameters() const
    {
      return static_cast<const JK40Parameters_t&>(*this);
    }
 
 
    /**
     * Set K40 parameters.
     *
     * \param  parameters   K40 parameters
     */
    void setK40Parameters(const JK40Parameters_t& parameters)
    {
      static_cast<JK40Parameters_t&>(*this) = parameters;
    }
 
 
    /**
     * Reset.
     */
    void reset()
    {
      R .set(0.0);
      p1.set(0.0);
      p2.set(0.0);
      p3.set(0.0);
      p4.set(0.0);
      cc.set(0.0);
      bc.set(0.0);
    }
 
 
    /**
     * Print model parameters to output stream conform include files.
     *
     * \param  out             output stream
     */
    void print(std::ostream& out) const
    {
      using namespace std;
 
      out << "JFitK40.hh" << endl;
      out << "parameters.R .set(" << FIXED(9,6) << this->R () << ");" << endl;
      out << "parameters.p1.set(" << FIXED(9,6) << this->p1() << ");" << endl;
      out << "parameters.p2.set(" << FIXED(9,6) << this->p2() << ");" << endl;
      out << "parameters.p3.set(" << FIXED(9,6) << this->p3() << ");" << endl;
      out << "parameters.p4.set(" << FIXED(9,6) << this->p4() << ");" << endl;
      out << "cc                " << FIXED(9,6) << this->cc() << endl;
      out << "bc                " << FIXED(9,6) << this->bc() << endl;
      out << endl;
 
      out << "JK40DefaultSimulator.hh" << endl;
      out << "static constexpr double p1  =  " << FIXED(9,6) << this->p1() << ";" << endl;
      out << "static constexpr double p2  =  " << FIXED(9,6) << this->p2() << ";" << endl;
      out << "static constexpr double p3  =  " << FIXED(9,6) << this->p3() << ";" << endl;
      out << "static constexpr double p4  =  " << FIXED(9,6) << this->p4() << ";" << endl;
      out << endl;
    }
 
 
    /**
     * Write model parameters to output stream.
     *
     * \param  out             output stream
     * \param  object          model parameters
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JK40Parameters_t& object)
    {
      using namespace std;
 
      out << "Rate [Hz] " << FIXED(12,6) << object.R  << endl;
      out << "p1        " << FIXED(12,6) << object.p1 << endl;
      out << "p2        " << FIXED(12,6) << object.p2 << endl;
      out << "p3        " << FIXED(12,6) << object.p3 << endl;
      out << "p4        " << FIXED(12,6) << object.p4 << endl;
      out << "cc        " << FIXED(12,6) << object.cc << endl;
      out << "bc        " << FIXED(12,6) << object.bc << endl;
      out << endl;
 
      return out;
    }
 
    JParameter_t R;                  //!< maximal coincidence rate [Hz]
    JParameter_t p1;                 //!< 1st order angle dependence coincidence rate
    JParameter_t p2;                 //!< 2nd order angle dependence coincidence rate
    JParameter_t p3;                 //!< 3rd order angle dependence coincidence rate
    JParameter_t p4;                 //!< 4th order angle dependence coincidence rate
    JParameter_t cc;                 //!< fraction of signal correlated background
    JParameter_t bc;                 //!< constant background
  };
 
 
  /**
   * Fit parameters for two-fold coincidence rate due to K40.
   */
  struct JK40Parameters :
    JK40Parameters_t
  {
    /**
     * Default constructor.
     */
    JK40Parameters() 
    {}
 
 
    /**
     * Get default values.
     *
     * Values obtained with $JPP_DIR/examples/JCalibrate/JOMGsim.sh type B (see $JPP_DIR/examples/JCalibrate/README.md).
     * If you change these values, you may also want to change the corresponding values in JK40DefaultSimulator.hh.
     *
     * \return              parameters
     */
    static const JK40Parameters& getInstance()
    {
      static JK40Parameters parameters;
      parameters.R .set(18.430675);
      parameters.p1.set( 2.919895);
      parameters.p2.set(-0.831970);
      parameters.p3.set( 1.407887);
      parameters.p4.set( 0.170510);
      parameters.cc.set( 0.0);
      parameters.bc.set( 0.0);
 
      return parameters;
    }
 
 
    /**
     * Get number of fit parameters.
     *
     * \return               number of parameters
     */
    inline size_t getN() const
    {
      return ((R .isFree() ? 1 : 0) +
              (p1.isFree() ? 1 : 0) +
              (p2.isFree() ? 1 : 0) +
              (p3.isFree() ? 1 : 0) +
              (p4.isFree() ? 1 : 0) +
              (cc.isFree() ? 1 : 0) +
              (bc.isFree() ? 1 : 0));
    }
 
 
    /**
     * Get index of parameter.
     *
     * \param  p            pointer to data member
     * \return              index
     */
    int getIndex(JParameter_t JK40Parameters::*p) const
    {
      if (!(this->*p).isFree()) {
        return INVALID_INDEX;
      }
 
      int N = 0;
 
      if (p == &JK40Parameters::R)  { return N; } if (R .isFree()) { ++N; }
      if (p == &JK40Parameters::p1) { return N; } if (p1.isFree()) { ++N; }
      if (p == &JK40Parameters::p2) { return N; } if (p2.isFree()) { ++N; }
      if (p == &JK40Parameters::p3) { return N; } if (p3.isFree()) { ++N; }
      if (p == &JK40Parameters::p4) { return N; } if (p4.isFree()) { ++N; }
      if (p == &JK40Parameters::cc) { return N; } if (cc.isFree()) { ++N; }
      if (p == &JK40Parameters::bc) { return N; } if (bc.isFree()) { ++N; }
 
      return INVALID_INDEX;
    }
 
 
    /**
     * Get K40 coincidence rate as a function of cosine angle between PMT axes.
     *
     * \param  ct           cosine angle between PMT axes
     * \return              rate [Hz]
     */
    double getValue(const double ct) const
    {
      return R * exp(ct*(p1+ct*(p2+ct*(p3+ct*p4))) - (p1+p2+p3+p4));
    }
 
 
    /**
     * Get gradient.
     *
     * \param  ct            cosine angle between PMT axes
     * \return               gradient
     */
    const JK40Parameters_t& getGradient(const double ct) const
    {
      gradient.reset();
 
      const double rate  = getValue(ct);
      const double ct2   = ct * ct;
 
      if (R .isFree()) { gradient.R  = rate / R; }
      if (p1.isFree()) { gradient.p1 = rate * ct         -  rate; }
      if (p2.isFree()) { gradient.p2 = rate * ct2        -  rate; }
      if (p3.isFree()) { gradient.p3 = rate * ct2 * ct   -  rate; }
      if (p4.isFree()) { gradient.p4 = rate * ct2 * ct2  -  rate; }
      if (cc.isFree()) { gradient.cc = rate; }
      if (bc.isFree()) { gradient.bc = 1.0; }
 
      return gradient;
    }
 
  private: 
    mutable JK40Parameters_t gradient;
  };
 
 
  /**
   * Auxiliary data structure for fraction of water to total. 
   */
  struct JWater :
    public JK40Parameters_t
  {
    /**
     * Default constructor.
     */
    JWater() 
    {
      //simple for Water
      this->R .set( 0.595927);
      this->p1.set(-0.032784);
      this->p2.set( 0.075297);
      this->p3.set(-0.076032);
      this->p4.set(-0.453102);
 
      /*
      //simple DOM low stat for Water
      this->R .set( 0.589235);
      this->p1.set(-0.018508);
      this->p2.set( 0.105162);
      this->p3.set(-0.107452);
      this->p4.set(-0.475067);
      */
 
      /*
      //full DOM for Water
      this->R .set( 0.649367);
      this->p1.set(-0.011764);
      this->p2.set(-0.282949);
      this->p3.set(-0.119943);
      this->p4.set(-0.019789);
      */
    }
 
 
    /**
     * Get fraction of water as a function of cosine angle between PMT axes.
     *
     * \param  ct           cosine angle between PMT axes
     * \return              fraction
     */
    double operator()(const double ct) const
    {
      return R * exp(ct*(p1+ct*(p2+ct*(p3+ct*p4))) - (p1+p2+p3+p4));
    }
  };
 
 
  /**
   * Function object for fraction of water to total. 
   */
  const JWater getWater;
 
 
  /**
   * Auxiliary data structure to define ring.
   * The ring values run from 'A' (down) to 'F' (up).
   */
  struct ring_type {
    /**
     * Constructor.
     *
     * \param  c            ring
     */
    ring_type(const char c) :
      c(c)
    {}
 
 
    /**
     * Get index.
     *
     * \return              index
     */
    int getIndex() const
    {
      return (c - 'A');
    }
 
 
    /**
     * Type conversion operator.
     *
     * \return              index
     */
    operator int() const
    {
      return getIndex();
    }
 
 
    /**
     * Get ring.
     *
     * \param  index        index
     * \return              ring
     */
    static inline ring_type getRing(const int index)
    {
      return { (char) ('A' + index) };
    }
 
 
    /**
     * Read ring from input stream.
     *
     * \param  in           input stream
     * \param  ring         ring
     * \return              input stream
     */
    friend inline std::istream& operator>>(std::istream& in, ring_type& ring)
    {
      return in >> ring.c;
    }
 
 
    /**
     * Write ring to output stream.
     *
     * \param  out          output stream
     * \param  ring         ring
     * \return              output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const ring_type& ring)
    {
      return out << ring.c;
    }
 
    char c;               // ring
  };
 
 
  /**
   * Type definition of indices of pair of rings.
   */
  typedef std::pair<int, int>   ring_pair;
 
 
  /**
   * Get ring.
   *
   * \param  dz           cosine zenith angle of PMT axis
   * \return              ring
   */
  inline ring_type getRing(const double dz)
  {
    const double rz[] = { // cosine zenith angle of PMT axes 
      -1.000, 
      -0.850, 
      -0.555,
      -0.295,
      +0.295,
      +0.555
    };
 
    const int N = sizeof(rz) / sizeof(rz[0]);
 
    for (int i = 1; i != N; ++i) {
      if (dz <= 0.5 * (rz[i-1] + rz[i])) {
        return ring_type::getRing(i-1);
      }
    }
 
    return ring_type::getRing(N-1);
  }
 
 
  /**
   * Auxiliary data structure to handle transmittance of glass sphere due to sedimentation.
   */
  struct JTransmittance_t :
    std::array<JParameter_t, NUMBER_OF_RINGS>
  {
    /**
     * Default constructor.
     */
    JTransmittance_t()
    {
      reset();
    }
 
 
    /**
     * Reset.
     */
    void reset()
    {
      for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
        (*this)[i].set(0.0);
      }
    }
 
 
    /**
     * Set parameters that are free to given values.
     *
     * \param  parameters   parameters
     */
    void set(const JTransmittance_t& parameters)
    {
      for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
        if ((*this)[i].isFree()) { (*this)[i].set(parameters[i].get()); }
      }
    }
  };
 
 
  /**
   * Auxiliary data structure to handle transmittance of glass sphere due to sedimentation.
   */
  struct JTransmittance :
    public JTransmittance_t
  {
    /**
     * Default constructor.
     */
    JTransmittance()
    {}
 
 
    /**
     * Get default values.
     *
     * \return              transmittance
     */
    static const JTransmittance& getInstance()
    {
      static JTransmittance transmittance;
 
      for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
        transmittance[i].set(1.0);
      }
 
      return transmittance;
    }
 
 
    /**
     * Get number of fit parameters.
     *
     * \return               number of parameters
     */
    inline size_t getN() const
    {
      int N = 0;
 
      for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
        if ((*this)[i].isFree()) {
          N += 1;
        }
      }
 
      return N;
    }
 
 
    /**
     * Get weighed contribution of water and glass.
     *
     * \param  ct            cosine angle between PMT axes
     * \param  pair          PMT ring pair
     * \return               fraction
     */
    double getValue(const double ct, const ring_pair pair) const
    {
      const double w = getWater(ct);
      
      return w * (*this)[pair.first].get() * (*this)[pair.second].get()  +  (1.0 - w);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  ct            cosine angle between PMT axes
     * \param  pair          PMT ring pair
     * \return               gradient
     */
    const JTransmittance_t& getGradient(const double ct, const ring_pair pair) const
    {
      gradient.reset();
 
      const double w = getWater(ct);
 
      gradient[pair.first]   +=  w * (*this)[pair.second].get();
      gradient[pair.second]  +=  w * (*this)[pair.first] .get();
 
      return gradient;
    }
 
 
    /**
     * Write transmittances to output stream.
     *
     * \param  out             output stream
     * \param  object          transmittances
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JTransmittance& object)
    {
      using namespace std;
 
      for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
        out << "transmittance[" << ring_type::getRing(i) << "] " << object[i].get() << endl;
      }
 
      return out;
    }
 
  private:
    mutable JTransmittance_t gradient;
  };
 
 
  /**
   * Fit model.
   */
  struct JModel_t {
 
    JK40Parameters    model;
    JTransmittance    transmittance;
    JPMTParameters_t  parameters[NUMBER_OF_PMTS];
 
 
    /**
     * Reset.
     */
    void reset()
    {
      model        .reset();
      transmittance.reset();
 
      for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
        parameters[i].reset();
      }
    }
 
 
    /**
     * Write model parameters to output stream.
     *
     * \param  out             output stream
     * \param  object          model parameters
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JModel_t& object)
    {
      using namespace std;
 
      out << object.model;
      out << object.transmittance;
 
      for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
        out << "PMT[" << FILL(2,'0') << i << FILL() << "]." << object.parameters[i].status << endl << object.parameters[i];
      }
 
      return out;
    }
  };
 
 
  /**
   * Fit model.
   *
   * In the absence of TDC constraints, the average time offset is fixed to zero.
   */
  struct JModel :
    public JModel_t,
    public JModule,
    public JCombinatorics_t
  {
    using JCombinatorics_t::getIndex;
 
    /**
     * Auxiliary data structure for derived quantities of a given PMT pair.
     */
    struct real_type {
      ring_pair pair;                      //!< PMT ring pair
      double ct;                           //!< cosine angle between PMT axes
      double t0;                           //!< time offset [ns]
      double sigma;                        //!< total width [ns]
      double signal;                       //!< combined signal
      double background;                   //!< combined background
      double cc;                           //!< correlated background
      double bc;                           //!< uncorrelated background
    };
 
 
    /**
     * Default constructor.
     */
    JModel()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  module       detector module
     * \param  parameters   K40 parameters
     * \param  TDC          TDC constraints
     * \param  option       option
     */
    JModel(const JModule&            module, 
           const JK40Parameters&     parameters,
           const JTDC_t::range_type& TDC,
           const int                 option) :
      JModule         (module),
      JCombinatorics_t(module)
    {
      model         = parameters;
      transmittance = JTransmittance::getInstance();
 
      for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
        this->parameters[i] = JPMTParameters_t::getInstance();
      }
 
      for (JTDC_t::const_iterator i = TDC.first; i != TDC.second; ++i) {
 
        if (i->second != JTDC_t::WILDCARD) {
 
          this->parameters[i->second].t0.fix();
 
        } else {
 
          for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
            this->parameters[i].t0.fix();
          }
        }
      }
 
      this->index = (TDC.first == TDC.second ? 0 : INVALID_INDEX);
 
      setOption(option);
    }
 
 
    /**
     * Constructor.
     *
     * \param  module       detector module
     * \param  parameters   K40 parameters
     */
    JModel(const JModule&            module, 
           const JK40Parameters&     parameters) :
      JModule         (module),
      JCombinatorics_t(module)
    {
      model         = parameters;
      transmittance = JTransmittance::getInstance();
 
      for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
        this->parameters[i] = JPMTParameters_t::getInstance();
      }
    }
 
 
    /**
     * Get fit option.
     *
     * \return              option
     */
    JOption_t getOption() const
    {
      return option;
    }
 
 
    /**
     * Set fit option.
     *
     * \param  option       option
     */
    inline void setOption(const int option)
    {
      switch (option) {
 
      case FIT_PMTS_t:
 
        model.R .fix();
        model.p1.fix();
        model.p2.fix();
        model.p3.fix();
        model.p4.fix();
 
        for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
          transmittance[i].fix();
        }
 
        for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
          parameters[i].bg.fix();
        }
 
        break;
 
      case FIT_PMTS_AND_ANGULAR_DEPENDENCE_t:
 
        model.R .fix();
 
        for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
          transmittance[i].fix();
        }
 
        for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
          parameters[i].bg.fix();
        }
 
        break;
 
      case FIT_PMTS_AND_BACKGROUND_t:
 
        model.R .fix();
        model.p1.fix();
        model.p2.fix();
        model.p3.fix();
        model.p4.fix();
 
        for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
          transmittance[i].fix();
        }
 
        break;
 
      case FIT_PMTS_QE_FIXED_t:
 
        model.R .fix();
        model.p1.fix();
        model.p2.fix();
        model.p3.fix();
        model.p4.fix();
 
        for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
          transmittance[i].fix();
        }
 
        for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
          parameters[i].QE.fix();
          parameters[i].bg.fix();
        }
 
        break;
 
      case FIT_MODEL_t:
 
        index = INVALID_INDEX;
 
        //cc.fix();
        //bc.fix();
 
        for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
          transmittance[i].fix();
        }
 
        for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
          parameters[i].QE.fix();
          parameters[i].t0.fix();
          parameters[i].bg.fix();
        }
 
        break;
 
      default:
 
        THROW(JValueOutOfRange, "Invalid option " << option);
      }
 
      this->option = static_cast<JOption_t>(option);
    }
 
 
    /**
     * Check if time offset is fixed.
     *
     * \return              true if time offset fixed; else false
     */
    bool hasFixedTimeOffset() const
    {
      return index != INVALID_INDEX;
    }
 
 
    /**
     * Get time offset.
     *
     * \return              time offset
     */
    double getFixedTimeOffset() const
    {
      double t0 = 0.0;
      size_t N  = 0;
 
      for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
        if (parameters[i].t0.isFree()) {
          t0 += parameters[i].t0;
          N  += 1;
        }
      }
      
      return t0 /= N;
    }
 
    
    /**
     * Get index of PMT used for fixed time offset.
     *
     * \return              index
     */
    int getIndex() const
    {
      return index;
    }
 
 
    /**
     * Set index of PMT used for fixed time offset.
     */
    void setIndex()
    {
      if (index != INVALID_INDEX) {
 
        for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
 
          if (parameters[i].status) {
 
            index = i;
 
            parameters[index].t0.fix();
 
            return;
          }
        }
 
        THROW(JValueOutOfRange, "No valid index.");
      }
    }
 
 
    /**
     * Get number of fit parameters.
     *
     * \return              number of parameters
     */
    inline size_t getN() const
    {
      size_t N = 0;
 
      N += model        .getN();
      N += transmittance.getN();
 
      for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
        N += parameters[i].getN();
      }
 
      return N;
    }
 
 
    /**
     * Get intrinsic K40 arrival time spread.
     *
     * \return              sigma [ns]
     */
    double getSigmaK40() const
    {
      return this->sigmaK40_ns;
    }
 
 
    /**
     * Set intrinsic K40 arrival time spread.
     *
     * \param  sigma        sigma [ns]
     */
    void setSigmaK40(const double sigma)
    {
      this->sigmaK40_ns = sigma;
    }
 
 
    /**
     * Get derived quantities.
     *
     * \param  pair         PMT pair
     * \return              quantities
     */
    const real_type& getReal(const pair_type& pair) const
    {
      real.pair.first  =  getRing((*this)[pair.first] .getDZ()).getIndex();
      real.pair.second =  getRing((*this)[pair.second].getDZ()).getIndex();
 
      real.ct         =  JPP::getDot((*this)[pair.first].getDirection(), (*this)[pair.second].getDirection());
 
      real.t0         = (pair.first  == this->index  ?  -this->parameters[pair.second].t0()  :
                         pair.second == this->index  ?  +this->parameters[pair.first ].t0()  :
                         this->parameters[pair.first].t0() - this->parameters[pair.second].t0());
      
      real.sigma      =  sqrt(this->parameters[pair.first ].TTS() * this->parameters[pair.first ].TTS()  +
                              this->parameters[pair.second].TTS() * this->parameters[pair.second].TTS()  +
                              this->getSigmaK40()                 * this->getSigmaK40());
      
      real.signal     =  this->parameters[pair.first].QE() * this->parameters[pair.second].QE();
 
      real.background =  this->parameters[pair.first].bg() + this->parameters[pair.second].bg();
 
      real.cc         =  real.signal * model.cc();
      real.bc         =  model.bc();
 
      const double z1 = (*this)[pair.first ].getDirection().getDZ();
      const double z2 = (*this)[pair.second].getDirection().getDZ();
 
      if (fabs(z1) <= TEROSTAT_DZ &&
          fabs(z2) <= TEROSTAT_DZ &&
          signbit(z1) != signbit(z2)) {
 
        real.signal *= TEROSTAT_R1;
      }
 
      return real;
    }
 
 
    /**
     * Get K40 coincidence rate as a function of cosine angle between PMT axes.
     *
     * \param  ct           cosine angle between PMT axes
     * \return              rate [Hz]
     */
    double getValue(const double ct) const
    {
      return model.getValue(ct);
    }
 
 
    /**
     * Get K40 coincidence rate.
     *
     * \param  pair         PMT pair
     * \param  dt_ns        time difference [ns]
     * \return              rate [Hz/ns]
     */
    double getValue(const pair_type& pair, const double dt_ns) const
    {
      using namespace std;
      using namespace JPP;
 
      const real_type& real = getReal(pair);
 
      const JBell bell(real.t0, real.sigma, real.signal, 0.0, BELL_SHAPE);
 
      const double R1  = model.getValue(real.ct);
      const double R2  = bell .getValue(dt_ns);
 
      return real.bc  +  real.background  +  R1 * (real.cc + R2);
    }
 
 
    /**
     * Write model parameters to output stream.
     *
     * \param  out             output stream
     * \param  object          model parameters
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JModel& object)
    {
      using namespace std;
 
      out << "Module    " << setw(10)  << object.getID() << endl;
      out << "option    " << object.option << endl;
      out << "index     " << object.index  << endl;
 
      out << static_cast<const JModel_t&>(object);
 
      return out;
    }
 
  private:
    int               index;                        //!< index of PMT used for fixed time offset
    double            sigmaK40_ns  =  0.54;         //!< intrinsic K40 arrival time spread [ns]
    JOption_t         option;                       //!< fit option (see JCALIBRATE::JOption_t)
    mutable real_type real;
  };
 
 
  /**
   * Fit.
   */
  class JFit
  {
  public:
    /**
     * Result type.
     */
    struct result_type {
      double chi2;
      int    ndf;
    };
 
    typedef std::shared_ptr<JMEstimator>  estimator_type;
  
 
    /**
     * Constructor
     *
     * \param  option       M-estimator
     * \param  debug        debug
     */
    JFit(const int option, const int debug) :
      debug(debug)
    {
      using namespace JPP;
 
      estimator.reset(getMEstimator(option));
    }
    
 
    /**
     * Fit.
     *
     * \param  data         data
     * \return              chi2, NDF
     */
    result_type operator()(const data_type& data)
    {
      using namespace std;
      using namespace JPP;
 
 
      value.setIndex();
 
      const size_t N = value.getN();
 
      V.resize(N);
      Y.resize(N);
      h.resize(N);
 
      double xmax = numeric_limits<double>::lowest();
      double xmin = numeric_limits<double>::max();
 
      int ndf = 0;
 
      for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
 
        const pair_type& pair = ix->first;
 
        if (value.parameters[pair.first ].status &&
            value.parameters[pair.second].status) {
 
          ndf += ix->second.size();
 
          for (const rate_type& iy : ix->second) {
            if (iy.dt_ns > xmax) { xmax = iy.dt_ns; }
            if (iy.dt_ns < xmin) { xmin = iy.dt_ns; }
          }
        }
      }
 
      ndf -= value.getN();
 
      if (ndf < 0) {
        return { 0.0, ndf };
      }
 
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        if (value.parameters[pmt].t0.isFree()) {
          value.parameters[pmt].t0.setLimits(xmin, xmax);
        }
      }
 
 
      lambda    = LAMBDA_MIN;
 
      double precessor = numeric_limits<double>::max();
 
      for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
 
        DEBUG("step:     " << numberOfIterations << endl);
 
        evaluate(data);
 
        DEBUG("lambda:   " << FIXED(12,5) << lambda    << endl);
        DEBUG("chi2:     " << FIXED(12,3) << successor << endl);
 
        if (successor < precessor) {
 
          if (numberOfIterations != 0) {
 
            if (fabs(precessor - successor) < EPSILON) {
 
              seterr(data);
 
              return { successor / estimator->getRho(1.0), ndf };
            }
 
            if (lambda > LAMBDA_MIN) {
              lambda /= LAMBDA_DOWN;
            }
          }
 
          precessor = successor;
          previous  = value;
 
        } else {
 
          value   = previous;
          lambda *= LAMBDA_UP;
 
          if (lambda > LAMBDA_MAX) {
            break;
          }
 
          evaluate(data);
        }
 
        if (debug >= debug_t) {
 
          size_t row = 0;
 
          if (value.model.R .isFree()) { cout << "R  " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.model.p1.isFree()) { cout << "p1 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.model.p2.isFree()) { cout << "p2 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.model.p3.isFree()) { cout << "p3 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.model.p4.isFree()) { cout << "p4 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.model.cc.isFree()) { cout << "cc " << FIXED(12,3) << Y[row] << endl; ++row; }
 
          for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
            if (value.transmittance[i].isFree()) { cout << "transmittance[" << ring_type::getRing(i) << "] " << FIXED(9,6) << Y[row] << endl; ++row; }
          }
 
          for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
            if (value.parameters[pmt].QE .isFree()) { cout << "PMT[" << setw(2) << pmt << "].QE  " << FIXED(12,5) << Y[row] << endl; ++row; }
            if (value.parameters[pmt].TTS.isFree()) { cout << "PMT[" << setw(2) << pmt << "].TTS " << FIXED(12,5) << Y[row] << endl; ++row; }
            if (value.parameters[pmt].t0 .isFree()) { cout << "PMT[" << setw(2) << pmt << "].t0  " << FIXED(12,5) << Y[row] << endl; ++row; }
            if (value.parameters[pmt].bg .isFree()) { cout << "PMT[" << setw(2) << pmt << "].bg  " << FIXED(12,5) << Y[row] << endl; ++row; }
          }
        }
 
        // force definite positiveness
 
        for (size_t i = 0; i != N; ++i) {
 
          if (V(i,i) < PIVOT) {
            V(i,i) = PIVOT;
          }
 
          h[i] = 1.0 / sqrt(V(i,i));
        }
 
        // normalisation
 
        for (size_t i = 0; i != N; ++i) {
          for (size_t j = 0; j != i; ++j) {
            V(j,i) *= h[i] * h[j];
            V(i,j)  = V(j,i);
          }
        }
 
        for (size_t i = 0; i != N; ++i) {
          V(i,i) = 1.0 + lambda;
        }
 
        // solve A x = b
 
        for (size_t col = 0; col != N; ++col) {
          Y[col] *= h[col];
        }
 
        try {
          V.solve(Y);
        }
        catch (const exception& error) {
 
          ERROR("JGandalf: " << error.what() << endl << V << endl);
 
          break;
        }
 
        // update value
 
        const double factor = 2.0;
 
        size_t row = 0;
 
        if (value.model.R .isFree()) { value.model.R  -= factor * h[row] * Y[row];  ++row; }
        if (value.model.p1.isFree()) { value.model.p1 -= factor * h[row] * Y[row];  ++row; }
        if (value.model.p2.isFree()) { value.model.p2 -= factor * h[row] * Y[row];  ++row; }
        if (value.model.p3.isFree()) { value.model.p3 -= factor * h[row] * Y[row];  ++row; }
        if (value.model.p4.isFree()) { value.model.p4 -= factor * h[row] * Y[row];  ++row; }
        if (value.model.cc.isFree()) { value.model.cc -= factor * h[row] * Y[row];  ++row; }
        if (value.model.bc.isFree()) { value.model.bc -= factor * h[row] * Y[row];  ++row; }
 
        for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
          if (value.transmittance[i].isFree()) { value.transmittance[i] -= factor * h[row] * Y[row];  ++row; }
        }
 
        for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
          if (value.parameters[pmt].QE .isFree()) { value.parameters[pmt].QE  -= factor * h[row] * Y[row];  ++row; }
          if (value.parameters[pmt].TTS.isFree()) { value.parameters[pmt].TTS -= factor * h[row] * Y[row];  ++row; }
          if (value.parameters[pmt].t0 .isFree()) { value.parameters[pmt].t0  -= factor * h[row] * Y[row];  ++row; }
          if (value.parameters[pmt].bg .isFree()) { value.parameters[pmt].bg  -= factor * h[row] * Y[row];  ++row; }
        }
      }
 
      seterr(data);
 
      return { precessor / estimator->getRho(1.0), ndf };
    }
    
 
    static constexpr int    MAXIMUM_ITERATIONS = 100000;                                  //!< maximal number of iterations.
    static constexpr double EPSILON            = 1.0e-3;                                  //!< maximal distance to minimum.
    static constexpr double LAMBDA_MIN         = 1.0e-2;                                  //!< minimal value control parameter
    static constexpr double LAMBDA_MAX         = 1.0e+4;                                  //!< maximal value control parameter
    static constexpr double LAMBDA_UP          =  10.0;                                   //!< multiplication factor control parameter
    static constexpr double LAMBDA_DOWN        =  10.0;                                   //!< multiplication factor control parameter
    static constexpr double PIVOT              = std::numeric_limits<double>::epsilon();  //!< minimal value diagonal element of matrix
 
    int                         debug;
    estimator_type              estimator;                                                //!< M-Estimator function
 
    double                      lambda;
    JModel                      value;
    JModel_t                    error;
    int                         numberOfIterations;
    JMATH::JMatrixNS            V;
 
    bool                        TEST = false;
 
  private:
    /**
     * Evaluation of fit.
     *
     * \param  data         data
     */
    void evaluate(const data_type& data)
    {
      using namespace std;
      using namespace JPP;
      
      typedef JModel::real_type real_type;
 
 
      successor = 0.0;
 
      V.reset();
      Y.reset();
 
 
      // model parameter indices
 
      const struct M_t {
        M_t(const JModel& value)
        {
          R  = value.model.getIndex(&JK40Parameters_t::R);
          p1 = value.model.getIndex(&JK40Parameters_t::p1);
          p2 = value.model.getIndex(&JK40Parameters_t::p2);
          p3 = value.model.getIndex(&JK40Parameters_t::p3);
          p4 = value.model.getIndex(&JK40Parameters_t::p4);
          cc = value.model.getIndex(&JK40Parameters_t::cc);
          bc = value.model.getIndex(&JK40Parameters_t::bc);
        }
 
        int R;
        int p1;
        int p2;
        int p3;
        int p4;
        int cc;
        int bc;
 
      } M(value);
 
 
      // transmittance indices
 
      const struct T_t : public std::array<int, NUMBER_OF_RINGS> {
        T_t(const JModel& value) 
        {
          for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
            (*this)[i] = INVALID_INDEX;
          }
 
          int N = value.model.getN();
 
          for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
            if (value.transmittance[i].isFree()) { (*this)[i] = N; ++N; } 
          }
        }
      } T(value);
 
 
      // PMT parameter indices
 
      struct i_t {
        i_t() :
          QE (INVALID_INDEX),
          TTS(INVALID_INDEX),
          t0 (INVALID_INDEX),
          bg (INVALID_INDEX)
        {}
 
        int QE;
        int TTS;
        int t0;
        int bg;
      };
 
      const struct I_t : public std::array<i_t, NUMBER_OF_PMTS> {
        I_t(const JModel& value)
        {
          int N = value.model.getN()  +  value.transmittance.getN();
 
          for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
            if (value.parameters[i].QE .isFree()) { (*this)[i].QE  = N; ++N; }
            if (value.parameters[i].TTS.isFree()) { (*this)[i].TTS = N; ++N; }
            if (value.parameters[i].t0 .isFree()) { (*this)[i].t0  = N; ++N; }
            if (value.parameters[i].bg .isFree()) { (*this)[i].bg  = N; ++N; }
          }
        }
          
      } I(value);
 
 
      struct buffer_type : public vector< pair<int, double> > {
        double operator[](const int index) const
        {
          for (const_iterator i = this->begin(); i != this->end(); ++i) {
            if (i->first == index) {
              return i->second;
            }
          }
 
          THROW(JValueOutOfRange, "Invalid index " << index);
        }
      };
 
      buffer_type  buffer;
 
#define PUSH_BACK(i,v) if (i != INVALID_INDEX) { buffer.push_back({i, v}); }
 
 
      size_t number_of_errors = 0;
 
      for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
 
        const pair_type& pair = ix->first;
 
        if (value.parameters[pair.first ].status &&
            value.parameters[pair.second].status) {
 
          const real_type& real = value.getReal(pair);
 
          const JBell bell(real.t0, real.sigma, real.signal, 0.0, BELL_SHAPE);
 
          const double            R1  = value.model        .getValue   (real.ct);
          const double            T1  = value.transmittance.getValue   (real.ct, real.pair);
          const JK40Parameters_t& R1p = value.model        .getGradient(real.ct);
          const JTransmittance_t& T1p = value.transmittance.getGradient(real.ct, real.pair);
 
          for (const rate_type& iy : ix->second) {
 
            const double   R2  = bell.getValue   (iy.dt_ns);
            const JBell_t& R2p = bell.getGradient(iy.dt_ns);        
 
            const double  R   = real.bc  +  real.background  +  T1 * R1 * (real.cc + R2);
            const double  u   =  (iy.value - R)       / iy.error;
            const double  w   = -estimator->getPsi(u) / iy.error;
 
            successor += estimator->getRho(u);
 
            buffer.clear();
 
            PUSH_BACK(M.R,  w * T1 * (real.cc + R2) * R1p.R () * value.model.R .getDerivative());
            PUSH_BACK(M.p1, w * T1 * (real.cc + R2) * R1p.p1() * value.model.p1.getDerivative());
            PUSH_BACK(M.p2, w * T1 * (real.cc + R2) * R1p.p2() * value.model.p2.getDerivative());
            PUSH_BACK(M.p3, w * T1 * (real.cc + R2) * R1p.p3() * value.model.p3.getDerivative());
            PUSH_BACK(M.p4, w * T1 * (real.cc + R2) * R1p.p4() * value.model.p4.getDerivative());
            PUSH_BACK(M.cc, w * T1 * real.signal    * R1p.cc() * value.model.cc.getDerivative());
            PUSH_BACK(M.bc, w *                       R1p.bc() * value.model.bc.getDerivative());
 
            for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
              PUSH_BACK(T[i], w * R1 * (real.cc + R2) * T1p[i] * value.transmittance[i].getDerivative());
            }
 
            PUSH_BACK(I[pair.first] .QE,  w * T1 * R1 * R2p.signal * value.parameters[pair.second].QE () * value.parameters[pair.first ].QE .getDerivative());
            PUSH_BACK(I[pair.second].QE,  w * T1 * R1 * R2p.signal * value.parameters[pair.first ].QE () * value.parameters[pair.second].QE .getDerivative());
            PUSH_BACK(I[pair.first] .TTS, w * T1 * R1 * R2p.sigma  * value.parameters[pair.first ].TTS() * value.parameters[pair.first ].TTS.getDerivative() / real.sigma);
            PUSH_BACK(I[pair.second].TTS, w * T1 * R1 * R2p.sigma  * value.parameters[pair.second].TTS() * value.parameters[pair.second].TTS.getDerivative() / real.sigma);
            PUSH_BACK(I[pair.first] .t0,  w * T1 * R1 * R2p.mean   * value.parameters[pair.first ].t0 .getDerivative() * +1.0);
            PUSH_BACK(I[pair.second].t0,  w * T1 * R1 * R2p.mean   * value.parameters[pair.second].t0 .getDerivative() * -1.0);
            PUSH_BACK(I[pair.first] .bg,  w * value.parameters[pair.first ].bg .getDerivative());
            PUSH_BACK(I[pair.second].bg,  w * value.parameters[pair.second].bg .getDerivative());
 
            if (TEST) {
 
              DEBUG("PMT pair(" << setw(2) << pair.first << "," << setw(2) << pair.second << ") " << FIXED(7,3) << iy.dt_ns << " [ns]" << endl);
 
              const double PRECISION = 1.0e-5;
 
#define MAKE_TEST(i,v) if (i != INVALID_INDEX) {                        \
                                                                        \
                const bool status  =  fabs(buffer[i] - v) <= PRECISION; \
                                                                        \
                DEBUG((status ? GREEN : RED)                            \
                      << setw(20)    << left << #i << right << ' '      \
                      << setw(3)     << i                   << ' '      \
                      << FIXED(12,5) << buffer[i]           << ' '      \
                      << FIXED(12,5) << v                   << ' '      \
                      << (!status ? "***" : "")                         \
                      << RESET << endl);                                \
                                                                        \
                if (!status) {                                          \
                  number_of_errors += 1;                                \
                }                                                       \
              }
 
              struct JTest_t : public JModel
              {
                JTest_t& operator=(const JModel& model)
                {
                  static_cast<JModel&>(*this) = model;
 
                  this->model.R .relax();
                  this->model.p1.relax();
                  this->model.p2.relax();
                  this->model.p3.relax();
                  this->model.p4.relax();
                  this->model.cc.relax();
                  this->model.bc.relax();
 
                  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
                    parameters[i].QE .relax();
                    parameters[i].TTS.relax();
                    parameters[i].t0 .relax();
                    parameters[i].bg .relax();
                  }
 
                  return *this;
                }
 
                double operator()(const pair_type& pair, const rate_type& iy) const
                {
                  return (iy.value - getValue(pair, iy.dt_ns)) / iy.error;
                }
              };
 
              const double DX =  1.0e-8;         // dx
              JTest_t      m1, m2;               // y1, y2
 
              // derivative
 
              auto fp = [&DX = DX, &m1 = m1, &m2 = m2, &estimator = estimator](const pair_type& pair, const rate_type& iy)
              {
                return (estimator->getRho(m2(pair, iy)) - estimator->getRho(m1(pair, iy))) / DX;
              };
 
              { m1 = m2 = value; m2.model.R  += 0.5*DX; m1.model.R  -= 0.5*DX;  MAKE_TEST(M.R,  fp(pair, iy) * value.model.R .getDerivative()); }
              { m1 = m2 = value; m2.model.p1 += 0.5*DX; m1.model.p1 -= 0.5*DX;  MAKE_TEST(M.p1, fp(pair, iy) * value.model.p1.getDerivative()); }
              { m1 = m2 = value; m2.model.p2 += 0.5*DX; m1.model.p2 -= 0.5*DX;  MAKE_TEST(M.p2, fp(pair, iy) * value.model.p2.getDerivative()); }
              { m1 = m2 = value; m2.model.p3 += 0.5*DX; m1.model.p3 -= 0.5*DX;  MAKE_TEST(M.p3, fp(pair, iy) * value.model.p3.getDerivative()); }
              { m1 = m2 = value; m2.model.p4 += 0.5*DX; m1.model.p4 -= 0.5*DX;  MAKE_TEST(M.p4, fp(pair, iy) * value.model.p4.getDerivative()); }
              { m1 = m2 = value; m2.model.cc += 0.5*DX; m1.model.cc -= 0.5*DX;  MAKE_TEST(M.cc, fp(pair, iy) * value.model.cc.getDerivative()); }
              { m1 = m2 = value; m2.model.bc += 0.5*DX; m1.model.bc -= 0.5*DX;  MAKE_TEST(M.bc, fp(pair, iy) * value.model.bc.getDerivative()); }
 
              { m1 = m2 = value; m2.parameters[pair.first] .QE  += 0.5*DX; m1.parameters[pair.first] .QE  -= 0.5*DX;  MAKE_TEST(I[pair.first] .QE,  fp(pair, iy) * value.parameters[pair.first] .QE .getDerivative()); }
              { m1 = m2 = value; m2.parameters[pair.second].QE  += 0.5*DX; m1.parameters[pair.second].QE  -= 0.5*DX;  MAKE_TEST(I[pair.second].QE,  fp(pair, iy) * value.parameters[pair.second].QE .getDerivative()); }
              { m1 = m2 = value; m2.parameters[pair.first] .TTS += 0.5*DX; m1.parameters[pair.first] .TTS -= 0.5*DX;  MAKE_TEST(I[pair.first] .TTS, fp(pair, iy) * value.parameters[pair.first] .TTS.getDerivative()); }
              { m1 = m2 = value; m2.parameters[pair.second].TTS += 0.5*DX; m1.parameters[pair.second].TTS -= 0.5*DX;  MAKE_TEST(I[pair.second].TTS, fp(pair, iy) * value.parameters[pair.second].TTS.getDerivative()); }
              if (pair.first  != value.getIndex()) {
                m1 = m2 = value; m2.parameters[pair.first] .t0  += 0.5*DX; m1.parameters[pair.first] .t0  -= 0.5*DX;  MAKE_TEST(I[pair.first] .t0,  fp(pair, iy) * value.parameters[pair.first] .t0 .getDerivative());
              }
              if (pair.second != value.getIndex()) {
                m1 = m2 = value; m2.parameters[pair.second].t0  += 0.5*DX; m1.parameters[pair.second].t0  -= 0.5*DX;  MAKE_TEST(I[pair.second].t0,  fp(pair, iy) * value.parameters[pair.second].t0 .getDerivative());
              }
              { m1 = m2 = value; m2.parameters[pair.first] .bg  += 0.5*DX; m1.parameters[pair.first] .bg  -= 0.5*DX;  MAKE_TEST(I[pair.first] .bg,  fp(pair, iy) * value.parameters[pair.first] .bg .getDerivative()); }
              { m1 = m2 = value; m2.parameters[pair.second].bg  += 0.5*DX; m1.parameters[pair.second].bg  -= 0.5*DX;  MAKE_TEST(I[pair.second].bg,  fp(pair, iy) * value.parameters[pair.second].bg .getDerivative()); }
 
              cout << endl;
            }
 
            for (buffer_type::const_iterator row = buffer.begin(); row != buffer.end(); ++row) {
 
              Y[row->first] += row->second;
 
              V[row->first][row->first] += row->second * row->second;
 
              for (buffer_type::const_iterator col = buffer.begin(); col != row; ++col) {
                V[row->first][col->first] += row->second * col->second;
                V[col->first][row->first]  = V[row->first][col->first];
              }
            }
          }
        }
      }
 
#undef PUSH_BACK
 
      if (TEST) {
 
        STATUS("Test finished with " << number_of_errors << " errors." << endl);
 
        exit(number_of_errors == 0 ? 0 : 1);
      }
    }
 
 
    /**
     * Set errors.
     *
     * \param  data         data
     */
    void seterr(const data_type& data)
    {
      using namespace std;
 
      error.reset();
 
      evaluate(data);
      
      try {
        V.invert();
      }
      catch (const exception& error) {}
 
#define SQRT(X) (X >= 0.0 ? sqrt(X) : std::numeric_limits<double>::max())
 
      size_t row = 0;
 
      if (value.model.R .isFree()) { error.model.R  = SQRT(V(row,row));  ++row; }
      if (value.model.p1.isFree()) { error.model.p1 = SQRT(V(row,row));  ++row; }
      if (value.model.p2.isFree()) { error.model.p2 = SQRT(V(row,row));  ++row; }
      if (value.model.p3.isFree()) { error.model.p3 = SQRT(V(row,row));  ++row; }
      if (value.model.p4.isFree()) { error.model.p4 = SQRT(V(row,row));  ++row; }
      if (value.model.cc.isFree()) { error.model.cc = SQRT(V(row,row));  ++row; }
      if (value.model.bc.isFree()) { error.model.bc = SQRT(V(row,row));  ++row; }
 
      for (int i = 0; i != NUMBER_OF_RINGS; ++i) {
        if (value.transmittance[i].isFree()) { error.transmittance[i] = SQRT(V(row,row));  ++row; }
      }
 
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        if (value.parameters[pmt].QE .isFree()) { error.parameters[pmt].QE  = SQRT(V(row,row));  ++row; }
        if (value.parameters[pmt].TTS.isFree()) { error.parameters[pmt].TTS = SQRT(V(row,row));  ++row; }
        if (value.parameters[pmt].t0 .isFree()) { error.parameters[pmt].t0  = SQRT(V(row,row));  ++row; }
        if (value.parameters[pmt].bg .isFree()) { error.parameters[pmt].bg  = SQRT(V(row,row));  ++row; }
      }
 
#undef SQRT
    }
 
    
    JMATH::JVectorND            Y;                             // gradient
    double                      successor;
    JModel                      previous;
    std::vector<double>         h;                             // normalisation vector
  };
}
 
#endif