JCALIBRATE::JFit Class Reference

Fit. More...

#include <JFitK40.hh>

Classes
struct	result_type
	Result type. More...

Public Types
typedef std::shared_ptr< JMEstimator >	estimator_type

Public Member Functions
	JFit (const int option, const int debug)
	Constructor.
result_type	operator() (const data_type &data)
	Fit.

Public Attributes
int	debug
estimator_type	estimator
	M-Estimator function.
double	lambda
JModel	value
JModel_t	error
int	numberOfIterations
JMATH::JMatrixNS	V

Static Public Attributes
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

Private Member Functions
void	evaluate (const data_type &data)
	Evaluation of fit.
void	seterr (const data_type &data)
	Set errors.

Private Attributes
JMATH::JVectorND	Y
double	successor
JModel	previous
std::vector< double >	h

Detailed Description

Fit.

Definition at line 1290 of file JFitK40.hh.

Member Typedef Documentation

estimator_type

std::shared_ptr<JMEstimator> JCALIBRATE::JFit::estimator_type

Definition at line 1301 of file JFitK40.hh.

Constructor & Destructor Documentation

JFit()

JCALIBRATE::JFit::JFit ( const int option, const int debug )

inline

Constructor.

Parameters

option	M-estimator
debug	debug

Definition at line 1310 of file JFitK40.hh.

                                            :
      debug(debug)
    {
      using namespace JPP;
 
      estimator.reset(getMEstimator(option));
    }

Member Function Documentation

operator()()

result_type JCALIBRATE::JFit::operator() ( const data_type & data )

inline

Fit.

Parameters

data

data

Returns

chi2, NDF

Definition at line 1325 of file JFitK40.hh.

    {
      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.setRange(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.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.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; }
            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.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.cc.isFree()) { value.cc -= factor * h[row] * Y[row];  ++row; }
        if (value.bc.isFree()) { value.bc -= 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 };
    }

evaluate()

void JCALIBRATE::JFit::evaluate ( const data_type & data )

inlineprivate

Evaluation of fit.

Parameters

data

data

Definition at line 1527 of file JFitK40.hh.

    {
      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);
          cc = model.getIndex(&JK40Parameters_t::cc);
          bc = model.getIndex(&JK40Parameters_t::bc);
        }
 
        int R;
        int p1;
        int p2;
        int p3;
        int p4;
        int cc;
        int bc;
 
      } 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),
          bg (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; }
          if (model.parameters[pmt].bg .isFree()) { bg  = index + N; ++N; }
        }
 
        int QE;
        int TTS;
        int t0;
        int bg;
      };
 
 
      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 JBell bell(real.t0, real.sigma, real.signal, 0.0, BELL_SHAPE);
 
          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  = bell.getValue   (iy.dt_ns);
            const JBell_t& R2p = bell.getGradient(iy.dt_ns);        
 
            const double  R   = real.bc  +  real.background  +  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();
 
            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.cc != INVALID_INDEX) { buffer.push_back({M.cc, w *  R1 * R1p.cc() * value.cc.getDerivative()}); }
            if (M.bc != INVALID_INDEX) { buffer.push_back({M.bc, w *       R1p.bc() * value.bc.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}); }
            if (PMT.first .bg  != INVALID_INDEX) { buffer.push_back({PMT.first .bg,  w * value.parameters[pair.first ].bg .getDerivative()}); }
            if (PMT.second.bg  != INVALID_INDEX) { buffer.push_back({PMT.second.bg,  w * value.parameters[pair.second].bg .getDerivative()}); }
                
            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];
              }
            }
          }
        }
      }
    }

seterr()

void JCALIBRATE::JFit::seterr ( const data_type & data )

inlineprivate

Set errors.

Parameters

data

data

