JFitK40.hh

Go to the documentation of this file.

#ifndef __JCALIBRATE_JGANDALFK40__
#define __JCALIBRATE_JGANDALFK40__

#include <vector>
#include <map>
#include <memory>
#include <limits>
#include <cmath>

#include "km3net-dataformat/online/JDAQ.hh"

#include "JLang/JException.hh"
#include "JLang/JManip.hh"

#include "JTools/JRange.hh"

#include "JDetector/JModule.hh"

#include "JMath/JVectorND.hh"
#include "JMath/JMatrixNS.hh"
#include "JMath/JMath.hh"
#include "JMath/JGauss.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_MODEL_t                        =   2,     //!< fit parameters of K40 rate and TTSs of PMTs
    FIT_PMTS_AND_ANGULAR_DEPENDENCE_t  =   3      //!< fit parameters of PMTs and angular dependence of K40 rate
  };

  static const int INVALID_INDEX       =  -1;     //!< invalid index


  /**
   * 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;
    double  value;
    double  error;
  };


  /**
   * 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();
    }


    /**
     * 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 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.value                 << ' ';
        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();
    }


    /**
     * Reset.
     */
    void reset()
    {
      status = true;
      QE     = 0.0;
      TTS    = 0.0;
      t0     = 0.0;
    }


    /**
     * 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       = 1.0;
      parameters.TTS      = TTS_NS;
      parameters.t0       = 0.0;

      return parameters;
    }

    
    /**
     * 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));
    }


    /**
     * Get index of parameter.
     *
     * \param  p            pointer to data member
     * \return              index
     */
    int getIndex(JParameter_t JPMTParameters_t::*p) const
    {
      if (!(this->*p).isFree()) {
        return INVALID_INDEX;
      }

      int N = 0;

      if (p == &JPMTParameters_t::QE ) { return N; }; if (QE .isFree()) { ++N; }
      if (p == &JPMTParameters_t::TTS) { return N; }; if (TTS.isFree()) { ++N; }
      if (p == &JPMTParameters_t::t0 ) { return N; }; if (t0 .isFree()) { ++N; }

      return INVALID_INDEX;
    }


    /**
     * Disable PMT.
     */
    void disable()
    {
      status = false;

      QE .fix(0.0);
      TTS.fix(TTS);
      t0 .fix(0.0);
    }


    /**
     * Enable PMT.
     */
    void enable()
    {
      status = true;

      QE .set();
      TTS.set();
      t0 .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  " << object.QE  << endl;
      out << "TTS " << object.TTS << endl;
      out << "t0  " << object.t0  << 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]
  };


  /**
   * Fit parameters for two-fold coincidence rate due to K40.
   */
  struct JK40Parameters_t {
    /**
     * Default constructor.
     */
    JK40Parameters_t()
    {
      reset();
    }


    /**
     * Reset.
     */
    void reset()
    {
      R  = 0.0;
      p1 = 0.0;
      p2 = 0.0;
      p3 = 0.0;
      p4 = 0.0;
      bg = 0.0;
      cc = 0.0;
    }


    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 bg;                 //!< constant background [Hz]
    JParameter_t cc;                 //!< fraction of signal correlated background
  };


  /**
   * Fit parameters for two-fold coincidence rate due to K40.
   */
  struct JK40Parameters :
    JK40Parameters_t
  {
    /**
     * Default constructor.
     */
    JK40Parameters() 
    {}


    /**
     * Get default values.
     *
     * The default parameter values are set to those obtained from a designated simulation
     * of K40 decays (see http://wiki.km3net.de/index.php/OMGsim_simulations_for_K40_fit).\n
     * 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  = 18.460546;
      parameters.p1 =  3.0767;
      parameters.p2 = -1.2078;
      parameters.p3 =  0.9905;
      parameters.p4 =  0.9379;
      parameters.bg =  0.0;
      parameters.cc =  0.0;

      return parameters;
    }


    /**
     * Get K40 parameters.
     *
     * \return              K40 parameters
     */
    const JK40Parameters& getK40Parameters() const
    {
      return static_cast<const JK40Parameters&>(*this);
    }


    /**
     * Set K40 parameters.
     *
     * \param  parameters   K40 parameters
     */
    void setK40Parameters(const JK40Parameters_t& parameters)
    {
      static_cast<JK40Parameters_t&>(*this) = 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) +
              (bg.isFree() ? 1 : 0) +
              (cc.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::bg) { return N; } if (bg.isFree()) { ++N; }
      if (p == &JK40Parameters::cc) { return N; } if (cc.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 (bg.isFree()) { gradient.bg = 1.0; }
      if (cc.isFree()) { gradient.cc = rate; }

      return gradient;
    }

  private: 
    mutable JK40Parameters_t gradient;
  };


  /**
   * Fit model.
   *
   * In the absence of TDC constraints, the average time offset is fixed to zero.
   */
  struct JModel :
    public JK40Parameters,
    public JModule,
    public JCombinatorics_t
  {
    using JK40Parameters::getIndex;
    using JK40Parameters::getValue;
    using JCombinatorics_t::getIndex;


    /**
     * Auxiliary data structure for derived quantities of a given PMT pair.
     */
    struct real_type {
      double ct;                           //!< cosine angle between PMT axes
      double t0;                           //!< time offset [ns]
      double sigma;                        //!< total width [ns]
      double signal;                       //!< combined signal
    };


    /**
     * 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)
    {
      setK40Parameters(parameters);

      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) {
        this->parameters[i->second].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)
    {
      setK40Parameters(parameters);

      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)
    {
      if        (option == FIT_PMTS_t) {

        R .fix();
        p1.fix();
        p2.fix();
        p3.fix();
        p4.fix();

        parameters[index].t0 .fix();

      } else if (option == FIT_MODEL_t) {

        bg.fix();
        cc.fix();

        for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
          parameters[i].QE .fix();
          parameters[i].t0 .fix();
        }

      } else if (option == FIT_PMTS_AND_ANGULAR_DEPENDENCE_t) {

        R .fix();

      } else {

        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 = JK40Parameters::getN();

      for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
        N += parameters[i].getN();
      }

      return N;
    }


    /**
     * Get index of parameter.
     *
     * \param  pmt          PMT
     * \return              index
     */
    int getIndex(int pmt) const
    {
      int N = JK40Parameters::getN();

      for (int i = 0; i != pmt; ++i) {
        N += parameters[i].getN();
      }

      return N;
    }


    /**
     * Get index of parameter.
     *
     * \param  pmt          PMT
     * \param  p            pointer to data member
     * \return              index
     */
    int getIndex(int pmt, JParameter_t JPMTParameters_t::*p) const
    {
      return getIndex(pmt) + parameters[pmt].getIndex(p);
    }


    /**
     * 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 parameters.
     *
     * \param  pair         PMT pair
     * \return              parameters
     */
    const real_type& getReal(const pair_type& pair) const
    {
      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();

      return real;
    }


    /**
     * 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 JGauss gauss(real.t0, real.sigma, real.signal);

      const double R1  = this->getValue(real.ct);
      const double R2  = gauss.getValue(dt_ns);

      return bg() + R1 * (cc() + R2);
    }


    /**
     * Get K40 coincidence rate.
     *
     * \param  pair         PMT pair
     * \return              rate [Hz]
     */
    double getValue(const pair_type& pair) const
    {
      using namespace std;
      using namespace JPP;

      const real_type& real = getReal(pair);

      const double R1  = this->getValue(real.ct);
      const double R2  = real.signal;

      return bg() + R1 * (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 << "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 << "bg        " << FIXED(12,6) << object.bg << endl;
      out << "cc        " << FIXED(12,6) << object.cc << endl;

      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;
    }

    JPMTParameters_t  parameters[NUMBER_OF_PMTS];

  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;
    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);

      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();
        }
      }

      ndf -= value.getN();


      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*fabs(precessor)) {
              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) {
            return { precessor / estimator->getRho(1.0), ndf }; // no improvement found
          }

          evaluate(data);
        }

        if (debug >= debug_t) {

          size_t row = 0;

          if (value.R .isFree()) { cout << "R  " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.p1.isFree()) { cout << "p1 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.p2.isFree()) { cout << "p2 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.p3.isFree()) { cout << "p3 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.p4.isFree()) { cout << "p4 " << FIXED(12,5) << Y[row] << endl; ++row; }
          if (value.bg.isFree()) { cout << "bg " << FIXED(12,3) << Y[row] << endl; ++row; }
          if (value.cc.isFree()) { cout << "cc " << FIXED(12,3) << 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; }
          }
        }

        // 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.R .isFree()) { value.R  -= factor * h[row] * Y[row];  ++row; }
        if (value.p1.isFree()) { value.p1 -= factor * h[row] * Y[row];  ++row; }
        if (value.p2.isFree()) { value.p2 -= factor * h[row] * Y[row];  ++row; }
        if (value.p3.isFree()) { value.p3 -= factor * h[row] * Y[row];  ++row; }
        if (value.p4.isFree()) { value.p4 -= factor * h[row] * Y[row];  ++row; }
        if (value.bg.isFree()) { value.bg -= factor * h[row] * Y[row];  ++row; }
        if (value.cc.isFree()) { value.cc -= 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; }
        }
      }

      return { precessor / estimator->getRho(1.0), ndf };
    }
    

    static constexpr int    MAXIMUM_ITERATIONS =  10000;                                  //!< maximal number of iterations.
    static constexpr double EPSILON            = 1.0e-4;                                  //!< maximal distance to minimum.
    static constexpr double LAMBDA_MIN         =   0.01;                                  //!< minimal value control parameter
    static constexpr double LAMBDA_MAX         = 100.0;                                   //!< 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;
    int                         numberOfIterations;
    JMATH::JMatrixNS            V;

  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& model)
        {
          R  = model.getIndex(&JK40Parameters_t::R);
          p1 = model.getIndex(&JK40Parameters_t::p1);
          p2 = model.getIndex(&JK40Parameters_t::p2);
          p3 = model.getIndex(&JK40Parameters_t::p3);
          p4 = model.getIndex(&JK40Parameters_t::p4);
          bg = model.getIndex(&JK40Parameters_t::bg);
          cc = model.getIndex(&JK40Parameters_t::cc);
        }

        int R;
        int p1;
        int p2;
        int p3;
        int p4;
        int bg;
        int cc;

      } M(value);


      // PMT parameter indices

      struct I_t {
        I_t(const JModel& model, const int pmt) :
          QE (INVALID_INDEX),
          TTS(INVALID_INDEX),
          t0 (INVALID_INDEX)
        {
          const int index = model.getIndex(pmt);

          int N = 0;

          if (model.parameters[pmt].QE .isFree()) { QE  = index + N; ++N; }
          if (model.parameters[pmt].TTS.isFree()) { TTS = index + N; ++N; }
          if (model.parameters[pmt].t0 .isFree()) { t0  = index + N; ++N; }
        }

        int QE;
        int TTS;
        int t0;
      };


      typedef vector< pair<int, double> >  buffer_type;
 
      buffer_type  buffer;

      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 JGauss gauss(real.t0, real.sigma, real.signal);

          const double            R1  = value.getValue   (real.ct);
          const JK40Parameters_t& R1p = value.getGradient(real.ct);

          const std::pair<I_t, I_t> PMT(I_t(value, pair.first),
                                        I_t(value, pair.second));

          for (const rate_type& iy : ix->second) {

            const double  R2  = gauss.getValue   (iy.dt_ns);
            const JGauss& R2p = gauss.getGradient(iy.dt_ns);        

            const double  R   = value.bg()  +  R1 * (value.cc() + R2);
            const double  u   =  (iy.value - R)       / iy.error;
            const double  w   = -estimator->getPsi(u) / iy.error;

            successor += estimator->getRho(u);

            buffer.clear();

            if (M.R  != INVALID_INDEX) { buffer.push_back({M.R,  w * (value.cc() + R2) * R1p.R () * value.R .getDerivative()}); }
            if (M.p1 != INVALID_INDEX) { buffer.push_back({M.p1, w * (value.cc() + R2) * R1p.p1() * value.p1.getDerivative()}); }
            if (M.p2 != INVALID_INDEX) { buffer.push_back({M.p2, w * (value.cc() + R2) * R1p.p2() * value.p2.getDerivative()}); }
            if (M.p3 != INVALID_INDEX) { buffer.push_back({M.p3, w * (value.cc() + R2) * R1p.p3() * value.p3.getDerivative()}); }
            if (M.p4 != INVALID_INDEX) { buffer.push_back({M.p4, w * (value.cc() + R2) * R1p.p4() * value.p4.getDerivative()}); }
            if (M.bg != INVALID_INDEX) { buffer.push_back({M.bg, w *  value.bg.getDerivative()}); }
            if (M.cc != INVALID_INDEX) { buffer.push_back({M.cc, w *  R1 * R1p.cc() * value.cc.getDerivative()}); }

            if (PMT.first .QE  != INVALID_INDEX) { buffer.push_back({PMT.first .QE , w * R1 * R2p.signal * value.parameters[pair.second].QE () * value.parameters[pair.first ].QE .getDerivative()}); }
            if (PMT.second.QE  != INVALID_INDEX) { buffer.push_back({PMT.second.QE , w * R1 * R2p.signal * value.parameters[pair.first ].QE () * value.parameters[pair.second].QE .getDerivative()}); }
            if (PMT.first .TTS != INVALID_INDEX) { buffer.push_back({PMT.first .TTS, w * R1 * R2p.sigma  * value.parameters[pair.first ].TTS() * value.parameters[pair.first ].TTS.getDerivative() / real.sigma}); }
            if (PMT.second.TTS != INVALID_INDEX) { buffer.push_back({PMT.second.TTS, w * R1 * R2p.sigma  * value.parameters[pair.second].TTS() * value.parameters[pair.second].TTS.getDerivative() / real.sigma}); }
            if (PMT.first .t0  != INVALID_INDEX) { buffer.push_back({PMT.first .t0,  w * R1 * R2p.mean   * value.parameters[pair.first ].t0 .getDerivative() * +1.0}); }
            if (PMT.second.t0  != INVALID_INDEX) { buffer.push_back({PMT.second.t0,  w * R1 * R2p.mean   * value.parameters[pair.second].t0 .getDerivative() * -1.0}); }

            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];
              }
            }
          }
        }
      }
    }

    JMATH::JVectorND            Y;                             // gradient
    double                      successor;
    JModel                      previous;
    std::vector<double>         h;                             // normalisation vector
  };
}

#endif