Jpp  18.2.0-rc.1
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

std::vector< double > 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 735 of file JKatoomba_t.hh.

Member Typedef Documentation

Definition at line 738 of file JKatoomba_t.hh.

Type definition internal data structure.

Definition at line 743 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 753 of file JKatoomba_t.hh.

755  :
756  JKatoomba<>(geometry, velocity, option)
757  {};

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

769  {
770  using namespace std;
771  using namespace JPP;
772 
773  value.setOption(this->option);
774 
775  const int N = value.getN();
776 
777  V.resize(N);
778  Y.resize(N);
779  h.resize(N);
780 
781  data_type data;
782 
783  for (T hit = __begin; hit != __end; ++hit) {
784  data[hit->getLocation()][hit->getEmitter()].push_back(*hit);
785  }
786 
787  lambda = LAMBDA_MIN;
788 
789  result_type precessor = numeric_limits<double>::max();
790 
792 
793  DEBUG("step: " << numberOfIterations << endl);
794 
795  evaluate(data);
796 
797  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
798  DEBUG("chi2: " << FIXED(12,5) << successor << endl);
799 
800  if (successor < precessor) {
801 
802  if (numberOfIterations != 0) {
803 
804  if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
805  return successor;
806  }
807 
808  if (lambda > LAMBDA_MIN) {
809  lambda /= LAMBDA_DOWN;
810  }
811  }
812 
813  precessor = successor;
814  previous = value;
815 
816  } else {
817 
818  value = previous;
819  lambda *= LAMBDA_UP;
820 
821  if (lambda > LAMBDA_MAX) {
822  return precessor; // no improvement found
823  }
824 
825  evaluate(data);
826  }
827 
828  // force definite positiveness
829 
830  for (int i = 0; i != N; ++i) {
831 
832  if (V(i,i) < PIVOT) {
833  V(i,i) = PIVOT;
834  }
835 
836  h[i] = 1.0 / sqrt(V(i,i));
837  }
838 
839  // normalisation
840 
841  for (int i = 0; i != N; ++i) {
842  for (int j = 0; j != i; ++j) {
843  V(j,i) *= h[i] * h[j];
844  V(i,j) = V(j,i);
845  }
846  }
847 
848  for (int i = 0; i != N; ++i) {
849  V(i,i) = 1.0 + lambda;
850  }
851 
852 
853  // solve A x = b
854 
855  for (int col = 0; col != N; ++col) {
856  Y[col] *= h[col];
857  }
858 
859  try {
860  V.solve(Y);
861  }
862  catch (const exception& error) {
863 
864  ERROR("JGandalf: " << error.what() << endl << V << endl);
865 
866  break;
867  }
868 
869  // update value
870 
871  for (int row = 0; row != N; ++row) {
872 
873  DEBUG("u[" << noshowpos << setw(3) << row << "] = " << showpos << FIXED(20,5) << value[row]);
874 
875  value[row] -= h[row] * Y[row];
876 
877  DEBUG(" -> " << FIXED(20,10) << value[row] << noshowpos << endl);
878  }
879  }
880 
881  return precessor;
882  }
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:889
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:886
static double LAMBDA_MIN
minimal value control parameter
Definition: JKatoomba_t.hh:887
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:888
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JKatoomba_t.hh:885
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:890
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:98
static double PIVOT
minimal value diagonal element of matrix
Definition: JKatoomba_t.hh:891
int j
Definition: JPolint.hh:703
void evaluate(const data_type &data)
Evaluation of fit.
Definition: JKatoomba_t.hh:904
#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 904 of file JKatoomba_t.hh.

