Jpp  18.2.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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
< JMEstimator
estimator_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-4
 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 1172 of file JFitK40.hh.

Member Typedef Documentation

Definition at line 1183 of file JFitK40.hh.

Constructor & Destructor Documentation

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

Constructor.

Parameters
optionM-estimator
debugdebug

Definition at line 1192 of file JFitK40.hh.

1192  :
1193  debug(debug)
1194  {
1195  using namespace JPP;
1196 
1197  estimator.reset(getMEstimator(option));
1198  }
estimator_type estimator
M-Estimator function.
Definition: JFitK40.hh:1371
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:203

Member Function Documentation

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

Fit.

Parameters
datadata
Returns
chi2, NDF

Definition at line 1207 of file JFitK40.hh.

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

Evaluation of fit.

Parameters
datadata

Definition at line 1384 of file JFitK40.hh.

1385  {
1386  using namespace std;
1387  using namespace JPP;
1388 
1389  typedef JModel::real_type real_type;
1390 
1391 
1392  successor = 0.0;
1393 
1394  V.reset();
1395  Y.reset();
1396 
1397 
1398  // model parameter indices
1399 
1400  const struct M_t {
1401  M_t(const JModel& model)
1402  {
1403  R = model.getIndex(&JK40Parameters_t::R);
1404  p1 = model.getIndex(&JK40Parameters_t::p1);
1405  p2 = model.getIndex(&JK40Parameters_t::p2);
1406  p3 = model.getIndex(&JK40Parameters_t::p3);
1407  p4 = model.getIndex(&JK40Parameters_t::p4);
1408  cc = model.getIndex(&JK40Parameters_t::cc);
1409  }
1410 
1411  int R;
1412  int p1;
1413  int p2;
1414  int p3;
1415  int p4;
1416  int cc;
1417 
1418  } M(value);
1419 
1420 
1421  // PMT parameter indices
1422 
1423  struct I_t {
1424  I_t(const JModel& model, const int pmt) :
1425  QE (INVALID_INDEX),
1426  TTS(INVALID_INDEX),
1427  t0 (INVALID_INDEX),
1428  bg (INVALID_INDEX)
1429  {
1430  const int index = model.getIndex(pmt);
1431 
1432  int N = 0;
1433 
1434  if (model.parameters[pmt].QE .isFree()) { QE = index + N; ++N; }
1435  if (model.parameters[pmt].TTS.isFree()) { TTS = index + N; ++N; }
1436  if (model.parameters[pmt].t0 .isFree()) { t0 = index + N; ++N; }
1437  if (model.parameters[pmt].bg .isFree()) { bg = index + N; ++N; }
1438  }
1439 
1440  int QE;
1441  int TTS;
1442  int t0;
1443  int bg;
1444  };
1445 
1446 
1448 
1449  buffer_type buffer;
1450 
1451  for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1452 
1453  const pair_type& pair = ix->first;
1454 
1455  if (value.parameters[pair.first ].status &&
1456  value.parameters[pair.second].status) {
1457 
1458  const real_type& real = value.getReal(pair);
1459 
1460  const JGauss gauss(real.t0, real.sigma, real.signal);
1461 
1462  const double R1 = value.getValue (real.ct);
1463  const JK40Parameters_t& R1p = value.getGradient(real.ct);
1464 
1465  const std::pair<I_t, I_t> PMT(I_t(value, pair.first),
1466  I_t(value, pair.second));
1467 
1468  for (const rate_type& iy : ix->second) {
1469 
1470  const double R2 = gauss.getValue (iy.dt_ns);
1471  const JGauss& R2p = gauss.getGradient(iy.dt_ns);
1472 
1473  const double R = real.background + R1 * (value.cc() + R2);
1474  const double u = (iy.value - R) / iy.error;
1475  const double w = -estimator->getPsi(u) / iy.error;
1476 
1477  successor += estimator->getRho(u);
1478 
1479  buffer.clear();
1480 
1481  if (M.R != INVALID_INDEX) { buffer.push_back({M.R, w * (value.cc() + R2) * R1p.R () * value.R .getDerivative()}); }
1482  if (M.p1 != INVALID_INDEX) { buffer.push_back({M.p1, w * (value.cc() + R2) * R1p.p1() * value.p1.getDerivative()}); }
1483  if (M.p2 != INVALID_INDEX) { buffer.push_back({M.p2, w * (value.cc() + R2) * R1p.p2() * value.p2.getDerivative()}); }
1484  if (M.p3 != INVALID_INDEX) { buffer.push_back({M.p3, w * (value.cc() + R2) * R1p.p3() * value.p3.getDerivative()}); }
1485  if (M.p4 != INVALID_INDEX) { buffer.push_back({M.p4, w * (value.cc() + R2) * R1p.p4() * value.p4.getDerivative()}); }
1486  if (M.cc != INVALID_INDEX) { buffer.push_back({M.cc, w * R1 * R1p.cc() * value.cc.getDerivative()}); }
1487 
1488  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()}); }
1489  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()}); }
1490  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}); }
1491  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}); }
1492  if (PMT.first .t0 != INVALID_INDEX) { buffer.push_back({PMT.first .t0, w * R1 * R2p.mean * value.parameters[pair.first ].t0 .getDerivative() * +1.0}); }
1493  if (PMT.second.t0 != INVALID_INDEX) { buffer.push_back({PMT.second.t0, w * R1 * R2p.mean * value.parameters[pair.second].t0 .getDerivative() * -1.0}); }
1494  if (PMT.first .bg != INVALID_INDEX) { buffer.push_back({PMT.first .bg, w * value.parameters[pair.first ].bg .getDerivative()}); }
1495  if (PMT.second.bg != INVALID_INDEX) { buffer.push_back({PMT.second.bg, w * value.parameters[pair.second].bg .getDerivative()}); }
1496 
1497  for (buffer_type::const_iterator row = buffer.begin(); row != buffer.end(); ++row) {
1498 
1499  Y[row->first] += row->second;
1500 
1501  V[row->first][row->first] += row->second * row->second;
1502 
1503  for (buffer_type::const_iterator col = buffer.begin(); col != row; ++col) {
1504  V[row->first][col->first] += row->second * col->second;
1505  V[col->first][row->first] = V[row->first][col->first];
1506  }
1507  }
1508  }
1509  }
1510  }
1511  }
Data structure for measured coincidence rate of pair of PMTs.
Definition: JFitK40.hh:64
data_type w[N+1][M+1]
Definition: JPolint.hh:867
const real_type & getReal(const pair_type &pair) const
Get derived parameters.
Definition: JFitK40.hh:1068
JParameter_t t0
time offset [ns]
Definition: JFitK40.hh:562
TPaveText * p1
JParameter_t R
maximal coincidence rate [Hz]
Definition: JFitK40.hh:594
static const int INVALID_INDEX
invalid index
Definition: JFitK40.hh:58
const JK40Parameters_t & getGradient(const double ct) const
Get gradient.
Definition: JFitK40.hh:721
static const JPBS_t PMT(3, 4, 2, 3)
PBS of photo-multiplier tube (PMT)
double getDerivative() const
Get derivative of value.
Definition: JFitK40.hh:332
Data structure for a pair of indices.
void reset()
Reset.
Definition: JVectorND.hh:45
#define R1(x)
int getIndex() const
Get index of PMT used for fixed time offset.
Definition: JFitK40.hh:962
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:456
std::vector< JHitW0 > buffer_type
hits
Definition: JPerth.cc:67
double dt_ns
time difference [ns]
Definition: JFitK40.hh:90
JParameter_t bg
background [Hz/ns]
Definition: JFitK40.hh:563
JMATH::JMatrixNS V
Definition: JFitK40.hh:1376
JParameter_t QE
relative quantum efficiency [unit]
Definition: JFitK40.hh:560
JParameter_t p3
3rd order angle dependence coincidence rate
Definition: JFitK40.hh:597
double error
error of rate [Hz/ns]
Definition: JFitK40.hh:92
estimator_type estimator
M-Estimator function.
Definition: JFitK40.hh:1371
JParameter_t TTS
transition-time spread [ns]
Definition: JFitK40.hh:561
JParameter_t p4
4th order angle dependence coincidence rate
Definition: JFitK40.hh:598
p2
Definition: module-Z:fit.sh:74
JMATH::JVectorND Y
Definition: JFitK40.hh:1513
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
bool isFree() const
Check if parameter is free.
Definition: JFitK40.hh:238
JParameter_t p2
2nd order angle dependence coincidence rate
Definition: JFitK40.hh:596
Auxiliary data structure for derived quantities of a given PMT pair.
Definition: JFitK40.hh:761
then set_variable NUMBER_OF_TESTS else set_variable NUMBER_OF_TESTS fi function gauss()
double getValue(const pair_type &pair, const double dt_ns) const
Get K40 coincidence rate.
Definition: JFitK40.hh:1095
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitK40.hh:570
double u[N+1]
Definition: JPolint.hh:865
Fit model.
Definition: JFitK40.hh:748
double successor
Definition: JFitK40.hh:1514
JPMTParameters_t parameters[NUMBER_OF_PMTS]
Definition: JFitK40.hh:1159
p3
Definition: module-Z:fit.sh:74
double value
value of rate [Hz/ns]
Definition: JFitK40.hh:91
JParameter_t p1
1st order angle dependence coincidence rate
Definition: JFitK40.hh:595
JParameter_t cc
fraction of signal correlated background
Definition: JFitK40.hh:599

