Jpp  19.1.0
the software that should make you happy
Public Types | Public Member Functions | Public Attributes | Static Public Attributes | Private Member Functions | Private Attributes | List of all members
JACOUSTICS::JKatoomba< JGandalf > Struct Reference

Template specialisation of fit function of acoustic model based on JGandalf minimiser. More...

#include <JKatoomba_t.hh>

Inheritance diagram for JACOUSTICS::JKatoomba< JGandalf >:
JACOUSTICS::JKatoomba< JMinimiser_t >

Public Types

typedef double result_type
 
typedef std::map< JLocation, std::map< JEmitter, std::vector< JHit > > > data_type
 Type definition internal data structure. More...
 

Public Member Functions

 JKatoomba (const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
 Constructor The option corresponds to the use of fit parameters in the model of the detector geometry. More...
 
template<class T >
result_type operator() (T __begin, T __end)
 Fit. More...
 

Public Attributes

double lambda
 
JModel value
 
int numberOfIterations
 
JMATH::JMatrixNS V
 

Static Public Attributes

static int debug = 0
 debug level More...
 
static int MAXIMUM_ITERATIONS = 1000
 maximal number of iterations More...
 
static double EPSILON = 1.0e-3
 maximal distance to minimum More...
 
static double LAMBDA_MIN = 0.01
 minimal value control parameter More...
 
static double LAMBDA_MAX = 100.0
 maximal value control parameter More...
 
static double LAMBDA_UP = 10.0
 multiplication factor control parameter More...
 
static double LAMBDA_DOWN = 10.0
 multiplication factor control parameter More...
 
static double PIVOT = std::numeric_limits<double>::epsilon()
 minimal value diagonal element of matrix More...
 

Private Member Functions

void evaluate (const data_type &data)
 Evaluation of fit. More...
 

Private Attributes

JMATH::JVectorND Y
 
result_type successor
 
JModel previous
 
std::vector< double > h
 

Detailed Description

Template specialisation of fit function of acoustic model based on JGandalf minimiser.

Definition at line 678 of file JKatoomba_t.hh.

Member Typedef Documentation

◆ result_type

typedef double JACOUSTICS::JKatoomba< JGandalf >::result_type

Definition at line 681 of file JKatoomba_t.hh.

◆ data_type

typedef std::map<JLocation, std::map<JEmitter, std::vector<JHit> > > JACOUSTICS::JKatoomba< JGandalf >::data_type

Type definition internal data structure.

Definition at line 686 of file JKatoomba_t.hh.

Constructor & Destructor Documentation

◆ JKatoomba()

JACOUSTICS::JKatoomba< JGandalf >::JKatoomba ( const JGeometry &  geometry,
const JSoundVelocity &  velocity,
const int  option 
)
inline

Constructor The option corresponds to the use of fit parameters in the model of the detector geometry.


A negative implies that all strings in the detector use common fit parameters.

Parameters
geometrydetector geometry
velocitysound velocity
optionoption

Definition at line 698 of file JKatoomba_t.hh.

700  :
701  JKatoomba<>(geometry, velocity, option)
702  {};
JACOUSTICS::JKatoomba
Template definition of fit function of acoustic model.
Definition: JKatoomba_t.hh:77

Member Function Documentation

◆ operator()()

template<class T >
result_type JACOUSTICS::JKatoomba< JGandalf >::operator() ( T  __begin,
T  __end 
)
inline

Fit.

Parameters
__beginbegin of hits
__endend of hits
Returns
chi2 and gradient

Definition at line 713 of file JKatoomba_t.hh.

714  {
715  using namespace std;
716  using namespace JPP;
717 
718  value.setOption(this->option);
719 
720  const int N = value.getN();
721 
722  V.resize(N);
723  Y.resize(N);
724  h.resize(N);
725 
726  data_type data;
727 
728  for (T hit = __begin; hit != __end; ++hit) {
729  data[hit->getLocation()][hit->getEmitter()].push_back(*hit);
730  }
731 
732  lambda = LAMBDA_MIN;
733 
734  result_type precessor = numeric_limits<double>::max();
735 
736  for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
737 
738  DEBUG("step: " << numberOfIterations << endl);
739 
740  evaluate(data);
741 
742  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
743  DEBUG("chi2: " << FIXED(12,5) << successor << endl);
744 
745  if (successor < precessor) {
746 
747  if (numberOfIterations != 0) {
748 
749  if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
750  return successor;
751  }
752 
753  if (lambda > LAMBDA_MIN) {
754  lambda /= LAMBDA_DOWN;
755  }
756  }
757 
758  precessor = successor;
759  previous = value;
760 
761  } else {
762 
763  value = previous;
764  lambda *= LAMBDA_UP;
765 
766  if (lambda > LAMBDA_MAX) {
767  return precessor; // no improvement found
768  }
769 
770  evaluate(data);
771  }
772 
773  // force definite positiveness
774 
775  for (int i = 0; i != N; ++i) {
776 
777  if (V(i,i) < PIVOT) {
778  V(i,i) = PIVOT;
779  }
780 
781  h[i] = 1.0 / sqrt(V(i,i));
782  }
783 
784  // normalisation
785 
786  for (int i = 0; i != N; ++i) {
787  for (int j = 0; j != i; ++j) {
788  V(j,i) *= h[i] * h[j];
789  V(i,j) = V(j,i);
790  }
791  }
792 
793  for (int i = 0; i != N; ++i) {
794  V(i,i) = 1.0 + lambda;
795  }
796 
797 
798  // solve A x = b
799 
800  for (int col = 0; col != N; ++col) {
801  Y[col] *= h[col];
802  }
803 
804  try {
805  V.solve(Y);
806  }
807  catch (const exception& error) {
808 
809  ERROR("JGandalf: " << error.what() << endl << V << endl);
810 
811  break;
812  }
813 
814  // update value
815 
816  for (int row = 0; row != N; ++row) {
817 
818  DEBUG("u[" << noshowpos << setw(3) << row << "] = " << showpos << FIXED(20,5) << value[row]);
819 
820  value[row] -= h[row] * Y[row];
821 
822  DEBUG(" -> " << FIXED(20,10) << value[row] << noshowpos << endl);
823  }
824  }
825 
826  return precessor;
827  }
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
JTOOLS::j
int j
Definition: JPolint.hh:792
std
Definition: JSTDTypes.hh:14
FIXED
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
JACOUSTICS::JKatoomba< JGandalf >::LAMBDA_DOWN
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JKatoomba_t.hh:835
JACOUSTICS::JKatoomba< JGandalf >::PIVOT
static double PIVOT
minimal value diagonal element of matrix
Definition: JKatoomba_t.hh:836
JACOUSTICS::JKatoomba< JGandalf >::LAMBDA_UP
static double LAMBDA_UP
multiplication factor control parameter
Definition: JKatoomba_t.hh:834
JACOUSTICS::JKatoomba< JGandalf >::MAXIMUM_ITERATIONS
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JKatoomba_t.hh:830
JACOUSTICS::JKatoomba< JGandalf >::LAMBDA_MAX
static double LAMBDA_MAX
maximal value control parameter
Definition: JKatoomba_t.hh:833
JACOUSTICS::JKatoomba< JGandalf >::EPSILON
static double EPSILON
maximal distance to minimum
Definition: JKatoomba_t.hh:831
JACOUSTICS::JKatoomba< JGandalf >::evaluate
void evaluate(const data_type &data)
Evaluation of fit.
Definition: JKatoomba_t.hh:849
JACOUSTICS::JKatoomba< JGandalf >::LAMBDA_MIN
static double LAMBDA_MIN
minimal value control parameter
Definition: JKatoomba_t.hh:832
JACOUSTICS::JModel::setOption
void setOption(const int option)
Set fit option.
Definition: JAcoustics/JModel.hh:553
JACOUSTICS::JModel::getN
size_t getN() const
Get number of fit parameters.
Definition: JAcoustics/JModel.hh:710
JCALIBRATE::data_type
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition: JFitK40.hh:103
JMATH::JMatrixND::resize
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:446
JMATH::JMatrixNS::solve
void solve(JVectorND_t &u)
Get solution of equation A x = b.
Definition: JMatrixNS.hh:308

