1 | /* ======================================================================== *\
|
---|
2 | ! $Name: not supported by cvs2svn $:$Id: MMath.cc,v 1.39 2008-06-18 20:39:48 tbretz Exp $
|
---|
3 | ! --------------------------------------------------------------------------
|
---|
4 | !
|
---|
5 | ! *
|
---|
6 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
|
---|
7 | ! * Software. It is distributed to you in the hope that it can be a useful
|
---|
8 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
|
---|
9 | ! * It is distributed WITHOUT ANY WARRANTY.
|
---|
10 | ! *
|
---|
11 | ! * Permission to use, copy, modify and distribute this software and its
|
---|
12 | ! * documentation for any purpose is hereby granted without fee,
|
---|
13 | ! * provided that the above copyright notice appear in all copies and
|
---|
14 | ! * that both that copyright notice and this permission notice appear
|
---|
15 | ! * in supporting documentation. It is provided "as is" without express
|
---|
16 | ! * or implied warranty.
|
---|
17 | ! *
|
---|
18 | !
|
---|
19 | !
|
---|
20 | ! Author(s): Thomas Bretz 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
|
---|
21 | !
|
---|
22 | ! Copyright: MAGIC Software Development, 2000-2005
|
---|
23 | !
|
---|
24 | !
|
---|
25 | \* ======================================================================== */
|
---|
26 |
|
---|
27 | /////////////////////////////////////////////////////////////////////////////
|
---|
28 | //
|
---|
29 | // MMath
|
---|
30 | //
|
---|
31 | // Mars - Math package (eg Significances, etc)
|
---|
32 | //
|
---|
33 | /////////////////////////////////////////////////////////////////////////////
|
---|
34 | #include "MMath.h"
|
---|
35 |
|
---|
36 | #ifndef ROOT_TVector2
|
---|
37 | #include <TVector2.h>
|
---|
38 | #endif
|
---|
39 |
|
---|
40 | #ifndef ROOT_TVector3
|
---|
41 | #include <TVector3.h>
|
---|
42 | #endif
|
---|
43 |
|
---|
44 | #ifndef ROOT_TArrayD
|
---|
45 | #include <TArrayD.h>
|
---|
46 | #endif
|
---|
47 |
|
---|
48 | #ifndef ROOT_TComplex
|
---|
49 | #include <TComplex.h>
|
---|
50 | #endif
|
---|
51 |
|
---|
52 | //NamespaceImp(MMath);
|
---|
53 |
|
---|
54 | // --------------------------------------------------------------------------
|
---|
55 | //
|
---|
56 | // Calculate Significance as
|
---|
57 | // significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
|
---|
58 | //
|
---|
59 | // s: total number of events in signal region
|
---|
60 | // b: number of background events in signal region
|
---|
61 | //
|
---|
62 | Double_t MMath::Significance(Double_t s, Double_t b)
|
---|
63 | {
|
---|
64 | const Double_t k = b==0 ? 0 : s/b;
|
---|
65 | const Double_t f = s+k*k*b;
|
---|
66 |
|
---|
67 | return f==0 ? 0 : (s-b)/TMath::Sqrt(f);
|
---|
68 | }
|
---|
69 |
|
---|
70 | // --------------------------------------------------------------------------
|
---|
71 | //
|
---|
72 | // Symmetrized significance - this is somehow analog to
|
---|
73 | // SignificanceLiMaSigned
|
---|
74 | //
|
---|
75 | // Returns Significance(s,b) if s>b otherwise -Significance(b, s);
|
---|
76 | //
|
---|
77 | Double_t MMath::SignificanceSym(Double_t s, Double_t b)
|
---|
78 | {
|
---|
79 | return s>b ? Significance(s, b) : -Significance(b, s);
|
---|
80 | }
|
---|
81 |
|
---|
82 | // --------------------------------------------------------------------------
|
---|
83 | //
|
---|
84 | // calculates the significance according to Li & Ma
|
---|
85 | // ApJ 272 (1983) 317, Formula 17
|
---|
86 | //
|
---|
87 | // s // s: number of on events
|
---|
88 | // b // b: number of off events
|
---|
89 | // alpha = t_on/t_off; // t: observation time
|
---|
90 | //
|
---|
91 | // The significance has the same (positive!) value for s>b and b>s.
|
---|
92 | //
|
---|
93 | // Returns -1 if s<0 or b<0 or alpha<0 or the argument of sqrt<0
|
---|
94 | //
|
---|
95 | // Here is some eMail written by Daniel Mazin about the meaning of the arguments:
|
---|
96 | //
|
---|
97 | // > Ok. Here is my understanding:
|
---|
98 | // > According to Li&Ma paper (correctly cited in MMath.cc) alpha is the
|
---|
99 | // > scaling factor. The mathematics behind the formula 17 (and/or 9) implies
|
---|
100 | // > exactly this. If you scale OFF to ON first (using time or using any other
|
---|
101 | // > method), then you cannot use formula 17 (9) anymore. You can just try
|
---|
102 | // > the formula before scaling (alpha!=1) and after scaling (alpha=1), you
|
---|
103 | // > will see the result will be different.
|
---|
104 | //
|
---|
105 | // > Here are less mathematical arguments:
|
---|
106 | //
|
---|
107 | // > 1) the better background determination you have (smaller alpha) the more
|
---|
108 | // > significant is your excess, thus your analysis is more sensitive. If you
|
---|
109 | // > normalize OFF to ON first, you loose this sensitivity.
|
---|
110 | //
|
---|
111 | // > 2) the normalization OFF to ON has an error, which naturally depends on
|
---|
112 | // > the OFF and ON. This error is propagating to the significance of your
|
---|
113 | // > excess if you use the Li&Ma formula 17 correctly. But if you normalize
|
---|
114 | // > first and use then alpha=1, the error gets lost completely, you loose
|
---|
115 | // > somehow the criteria of goodness of the normalization.
|
---|
116 | //
|
---|
117 | Double_t MMath::SignificanceLiMa(Double_t s, Double_t b, Double_t alpha)
|
---|
118 | {
|
---|
119 | const Double_t sum = s+b;
|
---|
120 |
|
---|
121 | if (s<0 || b<0 || alpha<=0)
|
---|
122 | return -1;
|
---|
123 |
|
---|
124 | const Double_t l = s==0 ? 0 : s*TMath::Log(s/sum*(alpha+1)/alpha);
|
---|
125 | const Double_t m = b==0 ? 0 : b*TMath::Log(b/sum*(alpha+1) );
|
---|
126 |
|
---|
127 | return l+m<0 ? -1 : TMath::Sqrt((l+m)*2);
|
---|
128 | }
|
---|
129 |
|
---|
130 | /*
|
---|
131 | Double_t MMath::SignificanceLiMaErr(Double_t s, Double_t b, Double_t alpha)
|
---|
132 | {
|
---|
133 | Double_t S = SignificanceLiMa(s, b, alpha);
|
---|
134 | if (S<0)
|
---|
135 | return -1;
|
---|
136 |
|
---|
137 | const Double_t sum = s+b;
|
---|
138 |
|
---|
139 |
|
---|
140 | Double_t l = TMath::Log(s/sum*(alpha+1)/alpha)/TMath::Sqrt(2*S);
|
---|
141 | Double_t m = TMath::Log(s/sum*(alpha+1)/alpha)/TMath::Sqrt(2*S);
|
---|
142 |
|
---|
143 |
|
---|
144 | const Double_t sum = s+b;
|
---|
145 |
|
---|
146 | if (s<0 || b<0 || alpha<=0)
|
---|
147 | return -1;
|
---|
148 |
|
---|
149 | const Double_t l = s==0 ? 0 : s*TMath::Log(s/sum*(alpha+1)/alpha);
|
---|
150 | const Double_t m = b==0 ? 0 : b*TMath::Log(b/sum*(alpha+1) );
|
---|
151 |
|
---|
152 | return l+m<0 ? -1 : TMath::Sqrt((l+m)*2);
|
---|
153 | }
|
---|
154 | */
|
---|
155 |
|
---|
156 | // --------------------------------------------------------------------------
|
---|
157 | //
|
---|
158 | // Calculates MMath::SignificanceLiMa(s, b, alpha). Returns 0 if the
|
---|
159 | // calculation has failed. Otherwise the Li/Ma significance which was
|
---|
160 | // calculated. If s<b a negative value is returned.
|
---|
161 | //
|
---|
162 | Double_t MMath::SignificanceLiMaSigned(Double_t s, Double_t b, Double_t alpha)
|
---|
163 | {
|
---|
164 | const Double_t sig = SignificanceLiMa(s, b, alpha);
|
---|
165 | if (sig<=0)
|
---|
166 | return 0;
|
---|
167 |
|
---|
168 | return TMath::Sign(sig, s-alpha*b);
|
---|
169 | }
|
---|
170 |
|
---|
171 | // --------------------------------------------------------------------------
|
---|
172 | //
|
---|
173 | // Return Li/Ma (5) for the error of the excess, under the assumption that
|
---|
174 | // the existance of a signal is already known.
|
---|
175 | //
|
---|
176 | Double_t MMath::SignificanceLiMaExc(Double_t s, Double_t b, Double_t alpha)
|
---|
177 | {
|
---|
178 | Double_t Ns = s - alpha*b;
|
---|
179 | Double_t sN = s + alpha*alpha*b;
|
---|
180 |
|
---|
181 | if (Ns<0 || sN<0)
|
---|
182 | return 0;
|
---|
183 |
|
---|
184 | if (Ns==0 && sN==0)
|
---|
185 | return 0;
|
---|
186 |
|
---|
187 | return Ns/TMath::Sqrt(sN);
|
---|
188 | }
|
---|
189 |
|
---|
190 | // --------------------------------------------------------------------------
|
---|
191 | //
|
---|
192 | // Returns: 2/(sigma*sqrt(2))*integral[0,x](exp(-(x-mu)^2/(2*sigma^2)))
|
---|
193 | //
|
---|
194 | Double_t MMath::GaussProb(Double_t x, Double_t sigma, Double_t mean)
|
---|
195 | {
|
---|
196 | if (x<mean)
|
---|
197 | return 0;
|
---|
198 |
|
---|
199 | static const Double_t sqrt2 = TMath::Sqrt(2.);
|
---|
200 |
|
---|
201 | const Double_t rc = TMath::Erf((x-mean)/(sigma*sqrt2));
|
---|
202 |
|
---|
203 | if (rc<0)
|
---|
204 | return 0;
|
---|
205 | if (rc>1)
|
---|
206 | return 1;
|
---|
207 |
|
---|
208 | return rc;
|
---|
209 | }
|
---|
210 |
|
---|
211 | // ------------------------------------------------------------------------
|
---|
212 | //
|
---|
213 | // Return the "median" (at 68.3%) value of the distribution of
|
---|
214 | // abs(a[i]-Median)
|
---|
215 | //
|
---|
216 | template <class Size, class Element>
|
---|
217 | Double_t MMath::MedianDevImp(Size n, const Element *a, Double_t &med)
|
---|
218 | {
|
---|
219 | static const Double_t prob = 0.682689477208650697; //MMath::GaussProb(1.0);
|
---|
220 |
|
---|
221 | // Sanity check
|
---|
222 | if (n <= 0 || !a)
|
---|
223 | return 0;
|
---|
224 |
|
---|
225 | // Get median of distribution
|
---|
226 | med = TMath::Median(n, a);
|
---|
227 |
|
---|
228 | // Create the abs(a[i]-med) distribution
|
---|
229 | Double_t arr[n];
|
---|
230 | for (int i=0; i<n; i++)
|
---|
231 | arr[i] = TMath::Abs(a[i]-med);
|
---|
232 |
|
---|
233 | //return TMath::Median(n, arr)/0.67449896936; //MMath::GaussProb(x)=0.5
|
---|
234 |
|
---|
235 | // Define where to divide (floor because the highest possible is n-1)
|
---|
236 | const Int_t div = TMath::FloorNint(n*prob);
|
---|
237 |
|
---|
238 | // Calculate result
|
---|
239 | Double_t dev = TMath::KOrdStat(n, arr, div);
|
---|
240 | if (n%2 == 0)
|
---|
241 | {
|
---|
242 | dev += TMath::KOrdStat(n, arr, div-1);
|
---|
243 | dev /= 2;
|
---|
244 | }
|
---|
245 |
|
---|
246 | return dev;
|
---|
247 | }
|
---|
248 |
|
---|
249 | // ------------------------------------------------------------------------
|
---|
250 | //
|
---|
251 | // Return the "median" (at 68.3%) value of the distribution of
|
---|
252 | // abs(a[i]-Median)
|
---|
253 | //
|
---|
254 | Double_t MMath::MedianDev(Long64_t n, const Short_t *a, Double_t &med)
|
---|
255 | {
|
---|
256 | return MedianDevImp(n, a, med);
|
---|
257 | }
|
---|
258 |
|
---|
259 | // ------------------------------------------------------------------------
|
---|
260 | //
|
---|
261 | // Return the "median" (at 68.3%) value of the distribution of
|
---|
262 | // abs(a[i]-Median)
|
---|
263 | //
|
---|
264 | Double_t MMath::MedianDev(Long64_t n, const Int_t *a, Double_t &med)
|
---|
265 | {
|
---|
266 | return MedianDevImp(n, a, med);
|
---|
267 | }
|
---|
268 |
|
---|
269 | // ------------------------------------------------------------------------
|
---|
270 | //
|
---|
271 | // Return the "median" (at 68.3%) value of the distribution of
|
---|
272 | // abs(a[i]-Median)
|
---|
273 | //
|
---|
274 | Double_t MMath::MedianDev(Long64_t n, const Float_t *a, Double_t &med)
|
---|
275 | {
|
---|
276 | return MedianDevImp(n, a, med);
|
---|
277 | }
|
---|
278 |
|
---|
279 | // ------------------------------------------------------------------------
|
---|
280 | //
|
---|
281 | // Return the "median" (at 68.3%) value of the distribution of
|
---|
282 | // abs(a[i]-Median)
|
---|
283 | //
|
---|
284 | Double_t MMath::MedianDev(Long64_t n, const Double_t *a, Double_t &med)
|
---|
285 | {
|
---|
286 | return MedianDevImp(n, a, med);
|
---|
287 | }
|
---|
288 |
|
---|
289 | // ------------------------------------------------------------------------
|
---|
290 | //
|
---|
291 | // Return the "median" (at 68.3%) value of the distribution of
|
---|
292 | // abs(a[i]-Median)
|
---|
293 | //
|
---|
294 | Double_t MMath::MedianDev(Long64_t n, const Long_t *a, Double_t &med)
|
---|
295 | {
|
---|
296 | return MedianDevImp(n, a, med);
|
---|
297 | }
|
---|
298 |
|
---|
299 | // ------------------------------------------------------------------------
|
---|
300 | //
|
---|
301 | // Return the "median" (at 68.3%) value of the distribution of
|
---|
302 | // abs(a[i]-Median)
|
---|
303 | //
|
---|
304 | Double_t MMath::MedianDev(Long64_t n, const Long64_t *a, Double_t &med)
|
---|
305 | {
|
---|
306 | return MedianDevImp(n, a, med);
|
---|
307 | }
|
---|
308 |
|
---|
309 | Double_t MMath::MedianDev(Long64_t n, const Short_t *a) { Double_t med; return MedianDevImp(n, a, med); }
|
---|
310 | Double_t MMath::MedianDev(Long64_t n, const Int_t *a) { Double_t med; return MedianDevImp(n, a, med); }
|
---|
311 | Double_t MMath::MedianDev(Long64_t n, const Float_t *a) { Double_t med; return MedianDevImp(n, a, med); }
|
---|
312 | Double_t MMath::MedianDev(Long64_t n, const Double_t *a) { Double_t med; return MedianDevImp(n, a, med); }
|
---|
313 | Double_t MMath::MedianDev(Long64_t n, const Long_t *a) { Double_t med; return MedianDevImp(n, a, med); }
|
---|
314 | Double_t MMath::MedianDev(Long64_t n, const Long64_t *a) { Double_t med; return MedianDevImp(n, a, med); }
|
---|
315 |
|
---|
316 | // --------------------------------------------------------------------------
|
---|
317 | //
|
---|
318 | // This function reduces the precision to roughly 0.5% of a Float_t by
|
---|
319 | // changing its bit-pattern (Be carefull, in rare cases this function must
|
---|
320 | // be adapted to different machines!). This is usefull to enforce better
|
---|
321 | // compression by eg. gzip.
|
---|
322 | //
|
---|
323 | void MMath::ReducePrecision(Float_t &val)
|
---|
324 | {
|
---|
325 | UInt_t &f = (UInt_t&)val;
|
---|
326 |
|
---|
327 | f += 0x00004000;
|
---|
328 | f &= 0xffff8000;
|
---|
329 | }
|
---|
330 |
|
---|
331 | // -------------------------------------------------------------------------
|
---|
332 | //
|
---|
333 | // Quadratic interpolation
|
---|
334 | //
|
---|
335 | // calculate the parameters of a parabula such that
|
---|
336 | // y(i) = a + b*x(i) + c*x(i)^2
|
---|
337 | //
|
---|
338 | // If the determinant==0 an empty TVector3 is returned.
|
---|
339 | //
|
---|
340 | TVector3 MMath::GetParab(const TVector3 &x, const TVector3 &y)
|
---|
341 | {
|
---|
342 | const Double_t x1 = x(0);
|
---|
343 | const Double_t x2 = x(1);
|
---|
344 | const Double_t x3 = x(2);
|
---|
345 |
|
---|
346 | const Double_t y1 = y(0);
|
---|
347 | const Double_t y2 = y(1);
|
---|
348 | const Double_t y3 = y(2);
|
---|
349 |
|
---|
350 | const double det =
|
---|
351 | + x2*x3*x3 + x1*x2*x2 + x3*x1*x1
|
---|
352 | - x2*x1*x1 - x3*x2*x2 - x1*x3*x3;
|
---|
353 |
|
---|
354 |
|
---|
355 | if (det==0)
|
---|
356 | return TVector3();
|
---|
357 |
|
---|
358 | const double det1 = 1.0/det;
|
---|
359 |
|
---|
360 | const double ai11 = x2*x3*x3 - x3*x2*x2;
|
---|
361 | const double ai12 = x3*x1*x1 - x1*x3*x3;
|
---|
362 | const double ai13 = x1*x2*x2 - x2*x1*x1;
|
---|
363 |
|
---|
364 | const double ai21 = x2*x2 - x3*x3;
|
---|
365 | const double ai22 = x3*x3 - x1*x1;
|
---|
366 | const double ai23 = x1*x1 - x2*x2;
|
---|
367 |
|
---|
368 | const double ai31 = x3 - x2;
|
---|
369 | const double ai32 = x1 - x3;
|
---|
370 | const double ai33 = x2 - x1;
|
---|
371 |
|
---|
372 | return TVector3((ai11*y1 + ai12*y2 + ai13*y3) * det1,
|
---|
373 | (ai21*y1 + ai22*y2 + ai23*y3) * det1,
|
---|
374 | (ai31*y1 + ai32*y2 + ai33*y3) * det1);
|
---|
375 | }
|
---|
376 |
|
---|
377 | // --------------------------------------------------------------------------
|
---|
378 | //
|
---|
379 | // Interpolate the points with x-coordinates vx and y-coordinates vy
|
---|
380 | // by a parabola (second order polynomial) and return the value at x.
|
---|
381 | //
|
---|
382 | Double_t MMath::InterpolParabLin(const TVector3 &vx, const TVector3 &vy, Double_t x)
|
---|
383 | {
|
---|
384 | const TVector3 c = GetParab(vx, vy);
|
---|
385 | return c(0) + c(1)*x + c(2)*x*x;
|
---|
386 | }
|
---|
387 |
|
---|
388 | // --------------------------------------------------------------------------
|
---|
389 | //
|
---|
390 | // Interpolate the points with x-coordinates vx=(-1,0,1) and
|
---|
391 | // y-coordinates vy by a parabola (second order polynomial) and return
|
---|
392 | // the value at x.
|
---|
393 | //
|
---|
394 | Double_t MMath::InterpolParabLin(const TVector3 &vy, Double_t x)
|
---|
395 | {
|
---|
396 | const TVector3 c(vy(1), (vy(2)-vy(0))/2, vy(0)/2 - vy(1) + vy(2)/2);
|
---|
397 | return c(0) + c(1)*x + c(2)*x*x;
|
---|
398 | }
|
---|
399 |
|
---|
400 | Double_t MMath::InterpolParabLog(const TVector3 &vx, const TVector3 &vy, Double_t x)
|
---|
401 | {
|
---|
402 | const Double_t l0 = TMath::Log10(vx(0));
|
---|
403 | const Double_t l1 = TMath::Log10(vx(1));
|
---|
404 | const Double_t l2 = TMath::Log10(vx(2));
|
---|
405 |
|
---|
406 | const TVector3 vx0(l0, l1, l2);
|
---|
407 | return InterpolParabLin(vx0, vy, TMath::Log10(x));
|
---|
408 | }
|
---|
409 |
|
---|
410 | Double_t MMath::InterpolParabCos(const TVector3 &vx, const TVector3 &vy, Double_t x)
|
---|
411 | {
|
---|
412 | const Double_t l0 = TMath::Cos(vx(0));
|
---|
413 | const Double_t l1 = TMath::Cos(vx(1));
|
---|
414 | const Double_t l2 = TMath::Cos(vx(2));
|
---|
415 |
|
---|
416 | const TVector3 vx0(l0, l1, l2);
|
---|
417 | return InterpolParabLin(vx0, vy, TMath::Cos(x));
|
---|
418 | }
|
---|
419 |
|
---|
420 | // --------------------------------------------------------------------------
|
---|
421 | //
|
---|
422 | // Analytically calculated result of a least square fit of:
|
---|
423 | // y = A*e^(B*x)
|
---|
424 | // Equal weights
|
---|
425 | //
|
---|
426 | // It returns TArrayD(2) = { A, B };
|
---|
427 | //
|
---|
428 | // see: http://mathworld.wolfram.com/LeastSquaresFittingExponential.html
|
---|
429 | //
|
---|
430 | TArrayD MMath::LeastSqFitExpW1(Int_t n, Double_t *x, Double_t *y)
|
---|
431 | {
|
---|
432 | Double_t sumxsqy = 0;
|
---|
433 | Double_t sumylny = 0;
|
---|
434 | Double_t sumxy = 0;
|
---|
435 | Double_t sumy = 0;
|
---|
436 | Double_t sumxylny = 0;
|
---|
437 | for (int i=0; i<n; i++)
|
---|
438 | {
|
---|
439 | sumylny += y[i]*TMath::Log(y[i]);
|
---|
440 | sumxy += x[i]*y[i];
|
---|
441 | sumxsqy += x[i]*x[i]*y[i];
|
---|
442 | sumxylny += x[i]*y[i]*TMath::Log(y[i]);
|
---|
443 | sumy += y[i];
|
---|
444 | }
|
---|
445 |
|
---|
446 | const Double_t dev = sumy*sumxsqy - sumxy*sumxy;
|
---|
447 |
|
---|
448 | const Double_t a = (sumxsqy*sumylny - sumxy*sumxylny)/dev;
|
---|
449 | const Double_t b = (sumy*sumxylny - sumxy*sumylny)/dev;
|
---|
450 |
|
---|
451 | TArrayD rc(2);
|
---|
452 | rc[0] = TMath::Exp(a);
|
---|
453 | rc[1] = b;
|
---|
454 | return rc;
|
---|
455 | }
|
---|
456 |
|
---|
457 | // --------------------------------------------------------------------------
|
---|
458 | //
|
---|
459 | // Analytically calculated result of a least square fit of:
|
---|
460 | // y = A*e^(B*x)
|
---|
461 | // Greater weights to smaller values
|
---|
462 | //
|
---|
463 | // It returns TArrayD(2) = { A, B };
|
---|
464 | //
|
---|
465 | // see: http://mathworld.wolfram.com/LeastSquaresFittingExponential.html
|
---|
466 | //
|
---|
467 | TArrayD MMath::LeastSqFitExp(Int_t n, Double_t *x, Double_t *y)
|
---|
468 | {
|
---|
469 | // -------- Greater weights to smaller values ---------
|
---|
470 | Double_t sumlny = 0;
|
---|
471 | Double_t sumxlny = 0;
|
---|
472 | Double_t sumxsq = 0;
|
---|
473 | Double_t sumx = 0;
|
---|
474 | for (int i=0; i<n; i++)
|
---|
475 | {
|
---|
476 | sumlny += TMath::Log(y[i]);
|
---|
477 | sumxlny += x[i]*TMath::Log(y[i]);
|
---|
478 |
|
---|
479 | sumxsq += x[i]*x[i];
|
---|
480 | sumx += x[i];
|
---|
481 | }
|
---|
482 |
|
---|
483 | const Double_t dev = n*sumxsq-sumx*sumx;
|
---|
484 |
|
---|
485 | const Double_t a = (sumlny*sumxsq - sumx*sumxlny)/dev;
|
---|
486 | const Double_t b = (n*sumxlny - sumx*sumlny)/dev;
|
---|
487 |
|
---|
488 | TArrayD rc(2);
|
---|
489 | rc[0] = TMath::Exp(a);
|
---|
490 | rc[1] = b;
|
---|
491 | return rc;
|
---|
492 | }
|
---|
493 |
|
---|
494 | // --------------------------------------------------------------------------
|
---|
495 | //
|
---|
496 | // Analytically calculated result of a least square fit of:
|
---|
497 | // y = A+B*ln(x)
|
---|
498 | //
|
---|
499 | // It returns TArrayD(2) = { A, B };
|
---|
500 | //
|
---|
501 | // see: http://mathworld.wolfram.com/LeastSquaresFittingLogarithmic.html
|
---|
502 | //
|
---|
503 | TArrayD MMath::LeastSqFitLog(Int_t n, Double_t *x, Double_t *y)
|
---|
504 | {
|
---|
505 | Double_t sumylnx = 0;
|
---|
506 | Double_t sumy = 0;
|
---|
507 | Double_t sumlnx = 0;
|
---|
508 | Double_t sumlnxsq = 0;
|
---|
509 | for (int i=0; i<n; i++)
|
---|
510 | {
|
---|
511 | sumylnx += y[i]*TMath::Log(x[i]);
|
---|
512 | sumy += y[i];
|
---|
513 | sumlnx += TMath::Log(x[i]);
|
---|
514 | sumlnxsq += TMath::Log(x[i])*TMath::Log(x[i]);
|
---|
515 | }
|
---|
516 |
|
---|
517 | const Double_t b = (n*sumylnx-sumy*sumlnx)/(n*sumlnxsq-sumlnx*sumlnx);
|
---|
518 | const Double_t a = (sumy-b*sumlnx)/n;
|
---|
519 |
|
---|
520 | TArrayD rc(2);
|
---|
521 | rc[0] = a;
|
---|
522 | rc[1] = b;
|
---|
523 | return rc;
|
---|
524 | }
|
---|
525 |
|
---|
526 | // --------------------------------------------------------------------------
|
---|
527 | //
|
---|
528 | // Analytically calculated result of a least square fit of:
|
---|
529 | // y = A*x^B
|
---|
530 | //
|
---|
531 | // It returns TArrayD(2) = { A, B };
|
---|
532 | //
|
---|
533 | // see: http://mathworld.wolfram.com/LeastSquaresFittingPowerLaw.html
|
---|
534 | //
|
---|
535 | TArrayD MMath::LeastSqFitPowerLaw(Int_t n, Double_t *x, Double_t *y)
|
---|
536 | {
|
---|
537 | Double_t sumlnxlny = 0;
|
---|
538 | Double_t sumlnx = 0;
|
---|
539 | Double_t sumlny = 0;
|
---|
540 | Double_t sumlnxsq = 0;
|
---|
541 | for (int i=0; i<n; i++)
|
---|
542 | {
|
---|
543 | sumlnxlny += TMath::Log(x[i])*TMath::Log(y[i]);
|
---|
544 | sumlnx += TMath::Log(x[i]);
|
---|
545 | sumlny += TMath::Log(y[i]);
|
---|
546 | sumlnxsq += TMath::Log(x[i])*TMath::Log(x[i]);
|
---|
547 | }
|
---|
548 |
|
---|
549 | const Double_t b = (n*sumlnxlny-sumlnx*sumlny)/(n*sumlnxsq-sumlnx*sumlnx);
|
---|
550 | const Double_t a = (sumlny-b*sumlnx)/n;
|
---|
551 |
|
---|
552 | TArrayD rc(2);
|
---|
553 | rc[0] = TMath::Exp(a);
|
---|
554 | rc[1] = b;
|
---|
555 | return rc;
|
---|
556 | }
|
---|
557 |
|
---|
558 | // --------------------------------------------------------------------------
|
---|
559 | //
|
---|
560 | // Calculate the intersection of two lines defined by (x1;y1) and (x2;x2)
|
---|
561 | // Returns the intersection point.
|
---|
562 | //
|
---|
563 | // It is assumed that the lines intersect. If there is no intersection
|
---|
564 | // TVector2() is returned (which is not destinguishable from
|
---|
565 | // TVector2(0,0) if the intersection is at the coordinate source)
|
---|
566 | //
|
---|
567 | // Formula from: http://mathworld.wolfram.com/Line-LineIntersection.html
|
---|
568 | //
|
---|
569 | TVector2 MMath::GetIntersectionPoint(const TVector2 &x1, const TVector2 &y1, const TVector2 &x2, const TVector2 &y2)
|
---|
570 | {
|
---|
571 | TMatrix d(2,2);
|
---|
572 | d[0][0] = x1.X()-y1.X();
|
---|
573 | d[0][1] = x2.X()-y2.X();
|
---|
574 | d[1][0] = x1.Y()-y1.Y();
|
---|
575 | d[1][1] = x2.Y()-y2.Y();
|
---|
576 |
|
---|
577 | const Double_t denom = d.Determinant();
|
---|
578 | if (denom==0)
|
---|
579 | return TVector2();
|
---|
580 |
|
---|
581 | TMatrix l1(2,2);
|
---|
582 | TMatrix l2(2,2);
|
---|
583 |
|
---|
584 | l1[0][0] = x1.X();
|
---|
585 | l1[0][1] = y1.X();
|
---|
586 | l2[0][0] = x2.X();
|
---|
587 | l2[0][1] = y2.X();
|
---|
588 |
|
---|
589 | l1[1][0] = x1.Y();
|
---|
590 | l1[1][1] = y1.Y();
|
---|
591 | l2[1][0] = x2.Y();
|
---|
592 | l2[1][1] = y2.Y();
|
---|
593 |
|
---|
594 | TMatrix a(2,2);
|
---|
595 | a[0][0] = l1.Determinant();
|
---|
596 | a[0][1] = l2.Determinant();
|
---|
597 | a[1][0] = x1.X()-y1.X();
|
---|
598 | a[1][1] = x2.X()-y2.X();
|
---|
599 |
|
---|
600 | const Double_t X = a.Determinant()/denom;
|
---|
601 |
|
---|
602 | a[1][0] = x1.Y()-y1.Y();
|
---|
603 | a[1][1] = x2.Y()-y2.Y();
|
---|
604 |
|
---|
605 | const Double_t Y = a.Determinant()/denom;
|
---|
606 |
|
---|
607 | return TVector2(X, Y);
|
---|
608 | }
|
---|
609 |
|
---|
610 | // --------------------------------------------------------------------------
|
---|
611 | //
|
---|
612 | // Solves: x^2 + ax + b = 0;
|
---|
613 | // Return number of solutions returned as x1, x2
|
---|
614 | //
|
---|
615 | Int_t MMath::SolvePol2(Double_t a, Double_t b, Double_t &x1, Double_t &x2)
|
---|
616 | {
|
---|
617 | const Double_t r = a*a - 4*b;
|
---|
618 | if (r<0)
|
---|
619 | return 0;
|
---|
620 |
|
---|
621 | if (r==0)
|
---|
622 | {
|
---|
623 | x1 = x2 = -a/2;
|
---|
624 | return 1;
|
---|
625 | }
|
---|
626 |
|
---|
627 | const Double_t s = TMath::Sqrt(r);
|
---|
628 |
|
---|
629 | x1 = (-a+s)/2;
|
---|
630 | x2 = (-a-s)/2;
|
---|
631 |
|
---|
632 | return 2;
|
---|
633 | }
|
---|
634 |
|
---|
635 | // --------------------------------------------------------------------------
|
---|
636 | //
|
---|
637 | // This is a helper function making the execution of SolverPol3 a bit faster
|
---|
638 | //
|
---|
639 | static inline Double_t ReMul(const TComplex &c1, const TComplex &th)
|
---|
640 | {
|
---|
641 | const TComplex c2 = TComplex::Cos(th/3.);
|
---|
642 | return c1.Re() * c2.Re() - c1.Im() * c2.Im();
|
---|
643 | }
|
---|
644 |
|
---|
645 | // --------------------------------------------------------------------------
|
---|
646 | //
|
---|
647 | // Solves: x^3 + ax^2 + bx + c = 0;
|
---|
648 | // Return number of the real solutions, returned as z1, z2, z3
|
---|
649 | //
|
---|
650 | // Algorithm adapted from http://home.att.net/~srschmitt/cubizen.heml
|
---|
651 | // Which is based on the solution given in
|
---|
652 | // http://mathworld.wolfram.com/CubicEquation.html
|
---|
653 | //
|
---|
654 | // -------------------------------------------------------------------------
|
---|
655 | //
|
---|
656 | // Exact solutions of cubic polynomial equations
|
---|
657 | // by Stephen R. Schmitt Algorithm
|
---|
658 | //
|
---|
659 | // An exact solution of the cubic polynomial equation:
|
---|
660 | //
|
---|
661 | // x^3 + a*x^2 + b*x + c = 0
|
---|
662 | //
|
---|
663 | // was first published by Gerolamo Cardano (1501-1576) in his treatise,
|
---|
664 | // Ars Magna. He did not discoverer of the solution; a professor of
|
---|
665 | // mathematics at the University of Bologna named Scipione del Ferro (ca.
|
---|
666 | // 1465-1526) is credited as the first to find an exact solution. In the
|
---|
667 | // years since, several improvements to the original solution have been
|
---|
668 | // discovered. Zeno source code
|
---|
669 | //
|
---|
670 | // http://home.att.net/~srschmitt/cubizen.html
|
---|
671 | //
|
---|
672 | // % compute real or complex roots of cubic polynomial
|
---|
673 | // function cubic( var z1, z2, z3 : real, a, b, c : real ) : real
|
---|
674 | //
|
---|
675 | // var Q, R, D, S, T : real
|
---|
676 | // var im, th : real
|
---|
677 | //
|
---|
678 | // Q := (3*b - a^2)/9
|
---|
679 | // R := (9*b*a - 27*c - 2*a^3)/54
|
---|
680 | // D := Q^3 + R^2 % polynomial discriminant
|
---|
681 | //
|
---|
682 | // if (D >= 0) then % complex or duplicate roots
|
---|
683 | //
|
---|
684 | // S := sgn(R + sqrt(D))*abs(R + sqrt(D))^(1/3)
|
---|
685 | // T := sgn(R - sqrt(D))*abs(R - sqrt(D))^(1/3)
|
---|
686 | //
|
---|
687 | // z1 := -a/3 + (S + T) % real root
|
---|
688 | // z2 := -a/3 - (S + T)/2 % real part of complex root
|
---|
689 | // z3 := -a/3 - (S + T)/2 % real part of complex root
|
---|
690 | // im := abs(sqrt(3)*(S - T)/2) % complex part of root pair
|
---|
691 | //
|
---|
692 | // else % distinct real roots
|
---|
693 | //
|
---|
694 | // th := arccos(R/sqrt( -Q^3))
|
---|
695 | //
|
---|
696 | // z1 := 2*sqrt(-Q)*cos(th/3) - a/3
|
---|
697 | // z2 := 2*sqrt(-Q)*cos((th + 2*pi)/3) - a/3
|
---|
698 | // z3 := 2*sqrt(-Q)*cos((th + 4*pi)/3) - a/3
|
---|
699 | // im := 0
|
---|
700 | //
|
---|
701 | // end if
|
---|
702 | //
|
---|
703 | // return im % imaginary part
|
---|
704 | //
|
---|
705 | // end function
|
---|
706 | //
|
---|
707 | // see also http://en.wikipedia.org/wiki/Cubic_equation
|
---|
708 | //
|
---|
709 | Int_t MMath::SolvePol3(Double_t a, Double_t b, Double_t c,
|
---|
710 | Double_t &x1, Double_t &x2, Double_t &x3)
|
---|
711 | {
|
---|
712 | // Double_t coeff[4] = { 1, a, b, c };
|
---|
713 | // return TMath::RootsCubic(coeff, x1, x2, x3) ? 1 : 3;
|
---|
714 |
|
---|
715 | const Double_t Q = (a*a - 3*b)/9;
|
---|
716 | const Double_t R = (9*b*a - 27*c - 2*a*a*a)/54;
|
---|
717 | const Double_t D = R*R - Q*Q*Q; // polynomial discriminant
|
---|
718 |
|
---|
719 | // ----- The single-real / duplicate-roots solution -----
|
---|
720 |
|
---|
721 | // D<0: three real roots
|
---|
722 | // D>0: one real root
|
---|
723 | // D==0: maximum two real roots (two identical roots)
|
---|
724 |
|
---|
725 | // R==0: only one unique root
|
---|
726 | // R!=0: two roots
|
---|
727 |
|
---|
728 | if (D==0)
|
---|
729 | {
|
---|
730 | const Double_t r = MMath::Sqrt3(R);
|
---|
731 |
|
---|
732 | x1 = r - a/3.; // real root
|
---|
733 | if (R==0)
|
---|
734 | return 1;
|
---|
735 |
|
---|
736 | x2 = 2*r - a/3.; // real root
|
---|
737 | return 2;
|
---|
738 | }
|
---|
739 |
|
---|
740 | if (D>0) // complex or duplicate roots
|
---|
741 | {
|
---|
742 | const Double_t sqrtd = TMath::Sqrt(D);
|
---|
743 |
|
---|
744 | const Double_t A = TMath::Sign(1., R)*MMath::Sqrt3(TMath::Abs(R)+sqrtd);
|
---|
745 |
|
---|
746 | // The case A==0 cannot happen. This would imply D==0
|
---|
747 | // if (A==0)
|
---|
748 | // {
|
---|
749 | // x1 = -a/3;
|
---|
750 | // return 1;
|
---|
751 | // }
|
---|
752 |
|
---|
753 | x1 = (A+Q/A)-a/3;
|
---|
754 |
|
---|
755 | //const Double_t S = MMath::Sqrt3(R + sqrtd);
|
---|
756 | //const Double_t T = MMath::Sqrt3(R - sqrtd);
|
---|
757 | //x1 = (S+T) - a/3.; // real root
|
---|
758 |
|
---|
759 | return 1;
|
---|
760 |
|
---|
761 | //z2 = (S + T)/2 - a/3.; // real part of complex root
|
---|
762 | //z3 = (S + T)/2 - a/3.; // real part of complex root
|
---|
763 | //im = fabs(sqrt(3)*(S - T)/2) // complex part of root pair
|
---|
764 | }
|
---|
765 |
|
---|
766 | // ----- The general solution with three roots ---
|
---|
767 |
|
---|
768 | if (Q==0)
|
---|
769 | return 0;
|
---|
770 |
|
---|
771 | if (Q>0) // This is here for speed reasons
|
---|
772 | {
|
---|
773 | const Double_t sqrtq = TMath::Sqrt(Q);
|
---|
774 | const Double_t rq = R/TMath::Abs(Q);
|
---|
775 |
|
---|
776 | const Double_t t = TMath::ACos(rq/sqrtq)/3;
|
---|
777 |
|
---|
778 | static const Double_t sqrt3 = TMath::Sqrt(3.);
|
---|
779 |
|
---|
780 | const Double_t sn = TMath::Sin(t)*sqrt3;
|
---|
781 | const Double_t cs = TMath::Cos(t);
|
---|
782 |
|
---|
783 | x1 = 2*sqrtq * cs - a/3;
|
---|
784 | x2 = -sqrtq * (sn + cs) - a/3;
|
---|
785 | x3 = sqrtq * (sn - cs) - a/3;
|
---|
786 |
|
---|
787 | /* --- Easier to understand but slower ---
|
---|
788 | const Double_t th1 = TMath::ACos(rq/sqrtq);
|
---|
789 | const Double_t th2 = th1 + TMath::TwoPi();
|
---|
790 | const Double_t th3 = th2 + TMath::TwoPi();
|
---|
791 |
|
---|
792 | x1 = 2.*sqrtq * TMath::Cos(th1/3.) - a/3.;
|
---|
793 | x2 = 2.*sqrtq * TMath::Cos(th2/3.) - a/3.;
|
---|
794 | x3 = 2.*sqrtq * TMath::Cos(th3/3.) - a/3.;
|
---|
795 | */
|
---|
796 | return 3;
|
---|
797 | }
|
---|
798 |
|
---|
799 | const TComplex sqrtq = TComplex::Sqrt(Q);
|
---|
800 | const Double_t rq = R/TMath::Abs(Q);
|
---|
801 |
|
---|
802 | const TComplex th1 = TComplex::ACos(rq/sqrtq);
|
---|
803 | const TComplex th2 = th1 + TMath::TwoPi();
|
---|
804 | const TComplex th3 = th2 + TMath::TwoPi();
|
---|
805 |
|
---|
806 | // For ReMul, see bove
|
---|
807 | x1 = ReMul(2.*sqrtq, th1) - a/3.;
|
---|
808 | x2 = ReMul(2.*sqrtq, th2) - a/3.;
|
---|
809 | x3 = ReMul(2.*sqrtq, th3) - a/3.;
|
---|
810 |
|
---|
811 | return 3;
|
---|
812 | }
|
---|
813 |
|
---|
814 | // --------------------------------------------------------------------------
|
---|
815 | //
|
---|
816 | // Format a value and its error corresponding to the rules (note
|
---|
817 | // this won't work if the error is more then eight orders smaller than
|
---|
818 | // the value)
|
---|
819 | //
|
---|
820 | void MMath::Format(Double_t &v, Double_t &e)
|
---|
821 | {
|
---|
822 | // Valid digits
|
---|
823 | Int_t i = TMath::FloorNint(TMath::Log10(v))-TMath::FloorNint(TMath::Log10(e));
|
---|
824 |
|
---|
825 | // Check if error starts with 1 or 2. In this case use one
|
---|
826 | // more valid digit
|
---|
827 | TString error = Form("%.0e", e);
|
---|
828 | if (error[0]=='1' || error[0]=='2')
|
---|
829 | {
|
---|
830 | i++;
|
---|
831 | error = Form("%.1e", e);
|
---|
832 | }
|
---|
833 |
|
---|
834 | const char *fmt = Form("%%.%de", i);
|
---|
835 |
|
---|
836 | v = atof(Form(fmt, v));
|
---|
837 | e = error.Atof();
|
---|
838 | }
|
---|