Jpp  18.3.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 1174 of file JFitK40.hh.

Member Typedef Documentation

Definition at line 1185 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 1194 of file JFitK40.hh.

1194  :
1195  debug(debug)
1196  {
1197  using namespace JPP;
1198 
1199  estimator.reset(getMEstimator(option));
1200  }
estimator_type estimator
M-Estimator function.
Definition: JFitK40.hh:1373
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 1209 of file JFitK40.hh.

1210  {
1211  using namespace std;
1212  using namespace JPP;
1213 
1214 
1215  value.setIndex();
1216 
1217  const size_t N = value.getN();
1218 
1219  V.resize(N);
1220  Y.resize(N);
1221  h.resize(N);
1222 
1223  int ndf = 0;
1224 
1225  for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1226 
1227  const pair_type& pair = ix->first;
1228 
1229  if (value.parameters[pair.first ].status &&
1230  value.parameters[pair.second].status) {
1231 
1232  ndf += ix->second.size();
1233  }
1234  }
1235 
1236  ndf -= value.getN();
1237 
1238 
1239  lambda = LAMBDA_MIN;
1240 
1241  double precessor = numeric_limits<double>::max();
1242 
1244 
1245  DEBUG("step: " << numberOfIterations << endl);
1246 
1247  evaluate(data);
1248 
1249  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
1250  DEBUG("chi2: " << FIXED(12,3) << successor << endl);
1251 
1252  if (successor < precessor) {
1253 
1254  if (numberOfIterations != 0) {
1255 
1256  if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
1257  return { successor / estimator->getRho(1.0), ndf };
1258  }
1259 
1260  if (lambda > LAMBDA_MIN) {
1261  lambda /= LAMBDA_DOWN;
1262  }
1263  }
1264 
1265  precessor = successor;
1266  previous = value;
1267 
1268  } else {
1269 
1270  value = previous;
1271  lambda *= LAMBDA_UP;
1272 
1273  if (lambda > LAMBDA_MAX) {
1274  return { precessor / estimator->getRho(1.0), ndf }; // no improvement found
1275  }
1276 
1277  evaluate(data);
1278  }
1279 
1280  if (debug >= debug_t) {
1281 
1282  size_t row = 0;
1283 
1284  if (value.R .isFree()) { cout << "R " << FIXED(12,5) << Y[row] << endl; ++row; }
1285  if (value.p1.isFree()) { cout << "p1 " << FIXED(12,5) << Y[row] << endl; ++row; }
1286  if (value.p2.isFree()) { cout << "p2 " << FIXED(12,5) << Y[row] << endl; ++row; }
1287  if (value.p3.isFree()) { cout << "p3 " << FIXED(12,5) << Y[row] << endl; ++row; }
1288  if (value.p4.isFree()) { cout << "p4 " << FIXED(12,5) << Y[row] << endl; ++row; }
1289  if (value.cc.isFree()) { cout << "cc " << FIXED(12,3) << Y[row] << endl; ++row; }
1290 
1291  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1292  if (value.parameters[pmt].QE .isFree()) { cout << "PMT[" << setw(2) << pmt << "].QE " << FIXED(12,5) << Y[row] << endl; ++row; }
1293  if (value.parameters[pmt].TTS.isFree()) { cout << "PMT[" << setw(2) << pmt << "].TTS " << FIXED(12,5) << Y[row] << endl; ++row; }
1294  if (value.parameters[pmt].t0 .isFree()) { cout << "PMT[" << setw(2) << pmt << "].t0 " << FIXED(12,5) << Y[row] << endl; ++row; }
1295  if (value.parameters[pmt].bg .isFree()) { cout << "PMT[" << setw(2) << pmt << "].bg " << FIXED(12,5) << Y[row] << endl; ++row; }
1296  }
1297  }
1298 
1299  // force definite positiveness
1300 
1301  for (size_t i = 0; i != N; ++i) {
1302 
1303  if (V(i,i) < PIVOT) {
1304  V(i,i) = PIVOT;
1305  }
1306 
1307  h[i] = 1.0 / sqrt(V(i,i));
1308  }
1309 
1310  // normalisation
1311 
1312  for (size_t i = 0; i != N; ++i) {
1313  for (size_t j = 0; j != i; ++j) {
1314  V(j,i) *= h[i] * h[j];
1315  V(i,j) = V(j,i);
1316  }
1317  }
1318 
1319  for (size_t i = 0; i != N; ++i) {
1320  V(i,i) = 1.0 + lambda;
1321  }
1322 
1323  // solve A x = b
1324 
1325  for (size_t col = 0; col != N; ++col) {
1326  Y[col] *= h[col];
1327  }
1328 
1329  try {
1330  V.solve(Y);
1331  }
1332  catch (const exception& error) {
1333 
1334  ERROR("JGandalf: " << error.what() << endl << V << endl);
1335 
1336  break;
1337  }
1338 
1339  // update value
1340 
1341  const double factor = 2.0;
1342 
1343  size_t row = 0;
1344 
1345  if (value.R .isFree()) { value.R -= factor * h[row] * Y[row]; ++row; }
1346  if (value.p1.isFree()) { value.p1 -= factor * h[row] * Y[row]; ++row; }
1347  if (value.p2.isFree()) { value.p2 -= factor * h[row] * Y[row]; ++row; }
1348  if (value.p3.isFree()) { value.p3 -= factor * h[row] * Y[row]; ++row; }
1349  if (value.p4.isFree()) { value.p4 -= factor * h[row] * Y[row]; ++row; }
1350  if (value.cc.isFree()) { value.cc -= factor * h[row] * Y[row]; ++row; }
1351 
1352  for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1353  if (value.parameters[pmt].QE .isFree()) { value.parameters[pmt].QE -= factor * h[row] * Y[row]; ++row; }
1354  if (value.parameters[pmt].TTS.isFree()) { value.parameters[pmt].TTS -= factor * h[row] * Y[row]; ++row; }
1355  if (value.parameters[pmt].t0 .isFree()) { value.parameters[pmt].t0 -= factor * h[row] * Y[row]; ++row; }
1356  if (value.parameters[pmt].bg .isFree()) { value.parameters[pmt].bg -= factor * h[row] * Y[row]; ++row; }
1357  }
1358  }
1359 
1360  return { precessor / estimator->getRho(1.0), ndf };
1361  }
std::vector< double > h
Definition: JFitK40.hh:1518
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:973
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:1367
JParameter_t bg
background [Hz/ns]
Definition: JFitK40.hh:563
JMATH::JMatrixNS V
Definition: JFitK40.hh:1378
static constexpr double EPSILON
maximal distance to minimum.
Definition: JFitK40.hh:1365
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:1377
void evaluate(const data_type &data)
Evaluation of fit.
Definition: JFitK40.hh:1386
estimator_type estimator
M-Estimator function.
Definition: JFitK40.hh:1373
#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:1366
JParameter_t p4
4th order angle dependence coincidence rate
Definition: JFitK40.hh:598
JMATH::JVectorND Y
Definition: JFitK40.hh:1515
size_t getN() const
Get number of fit parameters.
Definition: JFitK40.hh:999
static constexpr int MAXIMUM_ITERATIONS
maximal number of iterations.
Definition: JFitK40.hh:1364
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:1368
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
double successor
Definition: JFitK40.hh:1516
JPMTParameters_t parameters[NUMBER_OF_PMTS]
Definition: JFitK40.hh:1161
static constexpr double LAMBDA_DOWN
multiplication factor control parameter
Definition: JFitK40.hh:1369
static constexpr double PIVOT
minimal value diagonal element of matrix
Definition: JFitK40.hh:1370
#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 1386 of file JFitK40.hh.

