1 #ifndef __JCALIBRATE_JGANDALFK40__
2 #define __JCALIBRATE_JGANDALFK40__
35 namespace JCALIBRATE {}
36 namespace JPP {
using namespace JCALIBRATE; }
38 namespace JCALIBRATE {
100 public std::map<pair_type, std::vector<rate_type> >
108 public JMath<JParameter_t>
170 set(
get() + parameter.
get());
184 set(
get() - parameter.
get());
303 void set(
const double value)
319 void fix(
const double value)
357 operator double()
const
386 return in >>
object.value;
401 out <<
FIXED(12,6) <<
object.get() <<
' '
402 << setw(5) << (
object.isFixed() ?
"fixed" :
" ") <<
' ';
405 out <<
FIXED(12,6) <<
object.value <<
' ';
406 out <<
FIXED(12,6) <<
object.range.getLowerLimit() <<
' '
407 <<
FIXED(12,6) <<
object.range.getUpperLimit();
425 static constexpr
double QE_MIN = 0.0;
426 static constexpr
double QE_MAX = 2.0;
427 static constexpr
double TTS_NS = 2.0;
465 parameters.
QE .
set(1.0);
466 parameters.
TTS.
set(TTS_NS);
467 parameters.
t0 .
set(0.0);
468 parameters.
bg .
set(0.0);
481 return ((QE.
isFree() ? 1 : 0) +
482 (TTS.isFree() ? 1 : 0) +
496 if (!(this->*p).
isFree()) {
550 out <<
"QE " <<
FIXED(7,3) <<
object.QE << endl;
551 out <<
"TTS " <<
FIXED(7,3) <<
object.TTS << endl;
552 out <<
"t0 " <<
FIXED(7,3) <<
object.t0 << endl;
553 out <<
"bg " <<
FIXED(7,3) <<
object.bg << endl;
629 parameters.
R .
set(18.460546);
630 parameters.
p1.
set( 3.0767);
631 parameters.
p2.
set(-1.2078);
632 parameters.
p3.
set( 0.9905);
633 parameters.
p4.
set( 0.9379);
634 parameters.
cc.
set( 0.0);
669 return ((
R .
isFree() ? 1 : 0) +
670 (
p1.isFree() ? 1 : 0) +
671 (
p2.isFree() ? 1 : 0) +
672 (
p3.isFree() ? 1 : 0) +
673 (
p4.isFree() ? 1 : 0) +
674 (cc.isFree() ? 1 : 0));
686 if (!(this->*p).
isFree()) {
726 const double ct2 = ct * ct;
728 if (
R .
isFree()) { gradient.R = rate /
R; }
729 if (
p1.isFree()) { gradient.p1 = rate * ct -
rate; }
730 if (
p2.isFree()) { gradient.p2 = rate * ct2 -
rate; }
731 if (
p3.isFree()) { gradient.p3 = rate * ct2 * ct -
rate; }
732 if (
p4.isFree()) { gradient.p4 = rate * ct2 * ct2 -
rate; }
733 if (cc.isFree()) { gradient.cc =
rate; }
792 setK40Parameters(parameters);
798 for (JTDC_t::const_iterator
i = TDC.first;
i != TDC.second; ++
i) {
799 this->parameters[
i->second].t0.fix();
819 setK40Parameters(parameters);
1019 for (
int i = 0;
i != pmt; ++
i) {
1047 return this->sigmaK40_ns;
1058 this->sigmaK40_ns =
sigma;
1070 real.
ct =
JPP::getDot((*
this)[pair.first].getDirection(), (*this)[pair.second].getDirection());
1072 real.t0 = (pair.first == this->index ? -this->
parameters[pair.second].t0() :
1073 pair.second == this->index ? +this->
parameters[pair.first ].t0() :
1078 this->getSigmaK40() * this->getSigmaK40());
1097 using namespace std;
1098 using namespace JPP;
1102 const JGauss
gauss(real.t0, real.sigma, real.signal);
1105 const double R2 = gauss.getValue(dt_ns);
1107 return real.background +
R1 * (cc() + R2);
1119 using namespace std;
1120 using namespace JPP;
1125 const double R2 = real.signal;
1127 return real.background +
R1 * (cc() + R2);
1140 using namespace std;
1142 out <<
"Module " << setw(10) <<
object.getID() << endl;
1143 out <<
"option " <<
object.option << endl;
1144 out <<
"index " <<
object.index << endl;
1145 out <<
"Rate [Hz] " <<
FIXED(12,6) <<
object.R << endl;
1146 out <<
"p1 " <<
FIXED(12,6) <<
object.p1 << endl;
1147 out <<
"p2 " <<
FIXED(12,6) <<
object.p2 << endl;
1148 out <<
"p3 " <<
FIXED(12,6) <<
object.p3 << endl;
1149 out <<
"p4 " <<
FIXED(12,6) <<
object.p4 << endl;
1150 out <<
"cc " <<
FIXED(12,6) <<
object.cc << endl;
1153 out <<
"PMT[" <<
FILL(2,
'0') <<
i <<
FILL() <<
"]." <<
object.parameters[
i].status << endl <<
object.parameters[
i];
1163 double sigmaK40_ns = 0.54;
1195 using namespace JPP;
1209 using namespace std;
1210 using namespace JPP;
1215 const size_t N =
value.getN();
1223 for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1227 if (
value.parameters[pair.first ].status &&
1228 value.parameters[pair.second].status) {
1230 ndf += ix->second.size();
1234 ndf -=
value.getN();
1237 lambda = LAMBDA_MIN;
1239 double precessor = numeric_limits<double>::max();
1241 for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
1243 DEBUG(
"step: " << numberOfIterations << endl);
1247 DEBUG(
"lambda: " <<
FIXED(12,5) << lambda << endl);
1248 DEBUG(
"chi2: " <<
FIXED(12,3) << successor << endl);
1250 if (successor < precessor) {
1252 if (numberOfIterations != 0) {
1254 if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
1255 return { successor / estimator->getRho(1.0), ndf };
1258 if (lambda > LAMBDA_MIN) {
1259 lambda /= LAMBDA_DOWN;
1263 precessor = successor;
1269 lambda *= LAMBDA_UP;
1271 if (lambda > LAMBDA_MAX) {
1272 return { precessor / estimator->getRho(1.0), ndf };
1282 if (
value.R .isFree()) { cout <<
"R " <<
FIXED(12,5) <<
Y[row] << endl; ++row; }
1283 if (
value.p1.isFree()) { cout <<
"p1 " <<
FIXED(12,5) <<
Y[row] << endl; ++row; }
1284 if (
value.p2.isFree()) { cout <<
"p2 " <<
FIXED(12,5) <<
Y[row] << endl; ++row; }
1285 if (
value.p3.isFree()) { cout <<
"p3 " <<
FIXED(12,5) <<
Y[row] << endl; ++row; }
1286 if (
value.p4.isFree()) { cout <<
"p4 " <<
FIXED(12,5) <<
Y[row] << endl; ++row; }
1287 if (
value.cc.isFree()) { cout <<
"cc " <<
FIXED(12,3) <<
Y[row] << endl; ++row; }
1290 if (
value.parameters[pmt].QE .isFree()) { cout <<
"PMT[" << setw(2) << pmt <<
"].QE " <<
FIXED(12,5) <<
Y[row] << endl; ++row; }
1291 if (
value.parameters[pmt].TTS.isFree()) { cout <<
"PMT[" << setw(2) << pmt <<
"].TTS " <<
FIXED(12,5) <<
Y[row] << endl; ++row; }
1292 if (
value.parameters[pmt].t0 .isFree()) { cout <<
"PMT[" << setw(2) << pmt <<
"].t0 " <<
FIXED(12,5) <<
Y[row] << endl; ++row; }
1293 if (
value.parameters[pmt].bg .isFree()) { cout <<
"PMT[" << setw(2) << pmt <<
"].bg " <<
FIXED(12,5) <<
Y[row] << endl; ++row; }
1299 for (
size_t i = 0;
i !=
N; ++
i) {
1301 if (
V(
i,
i) < PIVOT) {
1305 h[
i] = 1.0 / sqrt(
V(
i,
i));
1310 for (
size_t i = 0;
i !=
N; ++
i) {
1311 for (
size_t j = 0;
j !=
i; ++
j) {
1312 V(
j,
i) *= h[
i] * h[
j];
1317 for (
size_t i = 0;
i !=
N; ++
i) {
1318 V(
i,
i) = 1.0 + lambda;
1323 for (
size_t col = 0; col !=
N; ++col) {
1330 catch (
const exception& error) {
1332 ERROR(
"JGandalf: " << error.what() << endl <<
V << endl);
1339 const double factor = 2.0;
1343 if (
value.R .isFree()) {
value.R -= factor * h[row] *
Y[row]; ++row; }
1344 if (
value.p1.isFree()) {
value.p1 -= factor * h[row] *
Y[row]; ++row; }
1345 if (
value.p2.isFree()) {
value.p2 -= factor * h[row] *
Y[row]; ++row; }
1346 if (
value.p3.isFree()) {
value.p3 -= factor * h[row] *
Y[row]; ++row; }
1347 if (
value.p4.isFree()) {
value.p4 -= factor * h[row] *
Y[row]; ++row; }
1348 if (
value.cc.isFree()) {
value.cc -= factor * h[row] *
Y[row]; ++row; }
1351 if (
value.parameters[pmt].QE .isFree()) {
value.parameters[pmt].QE -= factor * h[row] *
Y[row]; ++row; }
1352 if (
value.parameters[pmt].TTS.isFree()) {
value.parameters[pmt].TTS -= factor * h[row] *
Y[row]; ++row; }
1353 if (
value.parameters[pmt].t0 .isFree()) {
value.parameters[pmt].t0 -= factor * h[row] *
Y[row]; ++row; }
1354 if (
value.parameters[pmt].bg .isFree()) {
value.parameters[pmt].bg -= factor * h[row] *
Y[row]; ++row; }
1358 return { precessor / estimator->getRho(1.0), ndf };
1362 static constexpr
int MAXIMUM_ITERATIONS = 10000;
1363 static constexpr
double EPSILON = 1.0e-4;
1364 static constexpr
double LAMBDA_MIN = 0.01;
1365 static constexpr
double LAMBDA_MAX = 100.0;
1366 static constexpr
double LAMBDA_UP = 10.0;
1367 static constexpr
double LAMBDA_DOWN = 10.0;
1386 using namespace std;
1387 using namespace JPP;
1430 const int index = model.
getIndex(pmt);
1451 for (data_type::const_iterator ix = data.begin(); ix != data.end(); ++ix) {
1455 if (
value.parameters[pair.first ].status &&
1456 value.parameters[pair.second].status) {
1458 const real_type& real =
value.getReal(pair);
1460 const JGauss
gauss(real.t0, real.sigma, real.signal);
1462 const double R1 =
value.getValue (real.ct);
1466 I_t(
value, pair.second));
1468 for (
const rate_type& iy : ix->second) {
1470 const double R2 = gauss.getValue (iy.
dt_ns);
1471 const JGauss& R2p = gauss.getGradient(iy.
dt_ns);
1473 const double R = real.background + R1 * (
value.cc() + R2);
1475 const double w = -estimator->getPsi(u) / iy.
error;
1477 successor += estimator->getRho(u);
1486 if (M.cc !=
INVALID_INDEX) { buffer.push_back({M.cc,
w * R1 * R1p.
cc() *
value.cc.getDerivative()}); }
1488 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()}); }
1489 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()}); }
1490 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}); }
1491 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}); }
1492 if (PMT.first .t0 !=
INVALID_INDEX) { buffer.push_back({PMT.first .t0,
w * R1 * R2p.mean *
value.parameters[pair.first ].t0 .getDerivative() * +1.0}); }
1493 if (PMT.second.t0 !=
INVALID_INDEX) { buffer.push_back({PMT.second.t0,
w * R1 * R2p.mean *
value.parameters[pair.second].t0 .getDerivative() * -1.0}); }
1494 if (PMT.first .bg !=
INVALID_INDEX) { buffer.push_back({PMT.first .bg,
w *
value.parameters[pair.first ].bg .getDerivative()}); }
1495 if (PMT.second.bg !=
INVALID_INDEX) { buffer.push_back({PMT.second.bg,
w *
value.parameters[pair.second].bg .getDerivative()}); }
1497 for (buffer_type::const_iterator row = buffer.begin(); row != buffer.end(); ++row) {
1499 Y[row->first] += row->second;
1501 V[row->first][row->first] += row->second * row->second;
1503 for (buffer_type::const_iterator col = buffer.begin(); col != row; ++col) {
1504 V[row->first][col->first] += row->second * col->second;
1505 V[col->first][row->first] =
V[row->first][col->first];
Data structure for measured coincidence rate of pair of PMTs.
result_type operator()(const data_type &data)
Fit.
JTOOLS::JRange< double > range_type
Type definition for range of parameter values.
size_t getN() const
Get number of fit parameters.
void set()
Set current value.
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
int getIndex(JParameter_t JPMTParameters_t::*p) const
Get index of parameter.
const real_type & getReal(const pair_type &pair) const
Get derived parameters.
int getIndex(JParameter_t JK40Parameters::*p) const
Get index of parameter.
JParameter_t t0
time offset [ns]
Auxiliary base class for aritmetic operations of derived class types.
JK40Parameters()
Default constructor.
Data structure for a composite optical module.
void setSigmaK40(const double sigma)
Set intrinsic K40 arrival time spread.
JParameter_t R
maximal coincidence rate [Hz]
bool hasFixedTimeOffset() const
Check if time offset is fixed.
static const int INVALID_INDEX
invalid index
size_t getN() const
Get number of fit parameters.
Interface for maximum likelihood estimator (M-estimator).
JParameter_t & add(const JParameter_t ¶meter)
Add parameter.
JParameter_t & sub(const JParameter_t ¶meter)
Subtract parameter.
JK40Parameters_t()
Default constructor.
JParameter_t & negate()
Negate parameter.
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
const JK40Parameters_t & getGradient(const double ct) const
Get gradient.
static const JPBS_t PMT(3, 4, 2, 3)
PBS of photo-multiplier tube (PMT)
void setIndex()
Set index of PMT used for fixed time offset.
void setK40Parameters(const JK40Parameters_t ¶meters)
Set K40 parameters.
int index
index of PMT used for fixed time offset
void set(const double value)
Set value.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
then set_variable singlesRate set_variable doublesRate set_variable numberOfSlices echo Generating random background echo Singles rate
double getDerivative() const
Get derivative of value.
double getValue(const pair_type &pair) const
Get K40 coincidence rate.
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
int getIndex() const
Get index of PMT used for fixed time offset.
JK40Parameters_t gradient
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
Auxiliary data structure for floating point format specification.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
JModel()
Default constructor.
double background
combined background
JPMTParameters_t()
Default constructor.
std::vector< JHitW0 > buffer_type
hits
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
double getValue(const double ct) const
Get K40 coincidence rate as a function of cosine angle between PMT axes.
friend std::ostream & operator<<(std::ostream &out, const JModel &object)
Write model parameters to output stream.
double dt_ns
time difference [ns]
JFit(const int option, const int debug)
Constructor.
void setOption(const int option)
Set fit option.
JParameter_t bg
background [Hz/ns]
void fix()
Fix current value.
JParameter_t & operator=(double value)
Assignment operator.
Fit parameters for two-fold coincidence rate due to K40.
JParameter_t QE
relative quantum efficiency [unit]
JParameter_t p3
3rd order angle dependence coincidence rate
void evaluate(const data_type &data)
Evaluation of fit.
double ct
cosine angle between PMT axes
double error
error of rate [Hz/ns]
fit parameters of PMTs and background
estimator_type estimator
M-Estimator function.
JParameter_t()
Default constructor.
rate_type(double dt_ns, double value, double error)
Constructor.
fit parameters of PMTs and angular dependence of K40 rate
double t0
time offset [ns]
const JK40Parameters & getK40Parameters() const
Get K40 parameters.
friend std::ostream & operator<<(std::ostream &out, const JParameter_t &object)
Write parameter to output stream.
int getIndex()
Get index for user I/O manipulation.
void fix(const double value)
Fix value.
double sigma
total width [ns]
JParameter_t TTS
transition-time spread [ns]
static const JPMTParameters_t & getInstance()
Get default values.
Auxiliary data structure for sequence of same character.
Fit parameters for single PMT.
JParameter_t p4
4th order angle dependence coincidence rate
JParameter_t & mul(const JParameter_t &first, const JParameter_t &second)
Scale parameter.
fit parameters of PMTs with QE fixed
JOption_t getOption() const
Get fit option.
then JCookie sh JDataQuality D $DETECTOR_ID R
size_t getN() const
Get number of fit parameters.
friend std::istream & operator>>(std::istream &in, JParameter_t &object)
Read parameter from input stream.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
bool isFree() const
Check if parameter is free.
Auxiliary class to define a range between two values.
JParameter_t p2
2nd order angle dependence coincidence rate
Auxiliary data structure for derived quantities of a given PMT pair.
fit parameters of K40 rate and TTSs of PMTs
then set_variable NUMBER_OF_TESTS else set_variable NUMBER_OF_TESTS fi function gauss()
double getFixedTimeOffset() const
Get time offset.
JParameter_t(const double value, const range_type &range=range_type::DEFAULT_RANGE())
Constructor.
double signal
combined signal
PMT combinatorics for optical module.
JParameter_t & mul(const double factor)
Scale parameter.
double getValue(const pair_type &pair, const double dt_ns) const
Get K40 coincidence rate.
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
static const JK40Parameters & getInstance()
Get default values.
JModel(const JModule &module, const JK40Parameters ¶meters, const JTDC_t::range_type &TDC, const int option)
Constructor.
JModel(const JModule &module, const JK40Parameters ¶meters)
Constructor.
Base class for data structures with artithmetic capabilities.
int getIndex(int pmt) const
Get index of parameter.
double getSigmaK40() const
Get intrinsic K40 arrival time spread.
Exception for accessing a value in a collection that is outside of its range.
int getIndex(int pmt, JParameter_t JPMTParameters_t::*p) const
Get index of parameter.
Fit parameters for two-fold coincidence rate due to K40.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
KM3NeT DAQ constants, bit handling, etc.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
bool isBound() const
Check if parameter is bound.
bool isFixed() const
Check if parameter is fixed.
std::shared_ptr< JMEstimator > estimator_type
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
rate_type()
Default constructor.
void disable()
Disable PMT.
JPMTParameters_t parameters[NUMBER_OF_PMTS]
JParameter_t & div(const double factor)
Scale parameter.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` source JAcousticsToolkit.sh typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
double get() const
Get value.
friend std::ostream & operator<<(std::ostream &out, const JPMTParameters_t &object)
Write PMT parameters to output stream.
Maximum likelihood estimator (M-estimators).
double operator()() const
Type conversion operator.
double value
value of rate [Hz/ns]
#define DEBUG(A)
Message macros.
JParameter_t p1
1st order angle dependence coincidence rate
JParameter_t cc
fraction of signal correlated background
Data structure for optical module.
Auxiliary class for fit parameter with optional limits.