source: trunk/MagicSoft/Mars/mbase/MMath.cc@ 8986

Last change on this file since 8986 was 8978, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 24.2 KB
Line 
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//
62Double_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//
77Double_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//
117Double_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/*
131Double_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//
162Double_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//
176Double_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//
194Double_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//
216template <class Size, class Element>
217Double_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//
254Double_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//
264Double_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//
274Double_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//
284Double_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//
294Double_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//
304Double_t MMath::MedianDev(Long64_t n, const Long64_t *a, Double_t &med)
305{
306 return MedianDevImp(n, a, med);
307}
308
309Double_t MMath::MedianDev(Long64_t n, const Short_t *a) { Double_t med; return MedianDevImp(n, a, med); }
310Double_t MMath::MedianDev(Long64_t n, const Int_t *a) { Double_t med; return MedianDevImp(n, a, med); }
311Double_t MMath::MedianDev(Long64_t n, const Float_t *a) { Double_t med; return MedianDevImp(n, a, med); }
312Double_t MMath::MedianDev(Long64_t n, const Double_t *a) { Double_t med; return MedianDevImp(n, a, med); }
313Double_t MMath::MedianDev(Long64_t n, const Long_t *a) { Double_t med; return MedianDevImp(n, a, med); }
314Double_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//
323void 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//
340TVector3 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//
382Double_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//
394Double_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
400Double_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
410Double_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//
430TArrayD 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//
467TArrayD 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//
503TArrayD 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//
535TArrayD 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//
569TVector2 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//
615Int_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//
639static 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//
709Int_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//
820void 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}
Note: See TracBrowser for help on using the repository browser.