Jpp master_rocky-44-g75b7c4f75
the software that should make you happy
Loading...
Searching...
No Matches
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.
 
result_type operator() (const data_type &data)
 Fit.
 

Public Attributes

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

Static Public Attributes

static constexpr int MAXIMUM_ITERATIONS = 10000
 maximal number of iterations.
 
static constexpr double EPSILON = 1.0e-5
 maximal distance to minimum.
 
static constexpr double LAMBDA_MIN = 0.01
 minimal value control parameter
 
static constexpr double LAMBDA_MAX = 100.0
 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.
 

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

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 }
JCALIBRATE::JFit::estimator
estimator_type estimator
M-Estimator function.
Definition JFitK40.hh:1375
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43

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
1245 for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
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 }
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
ERROR
#define ERROR(A)
Definition JMessage.hh:66
JCALIBRATE::JFit::h
std::vector< double > h
Definition JFitK40.hh:1520
JCALIBRATE::JFit::LAMBDA_MIN
static constexpr double LAMBDA_MIN
minimal value control parameter
Definition JFitK40.hh:1368
JCALIBRATE::JFit::LAMBDA_DOWN
static constexpr double LAMBDA_DOWN
multiplication factor control parameter
Definition JFitK40.hh:1371
JCALIBRATE::JFit::LAMBDA_MAX
static constexpr double LAMBDA_MAX
maximal value control parameter
Definition JFitK40.hh:1369
JCALIBRATE::JFit::LAMBDA_UP
static constexpr double LAMBDA_UP
multiplication factor control parameter
Definition JFitK40.hh:1370
JCALIBRATE::JFit::V
JMATH::JMatrixNS V
Definition JFitK40.hh:1380
JCALIBRATE::JFit::EPSILON
static constexpr double EPSILON
maximal distance to minimum.
Definition JFitK40.hh:1367
JCALIBRATE::JFit::evaluate
void evaluate(const data_type &data)
Evaluation of fit.
Definition JFitK40.hh:1388
JCALIBRATE::JFit::MAXIMUM_ITERATIONS
static constexpr int MAXIMUM_ITERATIONS
maximal number of iterations.
Definition JFitK40.hh:1366
JCALIBRATE::JFit::PIVOT
static constexpr double PIVOT
minimal value diagonal element of matrix
Definition JFitK40.hh:1372
JCALIBRATE::JFit::Y
JMATH::JVectorND Y
Definition JFitK40.hh:1517
JCALIBRATE::JParameter_t::isFree
bool isFree() const
Check if parameter is free.
Definition JFitK40.hh:240
JTOOLS::j
int j
Definition JPolint.hh:792
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JCALIBRATE::JK40Parameters_t::R
JParameter_t R
maximal coincidence rate [Hz]
Definition JFitK40.hh:596
JCALIBRATE::JK40Parameters_t::p1
JParameter_t p1
1st order angle dependence coincidence rate
Definition JFitK40.hh:597
JCALIBRATE::JK40Parameters_t::p2
JParameter_t p2
2nd order angle dependence coincidence rate
Definition JFitK40.hh:598
JCALIBRATE::JK40Parameters_t::p3
JParameter_t p3
3rd order angle dependence coincidence rate
Definition JFitK40.hh:599
JCALIBRATE::JK40Parameters_t::p4
JParameter_t p4
4th order angle dependence coincidence rate
Definition JFitK40.hh:600
JCALIBRATE::JK40Parameters_t::cc
JParameter_t cc
fraction of signal correlated background
Definition JFitK40.hh:601
JCALIBRATE::JModel::getN
size_t getN() const
Get number of fit parameters.
Definition JFitK40.hh:1001
JCALIBRATE::JModel::parameters
JPMTParameters_t parameters[NUMBER_OF_PMTS]
Definition JFitK40.hh:1163
JCALIBRATE::JModel::setIndex
void setIndex()
Set index of PMT used for fixed time offset.
Definition JFitK40.hh:975
JCALIBRATE::JPMTParameters_t::t0
JParameter_t t0
time offset [ns]
Definition JFitK40.hh:564
JCALIBRATE::JPMTParameters_t::TTS
JParameter_t TTS
transition-time spread [ns]
Definition JFitK40.hh:563
JCALIBRATE::JPMTParameters_t::bg
JParameter_t bg
background [Hz/ns]
Definition JFitK40.hh:565
JCALIBRATE::JPMTParameters_t::status
bool status
status
Definition JFitK40.hh:561
JCALIBRATE::JPMTParameters_t::QE
JParameter_t QE
relative quantum efficiency [unit]
Definition JFitK40.hh:562
JMATH::JMatrixND::resize
void resize(const size_t size)
Resize matrix.
Definition JMatrixND.hh:446
JMATH::JMatrixNS::solve
void solve(JVectorND_t &u)
Get solution of equation A x = b.
Definition JMatrixNS.hh:308
JTOOLS::JCombinatorics::pair_type
Data structure for a pair of indices.
Definition JCombinatorics.hh:29

◆ 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
1451 typedef vector< pair<int, double> > buffer_type;
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 }
p1
TPaveText * p1
Definition JDrawModule3D.cc:35
JCALIBRATE::INVALID_INDEX
static const int INVALID_INDEX
invalid index
Definition JFitK40.hh:60
JTOOLS::if
if((p==this->begin() &&this->getDistance(x,(p++) ->getX()) > distance_type::precision)||(p==this->end() &&this->getDistance((--p) ->getX(), x) > distance_type::precision))
Template base class for polynomial interpolations with polynomial result.
Definition JPolint.hh:770
JACOUSTICS::JModel
Model for fit to acoustics data.
Definition JAcoustics/JModel.hh:480
JCALIBRATE::JModel::real_type
Auxiliary data structure for derived quantities of a given PMT pair.
Definition JFitK40.hh:763
JMATH::JMatrixND::reset
JMatrixND & reset()
Set matrix to the null matrix.
Definition JMatrixND.hh:459
JMATH::JVectorND::reset
void reset()
Reset.
Definition JVectorND.hh:45

Member Data Documentation

◆ MAXIMUM_ITERATIONS

int JCALIBRATE::JFit::MAXIMUM_ITERATIONS = 10000
staticconstexpr

maximal number of iterations.

Definition at line 1366 of file JFitK40.hh.

◆ EPSILON

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

maximal distance to minimum.

Definition at line 1367 of file JFitK40.hh.

◆ LAMBDA_MIN

double JCALIBRATE::JFit::LAMBDA_MIN = 0.01
staticconstexpr

minimal value control parameter

Definition at line 1368 of file JFitK40.hh.

◆ LAMBDA_MAX

double JCALIBRATE::JFit::LAMBDA_MAX = 100.0
staticconstexpr

maximal value control parameter

Definition at line 1369 of file JFitK40.hh.

◆ LAMBDA_UP

double JCALIBRATE::JFit::LAMBDA_UP = 10.0
staticconstexpr

multiplication factor control parameter

Definition at line 1370 of file JFitK40.hh.

◆ LAMBDA_DOWN

double JCALIBRATE::JFit::LAMBDA_DOWN = 10.0
staticconstexpr

multiplication factor control parameter

Definition at line 1371 of file JFitK40.hh.

◆ PIVOT

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:
Generated by doxygen 1.12.0