1 #ifndef __JACOUSTICS__JKATOOMBA__
2 #define __JACOUSTICS__JKATOOMBA__
9 #include "TMatrixTSym.h"
10 #include "TDecompSVD.h"
37 namespace JACOUSTICS {}
38 namespace JPP {
using namespace JACOUSTICS; }
40 namespace JACOUSTICS {
62 for (
T i = __begin; i != __end; ++i) {
114 template<
template<
class T>
class JMinimiser_t =
JType>
151 template<
class JPDF_t>
159 const double Vi = velocity.getInverseVelocity(D, hit.
getZ(), position.
getZ());
172 template<
class JPDF_t>
180 const double Vi = velocity.getInverseVelocity(D, hit.
getZ(), position.
getZ());
182 return hit.getValue() - D * Vi;
228 H_t&
mul(
const double factor)
300 template<
class JPDF_t>
303 const double toa_s = this->getToA(model, hit);
304 const double u = (toa_s - hit.getValue()) / hit.sigma;
306 return estimator->getRho(u);
335 static constexpr
double TOLERANCE = 1.0e-8;
364 using namespace JGEOMETRY;
369 unique_lock<mutex> lock(mtx);
371 value =
JModel(__begin, __end);
381 const size_t N = value.getN();
386 for (
T hit = __begin; hit != __end; ++hit) {
390 const double Vi = velocity.getInverseVelocity(hit->getDistance(position), hit->getZ(), position.
getZ());
392 const double h1 =
string.getHeight(hit->getFloor());
393 const JPosition3D p1 =
string.getPosition() - hit->getPosition();
395 const double y = hit->getValue() - Vi*ds;
396 const double W = sqrt(hit->getWeight());
399 H.tx = W * Vi * p1.
getX() * h1 / ds;
400 H.ty = W * Vi * p1.
getY() * h1 / ds;
402 i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
403 i.tx = value.getIndex(hit->getString(), &H_t::tx);
404 i.ty = value.getIndex(hit->getString(), &H_t::ty);
406 V(i.t1, i.t1) +=
H.t1 *
H.t1;
408 Y[i.t1] += W *
H.t1 * y;
410 if (hit->getFloor() != 0) {
414 V(i.t1, i.tx) +=
H.t1 *
H.tx;
V(i.t1, i.ty) +=
H.t1 *
H.ty;
415 V(i.tx, i.t1) +=
H.tx *
H.t1;
V(i.ty, i.t1) +=
H.ty *
H.t1;
417 V(i.tx, i.tx) +=
H.tx *
H.tx;
V(i.tx, i.ty) +=
H.tx *
H.ty;
418 V(i.ty, i.tx) +=
H.ty *
H.tx;
V(i.ty, i.ty) +=
H.ty *
H.ty;
420 Y[i.tx] += W *
H.tx * y;
421 Y[i.ty] += W *
H.ty * y;
428 TDecompSVD svd(V, TOLERANCE);
432 V = svd.Invert(status);
434 for (
size_t row = 0; row !=
N; ++row) {
435 for (
size_t col = 0; col !=
N; ++col) {
436 value[row] +=
V(row,col) * Y[col];
446 static std::mutex mtx;
487 template<
class JPDF_t>
490 const double toa_s = this->getToA(model, hit);
491 const double u = (toa_s - hit.getValue()) / hit.sigma;
493 return estimator->getRho(u);
513 switch (this->option) {
585 value.setOption(this->option);
587 const int N = value.getN();
595 result_type precessor = numeric_limits<double>::max();
597 for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
599 DEBUG(
"step: " << numberOfIterations << endl);
601 evaluate(__begin, __end);
603 DEBUG(
"lambda: " <<
FIXED(12,5) << lambda << endl);
604 DEBUG(
"chi2: " <<
FIXED(12,5) << successor << endl);
606 if (successor < precessor) {
608 if (numberOfIterations != 0) {
610 if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
614 if (lambda > LAMBDA_MIN) {
615 lambda /= LAMBDA_DOWN;
619 precessor = successor;
627 if (lambda > LAMBDA_MAX) {
631 evaluate(__begin, __end);
637 for (
int i = 0; i !=
N; ++i) {
639 if (
V(i,i) < PIVOT) {
643 h[i] = 1.0 / sqrt(
V(i,i));
649 for (
int i = 0; i !=
N; ++i) {
650 for (
int j = 0;
j != i; ++
j) {
651 V(
j,i) *= h[i] * h[
j];
656 for (
int i = 0; i !=
N; ++i) {
657 V(i,i) = 1.0 + lambda;
665 ERROR(
"JKatoomb<JGandalf>: " << error.
what() << endl);
670 for (
int i = 0; i !=
N; ++i) {
672 DEBUG(
"u[" << noshowpos << setw(3) << i <<
"] = " << showpos <<
FIXED(20,5) << value[i]);
676 for (
int j = 0;
j !=
N; ++
j) {
677 y +=
V(i,
j) *
Y[
j] * h[i] * h[
j];
682 DEBUG(
' ' <<
FIXED(20,10) << y << noshowpos << endl);
724 for (
T hit = __begin; hit != __end; ++hit) {
730 const double D = hit->getDistance(position);
731 const double Vi = velocity.getInverseVelocity(D, hit->getZ(), position.
getZ());
732 const double toa_s = value.emitter[hit->getEKey()].t1 + D * Vi;
734 const double u = (toa_s - hit->getValue()) / hit->sigma;
735 const double W = sqrt(hit->getWeight());
737 successor += (W*W) * estimator->getRho(u);
739 H_t
H(1.0,
string.getGradient(parameters, hit->getPosition(), hit->getFloor()) * Vi);
741 H *= W * estimator->getPsi(u) / hit->sigma;
745 i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
746 i.tx = value.getIndex(hit->getString(), &H_t::tx);
747 i.ty = value.getIndex(hit->getString(), &H_t::ty);
748 i.tx2 = value.getIndex(hit->getString(), &H_t::tx2);
749 i.ty2 = value.getIndex(hit->getString(), &H_t::ty2);
750 i.vs = value.getIndex(hit->getString(), &H_t::vs);
752 V(i.t1, i.t1) += H.t1 * H.t1;
756 if (hit->getFloor() != 0) {
758 switch (this->option) {
761 V(i.t1, i.vs) += H.t1 * H.vs;
V(i.tx, i.vs) += H.tx * H.vs;
V(i.ty, i.vs) += H.ty * H.vs;
V(i.tx2, i.vs) += H.tx2 * H.vs;
V(i.ty2, i.vs) += H.ty2 * H.vs;
763 V(i.vs, i.t1) += H.vs * H.t1;
764 V(i.vs, i.tx) += H.vs * H.tx;
765 V(i.vs, i.ty) += H.vs * H.ty;
766 V(i.vs, i.tx2) += H.vs * H.tx2;
767 V(i.vs, i.ty2) += H.vs * H.ty2;
769 V(i.vs, i.vs) += H.vs * H.vs;
774 V(i.t1, i.tx2) += H.t1 * H.tx2;
V(i.tx, i.tx2) += H.tx * H.tx2;
V(i.ty, i.tx2) += H.ty * H.tx2;
776 V(i.tx2, i.t1) += H.tx2 * H.t1;
777 V(i.tx2, i.tx) += H.tx2 * H.tx;
778 V(i.tx2, i.ty) += H.tx2 * H.ty;
780 V(i.t1, i.ty2) += H.t1 * H.ty2;
V(i.tx, i.ty2) += H.tx * H.ty2;
V(i.ty, i.ty2) += H.ty * H.ty2;
782 V(i.ty2, i.t1) += H.ty2 * H.t1;
783 V(i.ty2, i.tx) += H.ty2 * H.tx;
784 V(i.ty2, i.ty) += H.ty2 * H.ty;
786 V(i.tx2, i.tx2) += H.tx2 * H.tx2;
V(i.tx2, i.ty2) += H.tx2 * H.ty2;
787 V(i.ty2, i.tx2) += H.ty2 * H.tx2;
V(i.ty2, i.ty2) += H.ty2 * H.ty2;
789 Y[i.tx2] += W * H.tx2;
790 Y[i.ty2] += W * H.ty2;
793 V(i.t1, i.tx) += H.t1 * H.tx;
V(i.t1, i.ty) += H.t1 * H.ty;
794 V(i.tx, i.t1) += H.tx * H.t1;
V(i.ty, i.t1) += H.ty * H.t1;
796 V(i.tx, i.tx) += H.tx * H.tx;
V(i.tx, i.ty) += H.tx * H.ty;
797 V(i.ty, i.tx) += H.ty * H.tx;
V(i.ty, i.ty) += H.ty * H.ty;
878 parameters(parameters),
879 estimator(geometry, V, parameters.option),
880 evaluator(geometry, V, parameters.option),
881 gandalf (geometry, V, parameters.option)
913 result = estimator(__begin, __end);
917 for (
double chi2 = evaluator(
result, __begin, __end), chi2_old =
chi2; ; ) {
923 for (
T hit = __begin; hit != __end; ++hit) {
925 const double x = fabs(hit->getValue() - estimator.getToA(
result, *hit)) / hit->sigma;
935 iter_swap(out, --__end);
937 result = estimator(__begin, __end);
948 iter_swap(out, __end++);
950 result = estimator(__begin, __end);
961 const double chi2 = gandalf (__begin, __end) / gandalf.estimator->getRho(1.0);
962 const double ndf =
getWeight(__begin, __end) - gandalf.value.getN();
964 return { gandalf.value,
chi2, ndf, __begin, __end };
Worker class for fit function of acoustic model.
fit times of emission of emitters and tilt angles of strings with second order correction and stretch...
JFit_t
Enumeration for fit algorithms.
I_t()
Default constructor.
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity, const int option)
Constructor.
General purpose data regression method.
JString & mul(const double factor)
Scale string.
Wrapper class around STL string class.
Auxiliary base class for aritmetic operations of derived class types.
double getQ() const
Get quality.
JKatoomba< JEstimator > estimator
double getWeight(T __begin, T __end)
Get total weight of data points.
int getFloor() const
Get floor number.
static double LAMBDA_UP
multiplication factor control parameter
Interface for maximum likelihood estimator (M-estimator).
JKatoomba< JAbstractMinimiser > evaluator
result_type operator()(const JFunction_t &fit, T __begin, T __end)
Get chi2.
Template specialisation of fit function of acoustic model based on linear approximation.
JKatoomba_t(const JGeometry &geometry, const JSoundVelocity &V, const JFitParameters ¶meters)
Constructor.
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
static const double H
Planck constant [eV s].
*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
JEmitter & mul(const double factor)
Scale emitter.
const JDetector & detector
detector
double operator()(const JModel &model, const JHit< JPDF_t > &hit) const
Fit function.
H_t & mul(const double factor)
Scale H-equation.
Auxiliary class for a type holder.
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.
JEKey getEKey() const
Get emitter hash key of this hit.
Template definition of linear fit.
H_t()
Default constructor.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
double getDistance(const JVector3D &pos) const
Get distance to point.
Auxiliary data structure for return value of fit.
fit times of emission of emitters and tilt angles of strings
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity, const int option)
Constructor.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
Model for fit to acoustics data.
JFitParameters parameters
Template specialisation of fit function of acoustic model based on JGandalf minimiser.
static double EPSILON
maximal distance to minimum
result_type operator()(T __begin, T __end)
Fit.
double getToE(const JModel &model, const JHit< JPDF_t > &hit) const
Get estimated time-of-emission for given hit.
static double LAMBDA_MIN
minimal value control parameter
void evaluate(T __begin, T __end)
Evaluation of fit.
static int debug
debug level
do set_variable OUTPUT_DIRECTORY $WORKDIR T
double operator()(T __begin, T __end)
Fit.
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity, const int option)
Constructor.
double getToA(const JModel &model, const JHit< JPDF_t > &hit) const
Get estimated time-of-arrival for given hit.
static struct JACOUSTICS::@4 compare
Auxiliary data structure to sort transmissions.
result_type operator()(const JModel &model, const JHit< JPDF_t > &hit) const
Fit function.
double getY() const
Get y position.
const JPosition3D & getPosition() const
Get position.
fit times of emission of emitters and tilt angles of strings with second order correction ...
double getLengthSquared() const
Get length squared.
General purpose messaging.
Implementation for depth dependend velocity of sound.
static double LAMBDA_MAX
maximal value control parameter
static int MAXIMUM_ITERATIONS
maximal number of iterations
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity, const int option)
Constructor.
double operator()(const JFunction_t &fit, T __begin, T __end)
Multi-dimensional fit.
JACOUSTICS::JModel::string_type string
JSoundVelocity velocity
sound velocity
std::shared_ptr< JMEstimator > estimator_type
JKatoomba(const JDetector &detector, const JSoundVelocity &velocity, const int option)
Constructor.
JACOUSTICS::JModel::emitter_type emitter
int getString() const
Get string number.
virtual const char * what() const override
Get error message.
const JModel & operator()(T __begin, T __end) const
Get start values of string parameters.
static double LAMBDA_DOWN
multiplication factor control parameter
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++...
Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser...
result_type operator()(const JModel &model, T __begin, T __end)
Fit.
double getToA() const
Get calibrated time of arrival.
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
static double PIVOT
minimal value diagonal element of matrix
int getID() const
Get identifier.
double getX() const
Get x position.
int mestimator
M-estimator.
Data structure for position in three dimensions.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
H_t(const JMODEL::JEmitter &emitter, const JMODEL::JString &string)
Constructor.
Model for fit to acoutsics data.
fit only times of emission of emitters
estimator_type estimator
M-Estimator function.
do echo Generating $dir eval D
double getZ() const
Get z position.
Maximum likelihood estimator (M-estimators).
Template definition of fit function of acoustic model.
#define DEBUG(A)
Message macros.
JKatoomba< JGandalf > gandalf
result_type< T > operator()(T __begin, T __end)
Fit.