Jpp  master_rocky
the software that should make you happy
JGlobalfit.hh
Go to the documentation of this file.
1 #ifndef __JACOUSTICS__JGLOBALFIT__
2 #define __JACOUSTICS__JGLOBALFIT__
3 
4 #include "JFit/JEstimator.hh"
5 
9 #include "JAcoustics/JModel.hh"
10 #include "JAcoustics/JHit.hh"
12 #include "JAcoustics/JTimeRange.hh"
13 
14 
15 /**
16  * \file
17  *
18  * Fit function of acoustic model.
19  * \author mdejong
20  */
21 namespace JACOUSTICS {}
22 namespace JPP { using namespace JACOUSTICS; }
23 
24 namespace JACOUSTICS {
25 
26  /**
27  * Global fit of prameterised detector geometry to acoustics data.
28  */
29  struct JGlobalfit {
30  /**
31  * Constructor.
32  *
33  * \param geometry detector geometry
34  * \param V sound velocity
35  * \param parameters parameters
36  */
37  JGlobalfit(const JGeometry& geometry, const JSoundVelocity& V, const JFitParameters& parameters) :
39  estimator(geometry, V, parameters.option),
40  evaluator(geometry, V, parameters.option),
41  gandalf (geometry, V, parameters.option)
42  {
43  using namespace JPP;
44 
45  evaluator.estimator.reset(getMEstimator(parameters.mestimator));
46  gandalf .estimator.reset(getMEstimator(parameters.mestimator));
47  }
48 
49 
50  /**
51  * Auxiliary data structure for return value of fit.
52  */
53  template<class T>
54  struct result_type {
56  double chi2;
57  double ndf;
58  T begin;
59  T end;
60 
61  /**
62  * Get number of hits used in fit (after outlier removal).
63  *
64  * \return size
65  */
66  int size() const
67  {
68  return std::distance(begin, end);
69  }
70 
71  /**
72  * Get time range of data.
73  *
74  * \return time range
75  */
77  {
79 
80  for (T i = begin; i != end; ++i) {
81  range.include(i->getValue());
82  }
83 
84  return range;
85  }
86  };
87 
88 
89  /**
90  * Fit.
91  *
92  * \param __begin begin of data
93  * \param __end end of data
94  * \return fit result
95  */
96  template<class T>
97  result_type<T> operator()(T __begin, T __end)
98  {
99  using namespace std;
100 
101  JModel model;
102  JModel& result = gandalf.value; // start value
103 
104  result = estimator(__begin, __end); // pre-fit
105 
106  if (parameters.stdev > 0.0) { // remove outliers
107 
108  double chi2[] = { numeric_limits<double>::max(),
109  numeric_limits<double>::max() };
110 
111  for ( ; ; ) {
112 
113  T out = __end;
114 
115  double xmax = 0.0;
116 
117  for (T hit = __begin; hit != __end; ++hit) {
118 
119  const double x = fabs(hit->getValue() - estimator.getToA(result, *hit)) / hit->getSigma();
120 
121  if (x > xmax) {
122  xmax = x;
123  out = hit;
124  }
125  }
126 
127  if (xmax > parameters.stdev) {
128 
129  if (chi2[0] == numeric_limits<double>::max()) {
130  chi2[0] = evaluator(result, __begin, __end);
131  }
132 
133  model = result;
134 
135  iter_swap(out, --__end);
136 
137  result = estimator(__begin, __end);
138  chi2[1] = evaluator(result, __begin, __end);
139 
140  if (chi2[0] - chi2[1] > parameters.stdev * parameters.stdev) {
141 
142  chi2[0] = chi2[1];
143 
144  continue;
145 
146  } else {
147 
148  iter_swap(out, __end++);
149 
150  result = model;
151  }
152  }
153 
154  break;
155  }
156  }
157 
158  const double chi2 = gandalf (__begin, __end) / gandalf.estimator->getRho(1.0);
159  const double ndf = getWeight(__begin, __end) - gandalf.value.getN();
160 
161  return { gandalf.value, chi2, ndf, __begin, __end }; // return values
162  }
163 
164 
169  };
170 }
171 
172 #endif
Acoustic hit.
Model for fit to acoutsics data.
Linear fit methods.
Acoustic fit parameters.
Acoustic geometries.
Fit functions of acoustic model.
Sound velocity.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
range_type & include(argument_type x)
Include given value to range.
Definition: JRange.hh:397
static JRange< double, std::less< double > > DEFAULT_RANGE()
Default range.
Definition: JRange.hh:555
const double xmax
Definition: JQuadrature.cc:24
Auxiliary classes and methods for acoustic position calibration.
double getWeight(T __begin, T __end)
Get total weight of data points.
Definition: JKatoomba_t.hh:61
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:203
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JSTDTypes.hh:14
double stdev
standard deviation for outlier removal
Auxiliary data structure for return value of fit.
Definition: JGlobalfit.hh:54
JTimeRange getTimeRange() const
Get time range of data.
Definition: JGlobalfit.hh:76
int size() const
Get number of hits used in fit (after outlier removal).
Definition: JGlobalfit.hh:66
Global fit of prameterised detector geometry to acoustics data.
Definition: JGlobalfit.hh:29
result_type< T > operator()(T __begin, T __end)
Fit.
Definition: JGlobalfit.hh:97
JKatoomba< JEstimator > estimator
Definition: JGlobalfit.hh:166
JKatoomba< JGandalf > gandalf
Definition: JGlobalfit.hh:168
JGlobalfit(const JGeometry &geometry, const JSoundVelocity &V, const JFitParameters &parameters)
Constructor.
Definition: JGlobalfit.hh:37
JKatoomba< JAbstractMinimiser > evaluator
Definition: JGlobalfit.hh:167
JFitParameters parameters
Definition: JGlobalfit.hh:165
Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser.
Definition: JKatoomba_t.hh:272
Template specialisation of fit function of acoustic model based on linear approximation.
Definition: JKatoomba_t.hh:439
Template specialisation of fit function of acoustic model based on JGandalf minimiser.
Definition: JKatoomba_t.hh:680
Model for fit to acoustics data.
size_t getN() const
Get number of fit parameters.
Implementation for depth dependend velocity of sound.