Jpp 21.0.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
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.
 
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
 
bool TEST = false
 

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 1257 of file JFitK40.hh.

Member Typedef Documentation

◆ estimator_type

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

1277 :
1278 debug(debug)
1279 {
1280 using namespace JPP;
1281
1282 estimator.reset(getMEstimator(option));
1283 }
estimator_type estimator
M-Estimator function.
Definition JFitK40.hh:1479
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 1292 of file JFitK40.hh.

1293 {
1294 using namespace std;
1295 using namespace JPP;
1296
1297
1298 value.setIndex();
1299
1300 const size_t N = value.getN();
1301
1302 V.resize(N);
1303 Y.resize(N);
1304 h.resize(N);
1305
1306 double xmax = numeric_limits<double>::lowest();
1307 double xmin = numeric_limits<double>::max();
1308
1309 int ndf = 0;
1310
1311 for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1312
1313 const pair_type& pair = ix->first;
1314
1315 if (value.parameters[pair.first ].status &&
1316 value.parameters[pair.second].status) {
1317
1318 ndf += ix->second.size();
1319
1320 for (const rate_type& iy : ix->second) {
1321 if (iy.dt_ns > xmax) { xmax = iy.dt_ns; }
1322 if (iy.dt_ns < xmin) { xmin = iy.dt_ns; }
1323 }
1324 }
1325 }
1326
1327 ndf -= value.getN();
1328
1329 if (ndf < 0) {
1330 return { 0.0, ndf };
1331 }
1332
1333 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1334 if (value.parameters[pmt].t0.isFree()) {
1335 value.parameters[pmt].t0.setLimits(xmin, xmax);
1336 }
1337 }
1338
1340
1341 double precessor = numeric_limits<double>::max();
1342
1344
1345 DEBUG("step: " << numberOfIterations << endl);
1346
1347 evaluate(data);
1348
1349 DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
1350 DEBUG("chi2: " << FIXED(12,3) << successor << endl);
1351
1352 if (successor < precessor) {
1353
1354 if (numberOfIterations != 0) {
1355
1356 if (fabs(precessor - successor) < EPSILON) {
1357
1358 seterr(data);
1359
1360 return { successor / estimator->getRho(1.0), ndf };
1361 }
1362
1363 if (lambda > LAMBDA_MIN) {
1365 }
1366 }
1367
1368 precessor = successor;
1369 previous = value;
1370
1371 } else {
1372
1373 value = previous;
1374 lambda *= LAMBDA_UP;
1375
1376 if (lambda > LAMBDA_MAX) {
1377 break;
1378 }
1379
1380 evaluate(data);
1381 }
1382
1383 if (debug >= debug_t) {
1384
1385 size_t row = 0;
1386
1387 if (value.R .isFree()) { cout << "R " << FIXED(12,5) << Y[row] << endl; ++row; }
1388 if (value.p1.isFree()) { cout << "p1 " << FIXED(12,5) << Y[row] << endl; ++row; }
1389 if (value.p2.isFree()) { cout << "p2 " << FIXED(12,5) << Y[row] << endl; ++row; }
1390 if (value.p3.isFree()) { cout << "p3 " << FIXED(12,5) << Y[row] << endl; ++row; }
1391 if (value.p4.isFree()) { cout << "p4 " << FIXED(12,5) << Y[row] << endl; ++row; }
1392 if (value.cc.isFree()) { cout << "cc " << FIXED(12,3) << Y[row] << endl; ++row; }
1393
1394 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1395 if (value.parameters[pmt].QE .isFree()) { cout << "PMT[" << setw(2) << pmt << "].QE " << FIXED(12,5) << Y[row] << endl; ++row; }
1396 if (value.parameters[pmt].TTS.isFree()) { cout << "PMT[" << setw(2) << pmt << "].TTS " << FIXED(12,5) << Y[row] << endl; ++row; }
1397 if (value.parameters[pmt].t0 .isFree()) { cout << "PMT[" << setw(2) << pmt << "].t0 " << FIXED(12,5) << Y[row] << endl; ++row; }
1398 if (value.parameters[pmt].bg .isFree()) { cout << "PMT[" << setw(2) << pmt << "].bg " << FIXED(12,5) << Y[row] << endl; ++row; }
1399 }
1400 }
1401
1402 // force definite positiveness
1403
1404 for (size_t i = 0; i != N; ++i) {
1405
1406 if (V(i,i) < PIVOT) {
1407 V(i,i) = PIVOT;
1408 }
1409
1410 h[i] = 1.0 / sqrt(V(i,i));
1411 }
1412
1413 // normalisation
1414
1415 for (size_t i = 0; i != N; ++i) {
1416 for (size_t j = 0; j != i; ++j) {
1417 V(j,i) *= h[i] * h[j];
1418 V(i,j) = V(j,i);
1419 }
1420 }
1421
1422 for (size_t i = 0; i != N; ++i) {
1423 V(i,i) = 1.0 + lambda;
1424 }
1425
1426 // solve A x = b
1427
1428 for (size_t col = 0; col != N; ++col) {
1429 Y[col] *= h[col];
1430 }
1431
1432 try {
1433 V.solve(Y);
1434 }
1435 catch (const exception& error) {
1436
1437 ERROR("JGandalf: " << error.what() << endl << V << endl);
1438
1439 break;
1440 }
1441
1442 // update value
1443
1444 const double factor = 2.0;
1445
1446 size_t row = 0;
1447
1448 if (value.R .isFree()) { value.R -= factor * h[row] * Y[row]; ++row; }
1449 if (value.p1.isFree()) { value.p1 -= factor * h[row] * Y[row]; ++row; }
1450 if (value.p2.isFree()) { value.p2 -= factor * h[row] * Y[row]; ++row; }
1451 if (value.p3.isFree()) { value.p3 -= factor * h[row] * Y[row]; ++row; }
1452 if (value.p4.isFree()) { value.p4 -= factor * h[row] * Y[row]; ++row; }
1453 if (value.cc.isFree()) { value.cc -= factor * h[row] * Y[row]; ++row; }
1454 if (value.bc.isFree()) { value.bc -= factor * h[row] * Y[row]; ++row; }
1455
1456 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1457 if (value.parameters[pmt].QE .isFree()) { value.parameters[pmt].QE -= factor * h[row] * Y[row]; ++row; }
1458 if (value.parameters[pmt].TTS.isFree()) { value.parameters[pmt].TTS -= factor * h[row] * Y[row]; ++row; }
1459 if (value.parameters[pmt].t0 .isFree()) { value.parameters[pmt].t0 -= factor * h[row] * Y[row]; ++row; }
1460 if (value.parameters[pmt].bg .isFree()) { value.parameters[pmt].bg -= factor * h[row] * Y[row]; ++row; }
1461 }
1462 }
1463
1464 seterr(data);
1465
1466 return { precessor / estimator->getRho(1.0), ndf };
1467 }
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define ERROR(A)
Definition JMessage.hh:66
std::vector< double > h
Definition JFitK40.hh:1788
static constexpr double LAMBDA_MIN
minimal value control parameter
Definition JFitK40.hh:1472
static constexpr double LAMBDA_DOWN
multiplication factor control parameter
Definition JFitK40.hh:1475
void seterr(const data_type &data)
Set errors.
Definition JFitK40.hh:1749
static constexpr double LAMBDA_MAX
maximal value control parameter
Definition JFitK40.hh:1473
static constexpr double LAMBDA_UP
multiplication factor control parameter
Definition JFitK40.hh:1474
JMATH::JMatrixNS V
Definition JFitK40.hh:1485
static constexpr double EPSILON
maximal distance to minimum.
Definition JFitK40.hh:1471
void evaluate(const data_type &data)
Evaluation of fit.
Definition JFitK40.hh:1495
static constexpr int MAXIMUM_ITERATIONS
maximal number of iterations.
Definition JFitK40.hh:1470
static constexpr double PIVOT
minimal value diagonal element of matrix
Definition JFitK40.hh:1476
JMATH::JVectorND Y
Definition JFitK40.hh:1785
bool isFree() const
Check if parameter is free.
Definition JFitK40.hh:244
void setLimits(const double xmin, const double xmax)
Set limits.
Definition JFitK40.hh:326
int j
Definition JPolint.hh:801
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JParameter_t bc
constant background
Definition JFitK40.hh:714
JParameter_t R
maximal coincidence rate [Hz]
Definition JFitK40.hh:708
JParameter_t p1
1st order angle dependence coincidence rate
Definition JFitK40.hh:709
JParameter_t p2
2nd order angle dependence coincidence rate
Definition JFitK40.hh:710
JParameter_t p3
3rd order angle dependence coincidence rate
Definition JFitK40.hh:711
JParameter_t p4
4th order angle dependence coincidence rate
Definition JFitK40.hh:712
JParameter_t cc
fraction of signal correlated background
Definition JFitK40.hh:713
JPMTParameters_t parameters[NUMBER_OF_PMTS]
Definition JFitK40.hh:845
size_t getN() const
Get number of fit parameters.
Definition JFitK40.hh:1130
void setIndex()
Set index of PMT used for fixed time offset.
Definition JFitK40.hh:1104
JParameter_t t0
time offset [ns]
Definition JFitK40.hh:602
JParameter_t TTS
transition-time spread [ns]
Definition JFitK40.hh:601
JParameter_t bg
background [Hz/ns]
Definition JFitK40.hh:603
JParameter_t QE
relative quantum efficiency [unit]
Definition JFitK40.hh:600
Data structure for measured coincidence rate of pair of PMTs.
Definition JFitK40.hh:70
double dt_ns
time difference [ns]
Definition JFitK40.hh:96
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 1495 of file JFitK40.hh.