Definition at line 1665 of file JFitK40.hh.

    {
      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 i = 0;
 
      if (value.R .isFree()) { error.R  = SQRT(V(i,i));  ++i; }
      if (value.p1.isFree()) { error.p1 = SQRT(V(i,i));  ++i; }
      if (value.p2.isFree()) { error.p2 = SQRT(V(i,i));  ++i; }
      if (value.p3.isFree()) { error.p3 = SQRT(V(i,i));  ++i; }
      if (value.p4.isFree()) { error.p4 = SQRT(V(i,i));  ++i; }
      if (value.cc.isFree()) { error.cc = SQRT(V(i,i));  ++i; }
      if (value.bc.isFree()) { error.bc = SQRT(V(i,i));  ++i; }
 
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        if (value.parameters[pmt].QE .isFree()) { error.parameters[pmt].QE  = SQRT(V(i,i));  ++i; }
        if (value.parameters[pmt].TTS.isFree()) { error.parameters[pmt].TTS = SQRT(V(i,i));  ++i; }
        if (value.parameters[pmt].t0 .isFree()) { error.parameters[pmt].t0  = SQRT(V(i,i));  ++i; }
        if (value.parameters[pmt].bg .isFree()) { error.parameters[pmt].bg  = SQRT(V(i,i));  ++i; }
      }
 
#undef SQRT
    }

Member Data Documentation

MAXIMUM_ITERATIONS

int JCALIBRATE::JFit::MAXIMUM_ITERATIONS = 100000

staticconstexpr

maximal number of iterations.

Definition at line 1504 of file JFitK40.hh.

EPSILON

double JCALIBRATE::JFit::EPSILON = 1.0e-3

staticconstexpr

maximal distance to minimum.

Definition at line 1505 of file JFitK40.hh.

LAMBDA_MIN

double JCALIBRATE::JFit::LAMBDA_MIN = 1.0e-2

staticconstexpr

minimal value control parameter

Definition at line 1506 of file JFitK40.hh.

LAMBDA_MAX

double JCALIBRATE::JFit::LAMBDA_MAX = 1.0e+4

staticconstexpr

maximal value control parameter

Definition at line 1507 of file JFitK40.hh.

LAMBDA_UP

double JCALIBRATE::JFit::LAMBDA_UP = 10.0

staticconstexpr

multiplication factor control parameter

Definition at line 1508 of file JFitK40.hh.

LAMBDA_DOWN

double JCALIBRATE::JFit::LAMBDA_DOWN = 10.0

staticconstexpr

multiplication factor control parameter

Definition at line 1509 of file JFitK40.hh.

PIVOT

double JCALIBRATE::JFit::PIVOT = std::numeric_limits<double>::epsilon()

staticconstexpr

minimal value diagonal element of matrix

Definition at line 1510 of file JFitK40.hh.

debug

int JCALIBRATE::JFit::debug

Definition at line 1512 of file JFitK40.hh.

estimator

estimator_type JCALIBRATE::JFit::estimator

M-Estimator function.

Definition at line 1513 of file JFitK40.hh.

lambda

double JCALIBRATE::JFit::lambda

Definition at line 1515 of file JFitK40.hh.

value

JModel JCALIBRATE::JFit::value

Definition at line 1516 of file JFitK40.hh.

error

JModel_t JCALIBRATE::JFit::error

Definition at line 1517 of file JFitK40.hh.

numberOfIterations

int JCALIBRATE::JFit::numberOfIterations

Definition at line 1518 of file JFitK40.hh.

V

JMATH::JMatrixNS JCALIBRATE::JFit::V

Definition at line 1519 of file JFitK40.hh.

Y

JMATH::JVectorND JCALIBRATE::JFit::Y

private

Definition at line 1701 of file JFitK40.hh.

successor

double JCALIBRATE::JFit::successor

private

Definition at line 1702 of file JFitK40.hh.

previous

JModel JCALIBRATE::JFit::previous

private

Definition at line 1703 of file JFitK40.hh.

h

std::vector<double> JCALIBRATE::JFit::h

private

Definition at line 1704 of file JFitK40.hh.

---

The documentation for this class was generated from the following file:

• JFitK40.hh