1 #ifndef __JACOUSTICS__JKATOOMBA__
2 #define __JACOUSTICS__JKATOOMBA__
9 #include "TMatrixTSym.h"
10 #include "TDecompSVD.h"
38 namespace JACOUSTICS {}
39 namespace JPP {
using namespace JACOUSTICS; }
41 namespace JACOUSTICS {
64 for (
T i = __begin; i != __end; ++i) {
75 template<
template<
class T>
class JMinimiser_t =
JType>
91 static const struct compare {
109 return first.
getQ() > second.
getQ();
149 const double Vi = velocity.getInverseVelocity(D, hit.
getZ(), position.
getZ());
169 const double Vi = velocity.getInverseVelocity(D, hit.
getZ(), position.
getZ());
206 JMODEL::JEmission(emission),
217 H_t&
mul(
const double factor)
291 const double toa_s = this->getToA(model, hit);
294 return estimator->getRho(u) * hit.
getWeight();
308 this->value.setOption(this->option);
327 return (*
this)(__begin, __end);
368 static constexpr
double TOLERANCE = 1.0e-8;
395 switch (MATRIX_INVERSION) {
421 template<
class T,
class JMatrixNS_t>
426 using namespace JGEOMETRY;
428 value = move(
JModel(__begin, __end));
438 const size_t N = value.getN();
443 for (
T hit = __begin; hit != __end; ++hit) {
445 const JString&
string = geometry[hit->getString()];
447 const double Vi = velocity.getInverseVelocity(hit->getDistance(position), hit->getZ(), position.
getZ());
449 const double h1 =
string.getHeight(hit->getFloor());
450 const JPosition3D p1 =
string.getPosition() - hit->getPosition();
452 const double y = hit->getValue() - Vi*ds;
453 const double W = sqrt(hit->getWeight());
456 H.tx = W * Vi * p1.
getX() * h1 / ds;
457 H.ty = W * Vi * p1.
getY() * h1 / ds;
459 i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
460 i.tx = value.getIndex(hit->getString(), &H_t::tx);
461 i.ty = value.getIndex(hit->getString(), &H_t::ty);
463 V(i.t1, i.t1) +=
H.t1 *
H.t1;
465 Y[i.t1] += W *
H.t1 *
y;
467 if (hit->getFloor() != 0) {
471 V(i.t1, i.tx) +=
H.t1 *
H.tx;
V(i.t1, i.ty) +=
H.t1 *
H.ty;
472 V(i.tx, i.t1) +=
H.tx *
H.t1;
V(i.ty, i.t1) +=
H.ty *
H.t1;
474 V(i.tx, i.tx) +=
H.tx *
H.tx;
V(i.tx, i.ty) +=
H.tx *
H.ty;
475 V(i.ty, i.tx) +=
H.ty *
H.tx;
V(i.ty, i.ty) +=
H.ty *
H.ty;
477 Y[i.tx] += W *
H.tx *
y;
478 Y[i.ty] += W *
H.ty *
y;
487 for (
size_t row = 0; row !=
N; ++row) {
488 for (
size_t col = 0; col !=
N; ++col) {
489 value[row] +=
V(row,col) * Y[col];
507 unique_lock<mutex> lock(mtx);
509 TDecompSVD svd(V, TOLERANCE);
513 V = svd.Invert(status);
530 static std::mutex mtx;
574 const double toa_s = this->getToA(model, hit);
577 return estimator->getRho(u);
597 switch (this->option) {
672 value.setOption(this->option);
674 const int N = value.getN();
682 result_type precessor = numeric_limits<double>::max();
684 for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
686 DEBUG(
"step: " << numberOfIterations << endl);
688 evaluate(__begin, __end);
690 DEBUG(
"lambda: " <<
FIXED(12,5) << lambda << endl);
691 DEBUG(
"chi2: " <<
FIXED(12,5) << successor << endl);
693 if (successor < precessor) {
695 if (numberOfIterations != 0) {
697 if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
701 if (lambda > LAMBDA_MIN) {
702 lambda /= LAMBDA_DOWN;
706 precessor = successor;
714 if (lambda > LAMBDA_MAX) {
718 evaluate(__begin, __end);
724 for (
int i = 0; i !=
N; ++i) {
726 if (
V(i,i) < PIVOT) {
730 h[i] = 1.0 / sqrt(
V(i,i));
736 for (
int i = 0; i !=
N; ++i) {
737 for (
int j = 0;
j != i; ++
j) {
738 V(
j,i) *= h[i] * h[
j];
743 for (
int i = 0; i !=
N; ++i) {
744 V(i,i) = 1.0 + lambda;
752 ERROR(
"JKatoomb<JGandalf>: " << error.
what() << endl);
757 for (
int i = 0; i !=
N; ++i) {
759 DEBUG(
"u[" << noshowpos << setw(3) << i <<
"] = " << showpos <<
FIXED(20,5) << value[i]);
763 for (
int j = 0;
j !=
N; ++
j) {
764 y +=
V(i,
j) *
Y[
j] * h[i] * h[
j];
769 DEBUG(
' ' <<
FIXED(20,10) << y << noshowpos << endl);
811 for (
T hit = __begin; hit != __end; ++hit) {
817 const double D = hit->getDistance(position);
818 const double Vi = velocity.getInverseVelocity(D, hit->getZ(), position.
getZ());
819 const double toa_s = value.emission[hit->getEKey()].t1 + D * Vi;
821 const double u = (toa_s - hit->getValue()) / hit->getSigma();
822 const double W = sqrt(hit->getWeight());
824 successor += (W*W) * estimator->getRho(u);
826 H_t
H(1.0,
string.getGradient(parameters, hit->getPosition(), hit->getFloor()) * Vi);
828 H *= W * estimator->getPsi(u) / hit->getSigma();
832 i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
833 i.tx = value.getIndex(hit->getString(), &H_t::tx);
834 i.ty = value.getIndex(hit->getString(), &H_t::ty);
835 i.tx2 = value.getIndex(hit->getString(), &H_t::tx2);
836 i.ty2 = value.getIndex(hit->getString(), &H_t::ty2);
837 i.vs = value.getIndex(hit->getString(), &
H_t::vs);
839 V(i.t1, i.t1) += H.t1 * H.t1;
843 if (hit->getFloor() != 0) {
845 switch (this->option) {
848 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;
850 V(i.vs, i.t1) += H.vs * H.t1;
851 V(i.vs, i.tx) += H.vs * H.tx;
852 V(i.vs, i.ty) += H.vs * H.ty;
853 V(i.vs, i.tx2) += H.vs * H.tx2;
854 V(i.vs, i.ty2) += H.vs * H.ty2;
856 V(i.vs, i.vs) += H.vs * H.vs;
861 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;
863 V(i.tx2, i.t1) += H.tx2 * H.t1;
864 V(i.tx2, i.tx) += H.tx2 * H.tx;
865 V(i.tx2, i.ty) += H.tx2 * H.ty;
867 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;
869 V(i.ty2, i.t1) += H.ty2 * H.t1;
870 V(i.ty2, i.tx) += H.ty2 * H.tx;
871 V(i.ty2, i.ty) += H.ty2 * H.ty;
873 V(i.tx2, i.tx2) += H.tx2 * H.tx2;
V(i.tx2, i.ty2) += H.tx2 * H.ty2;
874 V(i.ty2, i.tx2) += H.ty2 * H.tx2;
V(i.ty2, i.ty2) += H.ty2 * H.ty2;
876 Y[i.tx2] += W * H.tx2;
877 Y[i.ty2] += W * H.ty2;
880 V(i.t1, i.tx) += H.t1 * H.tx;
V(i.t1, i.ty) += H.t1 * H.ty;
881 V(i.tx, i.t1) += H.tx * H.t1;
V(i.ty, i.t1) += H.ty * H.t1;
883 V(i.tx, i.tx) += H.tx * H.tx;
V(i.tx, i.ty) += H.tx * H.ty;
884 V(i.ty, i.tx) += H.ty * H.tx;
V(i.ty, i.ty) += H.ty * H.ty;
fit times of emission of emitters and tilt angles of strings with second order correction and stretch...
I_t()
Default 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.
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).
double getWeight() const
Get weight.
const JGeometry & geometry
detector geometry
bool operator()(const JTransmission &first, const JTransmission &second) const
Sort transmissions in following order.
result_type operator()(const JFunction_t &fit, T __begin, T __end)
Get chi2.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
static void invert(TMatrixD &V)
Invert matrix using SVD algorithm.
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
H_t & mul(const double factor)
Scale H-equation.
JEKey getEKey() const
Get emitter hash key of this hit.
double getValue() const
Get expectation value of time-of-arrival.
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.
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.
fit times of emission of emitters and tilt angles of strings
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
Place holder for custom implementation.
Model for fit to acoustics data.
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor.
JMatrix_t
Algorithm for matrix inversion.
static double EPSILON
maximal distance to minimum
double operator()(const JModel &model, const JHit &hit) const
Fit function.
result_type operator()(T __begin, T __end)
Fit.
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor.
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.
result_type operator()(T __begin, T __end)
Fit.
JACOUSTICS::JModel::emission_type emission
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor.
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 void invert(JMatrixNS &V)
Invert matrix using LDU algorithm.
static int MAXIMUM_ITERATIONS
maximal number of iterations
static JMatrix_t MATRIX_INVERSION
JEmission & mul(const double factor)
Scale emission.
double operator()(const JFunction_t &fit, T __begin, T __end)
Multi-dimensional fit.
void invert()
Invert matrix according LDU decomposition.
Auxiliary data structure for compatibility of symmetric matrix.
double getToE(const JModel &model, const JHit &hit) const
Get estimated time-of-emission for given hit.
JACOUSTICS::JModel::string_type string
JSoundVelocity velocity
sound velocity
JModel & operator()(T __begin, T __end) const
Get start values of string parameters.
std::shared_ptr< JMEstimator > estimator_type
H_t(const JMODEL::JEmission &emission, const JMODEL::JString &string)
Constructor.
int getString() const
Get string number.
virtual const char * what() const override
Get error message.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
static double LAMBDA_DOWN
multiplication factor control parameter
result_type operator()(const JModel &model, const JHit &hit) const
Fit function.
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++...
result_type operator()(const JModel &model, T __begin, T __end)
Fit.
double getToA() const
Get calibrated time of arrival.
double getToA(const JModel &model, const JHit &hit) const
Get estimated time-of-arrival for given hit.
static double PIVOT
minimal value diagonal element of matrix
int getID() const
Get identifier.
double getX() const
Get x position.
JModel & evaluate(T __begin, T __end, const JType< JMatrixNS_t > &type) const
Get start values of string parameters.
Data structure for position in three dimensions.
Exception for accessing a value in a collection that is outside of its range.
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
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.
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor.
Maximum likelihood estimator (M-estimators).
Template definition of fit function of acoustic model.
#define DEBUG(A)
Message macros.
double getSigma() const
Get resolution of time-of-arrival.