1496 {
1497 using namespace std;
1498 using namespace JPP;
1499
1500 typedef JModel::real_type real_type;
1501
1502
1503 successor = 0.0;
1504
1505 V.reset();
1506 Y.reset();
1507
1508
1509 // model parameter indices
1510
1511 const struct M_t {
1512 M_t(const JModel& model)
1513 {
1514 R = model.getIndex(&JK40Parameters_t::R);
1515 p1 = model.getIndex(&JK40Parameters_t::p1);
1516 p2 = model.getIndex(&JK40Parameters_t::p2);
1517 p3 = model.getIndex(&JK40Parameters_t::p3);
1518 p4 = model.getIndex(&JK40Parameters_t::p4);
1519 cc = model.getIndex(&JK40Parameters_t::cc);
1520 bc = model.getIndex(&JK40Parameters_t::bc);
1521 }
1522
1523 int R;
1524 int p1;
1525 int p2;
1526 int p3;
1527 int p4;
1528 int cc;
1529 int bc;
1530
1531 } M(value);
1532
1533
1534 // PMT parameter indices
1535
1536 struct i_t {
1537 i_t() :
1538 QE (INVALID_INDEX),
1539 TTS(INVALID_INDEX),
1540 t0 (INVALID_INDEX),
1541 bg (INVALID_INDEX)
1542 {}
1543
1544 int QE;
1545 int TTS;
1546 int t0;
1547 int bg;
1548 };
1549
1550 const struct I_t : public std::array<i_t, NUMBER_OF_PMTS> {
1551 I_t(const JModel& value)
1552 {
1553 int N = value.JK40Parameters::getN();
1554
1555 for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
1556 if (value.parameters[i].QE .isFree()) { (*this)[i].QE = N; ++N; }
1557 if (value.parameters[i].TTS.isFree()) { (*this)[i].TTS = N; ++N; }
1558 if (value.parameters[i].t0 .isFree()) { (*this)[i].t0 = N; ++N; }
1559 if (value.parameters[i].bg .isFree()) { (*this)[i].bg = N; ++N; }
1560 }
1561 }
1562
1563 } I(value);
1564
1565
1566
1567 struct buffer_type : public vector< pair<int, double> > {
1568 double operator[](const int index) const
1569 {
1570 for (const_iterator i = this->begin(); i != this->end(); ++i) {
1571 if (i->first == index) {
1572 return i->second;
1573 }
1574 }
1575
1576 THROW(JValueOutOfRange, "Invalid index " << index);
1577 }
1578 };
1579
1580 buffer_type buffer;
1581
1582#define PUSH_BACK(i,v) if (i != INVALID_INDEX) { buffer.push_back({i, v}); }
1583
1584
1585 size_t number_of_errors = 0;
1586
1587 for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1588
1589 const pair_type& pair = ix->first;
1590
1591 if (value.parameters[pair.first ].status &&
1592 value.parameters[pair.second].status) {
1593
1594 const real_type& real = value.getReal(pair);
1595
1596 const JBell bell(real.t0, real.sigma, real.signal, 0.0, BELL_SHAPE);
1597
1598 const double R1 = value.getValue (real.ct);
1599 const JK40Parameters_t& R1p = value.getGradient(real.ct);
1600
1601 for (const rate_type& iy : ix->second) {
1602
1603 const double R2 = bell.getValue (iy.dt_ns);
1604 const JBell_t& R2p = bell.getGradient(iy.dt_ns);
1605
1606 const double R = real.bc + real.background + R1 * (real.cc + R2);
1607 const double u = (iy.value - R) / iy.error;
1608 const double w = -estimator->getPsi(u) / iy.error;
1609
1610 successor += estimator->getRho(u);
1611
1612 buffer.clear();
1613
1614 PUSH_BACK(M.R, w * (real.cc + R2) * R1p.R () * value.R .getDerivative());
1615 PUSH_BACK(M.p1, w * (real.cc + R2) * R1p.p1() * value.p1.getDerivative());
1616 PUSH_BACK(M.p2, w * (real.cc + R2) * R1p.p2() * value.p2.getDerivative());
1617 PUSH_BACK(M.p3, w * (real.cc + R2) * R1p.p3() * value.p3.getDerivative());
1618 PUSH_BACK(M.p4, w * (real.cc + R2) * R1p.p4() * value.p4.getDerivative());
1619 PUSH_BACK(M.cc, w * real.signal * R1p.cc() * value.cc.getDerivative());
1620 PUSH_BACK(M.bc, w * R1p.bc() * value.bc.getDerivative());
1621
1622 PUSH_BACK(I[pair.first] .QE, w * R1 * R2p.signal * value.parameters[pair.second].QE () * value.parameters[pair.first ].QE .getDerivative());
1623 PUSH_BACK(I[pair.second].QE, w * R1 * R2p.signal * value.parameters[pair.first ].QE () * value.parameters[pair.second].QE .getDerivative());
1624 PUSH_BACK(I[pair.first] .TTS, w * R1 * R2p.sigma * value.parameters[pair.first ].TTS() * value.parameters[pair.first ].TTS.getDerivative() / real.sigma);
1625 PUSH_BACK(I[pair.second].TTS, w * R1 * R2p.sigma * value.parameters[pair.second].TTS() * value.parameters[pair.second].TTS.getDerivative() / real.sigma);
1626 PUSH_BACK(I[pair.first] .t0, w * R1 * R2p.mean * value.parameters[pair.first ].t0 .getDerivative() * +1.0);
1627 PUSH_BACK(I[pair.second].t0, w * R1 * R2p.mean * value.parameters[pair.second].t0 .getDerivative() * -1.0);
1628 PUSH_BACK(I[pair.first] .bg, w * value.parameters[pair.first ].bg .getDerivative());
1629 PUSH_BACK(I[pair.second].bg, w * value.parameters[pair.second].bg .getDerivative());
1630
1631 if (TEST) {
1632
1633 DEBUG("PMT pair(" << setw(2) << pair.first << "," << setw(2) << pair.second << ") " << FIXED(7,3) << iy.dt_ns << " [ns]" << endl);
1634
1635 const double PRECISION = 1.0e-5;
1636
1637#define MAKE_TEST(i,v) if (i != INVALID_INDEX) { \
1638 \
1639 const bool status = fabs(buffer[i] - v) <= PRECISION; \
1640 \
1641 DEBUG((status ? GREEN : RED) \
1642 << setw(20) << left << #i << right << ' ' \
1643 << setw(3) << i << ' ' \
1644 << FIXED(12,5) << buffer[i] << ' ' \
1645 << FIXED(12,5) << v << ' ' \
1646 << (!status ? "***" : "") \
1647 << RESET << endl); \
1648 \
1649 if (!status) { \
1650 number_of_errors += 1; \
1651 } \
1652 }
1653
1654 struct JTest_t : public JModel
1655 {
1656 JTest_t& operator=(const JModel& model)
1657 {
1658 static_cast<JModel&>(*this) = model;
1659
1660 R .relax();
1661 p1.relax();
1662 p2.relax();
1663 p3.relax();
1664 p4.relax();
1665 cc.relax();
1666 bc.relax();
1667
1668 for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
1669 parameters[i].QE .relax();
1670 parameters[i].TTS.relax();
1671 parameters[i].t0 .relax();
1672 parameters[i].bg .relax();
1673 }
1674
1675 return *this;
1676 }
1677
1678 double operator()(const pair_type& pair, const rate_type& iy) const
1679 {
1680 return (iy.value - getValue(pair, iy.dt_ns)) / iy.error;
1681 }
1682 };
1683
1684 const double DX = 1.0e-8; // dx
1685 JTest_t m1, m2; // y1, y2
1686
1687 // derivative
1688
1689 auto fp = [&DX = DX, &m1 = m1, &m2 = m2, &estimator = estimator](const pair_type& pair, const rate_type& iy)
1690 {
1691 return (estimator->getRho(m2(pair, iy)) - estimator->getRho(m1(pair, iy))) / DX;
1692 };
1693
1694 { m1 = m2 = value; m2.R += 0.5*DX; m1.R -= 0.5*DX; MAKE_TEST(M.R, fp(pair, iy) * value.R .getDerivative()); }
1695 { m1 = m2 = value; m2.p1 += 0.5*DX; m1.p1 -= 0.5*DX; MAKE_TEST(M.p1, fp(pair, iy) * value.p1.getDerivative()); }
1696 { m1 = m2 = value; m2.p2 += 0.5*DX; m1.p2 -= 0.5*DX; MAKE_TEST(M.p2, fp(pair, iy) * value.p2.getDerivative()); }
1697 { m1 = m2 = value; m2.p3 += 0.5*DX; m1.p3 -= 0.5*DX; MAKE_TEST(M.p3, fp(pair, iy) * value.p3.getDerivative()); }
1698 { m1 = m2 = value; m2.p4 += 0.5*DX; m1.p4 -= 0.5*DX; MAKE_TEST(M.p4, fp(pair, iy) * value.p4.getDerivative()); }
1699 { m1 = m2 = value; m2.cc += 0.5*DX; m1.cc -= 0.5*DX; MAKE_TEST(M.cc, fp(pair, iy) * value.cc.getDerivative()); }
1700 { m1 = m2 = value; m2.bc += 0.5*DX; m1.bc -= 0.5*DX; MAKE_TEST(M.bc, fp(pair, iy) * value.bc.getDerivative()); }
1701
1702 { m1 = m2 = value; m2.parameters[pair.first] .QE += 0.5*DX; m1.parameters[pair.first] .QE -= 0.5*DX; MAKE_TEST(I[pair.first] .QE, fp(pair, iy) * value.parameters[pair.first] .QE .getDerivative()); }
1703 { m1 = m2 = value; m2.parameters[pair.second].QE += 0.5*DX; m1.parameters[pair.second].QE -= 0.5*DX; MAKE_TEST(I[pair.second].QE, fp(pair, iy) * value.parameters[pair.second].QE .getDerivative()); }
1704 { m1 = m2 = value; m2.parameters[pair.first] .TTS += 0.5*DX; m1.parameters[pair.first] .TTS -= 0.5*DX; MAKE_TEST(I[pair.first] .TTS, fp(pair, iy) * value.parameters[pair.first] .TTS.getDerivative()); }
1705 { m1 = m2 = value; m2.parameters[pair.second].TTS += 0.5*DX; m1.parameters[pair.second].TTS -= 0.5*DX; MAKE_TEST(I[pair.second].TTS, fp(pair, iy) * value.parameters[pair.second].TTS.getDerivative()); }
1706 if (pair.first != value.getIndex()) {
1707 m1 = m2 = value; m2.parameters[pair.first] .t0 += 0.5*DX; m1.parameters[pair.first] .t0 -= 0.5*DX; MAKE_TEST(I[pair.first] .t0, fp(pair, iy) * value.parameters[pair.first] .t0 .getDerivative());
1708 }
1709 if (pair.second != value.getIndex()) {
1710 m1 = m2 = value; m2.parameters[pair.second].t0 += 0.5*DX; m1.parameters[pair.second].t0 -= 0.5*DX; MAKE_TEST(I[pair.second].t0, fp(pair, iy) * value.parameters[pair.second].t0 .getDerivative());
1711 }
1712 { m1 = m2 = value; m2.parameters[pair.first] .bg += 0.5*DX; m1.parameters[pair.first] .bg -= 0.5*DX; MAKE_TEST(I[pair.first] .bg, fp(pair, iy) * value.parameters[pair.first] .bg .getDerivative()); }
1713 { m1 = m2 = value; m2.parameters[pair.second].bg += 0.5*DX; m1.parameters[pair.second].bg -= 0.5*DX; MAKE_TEST(I[pair.second].bg, fp(pair, iy) * value.parameters[pair.second].bg .getDerivative()); }
1714
1715 cout << endl;
1716 }
1717
1718 for (buffer_type::const_iterator row = buffer.begin(); row != buffer.end(); ++row) {
1719
1720 Y[row->first] += row->second;
1721
1722 V[row->first][row->first] += row->second * row->second;
1723
1724 for (buffer_type::const_iterator col = buffer.begin(); col != row; ++col) {
1725 V[row->first][col->first] += row->second * col->second;
1726 V[col->first][row->first] = V[row->first][col->first];
1727 }
1728 }
1729 }
1730 }
1731 }
1732
1733#undef PUSH_BACK
1734
1735 if (TEST) {
1736
1737 STATUS("Test finished with " << number_of_errors << " errors." << endl);
1738
1739 exit(number_of_errors == 0 ? 0 : 1);
1740 }
1741 }
TPaveText * p1
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
#define PUSH_BACK(i, v)
#define MAKE_TEST(i, v)
#define STATUS(A)
Definition JMessage.hh:63
double getDerivative() const
Get derivative of value.
Definition JFitK40.hh:383
Interface to read input and write output for TObject tests.
Definition JTest_t.hh:42
Exception for accessing a value in a collection that is outside of its range.
#define R1(x)
static const int INVALID_INDEX
invalid index
Definition JFitK40.hh:60
static double BELL_SHAPE
Bell shape.
Definition JFitK40.hh:64
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition JScale.hh:47
Model for fit to acoustics data.
Fit parameters for two-fold coincidence rate due to K40.
Definition JFitK40.hh:610
const JK40Parameters_t & getGradient(const double ct) const
Get gradient.
Definition JFitK40.hh:816
Auxiliary data structure for derived quantities of a given PMT pair.
Definition JFitK40.hh:887
int getIndex() const
Get index of PMT used for fixed time offset.
Definition JFitK40.hh:1095
const real_type & getReal(const pair_type &pair) const
Get derived quantities.
Definition JFitK40.hh:1170
double getValue(const pair_type &pair, const double dt_ns) const
Get K40 coincidence rate.
Definition JFitK40.hh:1210
double error
error of rate [Hz/ns]
Definition JFitK40.hh:98
double value
value of rate [Hz/ns]
Definition JFitK40.hh:97
Bell function object.
Definition JBell.hh:32
Gauss model.
Definition JGauss.hh:32
double background
Definition JGauss.hh:164
double signal
Definition JGauss.hh:163
JMatrixND & reset()
Set matrix to the null matrix.
Definition JMatrixND.hh:459
void reset()
Reset.
Definition JVectorND.hh:45

