Jpp  19.1.0
the software that should make you happy
Classes | Public Types | Public Member Functions | Public Attributes | Static Public Attributes | Private Member Functions | Private Attributes | List of all members
JCALIBRATE::JFit Class Reference

Fit. More...

#include <JFitK40.hh>

Classes

struct  result_type
 Result type. More...
 

Public Types

typedef std::shared_ptr< JMEstimatorestimator_type
 

Public Member Functions

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

Public Attributes

int debug
 
estimator_type estimator
 M-Estimator function. More...
 
double lambda
 
JModel value
 
int numberOfIterations
 
JMATH::JMatrixNS V
 

Static Public Attributes

static constexpr int MAXIMUM_ITERATIONS = 10000
 maximal number of iterations. More...
 
static constexpr double EPSILON = 1.0e-5
 maximal distance to minimum. More...
 
static constexpr double LAMBDA_MIN = 0.01
 minimal value control parameter More...
 
static constexpr double LAMBDA_MAX = 100.0
 maximal value control parameter More...
 
static constexpr double LAMBDA_UP = 10.0
 multiplication factor control parameter More...
 
static constexpr double LAMBDA_DOWN = 10.0
 multiplication factor control parameter More...
 
static constexpr double PIVOT = std::numeric_limits<double>::epsilon()
 minimal value diagonal element of matrix More...
 

Private Member Functions

void evaluate (const data_type &data)
 Evaluation of fit. More...
 

Private Attributes

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

Detailed Description

Fit.

Definition at line 1176 of file JFitK40.hh.

Member Typedef Documentation

◆ estimator_type

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

Definition at line 1187 of file JFitK40.hh.

Constructor & Destructor Documentation

◆ JFit()

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

Constructor.

Parameters
optionM-estimator
debugdebug

Definition at line 1196 of file JFitK40.hh.

1196  :
1197  debug(debug)
1198  {
1199  using namespace JPP;
1200 
1201  estimator.reset(getMEstimator(option));
1202  }
estimator_type estimator
M-Estimator function.
Definition: JFitK40.hh:1375
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:203
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).

Member Function Documentation

◆ operator()()

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

Fit.

Parameters
datadata
Returns
chi2, NDF

Definition at line 1211 of file JFitK40.hh.