◆ evaluate()

void JACOUSTICS::JKatoomba< JGandalf >::evaluate ( const data_type &  data)
inlineprivate

Evaluation of fit.

Parameters
datadata

Definition at line 849 of file JKatoomba_t.hh.

850  {
851  using namespace std;
852  using namespace JPP;
853 
854  successor = 0.0;
855 
856  V.reset();
857  Y.reset();
858 
859  for (data_type::const_iterator p = data.begin(); p != data.end(); ++p) {
860 
861  const JGEOMETRY::JString& string = geometry [p->first.getString()];
862  const JMODEL ::JString& parameters = value.string[p->first.getString()];
863  const JPosition3D position = string.getPosition(parameters, p->first.getFloor());
864 
865  for (data_type::mapped_type::const_iterator emitter = p->second.begin(); emitter != p->second.end(); ++emitter) {
866 
867  const double D = emitter->first.getDistance(position);
868  const double Vi = velocity.getInverseVelocity(D, emitter->first.getZ(), position.getZ());
869 
870  const H_t H0(1.0, string.getGradient(parameters, emitter->first.getPosition(), p->first.getFloor()) * Vi);
871 
872  for (data_type::mapped_type::mapped_type::const_iterator hit = emitter->second.begin(); hit != emitter->second.end(); ++hit) {
873 
874  const double toa_s = value.emission[hit->getEKey()].t1 + D * Vi;
875 
876  const double u = (toa_s - hit->getValue()) / hit->getSigma();
877  const double W = sqrt(hit->getWeight());
878 
879  successor += (W*W) * estimator->getRho(u);
880 
881  const H_t H = H0 * (W * estimator->getPsi(u) / hit->getSigma());
882 
883  I_t i;
884 
885  i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
886  i.tx = value.getIndex(hit->getString(), &H_t::tx);
887  i.ty = value.getIndex(hit->getString(), &H_t::ty);
888  i.tx2 = value.getIndex(hit->getString(), &H_t::tx2);
889  i.ty2 = value.getIndex(hit->getString(), &H_t::ty2);
890  i.vs = value.getIndex(hit->getString(), &H_t::vs);
891 
892  V(i.t1, i.t1) += H.t1 * H.t1;
893 
894  Y[i.t1] += W * H.t1;
895 
896  if (hit->getFloor() != 0) {
897 
898  switch (this->option) {
899 
900  case JMODEL::FIT_EMITTERS_AND_STRINGS_2nd_ORDER_AND_STRETCHING_t:
901  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;
902 
903  V(i.vs, i.t1) = V(i.t1, i.vs);
904  V(i.vs, i.tx) = V(i.tx, i.vs);
905  V(i.vs, i.ty) = V(i.ty, i.vs);
906  V(i.vs, i.tx2) = V(i.tx2, i.vs);
907  V(i.vs, i.ty2) = V(i.ty2, i.vs);
908 
909  V(i.vs, i.vs) += H.vs * H.vs;
910 
911  Y[i.vs] += W * H.vs;
912 
913  case JMODEL::FIT_EMITTERS_AND_STRINGS_2nd_ORDER_t:
914  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;
915 
916  V(i.tx2, i.t1) = V(i.t1, i.tx2);
917  V(i.tx2, i.tx) = V(i.tx, i.tx2);
918  V(i.tx2, i.ty) = V(i.ty, i.tx2);
919 
920  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;
921 
922  V(i.ty2, i.t1) = V(i.t1, i.ty2);
923  V(i.ty2, i.tx) = V(i.tx, i.ty2);
924  V(i.ty2, i.ty) = V(i.ty, i.ty2);
925 
926  V(i.tx2, i.tx2) += H.tx2 * H.tx2; V(i.tx2, i.ty2) += H.tx2 * H.ty2;
927  V(i.ty2, i.tx2) = V(i.tx2, i.ty2); V(i.ty2, i.ty2) += H.ty2 * H.ty2;
928 
929  Y[i.tx2] += W * H.tx2;
930  Y[i.ty2] += W * H.ty2;
931 
932  case JMODEL::FIT_EMITTERS_AND_STRINGS_1st_ORDER_t:
933  V(i.t1, i.tx) += H.t1 * H.tx; V(i.t1, i.ty) += H.t1 * H.ty;
934  V(i.tx, i.t1) = V(i.t1, i.tx); V(i.ty, i.t1) = V(i.t1, i.ty);
935 
936  V(i.tx, i.tx) += H.tx * H.tx; V(i.tx, i.ty) += H.tx * H.ty;
937  V(i.ty, i.tx) = V(i.tx, i.ty); V(i.ty, i.ty) += H.ty * H.ty;
938 
939  Y[i.tx] += W * H.tx;
940  Y[i.ty] += W * H.ty;
941  break;
942 
943  default:
944  break;
945  }
946  }
947  }
948  }
949  }
950  }
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
JGEOMETRY3D::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
JGEOMETRY3D::JVector3D::getDistance
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
JGEOMETRY3D::JVector3D::getZ
double getZ() const
Get z position.
Definition: JVector3D.hh:115
JACOUSTICS::JMODEL::FIT_EMITTERS_AND_STRINGS_1st_ORDER_t
@ FIT_EMITTERS_AND_STRINGS_1st_ORDER_t
fit times of emission of emitters and tilt angles of strings
Definition: JAcoustics/JModel.hh:50
JACOUSTICS::JMODEL::FIT_EMITTERS_AND_STRINGS_2nd_ORDER_t
@ FIT_EMITTERS_AND_STRINGS_2nd_ORDER_t
fit times of emission of emitters and tilt angles of strings with second order correction
Definition: JAcoustics/JModel.hh:51
JACOUSTICS::JMODEL::FIT_EMITTERS_AND_STRINGS_2nd_ORDER_AND_STRETCHING_t
@ FIT_EMITTERS_AND_STRINGS_2nd_ORDER_AND_STRETCHING_t
fit times of emission of emitters and tilt angles of strings with second order correction and stretch...
Definition: JAcoustics/JModel.hh:52
JPHYSICS::H
static const double H
Planck constant [eV s].
Definition: JPhysics/JConstants.hh:29
JSIRENE::vs
std::vector< double > vs
Definition: JSireneToolkit.hh:139
JTOOLS::u
double u[N+1]
Definition: JPolint.hh:865
JACOUSTICS::JModel::emission
JACOUSTICS::JModel::emission_type emission
JACOUSTICS::JModel::getIndex
size_t getIndex(int id, double JString::*p) const
Get index of fit parameter for given string.
Definition: JAcoustics/JModel.hh:723
JACOUSTICS::JModel::string
JACOUSTICS::JModel::string_type string
JMATH::JMatrixND::reset
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:459
JMATH::JVectorND::reset
void reset()
Reset.
Definition: JVectorND.hh:45