◆ seterr()

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

Set errors.

Parameters
datadata

Definition at line 1749 of file JFitK40.hh.

1750 {
1751 using namespace std;
1752
1753 error.reset();
1754
1755 evaluate(data);
1756
1757 try {
1758 V.invert();
1759 }
1760 catch (const exception& error) {}
1761
1762#define SQRT(X) (X >= 0.0 ? sqrt(X) : std::numeric_limits<double>::max())
1763
1764 size_t i = 0;
1765
1766 if (value.R .isFree()) { error.R = SQRT(V(i,i)); ++i; }
1767 if (value.p1.isFree()) { error.p1 = SQRT(V(i,i)); ++i; }
1768 if (value.p2.isFree()) { error.p2 = SQRT(V(i,i)); ++i; }
1769 if (value.p3.isFree()) { error.p3 = SQRT(V(i,i)); ++i; }
1770 if (value.p4.isFree()) { error.p4 = SQRT(V(i,i)); ++i; }
1771 if (value.cc.isFree()) { error.cc = SQRT(V(i,i)); ++i; }
1772 if (value.bc.isFree()) { error.bc = SQRT(V(i,i)); ++i; }
1773
1774 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1775 if (value.parameters[pmt].QE .isFree()) { error.parameters[pmt].QE = SQRT(V(i,i)); ++i; }
1776 if (value.parameters[pmt].TTS.isFree()) { error.parameters[pmt].TTS = SQRT(V(i,i)); ++i; }
1777 if (value.parameters[pmt].t0 .isFree()) { error.parameters[pmt].t0 = SQRT(V(i,i)); ++i; }
1778 if (value.parameters[pmt].bg .isFree()) { error.parameters[pmt].bg = SQRT(V(i,i)); ++i; }
1779 }
1780
1781#undef SQRT
1782 }
#define SQRT(X)
void invert()
Invert matrix according LDU decomposition.
Definition JMatrixNS.hh:75