Member Data Documentation

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

maximal number of iterations.

Definition at line 1362 of file JFitK40.hh.

constexpr double JCALIBRATE::JFit::EPSILON = 1.0e-4
static

maximal distance to minimum.

Definition at line 1363 of file JFitK40.hh.

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

minimal value control parameter

Definition at line 1364 of file JFitK40.hh.

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

maximal value control parameter

Definition at line 1365 of file JFitK40.hh.

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

multiplication factor control parameter

Definition at line 1366 of file JFitK40.hh.

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

multiplication factor control parameter

Definition at line 1367 of file JFitK40.hh.

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

minimal value diagonal element of matrix

Definition at line 1368 of file JFitK40.hh.

int JCALIBRATE::JFit::debug

Definition at line 1370 of file JFitK40.hh.

estimator_type JCALIBRATE::JFit::estimator

M-Estimator function.

Definition at line 1371 of file JFitK40.hh.

double JCALIBRATE::JFit::lambda

Definition at line 1373 of file JFitK40.hh.

JModel JCALIBRATE::JFit::value

Definition at line 1374 of file JFitK40.hh.

int JCALIBRATE::JFit::numberOfIterations

Definition at line 1375 of file JFitK40.hh.

JMATH::JMatrixNS JCALIBRATE::JFit::V

Definition at line 1376 of file JFitK40.hh.

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

Definition at line 1513 of file JFitK40.hh.

double JCALIBRATE::JFit::successor
private

Definition at line 1514 of file JFitK40.hh.

JModel JCALIBRATE::JFit::previous
private

Definition at line 1515 of file JFitK40.hh.

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

Definition at line 1516 of file JFitK40.hh.


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