Member Data Documentation

◆ debug

int JACOUSTICS::JKatoomba< JGandalf >::debug = 0
static

debug level

debug level.

Definition at line 829 of file JKatoomba_t.hh.

◆ MAXIMUM_ITERATIONS

int JACOUSTICS::JKatoomba< JGandalf >::MAXIMUM_ITERATIONS = 1000
static

maximal number of iterations

maximal number of iterations.

Definition at line 830 of file JKatoomba_t.hh.

◆ EPSILON

double JACOUSTICS::JKatoomba< JGandalf >::EPSILON = 1.0e-3
static

maximal distance to minimum

maximal distance to minimum.

Definition at line 831 of file JKatoomba_t.hh.

◆ LAMBDA_MIN

double JACOUSTICS::JKatoomba< JGandalf >::LAMBDA_MIN = 0.01
static

minimal value control parameter

Definition at line 832 of file JKatoomba_t.hh.

◆ LAMBDA_MAX

double JACOUSTICS::JKatoomba< JGandalf >::LAMBDA_MAX = 100.0
static

maximal value control parameter

Definition at line 833 of file JKatoomba_t.hh.

◆ LAMBDA_UP

double JACOUSTICS::JKatoomba< JGandalf >::LAMBDA_UP = 10.0
static

multiplication factor control parameter

