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