Jpp  18.2.1-ARCA-DF-PATCH
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Types | Public Member Functions | Public Attributes | Static Public Attributes | Private Member Functions | Private Attributes | List of all members
JACOUSTICS::JKatoomba< JGandalf > Struct Template 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<>

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. 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<>
struct JACOUSTICS::JKatoomba< JGandalf >

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

Definition at line 733 of file JKatoomba_t.hh.

Member Typedef Documentation

Definition at line 736 of file JKatoomba_t.hh.

Type definition internal data structure.

Definition at line 741 of file JKatoomba_t.hh.

Constructor & Destructor Documentation

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

Constructor.

Parameters
geometrydetector geometry
velocitysound velocity
optionoption

Definition at line 751 of file JKatoomba_t.hh.

753  :
754  JKatoomba<>(geometry, velocity, option)
755  {};

Member Function Documentation

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 766 of file JKatoomba_t.hh.

767  {
768  using namespace std;
769  using namespace JPP;
770 
771  value.setOption(this->option);
772 
773  const int N = value.getN();
774 
775  V.resize(N);
776  Y.resize(N);
777  h.resize(N);
778 
779  data_type data;
780 
781  for (T hit = __begin; hit != __end; ++hit) {
782  data[hit->getLocation()][hit->getEmitter()].push_back(*hit);
783  }
784 
785  lambda = LAMBDA_MIN;
786 
787  result_type precessor = numeric_limits<double>::max();
788 
790 
791  DEBUG("step: " << numberOfIterations << endl);
792 
793  evaluate(data);
794 
795  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
796  DEBUG("chi2: " << FIXED(12,5) << successor << endl);
797 
798  if (successor < precessor) {
799 
800  if (numberOfIterations != 0) {
801 
802  if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
803  return successor;
804  }
805 
806  if (lambda > LAMBDA_MIN) {
807  lambda /= LAMBDA_DOWN;
808  }
809  }
810 
811  precessor = successor;
812  previous = value;
813 
814  } else {
815 
816  value = previous;
817  lambda *= LAMBDA_UP;
818 
819  if (lambda > LAMBDA_MAX) {
820  return precessor; // no improvement found
821  }
822 
823  evaluate(data);
824  }
825 
826  // force definite positiveness
827 
828  for (int i = 0; i != N; ++i) {
829 
830  if (V(i,i) < PIVOT) {
831  V(i,i) = PIVOT;
832  }
833 
834  h[i] = 1.0 / sqrt(V(i,i));
835  }
836 
837  // normalisation
838 
839  for (int i = 0; i != N; ++i) {
840  for (int j = 0; j != i; ++j) {
841  V(j,i) *= h[i] * h[j];
842  V(i,j) = V(j,i);
843  }
844  }
845 
846  for (int i = 0; i != N; ++i) {
847  V(i,i) = 1.0 + lambda;
848  }
849 
850 
851  // solve A x = b
852 
853  for (int col = 0; col != N; ++col) {
854  Y[col] *= h[col];
855  }
856 
857  try {
858  V.solve(Y);
859  }
860  catch (const exception& error) {
861 
862  ERROR("JGandalf: " << error.what() << endl << V << endl);
863 
864  break;
865  }
866 
867  // update value
868 
869  for (int row = 0; row != N; ++row) {
870 
871  DEBUG("u[" << noshowpos << setw(3) << row << "] = " << showpos << FIXED(20,5) << value[row]);
872 
873  value[row] -= h[row] * Y[row];
874 
875  DEBUG(" -> " << FIXED(20,10) << value[row] << noshowpos << endl);
876  }
877  }
878 
879  return precessor;
880  }
size_t getN() const
Get number of fit parameters.
void setOption(const int option)
Set fit option.
static double LAMBDA_UP
multiplication factor control parameter
Definition: JKatoomba_t.hh:887
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:443
static double EPSILON
maximal distance to minimum
Definition: JKatoomba_t.hh:884
static double LAMBDA_MIN
minimal value control parameter
Definition: JKatoomba_t.hh:885
do set_variable OUTPUT_DIRECTORY $WORKDIR T
#define ERROR(A)
Definition: JMessage.hh:66
static double LAMBDA_MAX
maximal value control parameter
Definition: JKatoomba_t.hh:886
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JKatoomba_t.hh:883
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JKatoomba_t.hh:888
void solve(JVectorND_t &u)
Get solution of equation A x = b.
Definition: JMatrixNS.hh:308
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition: JFitK40.hh:99
static double PIVOT
minimal value diagonal element of matrix
Definition: JKatoomba_t.hh:889
int j
Definition: JPolint.hh:792
void evaluate(const data_type &data)
Evaluation of fit.
Definition: JKatoomba_t.hh:902
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
void JACOUSTICS::JKatoomba< JGandalf >::evaluate ( const data_type data)
inlineprivate

Evaluation of fit.

Parameters
datadata

Definition at line 902 of file JKatoomba_t.hh.

