source: trunk/Mars/mbase/MMath.cc@ 13249

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