Evaluation of fit.
1488 {
1489 using namespace std;
1490 using namespace JPP;
1491
1493
1494
1496
1499
1500
1501
1502
1503 const struct M_t {
1505 {
1506 R =
model.getIndex(&JK40Parameters_t::R);
1507 p1 =
model.getIndex(&JK40Parameters_t::p1);
1508 p2 =
model.getIndex(&JK40Parameters_t::p2);
1509 p3 =
model.getIndex(&JK40Parameters_t::p3);
1510 p4 =
model.getIndex(&JK40Parameters_t::p4);
1511 cc =
model.getIndex(&JK40Parameters_t::cc);
1512 }
1513
1514 int R;
1516 int p2;
1517 int p3;
1518 int p4;
1519 int cc;
1520
1522
1523
1524
1525
1526 struct I_t {
1527 I_t(
const JModel& model,
const int pmt) :
1532 {
1533 const int index =
model.getIndex(pmt);
1534
1535 int N = 0;
1536
1537 if (
model.parameters[pmt].QE .isFree()) { QE = index + N; ++N; }
1538 if (
model.parameters[pmt].TTS.isFree()) { TTS = index + N; ++N; }
1539 if (
model.parameters[pmt].t0 .isFree()) { t0 = index + N; ++N; }
1540 if (
model.parameters[pmt].bg .isFree()) { bg = index + N; ++N; }
1541 }
1542
1543 int QE;
1544 int TTS;
1545 int t0;
1546 int bg;
1547 };
1548
1549
1551
1553
1554 for (data_type::const_iterator ix =
data.begin(); ix !=
data.end(); ++ix) {
1555
1556 const pair_type& pair = ix->first;
1557
1558 if (value.parameters[pair.first ].status &&
1559 value.parameters[pair.second].status) {
1560
1561 const real_type& real = value.getReal(pair);
1562
1563 const JGauss gauss(real.t0, real.sigma, real.signal);
1564
1565 const double R1 = value.getValue (real.ct);
1566 const JK40Parameters_t& R1p = value.getGradient(real.ct);
1567
1568 const std::pair<I_t, I_t> PMT(I_t(value, pair.first),
1569 I_t(value, pair.second));
1570
1571 for (const rate_type& iy : ix->second) {
1572
1573 const double R2 = gauss.getValue (iy.dt_ns);
1574 const JGauss& R2p = gauss.getGradient(iy.dt_ns);
1575
1576 const double R = real.background + R1 * (value.cc() + R2);
1577 const double u = (iy.value - R) / iy.error;
1578 const double w = -estimator->getPsi(u) / iy.error;
1579
1580 successor += estimator->getRho(u);
1581
1582 buffer.clear();
1583
1584 if (M.R != INVALID_INDEX) { buffer.push_back({M.R, w * (value.cc() + R2) * R1p.R () * value.R .getDerivative()}); }
1585 if (M.p1 != INVALID_INDEX) { buffer.push_back({M.p1, w * (value.cc() + R2) * R1p.p1() * value.p1.getDerivative()}); }
1586 if (M.p2 != INVALID_INDEX) { buffer.push_back({M.p2, w * (value.cc() + R2) * R1p.p2() * value.p2.getDerivative()}); }
1587 if (M.p3 != INVALID_INDEX) { buffer.push_back({M.p3, w * (value.cc() + R2) * R1p.p3() * value.p3.getDerivative()}); }
1588 if (M.p4 != INVALID_INDEX) { buffer.push_back({M.p4, w * (value.cc() + R2) * R1p.p4() * value.p4.getDerivative()}); }
1589 if (M.cc != INVALID_INDEX) { buffer.push_back({M.cc, w * R1 * R1p.cc() * value.cc.getDerivative()}); }
1590
1591 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()}); }
1592 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()}); }
1593 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}); }
1594 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}); }
1595 if (PMT.first .t0 != INVALID_INDEX) { buffer.push_back({PMT.first .t0, w * R1 * R2p.mean * value.parameters[pair.first ].t0 .getDerivative() * +1.0}); }
1596 if (PMT.second.t0 != INVALID_INDEX) { buffer.push_back({PMT.second.t0, w * R1 * R2p.mean * value.parameters[pair.second].t0 .getDerivative() * -1.0}); }
1597 if (PMT.first .bg != INVALID_INDEX) { buffer.push_back({PMT.first .bg, w * value.parameters[pair.first ].bg .getDerivative()}); }
1598 if (PMT.second.bg != INVALID_INDEX) { buffer.push_back({PMT.second.bg, w * value.parameters[pair.second].bg .getDerivative()}); }
1599
1600 for (buffer_type::const_iterator row = buffer.begin(); row != buffer.end(); ++row) {
1601
1602 Y[row->first] += row->second;
1603
1604 V[row->first][row->first] += row->second * row->second;
1605
1606 for (buffer_type::const_iterator col = buffer.begin(); col != row; ++col) {
1607 V[row->first][col->first] += row->second * col->second;
1608 V[col->first][row->first] = V[row->first][col->first];
1609 }
1610 }
1611 }
1612 }
1613 }
1614 }
static const int INVALID_INDEX
invalid index
Model for fit to acoustics data.
Auxiliary data structure for derived quantities of a given PMT pair.
JMatrixND & reset()
Set matrix to the null matrix.