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) {
104 template<
template<
class T>
class JMinimiser_t =
JType>
147 const double Vi = velocity.getInverseVelocity(D, hit.
getZ(), position.
getZ());
167 const double Vi = velocity.getInverseVelocity(D, hit.
getZ(), position.
getZ());
204 JMODEL::JEmission(emission),
215 H_t&
mul(
const double factor)
289 const double toa_s = this->getToA(model, hit);
292 return estimator->getRho(u) * hit.
getWeight();
306 this->value.setOption(this->option);
325 return (*
this)(__begin, __end);
337 static constexpr
double TOLERANCE = 1.0e-8;
366 using namespace JGEOMETRY;
371 unique_lock<mutex> lock(mtx);
373 value =
JModel(__begin, __end);
383 const size_t N = value.getN();
388 for (
T hit = __begin; hit != __end; ++hit) {
390 const JString&
string = geometry[hit->getString()];
392 const double Vi = velocity.getInverseVelocity(hit->getDistance(position), hit->getZ(), position.
getZ());
394 const double h1 =
string.getHeight(hit->getFloor());
395 const JPosition3D p1 =
string.getPosition() - hit->getPosition();
397 const double y = hit->getValue() - Vi*ds;
398 const double W = sqrt(hit->getWeight());
401 H.tx = W * Vi * p1.
getX() * h1 / ds;
402 H.ty = W * Vi * p1.
getY() * h1 / ds;
404 i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
405 i.tx = value.getIndex(hit->getString(), &H_t::tx);
406 i.ty = value.getIndex(hit->getString(), &H_t::ty);
408 V(i.t1, i.t1) +=
H.t1 *
H.t1;
410 Y[i.t1] += W *
H.t1 * y;
412 if (hit->getFloor() != 0) {
416 V(i.t1, i.tx) +=
H.t1 *
H.tx;
V(i.t1, i.ty) +=
H.t1 *
H.ty;
417 V(i.tx, i.t1) +=
H.tx *
H.t1;
V(i.ty, i.t1) +=
H.ty *
H.t1;
419 V(i.tx, i.tx) +=
H.tx *
H.tx;
V(i.tx, i.ty) +=
H.tx *
H.ty;
420 V(i.ty, i.tx) +=
H.ty *
H.tx;
V(i.ty, i.ty) +=
H.ty *
H.ty;
422 Y[i.tx] += W *
H.tx * y;
423 Y[i.ty] += W *
H.ty * y;
430 TDecompSVD svd(V, TOLERANCE);
434 V = svd.Invert(status);
436 for (
size_t row = 0; row !=
N; ++row) {
437 for (
size_t col = 0; col !=
N; ++col) {
438 value[row] +=
V(row,col) * Y[col];
448 static std::mutex mtx;
491 const double toa_s = this->getToA(model, hit);
494 return estimator->getRho(u);
514 switch (this->option) {
589 value.setOption(this->option);
591 const int N = value.getN();
599 result_type precessor = numeric_limits<double>::max();
601 for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
603 DEBUG(
"step: " << numberOfIterations << endl);
605 evaluate(__begin, __end);
607 DEBUG(
"lambda: " <<
FIXED(12,5) << lambda << endl);
608 DEBUG(
"chi2: " <<
FIXED(12,5) << successor << endl);
610 if (successor < precessor) {
612 if (numberOfIterations != 0) {
614 if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
618 if (lambda > LAMBDA_MIN) {
619 lambda /= LAMBDA_DOWN;
623 precessor = successor;
631 if (lambda > LAMBDA_MAX) {
635 evaluate(__begin, __end);
641 for (
int i = 0; i !=
N; ++i) {
643 if (
V(i,i) < PIVOT) {
647 h[i] = 1.0 / sqrt(
V(i,i));
653 for (
int i = 0; i !=
N; ++i) {
654 for (
int j = 0;
j != i; ++
j) {
655 V(
j,i) *= h[i] * h[
j];
660 for (
int i = 0; i !=
N; ++i) {
661 V(i,i) = 1.0 + lambda;
669 ERROR(
"JKatoomb<JGandalf>: " << error.
what() << endl);
674 for (
int i = 0; i !=
N; ++i) {
676 DEBUG(
"u[" << noshowpos << setw(3) << i <<
"] = " << showpos <<
FIXED(20,5) << value[i]);
680 for (
int j = 0;
j !=
N; ++
j) {
681 y +=
V(i,
j) *
Y[
j] * h[i] * h[
j];
686 DEBUG(
' ' <<
FIXED(20,10) << y << noshowpos << endl);
728 for (
T hit = __begin; hit != __end; ++hit) {
734 const double D = hit->getDistance(position);
735 const double Vi = velocity.getInverseVelocity(D, hit->getZ(), position.
getZ());
736 const double toa_s = value.emission[hit->getEKey()].t1 + D * Vi;
738 const double u = (toa_s - hit->getValue()) / hit->getSigma();
739 const double W = sqrt(hit->getWeight());
741 successor += (W*W) * estimator->getRho(u);
743 H_t
H(1.0,
string.getGradient(parameters, hit->getPosition(), hit->getFloor()) * Vi);
745 H *= W * estimator->getPsi(u) / hit->getSigma();
749 i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
750 i.tx = value.getIndex(hit->getString(), &H_t::tx);
751 i.ty = value.getIndex(hit->getString(), &H_t::ty);
752 i.tx2 = value.getIndex(hit->getString(), &H_t::tx2);
753 i.ty2 = value.getIndex(hit->getString(), &H_t::ty2);
754 i.vs = value.getIndex(hit->getString(), &
H_t::vs);
756 V(i.t1, i.t1) += H.t1 * H.t1;
760 if (hit->getFloor() != 0) {
762 switch (this->option) {
765 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;
767 V(i.vs, i.t1) += H.vs * H.t1;
768 V(i.vs, i.tx) += H.vs * H.tx;
769 V(i.vs, i.ty) += H.vs * H.ty;
770 V(i.vs, i.tx2) += H.vs * H.tx2;
771 V(i.vs, i.ty2) += H.vs * H.ty2;
773 V(i.vs, i.vs) += H.vs * H.vs;
778 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;
780 V(i.tx2, i.t1) += H.tx2 * H.t1;
781 V(i.tx2, i.tx) += H.tx2 * H.tx;
782 V(i.tx2, i.ty) += H.tx2 * H.ty;
784 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;
786 V(i.ty2, i.t1) += H.ty2 * H.t1;
787 V(i.ty2, i.tx) += H.ty2 * H.tx;
788 V(i.ty2, i.ty) += H.ty2 * H.ty;
790 V(i.tx2, i.tx2) += H.tx2 * H.tx2;
V(i.tx2, i.ty2) += H.tx2 * H.ty2;
791 V(i.ty2, i.tx2) += H.ty2 * H.tx2;
V(i.ty2, i.ty2) += H.ty2 * H.ty2;
793 Y[i.tx2] += W * H.tx2;
794 Y[i.ty2] += W * H.ty2;
797 V(i.t1, i.tx) += H.t1 * H.tx;
V(i.t1, i.ty) += H.t1 * H.ty;
798 V(i.tx, i.t1) += H.tx * H.t1;
V(i.ty, i.t1) += H.ty * H.t1;
800 V(i.tx, i.tx) += H.tx * H.tx;
V(i.tx, i.ty) += H.tx * H.ty;
801 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
result_type operator()(const JFunction_t &fit, T __begin, T __end)
Get chi2.
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.
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
static struct JACOUSTICS::@4 compare
Auxiliary data structure to sort transmissions.
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 int MAXIMUM_ITERATIONS
maximal number of iterations
JEmission & mul(const double factor)
Scale emission.
double operator()(const JFunction_t &fit, T __begin, T __end)
Multi-dimensional fit.
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
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
const JModel & operator()(T __begin, T __end) const
Get start values of string parameters.
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.
Data structure for position in three dimensions.
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.