Jpp  18.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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) :
38  parameters(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;
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  model = result;
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
const double xmax
Definition: JQuadrature.cc:24
Linear fit methods.
Acoustic hit.
size_t getN() const
Get number of fit parameters.
double getWeight(T __begin, T __end)
Get total weight of data points.
Definition: JKatoomba_t.hh:59
Sound velocity.
range_type & include(argument_type x)
Include given value to range.
Definition: JRange.hh:397
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
int size() const
Get number of hits used in fit (after outlier removal).
Definition: JGlobalfit.hh:66
Template specialisation of fit function of acoustic model based on linear approximation.
Definition: JKatoomba_t.hh:364
Acoustic geometries.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
JTimeRange getTimeRange() const
Get time range of data.
Definition: JGlobalfit.hh:76
Global fit of prameterised detector geometry to acoustics data.
Definition: JGlobalfit.hh:29
Acoustic fit parameters.
double stdev
standard deviation for outlier removal
Model for fit to acoustics data.
Template specialisation of fit function of acoustic model based on JGandalf minimiser.
Definition: JKatoomba_t.hh:629
do set_variable OUTPUT_DIRECTORY $WORKDIR T
JGlobalfit(const JGeometry &geometry, const JSoundVelocity &V, const JFitParameters &parameters)
Constructor.
Definition: JGlobalfit.hh:37
JFitParameters parameters
Definition: JGlobalfit.hh:165
Auxiliary data structure for return value of fit.
Definition: JGlobalfit.hh:54
Implementation for depth dependend velocity of sound.
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
Template specialisation of fit function of acoustic model based on JAbstractMinimiser minimiser...
Definition: JKatoomba_t.hh:259
Fit functions of acoustic model.
result_type< T > operator()(T __begin, T __end)
Fit.
Definition: JGlobalfit.hh:97
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Definition: JMatrixNZ.sh:106
int mestimator
M-estimator.
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:203
static JRange< double, std::less< double > > DEFAULT_RANGE()
Default range.
Definition: JRange.hh:555
JKatoomba< JGandalf > gandalf
Definition: JGlobalfit.hh:168
JKatoomba< JEstimator > estimator
Definition: JGlobalfit.hh:166
Model for fit to acoutsics data.
JKatoomba< JAbstractMinimiser > evaluator
Definition: JGlobalfit.hh:167