1#ifndef __JCALIBRATE_JGANDALFK40__
2#define __JCALIBRATE_JGANDALFK40__
42 using KM3NETDAQ::NUMBER_OF_PMTS;
102 public std::map<pair_type, std::vector<rate_type> >
110 public JMath<JParameter_t>
322 void setRange(
const double xmin,
const double xmax)
324 const double x =
get();
391 operator double()
const
420 return in >>
object.value;
435 out <<
FIXED(12,6) <<
object.get() <<
' '
436 << setw(5) << (
object.isFixed() ?
"fixed" :
" ") <<
' ';
439 out <<
FIXED(12,6) <<
object.value <<
' ';
440 out <<
FIXED(12,6) <<
object.range.getLowerLimit() <<
' '
441 <<
FIXED(12,6) <<
object.range.getUpperLimit();
485 parameters.
QE .
set(1.0);
487 parameters.
t0 .
set(0.0);
488 parameters.
bg .
set(0.0);
515 if (
QE .isFree()) {
QE .
set(parameters.
QE); }
517 if (
t0 .isFree()) {
t0 .
set(parameters.
t0); }
518 if (
bg .isFree()) {
bg .
set(parameters.
bg); }
529 return ((
QE. isFree() ? 1 : 0) +
544 if (!(this->*p).isFree()) {
598 out <<
"QE " <<
FIXED(7,3) <<
object.QE << endl;
599 out <<
"TTS " <<
FIXED(7,3) <<
object.TTS << endl;
600 out <<
"t0 " <<
FIXED(7,3) <<
object.t0 << endl;
601 out <<
"bg " <<
FIXED(7,3) <<
object.bg << endl;
699 parameters.
R .
set(18.460546);
700 parameters.
p1.
set( 3.0767);
701 parameters.
p2.
set(-1.2078);
702 parameters.
p3.
set( 0.9905);
703 parameters.
p4.
set( 0.9379);
704 parameters.
cc.
set( 0.0);
717 return ((
R .isFree() ? 1 : 0) +
734 if (!(this->*p).isFree()) {
774 const double ct2 = ct * ct;
811 out <<
"Rate [Hz] " <<
FIXED(12,6) <<
object.R << endl;
812 out <<
"p1 " <<
FIXED(12,6) <<
object.p1 << endl;
813 out <<
"p2 " <<
FIXED(12,6) <<
object.p2 << endl;
814 out <<
"p3 " <<
FIXED(12,6) <<
object.p3 << endl;
815 out <<
"p4 " <<
FIXED(12,6) <<
object.p4 << endl;
816 out <<
"cc " <<
FIXED(12,6) <<
object.cc << endl;
818 for (
int i = 0; i != NUMBER_OF_PMTS; ++i) {
819 out <<
"PMT[" <<
FILL(2,
'0') << i <<
FILL() <<
"]." <<
object.parameters[i].status << endl <<
object.parameters[i];
879 for (
int i = 0; i != NUMBER_OF_PMTS; ++i) {
883 for (JTDC_t::const_iterator i = TDC.first; i != TDC.second; ++i) {
887 this->parameters[i->second].t0.fix();
891 for (
int i = 0; i != NUMBER_OF_PMTS; ++i) {
892 this->parameters[i].t0.fix();
916 for (
int i = 0; i != NUMBER_OF_PMTS; ++i) {
950 for (
int i = 0; i != NUMBER_OF_PMTS; ++i) {
960 for (
int i = 0; i != NUMBER_OF_PMTS; ++i) {
984 for (
int i = 0; i != NUMBER_OF_PMTS; ++i) {
995 for (
int i = 0; i != NUMBER_OF_PMTS; ++i) {
1033 for (
int i = 0; i != NUMBER_OF_PMTS; ++i) {
1062 for (
int i = 0; i != NUMBER_OF_PMTS; ++i) {
1088 for (
int i = 0; i != NUMBER_OF_PMTS; ++i) {
1106 for (
int i = 0; i != pmt; ++i) {
1160 pair.second == this->index ? +this->parameters[
pair.first ].t0() :
1161 this->parameters[
pair.first].t0() - this->parameters[
pair.second].t0());
1164 this->parameters[
pair.second].TTS() * this->parameters[
pair.second].TTS() +
1165 this->getSigmaK40() * this->getSigmaK40());
1186 using namespace std;
1187 using namespace JPP;
1194 const double R2 = gauss.getValue(dt_ns);
1208 using namespace std;
1209 using namespace JPP;
1229 using namespace std;
1231 out <<
"Module " << setw(10) <<
object.getID() << endl;
1232 out <<
"option " <<
object.option << endl;
1233 out <<
"index " <<
object.index << endl;
1235 out << static_cast<const JModel_t&>(
object);
1274 using namespace JPP;
1288 using namespace std;
1289 using namespace JPP;
1300 double xmax = numeric_limits<double>::lowest();
1301 double xmin = numeric_limits<double>::max();
1305 for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1312 ndf += ix->second.size();
1314 for (
const rate_type& iy : ix->second) {
1324 return { 0.0, ndf };
1327 for (
int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1336 double precessor = numeric_limits<double>::max();
1378 if (
debug >= debug_t) {
1389 for (
int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1399 for (
size_t i = 0; i != N; ++i) {
1405 h[i] = 1.0 / sqrt(
V(i,i));
1410 for (
size_t i = 0; i != N; ++i) {
1411 for (
size_t j = 0; j != i; ++j) {
1412 V(j,i) *=
h[i] *
h[j];
1417 for (
size_t i = 0; i != N; ++i) {
1423 for (
size_t col = 0; col != N; ++col) {
1430 catch (
const exception&
error) {
1432 ERROR(
"JGandalf: " <<
error.what() << endl <<
V << endl);
1439 const double factor = 2.0;
1450 for (
int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
1460 return { precessor /
estimator->getRho(1.0), ndf };
1470 static constexpr double PIVOT = std::numeric_limits<double>::epsilon();
1489 using namespace std;
1490 using namespace JPP;
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);
1527 I_t(
const JModel& model,
const int pmt) :
1533 const int index = model.getIndex(pmt);
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; }
1554 for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1563 const JGauss gauss(real.t0, real.sigma, real.signal);
1571 for (
const rate_type& iy : ix->second) {
1573 const double R2 = gauss.getValue (iy.
dt_ns);
1574 const JGauss& R2p = gauss.getGradient(iy.
dt_ns);
1576 const double R = real.background +
R1 * (
value.
cc() + R2);
1600 for (buffer_type::const_iterator row = buffer.begin(); row != buffer.end(); ++row) {
1602 Y[row->first] += row->second;
1604 V[row->first][row->first] += row->second * row->second;
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];
1624 using namespace std;
1633 catch (
const exception&
error) {}
1635#define SQRT(X) (X >= 0.0 ? sqrt(X) : std::numeric_limits<double>::max())
1646 for (
int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
JDAQPMTIdentifier PMT
Command line options.
KM3NeT DAQ constants, bit handling, etc.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Maximum likelihood estimator (M-estimators).
Base class for data structures with artithmetic capabilities.
General purpose messaging.
#define DEBUG(A)
Message macros.
Data structure for optical module.
Auxiliary class to define a range between two values.
std::shared_ptr< JMEstimator > estimator_type
static constexpr double LAMBDA_MIN
minimal value control parameter
static constexpr double LAMBDA_DOWN
multiplication factor control parameter
result_type operator()(const data_type &data)
Fit.
void seterr(const data_type &data)
Set errors.
static constexpr double LAMBDA_MAX
maximal value control parameter
static constexpr double LAMBDA_UP
multiplication factor control parameter
static constexpr double EPSILON
maximal distance to minimum.
JFit(const int option, const int debug)
Constructor.
void evaluate(const data_type &data)
Evaluation of fit.
static constexpr int MAXIMUM_ITERATIONS
maximal number of iterations.
static constexpr double PIVOT
minimal value diagonal element of matrix
estimator_type estimator
M-Estimator function.
Auxiliary class for fit parameter with optional limits.
JParameter_t & mul(const double factor)
Scale parameter.
void set(const double value)
Set value.
void fix()
Fix current value.
JParameter_t & sub(const JParameter_t ¶meter)
Subtract parameter.
JParameter_t & operator=(double value)
Assignment operator.
bool isFree() const
Check if parameter is free.
friend std::ostream & operator<<(std::ostream &out, const JParameter_t &object)
Write parameter to output stream.
friend std::istream & operator>>(std::istream &in, JParameter_t &object)
Read parameter from input stream.
JParameter_t & div(const double factor)
Scale parameter.
JParameter_t & mul(const JParameter_t &first, const JParameter_t &second)
Scale parameter.
double operator()() const
Type conversion operator.
void set()
Set current value.
JParameter_t(const double value, const range_type &range=range_type::DEFAULT_RANGE())
Constructor.
void setRange(const double xmin, const double xmax)
Set range.
JParameter_t & negate()
Negate parameter.
JParameter_t()
Default constructor.
bool atLimit(const double precision) const
Check if parameter is at limit;.
JTOOLS::JRange< double > range_type
Type definition for range of parameter values.
double getDerivative() const
Get derivative of value.
JParameter_t & add(const JParameter_t ¶meter)
Add parameter.
void fix(const double value)
Fix value.
double get() const
Get value.
bool isBound() const
Check if parameter is bound.
bool isFixed() const
Check if parameter is fixed.
Data structure for a composite optical module.
Exception for accessing a value in a collection that is outside of its range.
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
Auxiliary classes and methods for PMT calibration.
static const int INVALID_INDEX
invalid index
@ FIT_PMTS_QE_FIXED_t
fit parameters of PMTs with QE fixed
@ FIT_PMTS_AND_ANGULAR_DEPENDENCE_t
fit parameters of PMTs and angular dependence of K40 rate
@ FIT_MODEL_t
fit parameters of K40 rate and TTSs of PMTs
@ FIT_PMTS_AND_BACKGROUND_t
fit parameters of PMTs and background
@ FIT_PMTS_t
fit parameters of PMTs
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary data structure for sequence of same character.
Auxiliary data structure for floating point format specification.
PMT combinatorics for optical module.
Fit parameters for two-fold coincidence rate due to K40.
JParameter_t R
maximal coincidence rate [Hz]
JParameter_t p1
1st order angle dependence coincidence rate
JParameter_t p2
2nd order angle dependence coincidence rate
JParameter_t p3
3rd order angle dependence coincidence rate
const JK40Parameters_t & getK40Parameters() const
Get K40 parameters.
JParameter_t p4
4th order angle dependence coincidence rate
JParameter_t cc
fraction of signal correlated background
JK40Parameters_t()
Default constructor.
void setK40Parameters(const JK40Parameters_t ¶meters)
Set K40 parameters.
Fit parameters for two-fold coincidence rate due to K40.
size_t getN() const
Get number of fit parameters.
const JK40Parameters_t & getGradient(const double ct) const
Get gradient.
JK40Parameters_t gradient
static const JK40Parameters & getInstance()
Get default values.
int getIndex(JParameter_t JK40Parameters::*p) const
Get index of parameter.
double getValue(const double ct) const
Get K40 coincidence rate as a function of cosine angle between PMT axes.
JK40Parameters()
Default constructor.
Auxiliary data structure for derived quantities of a given PMT pair.
double signal
combined signal
double sigma
total width [ns]
double cc
correlated background
double background
combined background
double t0
time offset [ns]
double ct
cosine angle between PMT axes
friend std::ostream & operator<<(std::ostream &out, const JModel_t &object)
Write model parameters to output stream.
JPMTParameters_t parameters[NUMBER_OF_PMTS]
JModel()
Default constructor.
int getIndex(int pmt, JParameter_t JPMTParameters_t::*p) const
Get index of parameter.
friend std::ostream & operator<<(std::ostream &out, const JModel &object)
Write model parameters to output stream.
int getIndex(int pmt) const
Get index of parameter.
size_t getN() const
Get number of fit parameters.
double getValue(const pair_type &pair) const
Get K40 coincidence rate.
double sigmaK40_ns
intrinsic K40 arrival time spread [ns]
JOption_t getOption() const
Get fit option.
double getFixedTimeOffset() const
Get time offset.
void setSigmaK40(const double sigma)
Set intrinsic K40 arrival time spread.
int getIndex() const
Get index of PMT used for fixed time offset.
double getSigmaK40() const
Get intrinsic K40 arrival time spread.
void setOption(const int option)
Set fit option.
const real_type & getReal(const pair_type &pair) const
Get derived parameters.
JModel(const JModule &module, const JK40Parameters ¶meters)
Constructor.
double getValue(const pair_type &pair, const double dt_ns) const
Get K40 coincidence rate.
void setIndex()
Set index of PMT used for fixed time offset.
JModel(const JModule &module, const JK40Parameters ¶meters, const JTDC_t::range_type &TDC, const int option)
Constructor.
bool hasFixedTimeOffset() const
Check if time offset is fixed.
int index
index of PMT used for fixed time offset
Fit parameters for single PMT.
static constexpr double QE_MIN
minimal QE
friend std::ostream & operator<<(std::ostream &out, const JPMTParameters_t &object)
Write PMT parameters to output stream.
JParameter_t t0
time offset [ns]
static constexpr double TTS_NS
start value transition-time spread [ns]
JParameter_t TTS
transition-time spread [ns]
void disable()
Disable PMT.
size_t getN() const
Get number of fit parameters.
JPMTParameters_t()
Default constructor.
void set(const JPMTParameters_t ¶meters)
Set parameters that are free to given values.
JParameter_t bg
background [Hz/ns]
static constexpr double QE_MAX
maximal QE
int getIndex(JParameter_t JPMTParameters_t::*p) const
Get index of parameter.
static const JPMTParameters_t & getInstance()
Get default values.
JParameter_t QE
relative quantum efficiency [unit]
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Data structure for measured coincidence rate of pair of PMTs.
rate_type(double dt_ns, double value, double error)
Constructor.
double error
error of rate [Hz/ns]
double value
value of rate [Hz/ns]
rate_type()
Default constructor.
double dt_ns
time difference [ns]
Interface for maximum likelihood estimator (M-estimator).
Auxiliary base class for aritmetic operations of derived class types.
void resize(const size_t size)
Resize matrix.
JMatrixND & reset()
Set matrix to the null matrix.
void solve(JVectorND_t &u)
Get solution of equation A x = b.
void invert()
Invert matrix according LDU decomposition.