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 {
63 for (
T i = __begin; i != __end; ++i) {
115 template<
template<
class T>
class JMinimiser_t =
JType>
152 template<
class JPDF_t>
160 const double Vi = velocity.getInverseVelocity(D, hit.
getZ(), position.
getZ());
173 template<
class JPDF_t>
181 const double Vi = velocity.getInverseVelocity(D, hit.
getZ(), position.
getZ());
183 return hit.getValue() - D * Vi;
229 H_t&
mul(
const double factor)
301 template<
class JPDF_t>
304 const double toa_s = this->getToA(model, hit);
305 const double u = (toa_s - hit.getValue()) / hit.sigma;
307 return estimator->getRho(u);
336 static constexpr
double TOLERANCE = 1.0e-8;
365 using namespace JGEOMETRY;
370 unique_lock<mutex> lock(mtx);
372 value =
JModel(__begin, __end);
382 const size_t N = value.getN();
387 for (
T hit = __begin; hit != __end; ++hit) {
391 const double Vi = velocity.getInverseVelocity(hit->getDistance(position), hit->getZ(), position.
getZ());
393 const double h1 =
string.getHeight(hit->getFloor());
394 const JPosition3D p1 =
string.getPosition() - hit->getPosition();
396 const double y = hit->getValue() - Vi*ds;
397 const double W = sqrt(hit->getWeight());
400 H.tx = W * Vi * p1.
getX() * h1 / ds;
401 H.ty = W * Vi * p1.
getY() * h1 / ds;
403 i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
404 i.tx = value.getIndex(hit->getString(), &H_t::tx);
405 i.ty = value.getIndex(hit->getString(), &H_t::ty);
407 V(i.t1, i.t1) +=
H.t1 *
H.t1;
409 Y[i.t1] += W *
H.t1 * y;
411 if (hit->getFloor() != 0) {
415 V(i.t1, i.tx) +=
H.t1 *
H.tx;
V(i.t1, i.ty) +=
H.t1 *
H.ty;
416 V(i.tx, i.t1) +=
H.tx *
H.t1;
V(i.ty, i.t1) +=
H.ty *
H.t1;
418 V(i.tx, i.tx) +=
H.tx *
H.tx;
V(i.tx, i.ty) +=
H.tx *
H.ty;
419 V(i.ty, i.tx) +=
H.ty *
H.tx;
V(i.ty, i.ty) +=
H.ty *
H.ty;
421 Y[i.tx] += W *
H.tx * y;
422 Y[i.ty] += W *
H.ty * y;
429 TDecompSVD svd(V, TOLERANCE);
433 V = svd.Invert(status);
435 for (
size_t row = 0; row !=
N; ++row) {
436 for (
size_t col = 0; col !=
N; ++col) {
437 value[row] +=
V(row,col) * Y[col];
447 static std::mutex mtx;
488 template<
class JPDF_t>
491 const double toa_s = this->getToA(model, hit);
492 const double u = (toa_s - hit.getValue()) / hit.sigma;
494 return estimator->getRho(u);
514 switch (this->option) {
586 value.setOption(this->option);
588 const int N = value.getN();
596 result_type precessor = numeric_limits<double>::max();
598 for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
600 DEBUG(
"step: " << numberOfIterations << endl);
602 evaluate(__begin, __end);
604 DEBUG(
"lambda: " <<
FIXED(12,5) << lambda << endl);
605 DEBUG(
"chi2: " <<
FIXED(12,5) << successor << endl);
607 if (successor < precessor) {
609 if (numberOfIterations != 0) {
611 if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
615 if (lambda > LAMBDA_MIN) {
616 lambda /= LAMBDA_DOWN;
620 precessor = successor;
628 if (lambda > LAMBDA_MAX) {
632 evaluate(__begin, __end);
638 for (
int i = 0; i !=
N; ++i) {
640 if (
V(i,i) < PIVOT) {
644 h[i] = 1.0 / sqrt(
V(i,i));
650 for (
int i = 0; i !=
N; ++i) {
651 for (
int j = 0;
j != i; ++
j) {
652 V(
j,i) *= h[i] * h[
j];
657 for (
int i = 0; i !=
N; ++i) {
658 V(i,i) = 1.0 + lambda;
666 ERROR(
"JKatoomb<JGandalf>: " << error.
what() << endl);
671 for (
int i = 0; i !=
N; ++i) {
673 DEBUG(
"u[" << noshowpos << setw(3) << i <<
"] = " << showpos <<
FIXED(20,5) << value[i]);
677 for (
int j = 0;
j !=
N; ++
j) {
678 y +=
V(i,
j) *
Y[
j] * h[i] * h[
j];
683 DEBUG(
' ' <<
FIXED(20,10) << y << noshowpos << endl);
725 for (
T hit = __begin; hit != __end; ++hit) {
731 const double D = hit->getDistance(position);
732 const double Vi = velocity.getInverseVelocity(D, hit->getZ(), position.
getZ());
733 const double toa_s = value.emitter[hit->getEKey()].t1 + D * Vi;
735 const double u = (toa_s - hit->getValue()) / hit->sigma;
736 const double W = sqrt(hit->getWeight());
738 successor += (W*W) * estimator->getRho(u);
740 H_t
H(1.0,
string.getGradient(parameters, hit->getPosition(), hit->getFloor()) * Vi);
742 H *= W * estimator->getPsi(u) / hit->sigma;
746 i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
747 i.tx = value.getIndex(hit->getString(), &H_t::tx);
748 i.ty = value.getIndex(hit->getString(), &H_t::ty);
749 i.tx2 = value.getIndex(hit->getString(), &H_t::tx2);
750 i.ty2 = value.getIndex(hit->getString(), &H_t::ty2);
751 i.vs = value.getIndex(hit->getString(), &
H_t::vs);
753 V(i.t1, i.t1) += H.t1 * H.t1;
757 if (hit->getFloor() != 0) {
759 switch (this->option) {
762 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;
764 V(i.vs, i.t1) += H.vs * H.t1;
765 V(i.vs, i.tx) += H.vs * H.tx;
766 V(i.vs, i.ty) += H.vs * H.ty;
767 V(i.vs, i.tx2) += H.vs * H.tx2;
768 V(i.vs, i.ty2) += H.vs * H.ty2;
770 V(i.vs, i.vs) += H.vs * H.vs;
775 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;
777 V(i.tx2, i.t1) += H.tx2 * H.t1;
778 V(i.tx2, i.tx) += H.tx2 * H.tx;
779 V(i.tx2, i.ty) += H.tx2 * H.ty;
781 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;
783 V(i.ty2, i.t1) += H.ty2 * H.t1;
784 V(i.ty2, i.tx) += H.ty2 * H.tx;
785 V(i.ty2, i.ty) += H.ty2 * H.ty;
787 V(i.tx2, i.tx2) += H.tx2 * H.tx2;
V(i.tx2, i.ty2) += H.tx2 * H.ty2;
788 V(i.ty2, i.tx2) += H.ty2 * H.tx2;
V(i.ty2, i.ty2) += H.ty2 * H.ty2;
790 Y[i.tx2] += W * H.tx2;
791 Y[i.ty2] += W * H.ty2;
794 V(i.t1, i.tx) += H.t1 * H.tx;
V(i.t1, i.ty) += H.t1 * H.ty;
795 V(i.tx, i.t1) += H.tx * H.t1;
V(i.ty, i.t1) += H.ty * H.t1;
797 V(i.tx, i.tx) += H.tx * H.tx;
V(i.tx, i.ty) += H.tx * H.ty;
798 V(i.ty, i.tx) += H.ty * H.tx;
V(i.ty, i.ty) += H.ty * H.ty;
879 parameters(parameters),
880 estimator(geometry, V, parameters.option),
881 evaluator(geometry, V, parameters.option),
882 gandalf (geometry, V, parameters.option)
921 for (
T i = begin; i != end; ++i) {
940 result = estimator(__begin, __end);
944 for (
double chi2 = evaluator(
result, __begin, __end), chi2_old =
chi2; ; ) {
950 for (
T hit = __begin; hit != __end; ++hit) {
952 const double x = fabs(hit->getValue() - estimator.getToA(
result, *hit)) / hit->sigma;
962 iter_swap(out, --__end);
964 result = estimator(__begin, __end);
975 iter_swap(out, __end++);
977 result = estimator(__begin, __end);
988 const double chi2 = gandalf (__begin, __end) / gandalf.estimator->getRho(1.0);
989 const double ndf =
getWeight(__begin, __end) - gandalf.value.getN();
991 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.
int size() const
Get size of data.
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).
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
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.
JTimeRange getTimeRange() const
Get time range of data.
JACOUSTICS::JModel::string_type string
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
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.