1387  {
1388  using namespace std;
1389  using namespace JPP;
1390 
1391  typedef JModel::real_type real_type;
1392 
1393 
1394  successor = 0.0;
1395 
1396  V.reset();
1397  Y.reset();
1398 
1399 
1400  // model parameter indices
1401 
1402  const struct M_t {
1403  M_t(const JModel& model)
1404  {
1405  R = model.getIndex(&JK40Parameters_t::R);
1406  p1 = model.getIndex(&JK40Parameters_t::p1);
1407  p2 = model.getIndex(&JK40Parameters_t::p2);
1408  p3 = model.getIndex(&JK40Parameters_t::p3);
1409  p4 = model.getIndex(&JK40Parameters_t::p4);
1410  cc = model.getIndex(&JK40Parameters_t::cc);
1411  }
1412 
1413  int R;
1414  int p1;
1415  int p2;
1416  int p3;
1417  int p4;
1418  int cc;
1419 
1420  } M(value);
1421 
1422 
1423  // PMT parameter indices
1424 
1425  struct I_t {
1426  I_t(const JModel& model, const int pmt) :
1427  QE (INVALID_INDEX),
1428  TTS(INVALID_INDEX),
1429  t0 (INVALID_INDEX),
1430  bg (INVALID_INDEX)
1431  {
1432  const int index = model.getIndex(pmt);
1433 
1434  int N = 0;
1435 
1436  if (model.parameters[pmt].QE .isFree()) { QE = index + N; ++N; }
1437  if (model.parameters[pmt].TTS.isFree()) { TTS = index + N; ++N; }
1438  if (model.parameters[pmt].t0 .isFree()) { t0 = index + N; ++N; }
1439  if (model.parameters[pmt].bg .isFree()) { bg = index + N; ++N; }
1440  }
1441 
1442  int QE;
1443  int TTS;
1444  int t0;
1445  int bg;
1446  };
1447 
1448 
1450 
1451  buffer_type buffer;
1452 
1453  for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1454 
1455  const pair_type& pair = ix->first;
1456 
1457  if (value.parameters[pair.first ].status &&
1458  value.parameters[pair.second].status) {
1459 
1460  const real_type& real = value.getReal(pair);
1461 
1462  const JGauss gauss(real.t0, real.sigma, real.signal);
1463 
1464  const double R1 = value.getValue (real.ct);
1465  const JK40Parameters_t& R1p = value.getGradient(real.ct);
1466 
1467  const std::pair<I_t, I_t> PMT(I_t(value, pair.first),
1468  I_t(value, pair.second));
1469 
1470  for (const rate_type& iy : ix->second) {
1471 
1472  const double R2 = gauss.getValue (iy.dt_ns);
1473  const JGauss& R2p = gauss.getGradient(iy.dt_ns);
1474 
1475  const double R = real.background + R1 * (value.cc() + R2);
1476  const double u = (iy.value - R) / iy.error;
1477  const double w = -estimator->getPsi(u) / iy.error;
1478 
1479  successor += estimator->getRho(u);
1480 
1481  buffer.clear();
1482 
1483  if (M.R != INVALID_INDEX) { buffer.push_back({M.R, w * (value.cc() + R2) * R1p.R () * value.R .getDerivative()}); }
1484  if (M.p1 != INVALID_INDEX) { buffer.push_back({M.p1, w * (value.cc() + R2) * R1p.p1() * value.p1.getDerivative()}); }
1485  if (M.p2 != INVALID_INDEX) { buffer.push_back({M.p2, w * (value.cc() + R2) * R1p.p2() * value.p2.getDerivative()}); }
1486  if (M.p3 != INVALID_INDEX) { buffer.push_back({M.p3, w * (value.cc() + R2) * R1p.p3() * value.p3.getDerivative()}); }
1487  if (M.p4 != INVALID_INDEX) { buffer.push_back({M.p4, w * (value.cc() + R2) * R1p.p4() * value.p4.getDerivative()}); }
1488  if (M.cc != INVALID_INDEX) { buffer.push_back({M.cc, w * R1 * R1p.cc() * value.cc.getDerivative()}); }
1489 
1490  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()}); }
1491  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()}); }
1492  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}); }
1493  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}); }
1494  if (PMT.first .t0 != INVALID_INDEX) { buffer.push_back({PMT.first .t0, w * R1 * R2p.mean * value.parameters[pair.first ].t0 .getDerivative() * +1.0}); }
1495  if (PMT.second.t0 != INVALID_INDEX) { buffer.push_back({PMT.second.t0, w * R1 * R2p.mean * value.parameters[pair.second].t0 .getDerivative() * -1.0}); }
1496  if (PMT.first .bg != INVALID_INDEX) { buffer.push_back({PMT.first .bg, w * value.parameters[pair.first ].bg .getDerivative()}); }
1497  if (PMT.second.bg != INVALID_INDEX) { buffer.push_back({PMT.second.bg, w * value.parameters[pair.second].bg .getDerivative()}); }
1498 
1499  for (buffer_type::const_iterator row = buffer.begin(); row != buffer.end(); ++row) {
1500 
1501  Y[row->first] += row->second;
1502 
1503  V[row->first][row->first] += row->second * row->second;
1504 
1505  for (buffer_type::const_iterator col = buffer.begin(); col != row; ++col) {
1506  V[row->first][col->first] += row->second * col->second;
1507  V[col->first][row->first] = V[row->first][col->first];
1508  }
1509  }
1510  }
1511  }
1512  }
1513  }
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:1070
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:964
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:1378
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:1373
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:1515
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:1097
p4
Definition: JGraph2D.sh:83
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:1516
JPMTParameters_t parameters[NUMBER_OF_PMTS]
Definition: JFitK40.hh:1161
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 1364 of file JFitK40.hh.

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