1212  {
1213  using namespace std;
1214  using namespace JPP;
1215 
1216 
1217  value.setIndex();
1218 
1219  const size_t N = value.getN();
1220 
1221  V.resize(N);
1222  Y.resize(N);
1223  h.resize(N);
1224 
1225  int ndf = 0;
1226 
1227  for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1228 
1229  const pair_type& pair = ix->first;
1230 
1231  if (value.parameters[pair.first ].status &&
1232  value.parameters[pair.second].status) {
1233 
1234  ndf += ix->second.size();
1235  }
1236  }
1237 
1238  ndf -= value.getN();
1239 
1240 
1241  lambda = LAMBDA_MIN;
1242 
1243  double precessor = numeric_limits<double>::max();
1244 
1246 
1247  DEBUG("step: " << numberOfIterations << endl);
1248 
1249  evaluate(data);
1250 
1251  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
1252  DEBUG("chi2: " << FIXED(12,3) << successor << endl);
1253 
1254  if (successor < precessor) {
1255 
1256  if (numberOfIterations != 0) {
1257 
1258  if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
1259  return { successor / estimator->getRho(1.0), ndf };
1260  }
1261 
1262  if (lambda > LAMBDA_MIN) {
1263  lambda /= LAMBDA_DOWN;
1264  }
1265  }
1266 
1267  precessor = successor;
1268  previous = value;
1269 
1270  } else {
1271 
1272  value = previous;
1273  lambda *= LAMBDA_UP;
1274 
1275  if (lambda > LAMBDA_MAX) {
1276  return { precessor / estimator->getRho(1.0), ndf }; // no improvement found
1277  }
1278 
1279  evaluate(data);
1280  }
1281 
1282  if (debug >= debug_t) {
1283 
1284  size_t row = 0;
1285 
1286  if (value.R .isFree()) { cout << "R " << FIXED(12,5) << Y[row] << endl; ++row; }
1287  if (value.p1.isFree()) { cout << "p1 " << FIXED(12,5) << Y[row] << endl; ++row; }
1288  if (value.p2.isFree()) { cout << "p2 " << FIXED(12,5) << Y[row] << endl; ++row; }
1289  if (value.p3.isFree()) { cout << "p3 " << FIXED(12,5) << Y[row] << endl; ++row; }
1290  if (value.p4.isFree()) { cout << "p4 " << FIXED(12,5) << Y[row] << endl; ++row; }
1291  if (value.cc.isFree()) { cout << "cc " << FIXED(12,3) << Y[row] << endl; ++row; }
1292 
1293  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1294  if (value.parameters[pmt].QE .isFree()) { cout << "PMT[" << setw(2) << pmt << "].QE " << FIXED(12,5) << Y[row] << endl; ++row; }
1295  if (value.parameters[pmt].TTS.isFree()) { cout << "PMT[" << setw(2) << pmt << "].TTS " << FIXED(12,5) << Y[row] << endl; ++row; }
1296  if (value.parameters[pmt].t0 .isFree()) { cout << "PMT[" << setw(2) << pmt << "].t0 " << FIXED(12,5) << Y[row] << endl; ++row; }
1297  if (value.parameters[pmt].bg .isFree()) { cout << "PMT[" << setw(2) << pmt << "].bg " << FIXED(12,5) << Y[row] << endl; ++row; }
1298  }
1299  }
1300 
1301  // force definite positiveness
1302 
1303  for (size_t i = 0; i != N; ++i) {
1304 
1305  if (V(i,i) < PIVOT) {
1306  V(i,i) = PIVOT;
1307  }
1308 
1309  h[i] = 1.0 / sqrt(V(i,i));
1310  }
1311 
1312  // normalisation
1313 
1314  for (size_t i = 0; i != N; ++i) {
1315  for (size_t j = 0; j != i; ++j) {
1316  V(j,i) *= h[i] * h[j];
1317  V(i,j) = V(j,i);
1318  }
1319  }
1320 
1321  for (size_t i = 0; i != N; ++i) {
1322  V(i,i) = 1.0 + lambda;
1323  }
1324 
1325  // solve A x = b
1326 
1327  for (size_t col = 0; col != N; ++col) {
1328  Y[col] *= h[col];
1329  }
1330 
1331  try {
1332  V.solve(Y);
1333  }
1334  catch (const exception& error) {
1335 
1336  ERROR("JGandalf: " << error.what() << endl << V << endl);
1337 
1338  break;
1339  }
1340 
1341  // update value
1342 
1343  const double factor = 2.0;
1344 
1345  size_t row = 0;
1346 
1347  if (value.R .isFree()) { value.R -= factor * h[row] * Y[row]; ++row; }
1348  if (value.p1.isFree()) { value.p1 -= factor * h[row] * Y[row]; ++row; }
1349  if (value.p2.isFree()) { value.p2 -= factor * h[row] * Y[row]; ++row; }
1350  if (value.p3.isFree()) { value.p3 -= factor * h[row] * Y[row]; ++row; }
1351  if (value.p4.isFree()) { value.p4 -= factor * h[row] * Y[row]; ++row; }
1352  if (value.cc.isFree()) { value.cc -= factor * h[row] * Y[row]; ++row; }
1353 
1354  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1355  if (value.parameters[pmt].QE .isFree()) { value.parameters[pmt].QE -= factor * h[row] * Y[row]; ++row; }
1356  if (value.parameters[pmt].TTS.isFree()) { value.parameters[pmt].TTS -= factor * h[row] * Y[row]; ++row; }
1357  if (value.parameters[pmt].t0 .isFree()) { value.parameters[pmt].t0 -= factor * h[row] * Y[row]; ++row; }
1358  if (value.parameters[pmt].bg .isFree()) { value.parameters[pmt].bg -= factor * h[row] * Y[row]; ++row; }
1359  }
1360  }
1361 
1362  return { precessor / estimator->getRho(1.0), ndf };
1363  }
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
#define ERROR(A)
Definition: JMessage.hh:66
std::vector< double > h
Definition: JFitK40.hh:1520
static constexpr double LAMBDA_MIN
minimal value control parameter
Definition: JFitK40.hh:1368
static constexpr double LAMBDA_DOWN
multiplication factor control parameter
Definition: JFitK40.hh:1371
static constexpr double LAMBDA_MAX
maximal value control parameter
Definition: JFitK40.hh:1369
static constexpr double LAMBDA_UP
multiplication factor control parameter
Definition: JFitK40.hh:1370
int numberOfIterations
Definition: JFitK40.hh:1379
JMATH::JMatrixNS V
Definition: JFitK40.hh:1380
static constexpr double EPSILON
maximal distance to minimum.
Definition: JFitK40.hh:1367
void evaluate(const data_type &data)
Evaluation of fit.
Definition: JFitK40.hh:1388
static constexpr int MAXIMUM_ITERATIONS
maximal number of iterations.
Definition: JFitK40.hh:1366
static constexpr double PIVOT
minimal value diagonal element of matrix
Definition: JFitK40.hh:1372
double successor
Definition: JFitK40.hh:1518
JMATH::JVectorND Y
Definition: JFitK40.hh:1517
bool isFree() const
Check if parameter is free.
Definition: JFitK40.hh:240
@ debug_t
debug
Definition: JMessage.hh:29
int j
Definition: JPolint.hh:792
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
Definition: JSTDTypes.hh:14
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
JParameter_t R
maximal coincidence rate [Hz]
Definition: JFitK40.hh:596
JParameter_t p1
1st order angle dependence coincidence rate
Definition: JFitK40.hh:597
JParameter_t p2
2nd order angle dependence coincidence rate
Definition: JFitK40.hh:598
JParameter_t p3
3rd order angle dependence coincidence rate
Definition: JFitK40.hh:599
JParameter_t p4
4th order angle dependence coincidence rate
Definition: JFitK40.hh:600
JParameter_t cc
fraction of signal correlated background
Definition: JFitK40.hh:601
size_t getN() const
Get number of fit parameters.
Definition: JFitK40.hh:1001
JPMTParameters_t parameters[NUMBER_OF_PMTS]
Definition: JFitK40.hh:1163
void setIndex()
Set index of PMT used for fixed time offset.
Definition: JFitK40.hh:975
JParameter_t t0
time offset [ns]
Definition: JFitK40.hh:564
JParameter_t TTS
transition-time spread [ns]
Definition: JFitK40.hh:563
JParameter_t bg
background [Hz/ns]
Definition: JFitK40.hh:565
JParameter_t QE
relative quantum efficiency [unit]
Definition: JFitK40.hh:562
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:446
void solve(JVectorND_t &u)
Get solution of equation A x = b.
Definition: JMatrixNS.hh:308
Data structure for a pair of indices.

