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

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