905  {
906  using namespace std;
907  using namespace JPP;
908 
909  successor = 0.0;
910 
911  V.reset();
912 
913  for (std::vector<double>::iterator i = Y.begin(); i != Y.end(); ++i) {
914  *i = 0.0;
915  }
916 
917  for (data_type::const_iterator p = data.begin(); p != data.end(); ++p) {
918 
919  const JGEOMETRY::JString& string = geometry [p->first.getString()];
920  const JMODEL ::JString& parameters = value.string[p->first.getString()];
921  const JPosition3D position = string.getPosition(parameters, p->first.getFloor());
922 
923  for (data_type::mapped_type::const_iterator emitter = p->second.begin(); emitter != p->second.end(); ++emitter) {
924 
925  const double D = emitter->first.getDistance(position);
926  const double Vi = velocity.getInverseVelocity(D, emitter->first.getZ(), position.getZ());
927 
928  const H_t H0(1.0, string.getGradient(parameters, emitter->first.getPosition(), p->first.getFloor()) * Vi);
929 
930  for (data_type::mapped_type::mapped_type::const_iterator hit = emitter->second.begin(); hit != emitter->second.end(); ++hit) {
931 
932  const double toa_s = value.emission[hit->getEKey()].t1 + D * Vi;
933 
934  const double u = (toa_s - hit->getValue()) / hit->getSigma();
935  const double W = sqrt(hit->getWeight());
936 
937  successor += (W*W) * estimator->getRho(u);
938 
939  const H_t H = H0 * (W * estimator->getPsi(u) / hit->getSigma());
940 
941  I_t i;
942 
943  i.t1 = value.getIndex(hit->getEKey(), &H_t::t1);
944  i.tx = value.getIndex(hit->getString(), &H_t::tx);
945  i.ty = value.getIndex(hit->getString(), &H_t::ty);
946  i.tx2 = value.getIndex(hit->getString(), &H_t::tx2);
947  i.ty2 = value.getIndex(hit->getString(), &H_t::ty2);
948  i.vs = value.getIndex(hit->getString(), &H_t::vs);
949 
950  V(i.t1, i.t1) += H.t1 * H.t1;
951 
952  Y[i.t1] += W * H.t1;
953 
954  if (hit->getFloor() != 0) {
955 
956  switch (this->option) {
957 
959  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;
960 
961  V(i.vs, i.t1) = V(i.t1, i.vs);
962  V(i.vs, i.tx) = V(i.tx, i.vs);
963  V(i.vs, i.ty) = V(i.ty, i.vs);
964  V(i.vs, i.tx2) = V(i.tx2, i.vs);
965  V(i.vs, i.ty2) = V(i.ty2, i.vs);
966 
967  V(i.vs, i.vs) += H.vs * H.vs;
968 
969  Y[i.vs] += W * H.vs;
970 
972  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;
973 
974  V(i.tx2, i.t1) = V(i.t1, i.tx2);
975  V(i.tx2, i.tx) = V(i.tx, i.tx2);
976  V(i.tx2, i.ty) = V(i.ty, i.tx2);
977 
978  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;
979 
980  V(i.ty2, i.t1) = V(i.t1, i.ty2);
981  V(i.ty2, i.tx) = V(i.tx, i.ty2);
982  V(i.ty2, i.ty) = V(i.ty, i.ty2);
983 
984  V(i.tx2, i.tx2) += H.tx2 * H.tx2; V(i.tx2, i.ty2) += H.tx2 * H.ty2;
985  V(i.ty2, i.tx2) = V(i.tx2, i.ty2); V(i.ty2, i.ty2) += H.ty2 * H.ty2;
986 
987  Y[i.tx2] += W * H.tx2;
988  Y[i.ty2] += W * H.ty2;
989 
991  V(i.t1, i.tx) += H.t1 * H.tx; V(i.t1, i.ty) += H.t1 * H.ty;
992  V(i.tx, i.t1) = V(i.t1, i.tx); V(i.ty, i.t1) = V(i.t1, i.ty);
993 
994  V(i.tx, i.tx) += H.tx * H.tx; V(i.tx, i.ty) += H.tx * H.ty;
995  V(i.ty, i.tx) = V(i.tx, i.ty); V(i.ty, i.ty) += H.ty * H.ty;
996 
997  Y[i.tx] += W * H.tx;
998  Y[i.ty] += W * H.ty;
999  break;
1000 
1001  default:
1002  break;
1003  }
1004  }
1005  }
1006  }
1007  }
1008  }
fit times of emission of emitters and tilt angles of strings with second order correction ...
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:776
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 884 of file JKatoomba_t.hh.

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

maximal number of iterations

maximal number of iterations.

Definition at line 885 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 886 of file JKatoomba_t.hh.

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

minimal value control parameter

Definition at line 887 of file JKatoomba_t.hh.

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

maximal value control parameter

Definition at line 888 of file JKatoomba_t.hh.

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

multiplication factor control parameter

Definition at line 889 of file JKatoomba_t.hh.

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

multiplication factor control parameter

Definition at line 890 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 891 of file JKatoomba_t.hh.

double JACOUSTICS::JKatoomba< JGandalf >::lambda

Definition at line 893 of file JKatoomba_t.hh.

Definition at line 894 of file JKatoomba_t.hh.

int JACOUSTICS::JKatoomba< JGandalf >::numberOfIterations

Definition at line 895 of file JKatoomba_t.hh.

Definition at line 896 of file JKatoomba_t.hh.

Definition at line 1011 of file JKatoomba_t.hh.

Definition at line 1012 of file JKatoomba_t.hh.

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

Definition at line 1013 of file JKatoomba_t.hh.

Definition at line 1014 of file JKatoomba_t.hh.


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