◆ evaluate()

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

Evaluation of fit.

Parameters
datadata

Definition at line 1388 of file JFitK40.hh.

1389  {
1390  using namespace std;
1391  using namespace JPP;
1392 
1393  typedef JModel::real_type real_type;
1394 
1395 
1396  successor = 0.0;
1397 
1398  V.reset();
1399  Y.reset();
1400 
1401 
1402  // model parameter indices
1403 
1404  const struct M_t {
1405  M_t(const JModel& model)
1406  {
1407  R = model.getIndex(&JK40Parameters_t::R);
1408  p1 = model.getIndex(&JK40Parameters_t::p1);
1409  p2 = model.getIndex(&JK40Parameters_t::p2);
1410  p3 = model.getIndex(&JK40Parameters_t::p3);
1411  p4 = model.getIndex(&JK40Parameters_t::p4);
1412  cc = model.getIndex(&JK40Parameters_t::cc);
1413  }
1414 
1415  int R;
1416  int p1;
1417  int p2;
1418  int p3;
1419  int p4;
1420  int cc;
1421 
1422  } M(value);
1423 
1424 
1425  // PMT parameter indices
1426 
1427  struct I_t {
1428  I_t(const JModel& model, const int pmt) :
1429  QE (INVALID_INDEX),
1430  TTS(INVALID_INDEX),
1431  t0 (INVALID_INDEX),
1432  bg (INVALID_INDEX)
1433  {
1434  const int index = model.getIndex(pmt);
1435 
1436  int N = 0;
1437 
1438  if (model.parameters[pmt].QE .isFree()) { QE = index + N; ++N; }
1439  if (model.parameters[pmt].TTS.isFree()) { TTS = index + N; ++N; }
1440  if (model.parameters[pmt].t0 .isFree()) { t0 = index + N; ++N; }
1441  if (model.parameters[pmt].bg .isFree()) { bg = index + N; ++N; }
1442  }
1443 
1444  int QE;
1445  int TTS;
1446  int t0;
1447  int bg;
1448  };
1449 
1450 
1452 
1453  buffer_type buffer;
1454 
1455  for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1456 
1457  const pair_type& pair = ix->first;
1458 
1459  if (value.parameters[pair.first ].status &&
1460  value.parameters[pair.second].status) {
1461 
1462  const real_type& real = value.getReal(pair);
1463 
1464  const JGauss gauss(real.t0, real.sigma, real.signal);
1465 
1466  const double R1 = value.getValue (real.ct);
1467  const JK40Parameters_t& R1p = value.getGradient(real.ct);
1468 
1469  const std::pair<I_t, I_t> PMT(I_t(value, pair.first),
1470  I_t(value, pair.second));
1471 
1472  for (const rate_type& iy : ix->second) {
1473 
1474  const double R2 = gauss.getValue (iy.dt_ns);
1475  const JGauss& R2p = gauss.getGradient(iy.dt_ns);
1476 
1477  const double R = real.background + R1 * (value.cc() + R2);
1478  const double u = (iy.value - R) / iy.error;
1479  const double w = -estimator->getPsi(u) / iy.error;
1480 
1481  successor += estimator->getRho(u);
1482 
1483  buffer.clear();
1484 
1485  if (M.R != INVALID_INDEX) { buffer.push_back({M.R, w * (value.cc() + R2) * R1p.R () * value.R .getDerivative()}); }
1486  if (M.p1 != INVALID_INDEX) { buffer.push_back({M.p1, w * (value.cc() + R2) * R1p.p1() * value.p1.getDerivative()}); }
1487  if (M.p2 != INVALID_INDEX) { buffer.push_back({M.p2, w * (value.cc() + R2) * R1p.p2() * value.p2.getDerivative()}); }
1488  if (M.p3 != INVALID_INDEX) { buffer.push_back({M.p3, w * (value.cc() + R2) * R1p.p3() * value.p3.getDerivative()}); }
1489  if (M.p4 != INVALID_INDEX) { buffer.push_back({M.p4, w * (value.cc() + R2) * R1p.p4() * value.p4.getDerivative()}); }
1490  if (M.cc != INVALID_INDEX) { buffer.push_back({M.cc, w * R1 * R1p.cc() * value.cc.getDerivative()}); }
1491 
1492  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()}); }
1493  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()}); }
1494  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}); }
1495  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}); }
1496  if (PMT.first .t0 != INVALID_INDEX) { buffer.push_back({PMT.first .t0, w * R1 * R2p.mean * value.parameters[pair.first ].t0 .getDerivative() * +1.0}); }
1497  if (PMT.second.t0 != INVALID_INDEX) { buffer.push_back({PMT.second.t0, w * R1 * R2p.mean * value.parameters[pair.second].t0 .getDerivative() * -1.0}); }
1498  if (PMT.first .bg != INVALID_INDEX) { buffer.push_back({PMT.first .bg, w * value.parameters[pair.first ].bg .getDerivative()}); }
1499  if (PMT.second.bg != INVALID_INDEX) { buffer.push_back({PMT.second.bg, w * value.parameters[pair.second].bg .getDerivative()}); }
1500 
1501  for (buffer_type::const_iterator row = buffer.begin(); row != buffer.end(); ++row) {
1502 
1503  Y[row->first] += row->second;
1504 
1505  V[row->first][row->first] += row->second * row->second;
1506 
1507  for (buffer_type::const_iterator col = buffer.begin(); col != row; ++col) {
1508  V[row->first][col->first] += row->second * col->second;
1509  V[col->first][row->first] = V[row->first][col->first];
1510  }
1511  }
1512  }
1513  }
1514  }
1515  }
JDAQPMTIdentifier PMT
Command line options.
TPaveText * p1
double getDerivative() const
Get derivative of value.
Definition: JFitK40.hh:334
#define R1(x)
static const int INVALID_INDEX
invalid index
Definition: JFitK40.hh:60
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
std::vector< JHitW0 > buffer_type
hits
Definition: JPerth.cc:71
data_type w[N+1][M+1]
Definition: JPolint.hh:867
double u[N+1]
Definition: JPolint.hh:865
Model for fit to acoustics data.
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitK40.hh:572
const JK40Parameters_t & getGradient(const double ct) const
Get gradient.
Definition: JFitK40.hh:723
Auxiliary data structure for derived quantities of a given PMT pair.
Definition: JFitK40.hh:763
const real_type & getReal(const pair_type &pair) const
Get derived parameters.
Definition: JFitK40.hh:1072
double getValue(const pair_type &pair, const double dt_ns) const
Get K40 coincidence rate.
Definition: JFitK40.hh:1099
Data structure for measured coincidence rate of pair of PMTs.
Definition: JFitK40.hh:66
double error
error of rate [Hz/ns]
Definition: JFitK40.hh:94
double value
value of rate [Hz/ns]
Definition: JFitK40.hh:93
double dt_ns
time difference [ns]
Definition: JFitK40.hh:92
double signal
Definition: JGauss.hh:163
double mean
Definition: JGauss.hh:161
Gauss function object.
Definition: JGauss.hh:175
double sigma
sigma
Definition: JMathlib.hh:1665
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:459
void reset()
Reset.
Definition: JVectorND.hh:45