Definition at line 834 of file JKatoomba_t.hh.

◆ LAMBDA_DOWN

double JACOUSTICS::JKatoomba< JGandalf >::LAMBDA_DOWN = 10.0
static

multiplication factor control parameter

Definition at line 835 of file JKatoomba_t.hh.

◆ PIVOT

double JACOUSTICS::JKatoomba< JGandalf >::PIVOT = std::numeric_limits<double>::epsilon()
static

minimal value diagonal element of matrix

Definition at line 836 of file JKatoomba_t.hh.

◆ lambda

double JACOUSTICS::JKatoomba< JGandalf >::lambda

Definition at line 838 of file JKatoomba_t.hh.

◆ value

JModel JACOUSTICS::JKatoomba< JGandalf >::value

Definition at line 839 of file JKatoomba_t.hh.

◆ numberOfIterations

int JACOUSTICS::JKatoomba< JGandalf >::numberOfIterations

Definition at line 840 of file JKatoomba_t.hh.

◆ V

JMATH::JMatrixNS JACOUSTICS::JKatoomba< JGandalf >::V

Definition at line 841 of file JKatoomba_t.hh.

◆ Y

JMATH::JVectorND JACOUSTICS::JKatoomba< JGandalf >::Y
private

Definition at line 953 of file JKatoomba_t.hh.

◆ successor

result_type JACOUSTICS::JKatoomba< JGandalf >::successor
private

Definition at line 954 of file JKatoomba_t.hh.

◆ previous

JModel JACOUSTICS::JKatoomba< JGandalf >::previous
private

Definition at line 955 of file JKatoomba_t.hh.

◆ h

std::vector<double> JACOUSTICS::JKatoomba< JGandalf >::h
private

Definition at line 956 of file JKatoomba_t.hh.

The documentation for this struct was generated from the following files:
Generated by doxygen 1.9.1