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

Last change on this file since 11426 was 11426, checked in by tbretz, 9 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.