Member Data Documentation

◆ MAXIMUM_ITERATIONS

constexpr int JCALIBRATE::JFit::MAXIMUM_ITERATIONS = 10000
staticconstexpr

maximal number of iterations.

Definition at line 1366 of file JFitK40.hh.

◆ EPSILON

constexpr double JCALIBRATE::JFit::EPSILON = 1.0e-5
staticconstexpr

maximal distance to minimum.

Definition at line 1367 of file JFitK40.hh.

◆ LAMBDA_MIN

constexpr double JCALIBRATE::JFit::LAMBDA_MIN = 0.01
staticconstexpr

minimal value control parameter

Definition at line 1368 of file JFitK40.hh.

◆ LAMBDA_MAX

constexpr double JCALIBRATE::JFit::LAMBDA_MAX = 100.0
staticconstexpr

maximal value control parameter

Definition at line 1369 of file JFitK40.hh.

◆ LAMBDA_UP

constexpr double JCALIBRATE::JFit::LAMBDA_UP = 10.0
staticconstexpr

multiplication factor control parameter

Definition at line 1370 of file JFitK40.hh.

◆ LAMBDA_DOWN

constexpr double JCALIBRATE::JFit::LAMBDA_DOWN = 10.0
staticconstexpr

multiplication factor control parameter

Definition at line 1371 of file JFitK40.hh.

◆ PIVOT

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

minimal value diagonal element of matrix

Definition at line 1372 of file JFitK40.hh.

◆ debug

int JCALIBRATE::JFit::debug

Definition at line 1374 of file JFitK40.hh.

◆ estimator

estimator_type JCALIBRATE::JFit::estimator

M-Estimator function.

Definition at line 1375 of file JFitK40.hh.

◆ lambda

double JCALIBRATE::JFit::lambda

Definition at line 1377 of file JFitK40.hh.

◆ value

JModel JCALIBRATE::JFit::value

Definition at line 1378 of file JFitK40.hh.

◆ numberOfIterations

int JCALIBRATE::JFit::numberOfIterations

Definition at line 1379 of file JFitK40.hh.

◆ V

JMATH::JMatrixNS JCALIBRATE::JFit::V

Definition at line 1380 of file JFitK40.hh.

◆ Y

JMATH::JVectorND JCALIBRATE::JFit::Y
private

Definition at line 1517 of file JFitK40.hh.

◆ successor

double JCALIBRATE::JFit::successor
private

Definition at line 1518 of file JFitK40.hh.

◆ previous

JModel JCALIBRATE::JFit::previous
private

Definition at line 1519 of file JFitK40.hh.

◆ h

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

Definition at line 1520 of file JFitK40.hh.


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