Jpp  18.0.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JKatoomba.hh
Go to the documentation of this file.
1 #ifndef __JACOUSTICS__JKATOOMBA__
2 #define __JACOUSTICS__JKATOOMBA__
3 
4 #include <vector>
5 #include <memory>
6 #ifdef THREAD_SAFE
7 #include <mutex>
8 #endif
9 #include "TMatrixTSym.h"
10 #include "TDecompSVD.h"
11 
12 #include "JFit/JEstimator.hh"
13 #include "JFit/JRegressor.hh"
14 #include "JFit/JSimplex.hh"
15 #include "JFit/JMEstimator.hh"
16 #include "JMath/JMatrixNS.hh"
17 #include "JLang/JType.hh"
18 #include "JLang/JSharedPointer.hh"
19 #include "JMath/JMatrixNS.hh"
20 #include "Jeep/JMessage.hh"
21 
24 #include "JAcoustics/JEKey.hh"
25 #include "JAcoustics/JGeometry.hh"
26 #include "JAcoustics/JModel.hh"
27 #include "JAcoustics/JHit.hh"
29 
30 
31 /**
32  * \file
33  *
34  * Fit functions of acoustic model.
35  * \author mdejong
36  */
37 namespace JACOUSTICS {}
38 namespace JPP { using namespace JACOUSTICS; }
39 
40 namespace JACOUSTICS {
41 
42  using JLANG::JType;
43  using JMATH::JMath;
44  using JFIT::JEstimator;
46  using JFIT::JSimplex;
47  using JFIT::JMEstimator;
48 
49 
50  /**
51  * Get total weight of data points.
52  *
53  * \param __begin begin of data
54  * \param __end end of data
55  * \return total weight
56  */
57  template<class T>
58  inline double getWeight(T __begin, T __end)
59  {
60  double W = 0;
61 
62  for (T i = __begin; i != __end; ++i) {
63  W += i->getWeight();
64  }
65 
66  return W;
67  }
68 
69 
70  /**
71  * Auxiliary data structure to sort transmissions.
72  */
73  static const struct {
74  /**
75  * Sort transmissions in following order.
76  * - ascending identifier;
77  * - descending quality;
78  * - ascending time-of-arrival;
79  *
80  * \param first first transmission
81  * \param second second transmission
82  * \return true if first transmission before second; else false
83  */
84  inline bool operator()(const JTransmission& first, const JTransmission& second) const
85  {
86  if (first.getID() == second.getID()) {
87 
88  if (first.getQ() == second.getQ())
89  return first.getToA() < second.getToA();
90  else
91  return first.getQ() > second.getQ();
92 
93  } else {
94 
95  return first.getID() < second.getID();
96  }
97  }
98  } compare;
99 
100 
101  /**
102  * Template definition of fit function of acoustic model.
103  */
104  template<template<class T> class JMinimiser_t = JType>
105  struct JKatoomba;
106 
107 
108  /**
109  * Auxiliary base class for fit function of acoustic model.
110  */
111  template<>
112  struct JKatoomba<JType> {
113 
114  typedef std::shared_ptr<JMEstimator> estimator_type;
115 
116 
117  /**
118  * Constructor
119  *
120  * \param geometry detector geometry
121  * \param velocity sound velocity
122  * \param option option
123  */
124  JKatoomba(const JGeometry& geometry,
125  const JSoundVelocity& velocity,
126  const int option) :
127  geometry(geometry),
128  velocity(velocity),
129  option (option)
130  {};
131 
132 
133  /**
134  * Get estimated time-of-arrival for given hit.
135  *
136  * \param model model
137  * \param hit hit
138  * \return time-of-arrival
139  */
140  double getToA(const JModel& model, const JHit& hit) const
141  {
142  const JGEOMETRY::JString& string = geometry[hit.getString()];
143  const JMODEL ::JString& parameters = model.string[hit.getString()];
144  const JPosition3D position = string.getPosition(parameters, hit.getFloor());
145 
146  const double D = hit.getDistance(position);
147  const double Vi = velocity.getInverseVelocity(D, hit.getZ(), position.getZ());
148 
149  return model.emission[hit.getEKey()].t1 + D * Vi;
150  }
151 
152 
153  /**
154  * Get estimated time-of-emission for given hit.
155  *
156  * \param model model
157  * \param hit hit
158  * \return time-of-arrival
159  */
160  double getToE(const JModel& model, const JHit& hit) const
161  {
162  const JGEOMETRY::JString& string = geometry[hit.getString()];
163  const JMODEL ::JString& parameters = model.string[hit.getString()];
164  const JPosition3D position = string.getPosition(parameters, hit.getFloor());
165 
166  const double D = hit.getDistance(position);
167  const double Vi = velocity.getInverseVelocity(D, hit.getZ(), position.getZ());
168 
169  return hit.getValue() - D * Vi;
170  }
171 
172 
173  const JGeometry& geometry; //!< detector geometry
174  JSoundVelocity velocity; //!< sound velocity
175  int option; //!< fit option
176  estimator_type estimator; //!< M-Estimator function
177 
178  protected:
179  /**
180  * H-equation as per hit.
181  */
182  struct H_t :
183  public JMODEL::JEmission,
184  public JMODEL::JString,
185  public JMath<H_t>
186  {
187  /**
188  * Default constructor.
189  */
190  H_t() :
191  JMODEL::JEmission(),
192  JMODEL::JString ()
193  {}
194 
195 
196  /**
197  * Constructor.
198  *
199  * \param emission emission
200  * \param string string
201  */
202  H_t(const JMODEL::JEmission& emission,
203  const JMODEL::JString& string) :
204  JMODEL::JEmission(emission),
205  JMODEL::JString (string)
206  {}
207 
208 
209  /**
210  * Scale H-equation.
211  *
212  * \param factor multiplication factor
213  * \return this H-equation
214  */
215  H_t& mul(const double factor)
216  {
217  static_cast<JMODEL::JEmission&>(*this).mul(factor);
218  static_cast<JMODEL::JString&> (*this).mul(factor);
219 
220  return *this;
221  }
222  };
223 
224 
225  /**
226  * Indices of H-equation in global model.
227  */
228  struct I_t
229  {
230  /**
231  * Default constructor.
232  */
233  I_t() :
234  t1 (0),
235  tx (0),
236  ty (0),
237  tx2(0),
238  ty2(0),
239  vs (0)
240  {}
241 
242 
243  int t1;
244  int tx;
245  int ty;
246  int tx2;
247  int ty2;
248  int vs;
249  };
250  };
251 
252 
253  /**
254  * Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser.\n
255  * This class can be used to evaluate the chi2.
256  */
257  template<>
259  public JAbstractMinimiser<JModel>,
260  public JKatoomba<>
261  {
262  typedef double result_type;
263 
264 
265  /**
266  * Constructor
267  *
268  * \param geometry detector geometry
269  * \param velocity sound velocity
270  * \param option option
271  */
272  JKatoomba(const JGeometry& geometry,
273  const JSoundVelocity& velocity,
274  const int option) :
275  JKatoomba<>(geometry, velocity, option)
276  {};
277 
278 
279  /**
280  * Fit function.\n
281  * This method is used to determine the chi2 of given hit with respect to actual model.
282  *
283  * \param model model
284  * \param hit hit
285  * \return chi2
286  */
287  result_type operator()(const JModel& model, const JHit& hit) const
288  {
289  const double toa_s = this->getToA(model, hit);
290  const double u = (toa_s - hit.getValue()) / hit.getSigma();
291 
292  return estimator->getRho(u) * hit.getWeight();
293  }
294 
295 
296  /**
297  * Fit.
298  *
299  * \param __begin begin of hits
300  * \param __end end of hits
301  * \return chi2
302  */
303  template<class T>
304  result_type operator()(T __begin, T __end)
305  {
306  this->value.setOption(this->option);
307 
308  return JAbstractMinimiser<JModel>::operator()(*this, __begin, __end);
309  }
310 
311 
312  /**
313  * Fit.
314  *
315  * \param model model
316  * \param __begin begin of hits
317  * \param __end end of hits
318  * \return chi2
319  */
320  template<class T>
321  result_type operator()(const JModel& model, T __begin, T __end)
322  {
323  this->value = model;
324 
325  return (*this)(__begin, __end);
326  }
327  };
328 
329 
330  /**
331  * Template specialisation of fit function of acoustic model based on linear approximation.
332  */
333  template<>
335  public JKatoomba<>
336  {
337  static constexpr double TOLERANCE = 1.0e-8; //!< Tolerance for SVD
338 
339  /**
340  * Constructor
341  *
342  * \param geometry detector geometry
343  * \param velocity sound velocity
344  * \param option option
345  */
346  JKatoomba(const JGeometry& geometry,
347  const JSoundVelocity& velocity,
348  const int option) :
349  JKatoomba<>(geometry, velocity, option)
350  {};
351 
352 
353  /**
354  * Get start values of string parameters.\n
355  * Note that this method may throw an exception.
356  *
357  * \param __begin begin of hits
358  * \param __end end of hits
359  * \return model
360  */
361  template<class T>
362  const JModel& operator()(T __begin, T __end) const
363  {
364  using namespace std;
365  using namespace JPP;
366  using namespace JGEOMETRY;
367 
368 #ifdef THREAD_SAFE
369  static mutex mtx; // TDecompSVD
370 
371  unique_lock<mutex> lock(mtx);
372 #endif
373  value = JModel(__begin, __end); // set up global model according data
374 
375  if (this->option == JMODEL::fit_emitters_only_t)
376  value.setOption(JMODEL::fit_emitters_only_t);
377  else
379 
380  H_t H; // H-equation as per hit
381  I_t i; // indices of H-equation in global model
382 
383  const size_t N = value.getN();
384 
385  TMatrixD V(N,N); // V = H^T x H
386  vector<double> Y(N, 0.0); // Y = H^T x y
387 
388  for (T hit = __begin; hit != __end; ++hit) {
389 
390  const JString& string = geometry[hit->getString()];
391  const JPosition3D position = string.getPosition(hit->getFloor());
392  const double Vi = velocity.getInverseVelocity(hit->getDistance(position), hit->getZ(), position.getZ());
393 
394  const double h1 = string.getHeight(hit->getFloor());
395  const JPosition3D p1 = string.getPosition() - hit->getPosition();
396  const double ds = sqrt(p1.getLengthSquared() + h1*h1 + 2.0*p1.getZ()*h1);
397  const double y = hit->getValue() - Vi*ds;
398  const double W = sqrt(hit->getWeight());
399 
400  H.t1 = W * 1.0;
401  H.tx = W * Vi * p1.getX() * h1 / ds;
402  H.ty = W * Vi * p1.getY() * h1 / ds;
403 
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);
407 
408  V(i.t1, i.t1) += H.t1 * H.t1;
409 
410  Y[i.t1] += W * H.t1 * y;
411 
412  if (hit->getFloor() != 0) {
413 
414  if (this->option != JMODEL::fit_emitters_only_t) {
415 
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;
418 
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;
421 
422  Y[i.tx] += W * H.tx * y;
423  Y[i.ty] += W * H.ty * y;
424  }
425  }
426  }
427 
428  // evaluate (H^T x H)^-1 x H^T x y
429 
430  TDecompSVD svd(V, TOLERANCE);
431 
432  Bool_t status;
433 
434  V = svd.Invert(status);
435 
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];
439  }
440  }
441 
442  return value;
443  }
444 
445  private:
446  mutable JModel value;
447 #ifdef THREAD_SAFE
448  static std::mutex mtx; // TDecompSVD
449 #endif
450  };
451 
452 #ifdef THREAD_SAFE
453  /**
454  * Declaration of mutex.
455  */
456  std::mutex JKatoomba<JEstimator>::mtx;
457 #endif
458 
459  /**
460  * Template specialisation of fit function of acoustic model based on JSimplex minimiser.
461  */
462  template<>
463  struct JKatoomba<JSimplex> :
464  public JSimplex<JModel>,
465  public JKatoomba<>
466  {
467  /**
468  * Constructor
469  *
470  * \param geometry detector geometry
471  * \param velocity sound velocity
472  * \param option option
473  */
474  JKatoomba(const JGeometry& geometry,
475  const JSoundVelocity& velocity,
476  const int option) :
477  JKatoomba<>(geometry, velocity, option)
478  {};
479 
480 
481  /**
482  * Fit function.\n
483  * This method is used to determine the chi2 of given hit with respect to actual model.
484  *
485  * \param model model
486  * \param hit hit
487  * \return chi2
488  */
489  double operator()(const JModel& model, const JHit& hit) const
490  {
491  const double toa_s = this->getToA(model, hit);
492  const double u = (toa_s - hit.getValue()) / hit.getSigma();
493 
494  return estimator->getRho(u);
495  }
496 
497 
498  /**
499  * Fit.
500  *
501  * \param __begin begin of hits
502  * \param __end end of hits
503  * \return chi2
504  */
505  template<class T>
506  double operator()(T __begin, T __end)
507  {
508  this->step.clear();
509 
510  for (JModel::string_type::const_iterator i = this->value.string.begin(); i != this->value.string.end(); ++i) {
511 
512  JModel model;
513 
514  switch (this->option) {
515 
517  model.string[i->first] = JMODEL::JString(0.0, 0.0, 0.0, 0.0, 1.0e-6); this->step.push_back(model);
518 
520  model.string[i->first] = JMODEL::JString(0.0, 0.0, 1.0e-6, 0.0, 0.0); this->step.push_back(model);
521  model.string[i->first] = JMODEL::JString(0.0, 0.0, 0.0, 1.0e-6, 0.0); this->step.push_back(model);
522 
524  model.string[i->first] = JMODEL::JString(1.0e-3, 0.0, 0.0, 0.0, 0.0); this->step.push_back(model);
525  model.string[i->first] = JMODEL::JString(0.0, 1.0e-3, 0.0, 0.0, 0.0); this->step.push_back(model);
526  break;
527 
528  default:
529  break;
530  }
531  }
532 
533  for (JModel::emission_type::const_iterator i = this->value.emission.begin(); i != this->value.emission.end(); ++i) {
534 
535  JModel model;
536 
537  model.emission[i->first] = JMODEL::JEmission(5.0e-5); this->step.push_back(model);
538  }
539 
540  return JSimplex<JModel>::operator()(*this, __begin, __end);
541  }
542  };
543 
544 
545  /**
546  * Place holder for custom implementation.
547  */
548  template<class T>
549  class JGandalf {};
550 
551 
552  /**
553  * Template specialisation of fit function of acoustic model based on JGandalf minimiser.
554  */
555  template<>
556  struct JKatoomba<JGandalf> :
557  public JKatoomba<>
558  {
559  typedef double result_type;
560 
561 
562  /**
563  * Constructor
564  *
565  * \param geometry detector geometry
566  * \param velocity sound velocity
567  * \param option option
568  */
569  JKatoomba(const JGeometry& geometry,
570  const JSoundVelocity& velocity,
571  const int option) :
572  JKatoomba<>(geometry, velocity, option)
573  {};
574 
575 
576  /**
577  * Fit.
578  *
579  * \param __begin begin of hits
580  * \param __end end of hits
581  * \return chi2 and gradient
582  */
583  template<class T>
584  result_type operator()(T __begin, T __end)
585  {
586  using namespace std;
587  using namespace JPP;
588 
589  value.setOption(this->option);
590 
591  const int N = value.getN();
592 
593  V.resize(N);
594  Y.resize(N);
595  h.resize(N);
596 
597  lambda = LAMBDA_MIN;
598 
599  result_type precessor = numeric_limits<double>::max();
600 
601  for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
602 
603  DEBUG("step: " << numberOfIterations << endl);
604 
605  evaluate(__begin, __end);
606 
607  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
608  DEBUG("chi2: " << FIXED(12,5) << successor << endl);
609 
610  if (successor < precessor) {
611 
612  if (numberOfIterations != 0) {
613 
614  if (fabs(precessor - successor) < EPSILON*fabs(precessor)) {
615  return successor;
616  }
617 
618  if (lambda > LAMBDA_MIN) {
619  lambda /= LAMBDA_DOWN;
620  }
621  }
622 
623  precessor = successor;
624  previous = value;
625 
626  } else {
627 
628  value = previous;
629  lambda *= LAMBDA_UP;
630 
631  if (lambda > LAMBDA_MAX) {
632  return precessor; // no improvement found
633  }
634 
635  evaluate(__begin, __end);
636  }
637 
638 
639  // force definite positiveness
640 
641  for (int i = 0; i != N; ++i) {
642 
643  if (V(i,i) < PIVOT) {
644  V(i,i) = PIVOT;
645  }
646 
647  h[i] = 1.0 / sqrt(V(i,i));
648  }
649 
650 
651  // normalisation
652 
653  for (int i = 0; i != N; ++i) {
654  for (int j = 0; j != i; ++j) {
655  V(j,i) *= h[i] * h[j];
656  V(i,j) = V(j,i);
657  }
658  }
659 
660  for (int i = 0; i != N; ++i) {
661  V(i,i) = 1.0 + lambda;
662  }
663 
664 
665  try {
666  V.invert();
667  }
668  catch (const JException& error) {
669  ERROR("JKatoomb<JGandalf>: " << error.what() << endl);
670  return precessor;
671  }
672 
673 
674  for (int i = 0; i != N; ++i) {
675 
676  DEBUG("u[" << noshowpos << setw(3) << i << "] = " << showpos << FIXED(20,5) << value[i]);
677 
678  double y = 0.0;
679 
680  for (int j = 0; j != N; ++j) {
681  y += V(i,j) * Y[j] * h[i] * h[j];
682  }
683 
684  value[i] -= y;
685 
686  DEBUG(' ' << FIXED(20,10) << y << noshowpos << endl);
687  }
688  }
689 
690  return precessor;
691  }
692 
693  static int debug; //!< debug level
694  static int MAXIMUM_ITERATIONS; //!< maximal number of iterations
695  static double EPSILON; //!< maximal distance to minimum
696  static double LAMBDA_MIN; //!< minimal value control parameter
697  static double LAMBDA_MAX; //!< maximal value control parameter
698  static double LAMBDA_UP; //!< multiplication factor control parameter
699  static double LAMBDA_DOWN; //!< multiplication factor control parameter
700  static double PIVOT; //!< minimal value diagonal element of matrix
701 
702  double lambda;
706 
707  private:
708  /**
709  * Evaluation of fit.
710  *
711  * \param __begin begin of data
712  * \param __end end of data
713  */
714  template<class T>
715  inline void evaluate(T __begin, T __end)
716  {
717  using namespace std;
718  using namespace JPP;
719 
720  successor = 0.0;
721 
722  V.reset();
723 
724  for (std::vector<double>::iterator i = Y.begin(); i != Y.end(); ++i) {
725  *i = 0.0;
726  }
727 
728  for (T hit = __begin; hit != __end; ++hit) {
729 
730  const JGEOMETRY::JString& string = geometry[hit->getString()];
731  const JMODEL ::JString& parameters = value.string[hit->getString()];
732  const JPosition3D position = string.getPosition(parameters, hit->getFloor());
733 
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;
737 
738  const double u = (toa_s - hit->getValue()) / hit->getSigma();
739  const double W = sqrt(hit->getWeight());
740 
741  successor += (W*W) * estimator->getRho(u);
742 
743  H_t H(1.0, string.getGradient(parameters, hit->getPosition(), hit->getFloor()) * Vi);
744 
745  H *= W * estimator->getPsi(u) / hit->getSigma();
746 
747  I_t i;
748 
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);
755 
756  V(i.t1, i.t1) += H.t1 * H.t1;
757 
758  Y[i.t1] += W * H.t1;
759 
760  if (hit->getFloor() != 0) {
761 
762  switch (this->option) {
763 
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;
766 
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;
772 
773  V(i.vs, i.vs) += H.vs * H.vs;
774 
775  Y[i.vs] += W * H.vs;
776 
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;
779 
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;
783 
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;
785 
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;
789 
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;
792 
793  Y[i.tx2] += W * H.tx2;
794  Y[i.ty2] += W * H.ty2;
795 
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;
799 
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;
802 
803  Y[i.tx] += W * H.tx;
804  Y[i.ty] += W * H.ty;
805  break;
806 
807  default:
808  break;
809  }
810  }
811  }
812  }
813 
814 
815  std::vector<double> Y; // gradient
818  std::vector<double> h; // normalisation vector
819  };
820 
821 
822  /**
823  * debug level.
824  */
826 
827 
828  /**
829  * maximal number of iterations.
830  */
832 
833 
834  /**
835  * maximal distance to minimum.
836  */
837  double JKatoomba<JGandalf>::EPSILON = 1.0e-3;
838 
839 
840  /**
841  * minimal value control parameter
842  */
843  double JKatoomba<JGandalf>::LAMBDA_MIN = 0.01;
844 
845 
846  /**
847  * maximal value control parameter
848  */
849  double JKatoomba<JGandalf>::LAMBDA_MAX = 100.0;
850 
851 
852  /**
853  * multiplication factor control parameter
854  */
855  double JKatoomba<JGandalf>::LAMBDA_UP = 9.0;
856 
857 
858  /**
859  * multiplication factor control parameter
860  */
861  double JKatoomba<JGandalf>::LAMBDA_DOWN = 11.0;
862 
863 
864  /**
865  * minimal value diagonal element of matrix
866  */
867  double JKatoomba<JGandalf>::PIVOT = 1.0e-3;
868 }
869 
870 #endif
fit times of emission of emitters and tilt angles of strings with second order correction and stretch...
Linear fit methods.
Acoustic hit.
I_t()
Default constructor.
Definition: JKatoomba.hh:233
General exception.
Definition: JException.hh:23
General purpose data regression method.
JString & mul(const double factor)
Scale string.
Wrapper class around STL string class.
Definition: JString.hh:27
Auxiliary base class for aritmetic operations of derived class types.
Definition: JMath.hh:109
double getQ() const
Get quality.
TPaveText * p1
double getWeight(T __begin, T __end)
Get total weight of data points.
Definition: JKatoomba.hh:58
Sound velocity.
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
static double LAMBDA_UP
multiplication factor control parameter
Definition: JKatoomba.hh:698
Interface for maximum likelihood estimator (M-estimator).
Definition: JMEstimator.hh:20
double getWeight() const
Get weight.
const JGeometry & geometry
detector geometry
Definition: JKatoomba.hh:173
result_type operator()(const JFunction_t &fit, T __begin, T __end)
Get chi2.
Definition: JRegressor.hh:46
Acoustic geometries.
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
H_t & mul(const double factor)
Scale H-equation.
Definition: JKatoomba.hh:215
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.
Definition: JType.hh:19
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.
Definition: JManip.hh:446
Abstract minimiser.
Definition: JRegressor.hh:25
Template definition of linear fit.
Definition: JEstimator.hh:25
H_t()
Default constructor.
Definition: JKatoomba.hh:190
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
fit times of emission of emitters and tilt angles of strings
Acoustics hit.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
Place holder for custom implementation.
Definition: JKatoomba.hh:549
Acoustic fit parameters.
Model for fit to acoustics data.
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor.
Definition: JKatoomba.hh:124
Acoustic transmission.
static double EPSILON
maximal distance to minimum
Definition: JKatoomba.hh:695
double operator()(const JModel &model, const JHit &hit) const
Fit function.
Definition: JKatoomba.hh:489
result_type operator()(T __begin, T __end)
Fit.
Definition: JKatoomba.hh:584
std::vector< double > vs
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor.
Definition: JKatoomba.hh:474
static double LAMBDA_MIN
minimal value control parameter
Definition: JKatoomba.hh:696
void evaluate(T __begin, T __end)
Evaluation of fit.
Definition: JKatoomba.hh:715
static int debug
debug level
Definition: JKatoomba.hh:693
do set_variable OUTPUT_DIRECTORY $WORKDIR T
double operator()(T __begin, T __end)
Fit.
Definition: JKatoomba.hh:506
result_type operator()(T __begin, T __end)
Fit.
Definition: JKatoomba.hh:304
JACOUSTICS::JModel::emission_type emission
#define ERROR(A)
Definition: JMessage.hh:66
static struct JACOUSTICS::@4 compare
Auxiliary data structure to sort transmissions.
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor.
Definition: JKatoomba.hh:569
double getY() const
Get y position.
Definition: JVector3D.hh:104
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
N x N symmetric matrix.
Definition: JMatrixNS.hh:27
fit times of emission of emitters and tilt angles of strings with second order correction ...
double getLengthSquared() const
Get length squared.
Definition: JVector3D.hh:235
General purpose messaging.
Implementation for depth dependend velocity of sound.
static double LAMBDA_MAX
maximal value control parameter
Definition: JKatoomba.hh:697
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JKatoomba.hh:694
std::vector< double > h
Definition: JKatoomba.hh:818
JEmission & mul(const double factor)
Scale emission.
double operator()(const JFunction_t &fit, T __begin, T __end)
Multi-dimensional fit.
Definition: JSimplex.hh:71
Emitter hash key.
double getToE(const JModel &model, const JHit &hit) const
Get estimated time-of-emission for given hit.
Definition: JKatoomba.hh:160
JACOUSTICS::JModel::string_type string
JSoundVelocity velocity
sound velocity
Definition: JKatoomba.hh:174
std::shared_ptr< JMEstimator > estimator_type
Definition: JKatoomba.hh:114
H_t(const JMODEL::JEmission &emission, const JMODEL::JString &string)
Constructor.
Definition: JKatoomba.hh:202
int getString() const
Get string number.
Definition: JLocation.hh:134
virtual const char * what() const override
Get error message.
Definition: JException.hh:48
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:36
const JModel & operator()(T __begin, T __end) const
Get start values of string parameters.
Definition: JKatoomba.hh:362
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JKatoomba.hh:699
result_type operator()(const JModel &model, const JHit &hit) const
Fit function.
Definition: JKatoomba.hh:287
Simple fit method based on Powell&#39;s algorithm, see reference: Numerical Recipes in C++...
Definition: JSimplex.hh:42
result_type operator()(const JModel &model, T __begin, T __end)
Fit.
Definition: JKatoomba.hh:321
double getToA() const
Get calibrated time of arrival.
Acoustic transmission.
double getToA(const JModel &model, const JHit &hit) const
Get estimated time-of-arrival for given hit.
Definition: JKatoomba.hh:140
static double PIVOT
minimal value diagonal element of matrix
Definition: JKatoomba.hh:700
int getID() const
Get identifier.
double getX() const
Get x position.
Definition: JVector3D.hh:94
int j
Definition: JPolint.hh:703
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
container_type::const_iterator const_iterator
Definition: JHashMap.hh:86
double u[N+1]
Definition: JPolint.hh:776
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor.
Definition: JKatoomba.hh:272
Model for fit to acoutsics data.
fit only times of emission of emitters
estimator_type estimator
M-Estimator function.
Definition: JKatoomba.hh:176
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
double getZ() const
Get z position.
Definition: JVector3D.hh:115
JKatoomba(const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
Constructor.
Definition: JKatoomba.hh:346
Maximum likelihood estimator (M-estimators).
Template definition of fit function of acoustic model.
Definition: JKatoomba.hh:105
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
double getSigma() const
Get resolution of time-of-arrival.
std::vector< double > Y
Definition: JKatoomba.hh:815