Jpp
JFitToolkit.hh
Go to the documentation of this file.
1 #ifndef __JFIT__JFITTOOLKIT__
2 #define __JFIT__JFITTOOLKIT__
3 
4 #include <cmath>
5 
6 #include "TMath.h"
7 
8 #include "JMath/JMath.hh"
9 #include "JFit/JVectorNZ.hh"
10 #include "JFit/JMatrixNZ.hh"
11 #include "JFit/JEvt.hh"
12 
13 
14 /**
15  * \file
16  *
17  * Auxiliary methods to evaluate Poisson probabilities and chi2.
18  * \author mdejong
19  */
20 
21 /**
22  * Auxiliary classes and methods for linear and iterative data regression.
23  */
24 namespace JFIT {}
25 namespace JPP { using namespace JFIT; }
26 
27 namespace JFIT {
28 
29  /**
30  * Get Poisson probability to observe a hit or not
31  * for given expectation value for the number of hits.
32  *
33  * \param expval expectation value
34  * \param hit hit
35  * \return probability
36  */
37  inline double getP(const double expval, bool hit)
38  {
39  if (hit)
40  return 1.0 - exp(-expval); // 1 - P(0)
41  else
42  return exp(-expval); // P(0)
43  }
44 
45 
46  /**
47  * Get Poisson probability to observe given number of hits or more
48  * for the given expectation value for the number of hits.
49  *
50  * \param expval expectation value
51  * \param numberOfHits number of hits
52  * \param precision precision
53  * \return propability
54  */
55  inline double getP(const double expval, const unsigned int numberOfHits, const double precision = 1.0e-5)
56  {
57  double P = 0.0;
58 
59  for (int i = numberOfHits; ; ++i) {
60 
61  const double p = TMath::PoissonI(i, expval);
62 
63  P += p;
64 
65  if ((double) i > expval && p < precision) {
66  break;
67  }
68  }
69 
70  return P;
71  }
72 
73 
74  /**
75  * Get chi2 corresponding to given probability.
76  *
77  * \param P probability
78  * \return chi2
79  */
80  inline double getChi2(const double P)
81  {
82  return -log(P);
83  }
84 
85 
86  /**
87  * Get chi2 to observe a hit or not for given expectation value for the number of hits.
88  *
89  * \param expval expectation value
90  * \param hit hit
91  * \return probability
92  */
93  inline double getChi2(const double expval, bool hit)
94  {
95  if (hit)
96  return -log(1.0 - exp(-expval));
97  else
98  return expval;
99  }
100 
101 
102  /**
103  * Determine chi2 of a hit for a given model and time resolution.
104  *
105  * The template argument <tt>JModel_t</tt> refers to a data structure
106  * which should provide for the following method:
107  * <pre>
108  * double %getT(const JHit_t& hit) const; // get expected time of hit
109  * </pre>
110  *
111  * \param model model
112  * \param hit hit
113  * \param sigma sigma
114  * \return chi2
115  */
116  template<class JModel_t, class JHit_t>
117  inline double getChi2(const JModel_t& model, const JHit_t& hit, const double sigma)
118  {
119  const double ds = (hit.getT() - model.getT(hit)) / sigma;
120 
121  return ds * ds;
122  }
123 
124 
125  /**
126  * Determine chi2 of data for given model and time resolution.
127  *
128  * The template argument <tt>JModel_t</tt> refers to a data structure
129  * which should provide for the following method:
130  * <pre>
131  * double %getT(const value_type& hit) const; // expected time of hit
132  * </pre>
133  * where <tt>value_type</tt> corresponds to the value type if the input data.
134  *
135  * \param model model
136  * \param __begin begin of data
137  * \param __end end of data
138  * \param sigma time resolution [ns]
139  * \return chi2
140  */
141  template<class JModel_t, class T>
142  inline double getChi2(const JModel_t& model, T __begin, T __end, const double sigma)
143  {
144  double chi2 = 0.0;
145 
146  for (T hit = __begin; hit != __end; ++hit) {
147  chi2 += getChi2(model, *hit, sigma);
148  }
149 
150  return chi2;
151  }
152 
153 
154  /**
155  * Determine chi2 of data using full covariance matrix.
156  *
157  * \param __begin begin of residuals
158  * \param __end end of residuals
159  * \param V covariance matrix
160  * \return chi2
161  */
162  template<class T>
163  inline double getChi2(T __begin, T __end, const JMatrixNZ& V)
164  {
165  return JMATH::getDot(__begin, __end, V);
166  }
167 
168 
169  /**
170  * Determine chi2 using full covariance matrix.
171  *
172  * \param Y residuals
173  * \param V covariance matrix
174  * \return chi2
175  */
176  inline double getChi2(const JVectorNZ& Y,
177  const JMatrixNZ& V)
178  {
179  return getChi2(Y.begin(), Y.end(), V);
180  }
181 
182 
183  /**
184  * Determine chi2 of data for given track using full covariance matrix.
185  *
186  * \param track track
187  * \param __begin begin of data
188  * \param __end end of data
189  * \param V covariance matrix
190  * \return chi2
191  */
192  template<class T>
193  inline double getChi2(const JLine1Z& track,
194  T __begin,
195  T __end,
196  const JMatrixNZ& V)
197  {
198  const JVectorNZ Y(track, __begin, __end);
199 
200  return getChi2(Y, V);
201  }
202 
203 
204  /**
205  * Determine chi2 of data for given track and angular and time resolution.
206  *
207  * \param track track
208  * \param __begin begin of data
209  * \param __end end of data
210  * \param alpha angular resolution [deg]
211  * \param sigma time resolution [ns]
212  * \return chi2
213  */
214  template<class T>
215  inline double getChi2(const JLine1Z& track, T __begin, T __end, const double alpha, const double sigma)
216  {
217  JMatrixNZ V(track, __begin, __end, alpha, sigma);
218 
219  V.invert();
220 
221  return getChi2(track, __begin, __end, V);
222  }
223 
224 
225  /**
226  * Determine difference between chi2 with and without hit using full covariance matrix.
227  *
228  * \param __begin begin of residuals
229  * \param __end end of residuals
230  * \param V covariance matrix
231  * \param i index of excluded hit
232  * \return difference between chi2 w/o hit i
233  */
234  template<class T>
235  inline double getChi2(T __begin, T __end, const JMatrixNZ& V, const int i)
236  {
237  double chi2 = 0.0;
238 
239  T y = __begin;
240 
241  for (unsigned int j = 0; j != V.size(); ++j, ++y) {
242  chi2 += V(i,j) * (*y);
243  }
244 
245  return chi2*chi2 / V(i,i);
246  }
247 
248 
249  /**
250  * Determine difference between chi2 with and without hit using full covariance matrix.
251  *
252  * \param Y residuals
253  * \param V covariance matrix
254  * \param i index of excluded hit
255  * \return chi2
256  */
257  inline double getChi2(const JVectorNZ& Y,
258  const JMatrixNZ& V,
259  const int i)
260  {
261  return getChi2(Y.begin(), Y.end(), V, i);
262  }
263 
264 
265  /**
266  * Get chi2 of data for given model and fit function.
267  *
268  * The template argument <tt>JFit_t</tt> refers to a data structure
269  * which should provide for the function object operator:
270  * <pre>
271  * double operator()(const JModel_t& model, const value_type&) const; // chi2
272  * </pre>
273  * where
274  * <tt>JModel_t</tt> refers to the given model and
275  * <tt>value_type</tt> to the value type if the input data.
276  * The return value should correspond to the chi2 of the hit.
277  *
278  * \param model model
279  * \param fit fit function
280  * \param __begin begin of data
281  * \param __end end of data
282  */
283  template<class JModel_t, class JFit_t, class T>
284  inline double getChi2(const JModel_t& model, const JFit_t& fit, T __begin, T __end)
285  {
286  double chi2 = 0.0;
287 
288  for (T hit = __begin; hit != __end; ++hit) {
289  chi2 += fit(model, *hit);
290  }
291 
292  return chi2;
293  }
294 
295 
296 }
297 
298 #endif
JVectorNZ.hh
JAANET::JHit_t
Auxiliary class to set-up Hit.
Definition: JHit_t.hh:25
JFIT
Auxiliary classes and methods for linear and iterative data regression.
Definition: JEnergy.hh:15
JTOOLS::j
int j
Definition: JPolint.hh:634
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JEvt.hh
JFIT::getChi2
double getChi2(const JModel_t &model, const JFit_t &fit, T __begin, T __end)
Get chi2 of data for given model and fit function.
Definition: JFitToolkit.hh:284
JMatrixNZ.hh
JFIT::JLine1Z
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
JMath.hh
JMATH::JMatrixNS::invert
void invert()
Invert matrix.
Definition: JMatrixNS.hh:61
JFIT::JVectorNZ
Determination of the time residual vector of hits for a track along z-axis (JFIT::JLine1Z).
Definition: JVectorNZ.hh:25
JFIT::JMatrixNZ
Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z).
Definition: JMatrixNZ.hh:28
JFIT::getP
double getP(const double expval, const unsigned int numberOfHits, const double precision=1.0e-5)
Get Poisson probability to observe given number of hits or more for the given expectation value for t...
Definition: JFitToolkit.hh:55
JMATH::getDot
double getDot(const JFirst_t &first, const JSecond_t &second)
Get dot product of objects.
Definition: JMathToolkit.hh:131