maximal distance to minimum.

Definition at line 1365 of file JFitK40.hh.

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

minimal value control parameter

Definition at line 1366 of file JFitK40.hh.

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

maximal value control parameter

Definition at line 1367 of file JFitK40.hh.

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

multiplication factor control parameter

Definition at line 1368 of file JFitK40.hh.

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

multiplication factor control parameter

Definition at line 1369 of file JFitK40.hh.

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

minimal value diagonal element of matrix

Definition at line 1370 of file JFitK40.hh.

int JCALIBRATE::JFit::debug

Definition at line 1372 of file JFitK40.hh.

estimator_type JCALIBRATE::JFit::estimator

M-Estimator function.

Definition at line 1373 of file JFitK40.hh.

double JCALIBRATE::JFit::lambda

Definition at line 1375 of file JFitK40.hh.

JModel JCALIBRATE::JFit::value

Definition at line 1376 of file JFitK40.hh.

int JCALIBRATE::JFit::numberOfIterations

Definition at line 1377 of file JFitK40.hh.

JMATH::JMatrixNS JCALIBRATE::JFit::V

Definition at line 1378 of file JFitK40.hh.

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

Definition at line 1515 of file JFitK40.hh.

double JCALIBRATE::JFit::successor
private

Definition at line 1516 of file JFitK40.hh.

JModel JCALIBRATE::JFit::previous
private

Definition at line 1517 of file JFitK40.hh.

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

Definition at line 1518 of file JFitK40.hh.


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