Member Data Documentation

◆ MAXIMUM_ITERATIONS

int JCALIBRATE::JFit::MAXIMUM_ITERATIONS = 100000
staticconstexpr

maximal number of iterations.

Definition at line 1470 of file JFitK40.hh.

◆ EPSILON

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

maximal distance to minimum.

Definition at line 1471 of file JFitK40.hh.

◆ LAMBDA_MIN

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

minimal value control parameter

Definition at line 1472 of file JFitK40.hh.

◆ LAMBDA_MAX

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

maximal value control parameter

Definition at line 1473 of file JFitK40.hh.

◆ LAMBDA_UP

double JCALIBRATE::JFit::LAMBDA_UP = 10.0
staticconstexpr

multiplication factor control parameter

Definition at line 1474 of file JFitK40.hh.

◆ LAMBDA_DOWN

double JCALIBRATE::JFit::LAMBDA_DOWN = 10.0
staticconstexpr

multiplication factor control parameter

Definition at line 1475 of file JFitK40.hh.

◆ PIVOT

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

minimal value diagonal element of matrix

Definition at line 1476 of file JFitK40.hh.

◆ debug

int JCALIBRATE::JFit::debug

Definition at line 1478 of file JFitK40.hh.

◆ estimator

estimator_type JCALIBRATE::JFit::estimator

M-Estimator function.

Definition at line 1479 of file JFitK40.hh.

◆ lambda

double JCALIBRATE::JFit::lambda

Definition at line 1481 of file JFitK40.hh.

◆ value

JModel JCALIBRATE::JFit::value

Definition at line 1482 of file JFitK40.hh.

◆ error

JModel_t JCALIBRATE::JFit::error

Definition at line 1483 of file JFitK40.hh.

◆ numberOfIterations

int JCALIBRATE::JFit::numberOfIterations

Definition at line 1484 of file JFitK40.hh.

◆ V

JMATH::JMatrixNS JCALIBRATE::JFit::V

Definition at line 1485 of file JFitK40.hh.

◆ TEST

bool JCALIBRATE::JFit::TEST = false

Definition at line 1487 of file JFitK40.hh.

◆ Y

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

Definition at line 1785 of file JFitK40.hh.

◆ successor

double JCALIBRATE::JFit::successor
private

Definition at line 1786 of file JFitK40.hh.

◆ previous

JModel JCALIBRATE::JFit::previous
private

Definition at line 1787 of file JFitK40.hh.

◆ h

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

Definition at line 1788 of file JFitK40.hh.


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