903  {
904  using namespace std;
905  using namespace JPP;
906 
907  successor = 0.0;
908 
909  V.reset();
910  Y.reset();
911 
912  for (data_type::const_iterator p = data.begin(); p != data.end(); ++p) {
913 
914  const JGEOMETRY::JString& string = geometry [p->first.getString()];
915  const JMODEL ::JString& parameters = value.string[p->first.getString()];
916  const JPosition3D position = string.getPosition(parameters, p->first.getFloor());
917 
918  for (data_type::mapped_type::const_iterator emitter = p->second.begin(); emitter != p->second.end(); ++emitter) {
919 
920  const double D = emitter->first.getDistance(position);
921  const double Vi = velocity.getInverseVelocity(D, emitter->first.getZ(), position.getZ());
922 
923  const H_t H0(1.0, string.getGradient(parameters, emitter->first.getPosition(), p->first.getFloor()) * Vi);
924 
925  for (data_type::mapped_type::mapped_type::const_iterator hit = emitter->second.begin(); hit != emitter->second.end(); ++hit) {
926 
927  const double toa_s = value.emission[hit->getEKey()].t1 + D * Vi;
928 
929  const double u = (toa_s - hit->getValue()) / hit->getSigma();
930  const double W = sqrt(hit->getWeight());
931 
932  successor += (W*W) * estimator->getRho(u);
933 
934  const H_t H = H0 * (W * estimator->getPsi(u) / hit->getSigma());
935 
936  I_t i;
937 
938  i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
939  i.tx = value.getIndex(hit->getString(), &H_t::tx);
940  i.ty = value.getIndex(hit->getString(), &H_t::ty);
941  i.tx2 = value.getIndex(hit->getString(), &H_t::tx2);
942  i.ty2 = value.getIndex(hit->getString(), &H_t::ty2);
943  i.vs = value.getIndex(hit->getString(), &H_t::vs);
944 
945  V(i.t1, i.t1) += H.t1 * H.t1;
946 
947  Y[i.t1] += W * H.t1;
948 
949  if (hit->getFloor() != 0) {
950 
951  switch (this->option) {
952 
954  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;
955 
956  V(i.vs, i.t1) = V(i.t1, i.vs);
957  V(i.vs, i.tx) = V(i.tx, i.vs);
958  V(i.vs, i.ty) = V(i.ty, i.vs);
959  V(i.vs, i.tx2) = V(i.tx2, i.vs);
960  V(i.vs, i.ty2) = V(i.ty2, i.vs);
961 
962  V(i.vs, i.vs) += H.vs * H.vs;
963 
964  Y[i.vs] += W * H.vs;
965 
967  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;
968 
969  V(i.tx2, i.t1) = V(i.t1, i.tx2);
970  V(i.tx2, i.tx) = V(i.tx, i.tx2);
971  V(i.tx2, i.ty) = V(i.ty, i.tx2);
972 
973  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;
974 
975  V(i.ty2, i.t1) = V(i.t1, i.ty2);
976  V(i.ty2, i.tx) = V(i.tx, i.ty2);
977  V(i.ty2, i.ty) = V(i.ty, i.ty2);
978 
979  V(i.tx2, i.tx2) += H.tx2 * H.tx2; V(i.tx2, i.ty2) += H.tx2 * H.ty2;
980  V(i.ty2, i.tx2) = V(i.tx2, i.ty2); V(i.ty2, i.ty2) += H.ty2 * H.ty2;
981 
982  Y[i.tx2] += W * H.tx2;
983  Y[i.ty2] += W * H.ty2;
984 
986  V(i.t1, i.tx) += H.t1 * H.tx; V(i.t1, i.ty) += H.t1 * H.ty;
987  V(i.tx, i.t1) = V(i.t1, i.tx); V(i.ty, i.t1) = V(i.t1, i.ty);
988 
989  V(i.tx, i.tx) += H.tx * H.tx; V(i.tx, i.ty) += H.tx * H.ty;
990  V(i.ty, i.tx) = V(i.tx, i.ty); V(i.ty, i.ty) += H.ty * H.ty;
991 
992  Y[i.tx] += W * H.tx;
993  Y[i.ty] += W * H.ty;
994  break;
995 
996  default:
997  break;
998  }
999  }
1000  }
1001  }
1002  }
1003  }
fit times of emission of emitters and tilt angles of strings with second order correction ...
void reset()
Reset.
Definition: JVectorND.hh:45
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
Definition: diff-Tuna.sh:38
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:456
std::vector< double > vs
fit times of emission of emitters and tilt angles of strings
JACOUSTICS::JModel::emission_type emission
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
size_t getIndex(int id, double JString::*p) const
Get index of fit parameter for given string.
JACOUSTICS::JModel::string_type string
fit times of emission of emitters and tilt angles of strings with second order correction and stretch...
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
double u[N+1]
Definition: JPolint.hh:865
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
double getZ() const
Get z position.
Definition: JVector3D.hh:115

Member Data Documentation

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

debug level

debug level.

Definition at line 882 of file JKatoomba_t.hh.

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

maximal number of iterations

maximal number of iterations.

Definition at line 883 of file JKatoomba_t.hh.

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

maximal distance to minimum

maximal distance to minimum.

Definition at line 884 of file JKatoomba_t.hh.

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

minimal value control parameter

Definition at line 885 of file JKatoomba_t.hh.

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

maximal value control parameter

Definition at line 886 of file JKatoomba_t.hh.

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

multiplication factor control parameter

Definition at line 887 of file JKatoomba_t.hh.

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

multiplication factor control parameter

Definition at line 888 of file JKatoomba_t.hh.

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

minimal value diagonal element of matrix

Definition at line 889 of file JKatoomba_t.hh.

double JACOUSTICS::JKatoomba< JGandalf >::lambda

Definition at line 891 of file JKatoomba_t.hh.

Definition at line 892 of file JKatoomba_t.hh.

int JACOUSTICS::JKatoomba< JGandalf >::numberOfIterations

Definition at line 893 of file JKatoomba_t.hh.

Definition at line 894 of file JKatoomba_t.hh.

Definition at line 1006 of file JKatoomba_t.hh.

Definition at line 1007 of file JKatoomba_t.hh.

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

Definition at line 1008 of file JKatoomba_t.hh.

Definition at line 1009 of file JKatoomba_t.hh.


The documentation for this struct was generated from the following files: