Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TMath.h
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Authors: Rene Brun, Anna Kreshuk, Eddy Offermann, Fons Rademakers 29/07/95
3
4/*************************************************************************
5 * Copyright (C) 1995-2004, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12#ifndef ROOT_TMath
13#define ROOT_TMath
14
15#include "TMathBase.h"
16
17#include "TError.h"
18#include "ROOT/RSpan.hxx"
19
20#include <algorithm>
21#include <limits>
22#include <cmath>
23#include <vector>
24
25////////////////////////////////////////////////////////////////////////////////
26///
27/// TMath
28///
29/// Encapsulate most frequently used Math functions.
30/// NB. The basic functions Min, Max, Abs and Sign are defined
31/// in TMathBase.
32
33namespace TMath {
34
35////////////////////////////////////////////////////////////////////////////////
36// Fundamental constants
37
38////////////////////////////////////////////////////////////////////////////////
39/// \f$ \pi\f$
40constexpr Double_t Pi()
41{
42 return 3.14159265358979323846;
43}
44
45////////////////////////////////////////////////////////////////////////////////
46/// \f$ 2\pi\f$
47constexpr Double_t TwoPi()
48{
49 return 2.0 * Pi();
50}
51
52////////////////////////////////////////////////////////////////////////////////
53/// \f$ \frac{\pi}{2} \f$
54constexpr Double_t PiOver2()
55{
56 return Pi() / 2.0;
57}
58
59////////////////////////////////////////////////////////////////////////////////
60/// \f$ \frac{\pi}{4} \f$
61constexpr Double_t PiOver4()
62{
63 return Pi() / 4.0;
64}
65
66////////////////////////////////////////////////////////////////////////////////
67/// \f$ \frac{1.}{\pi}\f$
68constexpr Double_t InvPi()
69{
70 return 1.0 / Pi();
71}
72
73////////////////////////////////////////////////////////////////////////////////
74/// Conversion from radian to degree: \f$ \frac{180}{\pi} \f$
75constexpr Double_t RadToDeg()
76{
77 return 180.0 / Pi();
78}
79
80////////////////////////////////////////////////////////////////////////////////
81/// Conversion from degree to radian: \f$ \frac{\pi}{180} \f$
82constexpr Double_t DegToRad()
83{
84 return Pi() / 180.0;
85}
86
87////////////////////////////////////////////////////////////////////////////////
88/// \f$ \sqrt{2} \f$
89constexpr Double_t Sqrt2()
90{
91 return 1.4142135623730950488016887242097;
92}
93
94////////////////////////////////////////////////////////////////////////////////
95/// Base of natural log: \f$ e \f$
96constexpr Double_t E()
97{
98 return 2.71828182845904523536;
99}
100
101////////////////////////////////////////////////////////////////////////////////
102/// Natural log of 10 (to convert log to ln)
103constexpr Double_t Ln10()
104{
105 return 2.30258509299404568402;
106}
107
108////////////////////////////////////////////////////////////////////////////////
109/// Base-10 log of e (to convert ln to log)
110constexpr Double_t LogE()
111{
112 return 0.43429448190325182765;
113}
114
115////////////////////////////////////////////////////////////////////////////////
116/// Velocity of light in \f$ m s^{-1} \f$
117constexpr Double_t C()
118{
119 return 2.99792458e8;
120}
121
122////////////////////////////////////////////////////////////////////////////////
123/// \f$ cm s^{-1} \f$
124constexpr Double_t Ccgs()
125{
126 return 100.0 * C();
127}
128
129////////////////////////////////////////////////////////////////////////////////
130/// Speed of light uncertainty.
132{
133 return 0.0;
134}
135
136////////////////////////////////////////////////////////////////////////////////
137/// Gravitational constant in: \f$ m^{3} kg^{-1} s^{-2} \f$
138constexpr Double_t G()
139{
140 // use 2018 value from NIST (https://physics.nist.gov/cgi-bin/cuu/Value?bg|search_for=G)
141 return 6.67430e-11;
142}
143
144////////////////////////////////////////////////////////////////////////////////
145/// \f$ cm^{3} g^{-1} s^{-2} \f$
146constexpr Double_t Gcgs()
147{
148 return G() * 1000.0;
149}
150
151////////////////////////////////////////////////////////////////////////////////
152/// Gravitational constant uncertainty.
154{
155 // use 2018 value from NIST
156 return 0.00015e-11;
157}
158
159////////////////////////////////////////////////////////////////////////////////
160/// \f$ \frac{G}{\hbar C} \f$ in \f$ (GeV/c^{2})^{-2} \f$
161constexpr Double_t GhbarC()
162{
163 // use new value from NIST (2018)
164 return 6.70883e-39;
165}
166
167////////////////////////////////////////////////////////////////////////////////
168/// \f$ \frac{G}{\hbar C} \f$ uncertainty.
170{
171 // use new value from NIST (2018)
172 return 0.00015e-39;
173}
174
175////////////////////////////////////////////////////////////////////////////////
176/// Standard acceleration of gravity in \f$ m s^{-2} \f$
177constexpr Double_t Gn()
178{
179 return 9.80665;
180}
181
182////////////////////////////////////////////////////////////////////////////////
183/// Standard acceleration of gravity uncertainty.
185{
186 return 0.0;
187}
188
189////////////////////////////////////////////////////////////////////////////////
190/// Planck's constant in \f$ J s \f$: \f$ h \f$
191constexpr Double_t H()
192{
193 return 6.62607015e-34;
194}
195
196////////////////////////////////////////////////////////////////////////////////
197/// \f$ erg s \f$
198constexpr Double_t Hcgs()
199{
200 return 1.0e7 * H();
201}
202
203////////////////////////////////////////////////////////////////////////////////
204/// Planck's constant uncertainty.
206{
207 // Planck constant is exact according to 2019 redefinition
208 // (https://en.wikipedia.org/wiki/2019_redefinition_of_the_SI_base_units)
209 return 0.0;
210}
211
212////////////////////////////////////////////////////////////////////////////////
213/// \f$ \hbar \f$ in \f$ J s \f$: \f$ \hbar = \frac{h}{2\pi} \f$
214constexpr Double_t Hbar()
215{
216 return 1.054571817e-34;
217}
218
219////////////////////////////////////////////////////////////////////////////////
220/// \f$ erg s \f$
221constexpr Double_t Hbarcgs()
222{
223 return 1.0e7 * Hbar();
224}
225
226////////////////////////////////////////////////////////////////////////////////
227/// \f$ \hbar \f$ uncertainty.
229{
230 // hbar is an exact constant
231 return 0.0;
232}
233
234////////////////////////////////////////////////////////////////////////////////
235/// \f$ hc \f$ in \f$ J m \f$
236constexpr Double_t HC()
237{
238 return H() * C();
239}
240
241////////////////////////////////////////////////////////////////////////////////
242/// \f$ erg cm \f$
243constexpr Double_t HCcgs()
244{
245 return Hcgs() * Ccgs();
246}
247
248////////////////////////////////////////////////////////////////////////////////
249/// Boltzmann's constant in \f$ J K^{-1} \f$: \f$ k \f$
250constexpr Double_t K()
251{
252 return 1.380649e-23;
253}
254
255////////////////////////////////////////////////////////////////////////////////
256/// \f$ erg K^{-1} \f$
257constexpr Double_t Kcgs()
258{
259 return 1.0e7 * K();
260}
261
262////////////////////////////////////////////////////////////////////////////////
263/// Boltzmann's constant uncertainty.
265{
266 // constant is exact according to 2019 redefinition
267 // (https://en.wikipedia.org/wiki/2019_redefinition_of_the_SI_base_units)
268 return 0.0;
269}
270
271////////////////////////////////////////////////////////////////////////////////
272/// Stefan-Boltzmann constant in \f$ W m^{-2} K^{-4}\f$: \f$ \sigma \f$
273constexpr Double_t Sigma()
274{
275 return 5.670373e-8;
276}
277
278////////////////////////////////////////////////////////////////////////////////
279/// Stefan-Boltzmann constant uncertainty.
281{
282 return 0.0;
283}
284
285////////////////////////////////////////////////////////////////////////////////
286/// Avogadro constant (Avogadro's Number) in \f$ mol^{-1} \f$
287constexpr Double_t Na()
288{
289 return 6.02214076e+23;
290}
291
292////////////////////////////////////////////////////////////////////////////////
293/// Avogadro constant (Avogadro's Number) uncertainty.
295{
296 // constant is exact according to 2019 redefinition
297 // (https://en.wikipedia.org/wiki/2019_redefinition_of_the_SI_base_units)
298 return 0.0;
299}
300
301////////////////////////////////////////////////////////////////////////////////
302/// [Universal gas constant](http://scienceworld.wolfram.com/physics/UniversalGasConstant.html)
303/// (\f$ Na K \f$) in \f$ J K^{-1} mol^{-1} \f$
304//
305constexpr Double_t R()
306{
307 return K() * Na();
308}
309
310////////////////////////////////////////////////////////////////////////////////
311/// Universal gas constant uncertainty.
313{
314 return R() * ((KUncertainty() / K()) + (NaUncertainty() / Na()));
315}
316
317////////////////////////////////////////////////////////////////////////////////
318/// [Molecular weight of dry air 1976 US Standard Atmosphere](http://atmos.nmsu.edu/jsdap/encyclopediawork.html)
319/// in \f$ kg kmol^{-1} \f$ or \f$ gm mol^{-1} \f$
320constexpr Double_t MWair()
321{
322 return 28.9644;
323}
324
325////////////////////////////////////////////////////////////////////////////////
326/// [Dry Air Gas Constant (R / MWair)](http://atmos.nmsu.edu/education_and_outreach/encyclopedia/gas_constant.htm)
327/// in \f$ J kg^{-1} K^{-1} \f$
328constexpr Double_t Rgair()
329{
330 return (1000.0 * R()) / MWair();
331}
332
333////////////////////////////////////////////////////////////////////////////////
334/// Euler-Mascheroni Constant.
336{
337 return 0.577215664901532860606512090082402431042;
338}
339
340////////////////////////////////////////////////////////////////////////////////
341/// Elementary charge in \f$ C \f$ .
342constexpr Double_t Qe()
343{
344 return 1.602176634e-19;
345}
346
347////////////////////////////////////////////////////////////////////////////////
348/// Elementary charge uncertainty.
350{
351 // constant is exact according to 2019 redefinition
352 // (https://en.wikipedia.org/wiki/2019_redefinition_of_the_SI_base_units)
353 return 0.0;
354}
355
356////////////////////////////////////////////////////////////////////////////////
357// Mathematical Functions
358
359////////////////////////////////////////////////////////////////////////////////
360// Trigonometrical Functions
361
362inline Double_t Sin(Double_t);
363inline Double_t Cos(Double_t);
364inline Double_t Tan(Double_t);
365inline Double_t SinH(Double_t);
366inline Double_t CosH(Double_t);
367inline Double_t TanH(Double_t);
368inline Double_t ASin(Double_t);
369inline Double_t ACos(Double_t);
370inline Double_t ATan(Double_t);
376
377////////////////////////////////////////////////////////////////////////////////
378// Elementary Functions
379
380inline Double_t Ceil(Double_t x);
381inline Int_t CeilNint(Double_t x);
382inline Double_t Floor(Double_t x);
383inline Int_t FloorNint(Double_t x);
384template <typename T>
385inline Int_t Nint(T x);
386
387inline Double_t Sq(Double_t x);
388inline Double_t Sqrt(Double_t x);
389inline Double_t Exp(Double_t x);
390inline Double_t Ldexp(Double_t x, Int_t exp);
396inline Double_t Power(Double_t x, Int_t y);
397inline Double_t Log(Double_t x);
399inline Double_t Log10(Double_t x);
400inline Int_t Finite(Double_t x);
401inline Int_t Finite(Float_t x);
402inline Bool_t IsNaN(Double_t x);
403inline Bool_t IsNaN(Float_t x);
404
405inline Double_t QuietNaN();
406inline Double_t SignalingNaN();
407inline Double_t Infinity();
408
409template <typename T>
410struct Limits {
411 inline static T Min();
412 inline static T Max();
413 inline static T Epsilon();
414 };
415
416 // Some integer math
418
419 /// Comparing floating points.
420 /// Returns `kTRUE` if the absolute difference between `af` and `bf` is less than `epsilon`.
422 return af == bf || // shortcut for exact equality (necessary because otherwise (inf != inf))
423 TMath::Abs(af-bf) < epsilon ||
424 TMath::Abs(af - bf) < Limits<Double_t>::Min(); // handle 0 < 0 case
425
426 }
427 /// Comparing floating points.
428 /// Returns `kTRUE` if the relative difference between `af` and `bf` is less than `relPrec`.
430 return af == bf || // shortcut for exact equality (necessary because otherwise (inf != inf))
431 TMath::Abs(af - bf) <= 0.5 * relPrec * (TMath::Abs(af) + TMath::Abs(bf)) ||
432 TMath::Abs(af - bf) < Limits<Double_t>::Min(); // handle denormals
433 }
434
435 /////////////////////////////////////////////////////////////////////////////
436 // Array Algorithms
437
438 // Min, Max of an array
439 template <typename T> T MinElement(Long64_t n, const T *a);
440 template <typename T> T MaxElement(Long64_t n, const T *a);
441
442 // Locate Min, Max element number in an array
443 template <typename T> Long64_t LocMin(Long64_t n, const T *a);
444 template <typename Iterator> Iterator LocMin(Iterator first, Iterator last);
445 template <typename T> Long64_t LocMax(Long64_t n, const T *a);
446 template <typename Iterator> Iterator LocMax(Iterator first, Iterator last);
447
448 // Derivatives of an array
449 template <typename T> T *Gradient(Long64_t n, T *f, double h = 1);
450 template <typename T> T *Laplacian(Long64_t n, T *f, double h = 1);
451
452 // Hashing
453 ULong_t Hash(const void *txt, Int_t ntxt);
454 ULong_t Hash(const char *str);
455
458
459 Bool_t Permute(Int_t n, Int_t *a); // Find permutations
460
461 /////////////////////////////////////////////////////////////////////////////
462 // Geometrical Functions
463
464 //Sample quantiles
466 Bool_t isSorted=kTRUE, Int_t *index = nullptr, Int_t type=7);
467
468 // IsInside
469 template <typename T> Bool_t IsInside(T xp, T yp, Int_t np, T *x, T *y);
470
471 // Calculate the Cross Product of two vectors
472 template <typename T> T *Cross(const T v1[3],const T v2[3], T out[3]);
473
474 Float_t Normalize(Float_t v[3]); // Normalize a vector
475 Double_t Normalize(Double_t v[3]); // Normalize a vector
476
477 //Calculate the Normalized Cross Product of two vectors
478 template <typename T> inline T NormCross(const T v1[3],const T v2[3],T out[3]);
479
480 // Calculate a normal vector of a plane
481 template <typename T> T *Normal2Plane(const T v1[3],const T v2[3],const T v3[3], T normal[3]);
482
483 /////////////////////////////////////////////////////////////////////////////
484 // Polynomial Functions
485
486 Bool_t RootsCubic(const Double_t coef[4],Double_t &a, Double_t &b, Double_t &c);
487
488 /////////////////////////////////////////////////////////////////////////////
489 // Statistic Functions
490
491 Double_t Binomial(Int_t n,Int_t k); // Calculate the binomial coefficient n over k
517
518 /////////////////////////////////////////////////////////////////////////////
519 // Statistics over arrays
520
521 // Mean, Geometric Mean, Median, RMS(sigma), ModeHalfSample
522
523 template <typename T> Double_t Mean(Long64_t n, const T *a, const Double_t *w=nullptr);
524 template <typename Iterator> Double_t Mean(Iterator first, Iterator last);
525 template <typename Iterator, typename WeightIterator> Double_t Mean(Iterator first, Iterator last, WeightIterator wfirst);
526
527 template <typename T> Double_t GeomMean(Long64_t n, const T *a);
528 template <typename Iterator> Double_t GeomMean(Iterator first, Iterator last);
529
530 template <typename T> Double_t RMS(Long64_t n, const T *a, const Double_t *w=nullptr);
531 template <typename Iterator> Double_t RMS(Iterator first, Iterator last);
532 template <typename Iterator, typename WeightIterator> Double_t RMS(Iterator first, Iterator last, WeightIterator wfirst);
533
534 template <typename T> Double_t StdDev(Long64_t n, const T *a, const Double_t * w = nullptr) { return RMS<T>(n,a,w); } /// Same as RMS
535 template <typename Iterator> Double_t StdDev(Iterator first, Iterator last) { return RMS<Iterator>(first,last); } /// Same as RMS
536 template <typename Iterator, typename WeightIterator> Double_t StdDev(Iterator first, Iterator last, WeightIterator wfirst) { return RMS<Iterator,WeightIterator>(first,last,wfirst); } /// Same as RMS
537
538 template <typename T> Double_t Median(Long64_t n, const T *a, const Double_t *w=nullptr, Long64_t *work=nullptr);
539 template <typename T> Double_t ModeHalfSample(Long64_t n, const T *a, const Double_t *w = nullptr);
540
541 //k-th order statistic
542 template <class Element, typename Size> Element KOrdStat(Size n, const Element *a, Size k, Size *work = 0);
543
544 /////////////////////////////////////////////////////////////////////////////
545 // Special Functions
546
552
553 // Bessel functions
554 Double_t BesselI(Int_t n,Double_t x); /// Integer order modified Bessel function I_n(x)
555 Double_t BesselK(Int_t n,Double_t x); /// Integer order modified Bessel function K_n(x)
556 Double_t BesselI0(Double_t x); /// Modified Bessel function I_0(x)
557 Double_t BesselK0(Double_t x); /// Modified Bessel function K_0(x)
558 Double_t BesselI1(Double_t x); /// Modified Bessel function I_1(x)
559 Double_t BesselK1(Double_t x); /// Modified Bessel function K_1(x)
560 Double_t BesselJ0(Double_t x); /// Bessel function J0(x) for any real x
561 Double_t BesselJ1(Double_t x); /// Bessel function J1(x) for any real x
562 Double_t BesselY0(Double_t x); /// Bessel function Y0(x) for positive x
563 Double_t BesselY1(Double_t x); /// Bessel function Y1(x) for positive x
564 Double_t StruveH0(Double_t x); /// Struve functions of order 0
565 Double_t StruveH1(Double_t x); /// Struve functions of order 1
566 Double_t StruveL0(Double_t x); /// Modified Struve functions of order 0
567 Double_t StruveL1(Double_t x); /// Modified Struve functions of order 1
568
577 Double_t GammaDist(Double_t x, Double_t gamma, Double_t mu=0, Double_t beta=1);
579
580 void KNNDensity(std::span<const double> observations, std::span<const double> queries, std::span<double> result,
581 int k, double dmin = 0.0);
582}
583
584////////////////////////////////////////////////////////////////////////////////
585// Trig and other functions
586
587#include <cfloat>
588#include <cmath>
589
590#if defined(R__WIN32)
591# ifndef finite
592# define finite _finite
593# endif
594#endif
595
596////////////////////////////////////////////////////////////////////////////////
597/// Returns the sine of an angle of `x` radians.
598
600 { return sin(x); }
601
602////////////////////////////////////////////////////////////////////////////////
603/// Returns the cosine of an angle of `x` radians.
604
606 { return cos(x); }
607
608////////////////////////////////////////////////////////////////////////////////
609/// Returns the tangent of an angle of `x` radians.
610
612 { return tan(x); }
613
614////////////////////////////////////////////////////////////////////////////////
615/// Returns the hyperbolic sine of `x.
616
618 { return sinh(x); }
619
620////////////////////////////////////////////////////////////////////////////////
621/// Returns the hyperbolic cosine of `x`.
622
624 { return cosh(x); }
625
626////////////////////////////////////////////////////////////////////////////////
627/// Returns the hyperbolic tangent of `x`.
628
630 { return tanh(x); }
631
632////////////////////////////////////////////////////////////////////////////////
633/// Returns the principal value of the arc sine of `x`, expressed in radians.
634
636 {
637 return asin(x);
638 }
639
640////////////////////////////////////////////////////////////////////////////////
641/// Returns the principal value of the arc cosine of `x`, expressed in radians.
642
644 {
645 return acos(x);
646 }
647
648////////////////////////////////////////////////////////////////////////////////
649/// Returns the principal value of the arc tangent of `x`, expressed in radians.
650
652 { return atan(x); }
653
654////////////////////////////////////////////////////////////////////////////////
655/// Returns the principal value of the arc tangent of `y/x`, expressed in radians.
656
658 { if (x != 0) return atan2(y, x);
659 if (y == 0) return 0;
660 if (y > 0) return Pi()/2;
661 else return -Pi()/2;
662 }
663
664////////////////////////////////////////////////////////////////////////////////
665/// Returns `x*x`.
666
668 { return x*x; }
669
670////////////////////////////////////////////////////////////////////////////////
671/// Returns the square root of x.
672
674 { return sqrt(x); }
675
676////////////////////////////////////////////////////////////////////////////////
677/// Rounds `x` upward, returning the smallest integral value that is not less than `x`.
678
680 { return ceil(x); }
681
682////////////////////////////////////////////////////////////////////////////////
683/// Returns the nearest integer of `TMath::Ceil(x)`.
684
686 { return TMath::Nint(ceil(x)); }
687
688////////////////////////////////////////////////////////////////////////////////
689/// Rounds `x` downward, returning the largest integral value that is not greater than `x`.
690
692 { return floor(x); }
693
694////////////////////////////////////////////////////////////////////////////////
695/// Returns the nearest integer of `TMath::Floor(x)`.
696
698 { return TMath::Nint(floor(x)); }
699
700////////////////////////////////////////////////////////////////////////////////
701/// Round to nearest integer. Rounds half integers to the nearest even integer.
702
703template<typename T>
705{
706 int i;
707 if (x >= 0) {
708 i = int(x + 0.5);
709 if ( i & 1 && x + 0.5 == T(i) ) i--;
710 } else {
711 i = int(x - 0.5);
712 if ( i & 1 && x - 0.5 == T(i) ) i++;
713 }
714 return i;
715}
716
717////////////////////////////////////////////////////////////////////////////////
718/// Returns the base-e exponential function of x, which is e raised to the power `x`.
719
721 { return exp(x); }
722
723////////////////////////////////////////////////////////////////////////////////
724/// Returns the result of multiplying `x` (the significant) by 2 raised to the power of `exp` (the exponent).
725
727 { return ldexp(x, exp); }
728
729////////////////////////////////////////////////////////////////////////////////
730/// Returns `x` raised to the power `y`.
731
733 { return std::pow(x,y); }
734
735////////////////////////////////////////////////////////////////////////////////
736/// Returns `x` raised to the power `y`.
737
739 { return std::pow(x,(LongDouble_t)y); }
740
741////////////////////////////////////////////////////////////////////////////////
742/// Returns `x` raised to the power `y`.
743
745 { return std::pow(x,y); }
746
747////////////////////////////////////////////////////////////////////////////////
748/// Returns `x` raised to the power `y`.
749
751 { return pow(x, y); }
752
753////////////////////////////////////////////////////////////////////////////////
754/// Returns `x` raised to the power `y`.
755
757#ifdef R__ANSISTREAM
758 return std::pow(x, y);
759#else
760 return pow(x, (Double_t) y);
761#endif
762}
763
764////////////////////////////////////////////////////////////////////////////////
765/// Returns the natural logarithm of `x`.
766
768 { return log(x); }
769
770////////////////////////////////////////////////////////////////////////////////
771/// Returns the common (base-10) logarithm of `x`.
772
774 { return log10(x); }
775
776////////////////////////////////////////////////////////////////////////////////
777/// Check if it is finite with a mask in order to be consistent in presence of
778/// fast math.
779/// Inspired from the CMSSW FWCore/Utilities package
780
782#if defined(R__FAST_MATH)
783
784{
785 const unsigned long long mask = 0x7FF0000000000000LL;
786 union { unsigned long long l; double d;} v;
787 v.d =x;
788 return (v.l&mask)!=mask;
789}
790#else
791# if defined(R__HPUX11)
792 { return isfinite(x); }
793# elif defined(R__MACOSX)
794# ifdef isfinite
795 // from math.h
796 { return isfinite(x); }
797# else
798 // from cmath
799 { return std::isfinite(x); }
800# endif
801# else
802 { return finite(x); }
803# endif
804#endif
805
806////////////////////////////////////////////////////////////////////////////////
807/// Check if it is finite with a mask in order to be consistent in presence of
808/// fast math.
809/// Inspired from the CMSSW FWCore/Utilities package
810
812#if defined(R__FAST_MATH)
813
814{
815 const unsigned int mask = 0x7f800000;
816 union { unsigned int l; float d;} v;
817 v.d =x;
818 return (v.l&mask)!=mask;
819}
820#else
821{ return std::isfinite(x); }
822#endif
823
824// This namespace provides all the routines necessary for checking if a number
825// is a NaN also in presence of optimisations affecting the behaviour of the
826// floating point calculations.
827// Inspired from the CMSSW FWCore/Utilities package
828
829#if defined (R__FAST_MATH)
830namespace ROOT {
831namespace Internal {
832namespace Math {
833// abridged from GNU libc 2.6.1 - in detail from
834// math/math_private.h
835// sysdeps/ieee754/ldbl-96/math_ldbl.h
836
837// part of this file:
838 /*
839 * ====================================================
840 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
841 *
842 * Developed at SunPro, a Sun Microsystems, Inc. business.
843 * Permission to use, copy, modify, and distribute this
844 * software is freely granted, provided that this notice
845 * is preserved.
846 * ====================================================
847 */
848
849 // A union which permits us to convert between a double and two 32 bit ints.
850 typedef union {
852 struct {
853 UInt_t lsw;
854 UInt_t msw;
855 } parts;
857
858#define EXTRACT_WORDS(ix0,ix1,d) \
859 do { \
860 ieee_double_shape_type ew_u; \
861 ew_u.value = (d); \
862 (ix0) = ew_u.parts.msw; \
863 (ix1) = ew_u.parts.lsw; \
864 } while (0)
865
866 inline Bool_t IsNaN(Double_t x)
867 {
868 UInt_t hx, lx;
869
870 EXTRACT_WORDS(hx, lx, x);
871
872 lx |= hx & 0xfffff;
873 hx &= 0x7ff00000;
874 return (hx == 0x7ff00000) && (lx != 0);
875 }
876
877 typedef union {
879 UInt_t word;
881
882#define GET_FLOAT_WORD(i,d) \
883 do { \
884 ieee_float_shape_type gf_u; \
885 gf_u.value = (d); \
886 (i) = gf_u.word; \
887 } while (0)
888
889 inline Bool_t IsNaN(Float_t x)
890 {
891 UInt_t wx;
893 wx &= 0x7fffffff;
894 return (Bool_t)(wx > 0x7f800000);
895 }
896} } } // end NS ROOT::Internal::Math
897#endif // End R__FAST_MATH
898
899#if defined(R__FAST_MATH)
900 inline Bool_t TMath::IsNaN(Double_t x) { return ROOT::Internal::Math::IsNaN(x); }
901 inline Bool_t TMath::IsNaN(Float_t x) { return ROOT::Internal::Math::IsNaN(x); }
902#else
903 inline Bool_t TMath::IsNaN(Double_t x) { return std::isnan(x); }
904 inline Bool_t TMath::IsNaN(Float_t x) { return std::isnan(x); }
905#endif
906
907////////////////////////////////////////////////////////////////////////////////
908// Wrapper to numeric_limits
909
910////////////////////////////////////////////////////////////////////////////////
911/// Returns a quiet NaN as [defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Quiet_NaN).
912
914
915 return std::numeric_limits<Double_t>::quiet_NaN();
916}
917
918////////////////////////////////////////////////////////////////////////////////
919/// Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN).
920
922 return std::numeric_limits<Double_t>::signaling_NaN();
923}
924
925////////////////////////////////////////////////////////////////////////////////
926/// Returns an infinity as defined by the IEEE standard.
927
929 return std::numeric_limits<Double_t>::infinity();
930}
931
932////////////////////////////////////////////////////////////////////////////////
933/// Returns maximum representation for type T.
934
935template<typename T>
937 return (std::numeric_limits<T>::min)(); //N.B. use this signature to avoid class with macro min() on Windows
938}
939
940////////////////////////////////////////////////////////////////////////////////
941/// Returns minimum double representation.
942
943template<typename T>
945 return (std::numeric_limits<T>::max)(); //N.B. use this signature to avoid class with macro max() on Windows
946}
947
948////////////////////////////////////////////////////////////////////////////////
949/// Returns minimum double representation.
950
951template<typename T>
953 return std::numeric_limits<T>::epsilon();
954}
955
956////////////////////////////////////////////////////////////////////////////////
957// Advanced.
958
959////////////////////////////////////////////////////////////////////////////////
960/// Calculates the Normalized Cross Product of two vectors.
961
962template <typename T> inline T TMath::NormCross(const T v1[3],const T v2[3],T out[3])
963{
964 return Normalize(Cross(v1,v2,out));
965}
966
967////////////////////////////////////////////////////////////////////////////////
968/// Returns minimum of array a of length n.
969
970template <typename T>
972 return *std::min_element(a,a+n);
973}
974
975////////////////////////////////////////////////////////////////////////////////
976/// Returns maximum of array a of length n.
977
978template <typename T>
980 return *std::max_element(a,a+n);
981}
982
983////////////////////////////////////////////////////////////////////////////////
984/// Returns index of array with the minimum element.
985/// If more than one element is minimum returns first found.
986///
987/// Implement here since this one is found to be faster (mainly on 64 bit machines)
988/// than stl generic implementation.
989/// When performing the comparison, the STL implementation needs to de-reference both the array iterator
990/// and the iterator pointing to the resulting minimum location
991
992template <typename T>
994 if (n <= 0 || !a) return -1;
995 T xmin = a[0];
996 Long64_t loc = 0;
997 for (Long64_t i = 1; i < n; i++) {
998 if (xmin > a[i]) {
999 xmin = a[i];
1000 loc = i;
1001 }
1002 }
1003 return loc;
1004}
1005
1006////////////////////////////////////////////////////////////////////////////////
1007/// Returns index of array with the minimum element.
1008/// If more than one element is minimum returns first found.
1009
1010template <typename Iterator>
1011Iterator TMath::LocMin(Iterator first, Iterator last) {
1012
1013 return std::min_element(first, last);
1014}
1015
1016////////////////////////////////////////////////////////////////////////////////
1017/// \brief Calculate the one-dimensional gradient of an array with length n.
1018/// The first value in the returned array is a forward difference,
1019/// the next n-2 values are central differences, and the last is a backward difference.
1020///
1021/// \note Function leads to undefined behavior if n does not match the length of f
1022/// \param n the number of points in the array
1023/// \param f the array of points.
1024/// \param h the step size. The default step size is 1.
1025/// \return an array of size n with the gradient. Returns nullptr if n < 2 or f empty. Ownership is transferred to the
1026/// caller.
1027
1028template <typename T>
1029T *TMath::Gradient(const Long64_t n, T *f, const double h)
1030{
1031 if (!f) {
1032 ::Error("TMath::Gradient", "Input parameter f is empty.");
1033 return nullptr;
1034 } else if (n < 2) {
1035 ::Error("TMath::Gradient", "Input parameter n=%lld is smaller than 2.", n);
1036 return nullptr;
1037 }
1038 Long64_t i = 1;
1039 T *result = new T[n];
1040
1041 // Forward difference
1042 result[0] = (f[1] - f[0]) / h;
1043
1044 // Central difference
1045 while (i < n - 1) {
1046 result[i] = (f[i + 1] - f[i - 1]) / (2 * h);
1047 i++;
1048 }
1049 // Backward difference
1050 result[i] = (f[i] - f[i - 1]) / h;
1051 return result;
1052}
1053
1054////////////////////////////////////////////////////////////////////////////////
1055/// \brief Calculate the Laplacian of an array with length n.
1056/// The first value in the returned array is a forward difference,
1057/// the next n-2 values are central differences, and the last is a backward difference.
1058///
1059/// \note Function leads to undefined behavior if n does not match the length of f
1060/// \param n the number of points in the array
1061/// \param f the array of points.
1062/// \param h the step size. The default step size is 1.
1063/// \return an array of size n with the laplacian. Returns nullptr if n < 4 or f empty. Ownership is transferred to the
1064/// caller.
1065
1066template <typename T>
1067T *TMath::Laplacian(const Long64_t n, T *f, const double h)
1068{
1069 if (!f) {
1070 ::Error("TMath::Laplacian", "Input parameter f is empty.");
1071 return nullptr;
1072 } else if (n < 4) {
1073 ::Error("TMath::Laplacian", "Input parameter n=%lld is smaller than 4.", n);
1074 return nullptr;
1075 }
1076 Long64_t i = 1;
1077 T *result = new T[n];
1078
1079 // Forward difference
1080 result[0] = (4 * f[2] + 2 * f[0] - 5 * f[1] - f[3]) / (4 * h * h);
1081
1082 // Central difference
1083 while (i < n - 1) {
1084 result[i] = (f[i + 1] + f[i - 1] - 2 * f[i]) / (4 * h * h);
1085 i++;
1086 }
1087 // Backward difference
1088 result[i] = (2 * f[i] - 5 * f[i - 1] + 4 * f[i - 2] - f[i - 3]) / (4 * h * h);
1089 return result;
1090}
1091
1092////////////////////////////////////////////////////////////////////////////////
1093/// Returns index of array with the maximum element.
1094/// If more than one element is maximum returns first found.
1095///
1096/// Implement here since it is faster (see comment in LocMin function)
1097
1098template <typename T>
1100 if (n <= 0 || !a) return -1;
1101 T xmax = a[0];
1102 Long64_t loc = 0;
1103 for (Long64_t i = 1; i < n; i++) {
1104 if (xmax < a[i]) {
1105 xmax = a[i];
1106 loc = i;
1107 }
1108 }
1109 return loc;
1110}
1111
1112////////////////////////////////////////////////////////////////////////////////
1113/// Returns index of array with the maximum element.
1114/// If more than one element is maximum returns first found.
1115
1116template <typename Iterator>
1117Iterator TMath::LocMax(Iterator first, Iterator last)
1118{
1119
1120 return std::max_element(first, last);
1121}
1122
1123////////////////////////////////////////////////////////////////////////////////
1124/// Returns the weighted mean of an array defined by the iterators.
1125
1126template <typename Iterator>
1127Double_t TMath::Mean(Iterator first, Iterator last)
1128{
1129 Double_t sum = 0;
1130 Double_t sumw = 0;
1131 while ( first != last )
1132 {
1133 sum += *first;
1134 sumw += 1;
1135 first++;
1136 }
1137
1138 return sum/sumw;
1139}
1140
1141////////////////////////////////////////////////////////////////////////////////
1142/// Returns the weighted mean of an array defined by the first and
1143/// last iterators. The w iterator should point to the first element
1144/// of a vector of weights of the same size as the main array.
1145
1146template <typename Iterator, typename WeightIterator>
1147Double_t TMath::Mean(Iterator first, Iterator last, WeightIterator w)
1148{
1149
1150 Double_t sum = 0;
1151 Double_t sumw = 0;
1152 int i = 0;
1153 while ( first != last ) {
1154 if ( *w < 0) {
1155 ::Error("TMath::Mean","w[%d] = %.4e < 0 ?!",i,*w);
1156 return 0;
1157 }
1158 sum += (*w) * (*first);
1159 sumw += (*w) ;
1160 ++w;
1161 ++first;
1162 ++i;
1163 }
1164 if (sumw <= 0) {
1165 ::Error("TMath::Mean","sum of weights == 0 ?!");
1166 return 0;
1167 }
1168
1169 return sum/sumw;
1170}
1171
1172////////////////////////////////////////////////////////////////////////////////
1173/// Returns the weighted mean of an array a with length n.
1174
1175template <typename T>
1177{
1178 if (w) {
1179 return TMath::Mean(a, a+n, w);
1180 } else {
1181 return TMath::Mean(a, a+n);
1182 }
1183}
1184
1185////////////////////////////////////////////////////////////////////////////////
1186/// Returns the geometric mean of an array defined by the iterators.
1187/// \f[ GeomMean = (\prod_{i=0}^{n-1} |a[i]|)^{1/n} \f]
1188
1189template <typename Iterator>
1190Double_t TMath::GeomMean(Iterator first, Iterator last)
1191{
1192 Double_t logsum = 0.;
1193 Long64_t n = 0;
1194 while ( first != last ) {
1195 if (*first == 0) return 0.;
1196 Double_t absa = (Double_t) TMath::Abs(*first);
1198 ++first;
1199 ++n;
1200 }
1201
1202 return TMath::Exp(logsum/n);
1203}
1204
1205////////////////////////////////////////////////////////////////////////////////
1206/// Returns the geometric mean of an array a of size n.
1207/// \f[ GeomMean = (\prod_{i=0}^{n-1} |a[i]|)^{1/n} \f]
1208
1209template <typename T>
1211{
1212 return TMath::GeomMean(a, a+n);
1213}
1214
1215////////////////////////////////////////////////////////////////////////////////
1216/// Returns the Standard Deviation of an array defined by the iterators.
1217/// Note that this function returns the sigma(standard deviation) and
1218/// not the root mean square of the array.
1219///
1220/// Use the two pass algorithm, which is slower (! a factor of 2) but much more
1221/// precise. Since we have a vector the 2 pass algorithm is still faster than the
1222/// Welford algorithm. (See also ROOT-5545)
1223
1224template <typename Iterator>
1225Double_t TMath::RMS(Iterator first, Iterator last)
1226{
1227
1228 Double_t n = 0;
1229
1230 Double_t tot = 0;
1231 Double_t mean = TMath::Mean(first,last);
1232 while ( first != last ) {
1233 Double_t x = Double_t(*first);
1234 tot += (x - mean)*(x - mean);
1235 ++first;
1236 ++n;
1237 }
1238 Double_t rms = (n > 1) ? TMath::Sqrt(tot/(n-1)) : 0.0;
1239 return rms;
1240}
1241
1242////////////////////////////////////////////////////////////////////////////////
1243/// Returns the weighted Standard Deviation of an array defined by the iterators.
1244/// Note that this function returns the sigma(standard deviation) and
1245/// not the root mean square of the array.
1246///
1247/// As in the unweighted case use the two pass algorithm
1248
1249template <typename Iterator, typename WeightIterator>
1250Double_t TMath::RMS(Iterator first, Iterator last, WeightIterator w)
1251{
1252 Double_t tot = 0;
1253 Double_t sumw = 0;
1254 Double_t sumw2 = 0;
1255 Double_t mean = TMath::Mean(first,last,w);
1256 while ( first != last ) {
1257 Double_t x = Double_t(*first);
1258 sumw += *w;
1259 sumw2 += (*w) * (*w);
1260 tot += (*w) * (x - mean)*(x - mean);
1261 ++first;
1262 ++w;
1263 }
1264 // use the correction neff/(neff -1) for the unbiased formula
1266 return rms;
1267}
1268
1269////////////////////////////////////////////////////////////////////////////////
1270/// Returns the Standard Deviation of an array a with length n.
1271/// Note that this function returns the sigma(standard deviation) and
1272/// not the root mean square of the array.
1273
1274template <typename T>
1276{
1277 return (w) ? TMath::RMS(a, a+n, w) : TMath::RMS(a, a+n);
1278}
1279
1280////////////////////////////////////////////////////////////////////////////////
1281/// Calculates the Cross Product of two vectors:
1282/// out = [v1 x v2]
1283
1284template <typename T> T *TMath::Cross(const T v1[3],const T v2[3], T out[3])
1285{
1286 out[0] = v1[1] * v2[2] - v1[2] * v2[1];
1287 out[1] = v1[2] * v2[0] - v1[0] * v2[2];
1288 out[2] = v1[0] * v2[1] - v1[1] * v2[0];
1289
1290 return out;
1291}
1292
1293////////////////////////////////////////////////////////////////////////////////
1294/// Calculates a normal vector of a plane.
1295///
1296/// \param[in] p1, p2,p3 3 3D points belonged the plane to define it.
1297/// \param[out] normal Pointer to 3D normal vector (normalized)
1298
1299template <typename T> T * TMath::Normal2Plane(const T p1[3],const T p2[3],const T p3[3], T normal[3])
1300{
1301 T v1[3], v2[3];
1302
1303 v1[0] = p2[0] - p1[0];
1304 v1[1] = p2[1] - p1[1];
1305 v1[2] = p2[2] - p1[2];
1306
1307 v2[0] = p3[0] - p1[0];
1308 v2[1] = p3[1] - p1[1];
1309 v2[2] = p3[2] - p1[2];
1310
1312 return normal;
1313}
1314
1315////////////////////////////////////////////////////////////////////////////////
1316/// Function which returns kTRUE if point xp,yp lies inside the
1317/// polygon defined by the np points in arrays x and y, kFALSE otherwise.
1318/// Note that the polygon may be open or closed.
1319
1320template <typename T> Bool_t TMath::IsInside(T xp, T yp, Int_t np, T *x, T *y)
1321{
1322 Int_t i, j = np-1 ;
1324
1325 for (i=0; i<np; i++) {
1326 if ((y[i]<yp && y[j]>=yp) || (y[j]<yp && y[i]>=yp)) {
1327 if (x[i]+(yp-y[i])/(y[j]-y[i])*(x[j]-x[i])<xp) {
1328 oddNodes = !oddNodes;
1329 }
1330 }
1331 j=i;
1332 }
1333
1334 return oddNodes;
1335}
1336
1337////////////////////////////////////////////////////////////////////////////////
1338/// Returns the median of the array a where each entry i has weight w[i] .
1339/// Both arrays have a length of at least n . The median is a number obtained
1340/// from the sorted array a through
1341///
1342/// median = (a[jl]+a[jh])/2. where (using also the sorted index on the array w)
1343///
1344/// sum_i=0,jl w[i] <= sumTot/2
1345/// sum_i=0,jh w[i] >= sumTot/2
1346/// sumTot = sum_i=0,n w[i]
1347///
1348/// If w=0, the algorithm defaults to the median definition where it is
1349/// a number that divides the sorted sequence into 2 halves.
1350/// When n is odd or n > 1000, the median is kth element k = (n + 1) / 2.
1351/// when n is even and n < 1000the median is a mean of the elements k = n/2 and k = n/2 + 1.
1352///
1353/// If the weights are supplied (w not 0) all weights must be >= 0
1354///
1355/// If work is supplied, it is used to store the sorting index and assumed to be
1356/// >= n . If work=0, local storage is used, either on the stack if n < kWorkMax
1357/// or on the heap for n >= kWorkMax .
1358
1359template <typename T> Double_t TMath::Median(Long64_t n, const T *a, const Double_t *w, Long64_t *work)
1360{
1361
1362 const Int_t kWorkMax = 100;
1363
1364 if (n <= 0 || !a) return 0;
1367 Long64_t *ind;
1368 Long64_t workLocal[kWorkMax];
1369
1370 if (work) {
1371 ind = work;
1372 } else {
1373 ind = workLocal;
1374 if (n > kWorkMax) {
1376 ind = new Long64_t[n];
1377 }
1378 }
1379
1380 if (w) {
1381 Double_t sumTot2 = 0;
1382 for (Int_t j = 0; j < n; j++) {
1383 if (w[j] < 0) {
1384 ::Error("TMath::Median","w[%d] = %.4e < 0 ?!",j,w[j]);
1385 if (isAllocated) delete [] ind;
1386 return 0;
1387 }
1388 sumTot2 += w[j];
1389 }
1390
1391 sumTot2 /= 2.;
1392
1393 Sort(n, a, ind, kFALSE);
1394
1395 Double_t sum = 0.;
1396 Int_t jl;
1397 for (jl = 0; jl < n; jl++) {
1398 sum += w[ind[jl]];
1399 if (sum >= sumTot2) break;
1400 }
1401
1402 Int_t jh;
1403 sum = 2.*sumTot2;
1404 for (jh = n-1; jh >= 0; jh--) {
1405 sum -= w[ind[jh]];
1406 if (sum <= sumTot2) break;
1407 }
1408
1409 median = 0.5*(a[ind[jl]]+a[ind[jh]]);
1410
1411 } else {
1412
1413 if (n%2 == 1)
1414 median = KOrdStat(n, a,n/2, ind);
1415 else {
1416 median = 0.5*(KOrdStat(n, a, n/2 -1, ind)+KOrdStat(n, a, n/2, ind));
1417 }
1418 }
1419
1420 if (isAllocated)
1421 delete [] ind;
1422 return median;
1423}
1424
1425////////////////////////////////////////////////////////////////////////////////
1426/// Returns the half-sample mode of the array a where each entry i has weight w[i].
1427/// Both arrays must have a length of n. The mode is a number obtained
1428/// according to the paper http://arxiv.org/ftp/math/papers/0505/0505419.pdf (page 19).
1429/// The algorithm searches for the smallest range of sorted a values that
1430/// contains ~n/2 values. Then it re-applies the same algorithm to
1431/// that range until the range has <= 3 elements in it.
1432///
1433/// \note See David R. Bickel: "On a Fast, Robust Estimator of the Mode"
1434/// http://arxiv.org/ftp/math/papers/0505/0505419.pdf (page 19)
1435/// Initial implementation suggestion by Jean-Fran\c{c}ois Caron, 2014 on JIRA-6849
1436///
1437/// \tparam T data type of the values array `a`
1438/// \param n Number of elements in array `a`
1439/// \param a Pointer to array of values (owned by user, preallocated), must have length of n, and none of them should be NaN. The array must not be necessarily sorted nor contain unique values.
1440/// \param w Pointer to array of weights (owned by user, preallocated), must have length of n. If the weights are supplied (w not nullptr) all weights must be >= 0 to properly work.
1441/// \return NaN if n <= 0, a[0] if n == 1, TMath::Mean(n,a,w) if n == 2,
1442/// the (weighted) mean of the closest two if n == 3; for the rest of the cases,
1443/// if all values of `a` are unique or have equal weights (or `w` is nullptr), then the Bickel half-sample-mode algorithm is applied; otherwise the
1444/// value associated with the highest sum of weights (if `w` is not nullptr)
1445/// is returned.
1446/// \note If in a given iteration step, there are several equally small subranges, then the first occurrence is chosen.
1447template <typename T> Double_t TMath::ModeHalfSample(Long64_t n, const T *a, const Double_t *w)
1448{
1449 if (n <= 0)
1450 return TMath::QuietNaN();
1451 else if (n == 1)
1452 return a[0];
1453 else if (n == 2)
1454 return TMath::Mean(n, a, w);
1455
1456 if (w && std::adjacent_find(w, w + n, std::not_equal_to<>()) == w + n) // If all weights are equal, ignore those
1457 w = nullptr;
1458
1459 // Sort the array and remove duplicates (if w != nullptr)
1460 std::vector<T> values;
1461 std::vector<Double_t> weights;
1462 values.reserve(n);
1463 weights.reserve(n);
1464 for (Long64_t i = 0; i < n; ++i) {
1465 const T value = a[i];
1466 const Double_t weight = w ? w[i] : 1;
1467 const auto vstart = values.begin();
1468 const auto vstop = values.end();
1469 const auto viter = std::lower_bound(vstart, vstop, value);
1470 const auto vidx = std::distance(vstart, viter);
1471 if (viter != vstop && w) {
1472 weights[vidx] += weight;
1473 } else {
1474 values.insert(viter, value);
1475 weights.insert(weights.begin() + vidx, weight);
1476 }
1477 }
1478
1479 const size_t sn = values.size();
1480 if (sn <= 0)
1481 return TMath::QuietNaN();
1482 else if (sn == 1)
1483 return values[0];
1484 else if (sn == 2)
1485 return TMath::Mean(sn, values.data(), weights.data());
1486
1487 if (sn == 3) {
1488 const double d1_0 = values[1] - values[0];
1489 const double d2_1 = values[2] - values[1];
1490 if (d2_1 < d1_0)
1491 return TMath::Mean(2, values.data() + 1, weights.data() + 1);
1492 else if (d2_1 > d1_0)
1493 return TMath::Mean(2, values.data(), weights.data());
1494 else
1495 return values[1];
1496 }
1497
1498 const auto wstart = weights.begin();
1499 const auto wstop = weights.end();
1500 if (std::adjacent_find(wstart, wstop, std::not_equal_to<>()) != wstop) {
1501 // Get highest (sum of) weight element
1502 const auto wmaxiter = std::max_element(wstart, wstop);
1503 const auto wmaxidx = std::distance(wstart, wmaxiter);
1504 return values[wmaxidx];
1505 }
1506
1507 // All elements have equal weights (and are unique if w != nullptr)
1508
1509 // Initialize search
1510 n = sn;
1511 size_t jMin = 0;
1512
1513 // Do recursive calls dividing each time the interval by two
1514 while (n > 3) {
1515 Double_t min_v_range = values[jMin + n - 1] - values[jMin];
1516 const size_t N = std::ceil(n * 0.5);
1517 const size_t start = jMin;
1518 const size_t stop = start + n - N + 1; // +1 since we use < and not <=
1519 // Find sequentally what v_range is smallest by sliding the half-window
1520 for (size_t i = start; i < stop; i++) {
1521 Double_t range = values[i + N - 1] - values[i];
1522 if (range < min_v_range) {
1524 jMin = i;
1525 }
1526 }
1527 // assert(min_v_range == values[N - 1 + start] - values[start]);
1528 n = N;
1529 }
1530 if (n == 3) {
1531 const double d1_0 = values[jMin + 1] - values[jMin + 0];
1532 const double d2_1 = values[jMin + 2] - values[jMin + 1];
1533 if (d2_1 < d1_0)
1534 return (values[jMin + 1] + values[jMin + 2]) * 0.5;
1535 else if (d2_1 > d1_0)
1536 return (values[jMin] + values[jMin + 1]) * 0.5;
1537 else
1538 return values[jMin + 1];
1539 } else if (n == 2) {
1540 return (values[jMin] + values[jMin + 1]) * 0.5;
1541 } else {
1542 Error("ModeHalfSample", "Error in recursive algorithm, returning NaN"); // this should not happen
1543 return TMath::QuietNaN();
1544 }
1545}
1546
1547////////////////////////////////////////////////////////////////////////////////
1548/// Returns k_th order statistic of the array a of size n
1549/// (k_th smallest element out of n elements).
1550///
1551/// C-convention is used for array indexing, so if you want
1552/// the second smallest element, call KOrdStat(n, a, 1).
1553///
1554/// If work is supplied, it is used to store the sorting index and
1555/// assumed to be >= n. If work=0, local storage is used, either on
1556/// the stack if n < kWorkMax or on the heap for n >= kWorkMax.
1557/// Note that the work index array will not contain the sorted indices but
1558/// all indices of the smaller element in arbitrary order in work[0,...,k-1] and
1559/// all indices of the larger element in arbitrary order in work[k+1,..,n-1]
1560/// work[k] will contain instead the index of the returned element.
1561///
1562/// Taken from "Numerical Recipes in C++" without the index array
1563/// implemented by Anna Khreshuk.
1564///
1565/// See also the declarations at the top of this file
1566
1567template <class Element, typename Size>
1568Element TMath::KOrdStat(Size n, const Element *a, Size k, Size *work)
1569{
1570
1571 const Int_t kWorkMax = 100;
1572
1573 typedef Size Index;
1574
1576 Size i, ir, j, l, mid;
1577 Index arr;
1578 Index *ind;
1579 Index workLocal[kWorkMax];
1580 Index temp;
1581
1582 if (work) {
1583 ind = work;
1584 } else {
1585 ind = workLocal;
1586 if (n > kWorkMax) {
1588 ind = new Index[n];
1589 }
1590 }
1591
1592 for (Size ii=0; ii<n; ii++) {
1593 ind[ii]=ii;
1594 }
1595 Size rk = k;
1596 l=0;
1597 ir = n-1;
1598 for(;;) {
1599 if (ir<=l+1) { //active partition contains 1 or 2 elements
1600 if (ir == l+1 && a[ind[ir]]<a[ind[l]])
1601 {temp = ind[l]; ind[l]=ind[ir]; ind[ir]=temp;}
1602 Element tmp = a[ind[rk]];
1603 if (isAllocated)
1604 delete [] ind;
1605 return tmp;
1606 } else {
1607 mid = (l+ir) >> 1; //choose median of left, center and right
1608 {temp = ind[mid]; ind[mid]=ind[l+1]; ind[l+1]=temp;}//elements as partitioning element arr.
1609 if (a[ind[l]]>a[ind[ir]]) //also rearrange so that a[l]<=a[l+1]
1610 {temp = ind[l]; ind[l]=ind[ir]; ind[ir]=temp;}
1611
1612 if (a[ind[l+1]]>a[ind[ir]])
1613 {temp=ind[l+1]; ind[l+1]=ind[ir]; ind[ir]=temp;}
1614
1615 if (a[ind[l]]>a[ind[l+1]])
1616 {temp = ind[l]; ind[l]=ind[l+1]; ind[l+1]=temp;}
1617
1618 i=l+1; //initialize pointers for partitioning
1619 j=ir;
1620 arr = ind[l+1];
1621 for (;;){
1622 do i++; while (a[ind[i]]<a[arr]);
1623 do j--; while (a[ind[j]]>a[arr]);
1624 if (j<i) break; //pointers crossed, partitioning complete
1625 {temp=ind[i]; ind[i]=ind[j]; ind[j]=temp;}
1626 }
1627 ind[l+1]=ind[j];
1628 ind[j]=arr;
1629 if (j>=rk) ir = j-1; //keep active the partition that
1630 if (j<=rk) l=i; //contains the k_th element
1631 }
1632 }
1633}
1634
1635#endif
#define IsNaN(a)
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define h(i)
Definition RSha256.hxx:106
bool Bool_t
Boolean (0=false, 1=true) (bool)
Definition RtypesCore.h:77
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
unsigned long ULong_t
Unsigned long integer 4 bytes (unsigned long). Size depends on architecture.
Definition RtypesCore.h:69
long Long_t
Signed long integer 4 bytes (long). Size depends on architecture.
Definition RtypesCore.h:68
unsigned int UInt_t
Unsigned integer 4 bytes (unsigned int)
Definition RtypesCore.h:60
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:71
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
long double LongDouble_t
Long Double (not portable)
Definition RtypesCore.h:75
long long Long64_t
Portable signed long integer 8 bytes.
Definition RtypesCore.h:83
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:208
#define N
TGLVector3 Cross(const TGLVector3 &v1, const TGLVector3 &v2)
Definition TGLUtil.h:323
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t mask
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
void Sort()
float xmin
float * q
float xmax
const Double_t sigma
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Namespace for new Math classes and functions.
TMath.
Definition TMathBase.h:35
Double_t FDistI(Double_t F, Double_t N, Double_t M)
Calculates the cumulative distribution function of F-distribution (see ROOT::Math::fdistribution_cdf)...
Definition TMath.cxx:2300
Double_t GeomMean(Long64_t n, const T *a)
Returns the geometric mean of an array a of size n.
Definition TMath.h:1210
Double_t LogNormal(Double_t x, Double_t sigma, Double_t theta=0, Double_t m=1)
Computes the density of LogNormal distribution at point x.
Definition TMath.cxx:2440
T * Laplacian(Long64_t n, T *f, double h=1)
Calculate the Laplacian of an array with length n.
Definition TMath.h:1067
constexpr Double_t G()
Gravitational constant in: .
Definition TMath.h:138
Double_t CosH(Double_t)
Returns the hyperbolic cosine of x.
Definition TMath.h:623
T * Normal2Plane(const T v1[3], const T v2[3], const T v3[3], T normal[3])
Calculates a normal vector of a plane.
Definition TMath.h:1299
Double_t DiLog(Double_t x)
Modified Struve functions of order 1.
Definition TMath.cxx:116
Double_t BetaDist(Double_t x, Double_t p, Double_t q)
Computes the probability density function of the Beta distribution (the cumulative distribution funct...
Definition TMath.cxx:2071
void KNNDensity(std::span< const double > observations, std::span< const double > queries, std::span< double > result, int k, double dmin=0.0)
Computes a 1D k-nearest neighbor density estimate for a set of query points.
Definition TMath.cxx:3227
Double_t ACos(Double_t)
Returns the principal value of the arc cosine of x, expressed in radians.
Definition TMath.h:643
Double_t BesselI(Int_t n, Double_t x)
Computes the Integer Order Modified Bessel function I_n(x) for n=0,1,2,... and any real x.
Definition TMath.cxx:1590
Element KOrdStat(Size n, const Element *a, Size k, Size *work=0)
Returns k_th order statistic of the array a of size n (k_th smallest element out of n elements).
Definition TMath.h:1568
T * Gradient(Long64_t n, T *f, double h=1)
Calculate the one-dimensional gradient of an array with length n.
Definition TMath.h:1029
Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE)
Calculates a gaussian function with mean and sigma.
Definition TMath.cxx:471
Bool_t IsNaN(Double_t x)
Definition TMath.h:903
constexpr Double_t GUncertainty()
Gravitational constant uncertainty.
Definition TMath.h:153
constexpr Double_t C()
Velocity of light in .
Definition TMath.h:117
Double_t Factorial(Int_t i)
Computes factorial(n).
Definition TMath.cxx:252
Double_t RMS(Long64_t n, const T *a, const Double_t *w=nullptr)
Returns the Standard Deviation of an array a with length n.
Definition TMath.h:1275
Double_t KolmogorovTest(Int_t na, const Double_t *a, Int_t nb, const Double_t *b, Option_t *option)
Statistical test whether two one-dimensional sets of points are compatible with coming from the same ...
Definition TMath.cxx:805
constexpr Double_t GhbarCUncertainty()
uncertainty.
Definition TMath.h:169
Long64_t LocMin(Long64_t n, const T *a)
Returns index of array with the minimum element.
Definition TMath.h:993
constexpr Double_t Ccgs()
Definition TMath.h:124
constexpr Double_t SigmaUncertainty()
Stefan-Boltzmann constant uncertainty.
Definition TMath.h:280
Int_t Nint(T x)
Round to nearest integer. Rounds half integers to the nearest even integer.
Definition TMath.h:704
Double_t ModeHalfSample(Long64_t n, const T *a, const Double_t *w=nullptr)
Returns the half-sample mode of the array a where each entry i has weight w[i].
Definition TMath.h:1447
Double_t BinomialI(Double_t p, Int_t n, Int_t k)
Suppose an event occurs with probability p per trial Then the probability P of its occurring k or mor...
Definition TMath.cxx:2144
Double_t Vavilov(Double_t x, Double_t kappa, Double_t beta2)
Returns the value of the Vavilov probability density function.
Definition TMath.cxx:2780
Double_t Binomial(Int_t n, Int_t k)
Calculates the binomial coefficient n over k.
Definition TMath.cxx:2112
Float_t Normalize(Float_t v[3])
Normalize a vector v in place.
Definition TMath.cxx:518
constexpr Double_t NaUncertainty()
Avogadro constant (Avogadro's Number) uncertainty.
Definition TMath.h:294
Double_t Prob(Double_t chi2, Int_t ndf)
Computation of the probability for a certain Chi-squared (chi2) and number of degrees of freedom (ndf...
Definition TMath.cxx:637
Bool_t IsInside(T xp, T yp, Int_t np, T *x, T *y)
Function which returns kTRUE if point xp,yp lies inside the polygon defined by the np points in array...
Definition TMath.h:1320
Double_t ASin(Double_t)
Returns the principal value of the arc sine of x, expressed in radians.
Definition TMath.h:635
Double_t Log2(Double_t x)
Returns the binary (base-2) logarithm of x.
Definition TMath.cxx:107
Double_t BesselK1(Double_t x)
Modified Bessel function I_1(x)
Definition TMath.cxx:1529
Double_t Median(Long64_t n, const T *a, const Double_t *w=nullptr, Long64_t *work=nullptr)
Same as RMS.
Definition TMath.h:1359
void BubbleHigh(Int_t Narr, Double_t *arr1, Int_t *arr2)
Bubble sort variant to obtain the order of an array's elements into an index in order to do more usef...
Definition TMath.cxx:1314
Double_t Exp(Double_t x)
Returns the base-e exponential function of x, which is e raised to the power x.
Definition TMath.h:720
Double_t BesselI1(Double_t x)
Modified Bessel function K_0(x)
Definition TMath.cxx:1494
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition TMath.cxx:190
Bool_t Permute(Int_t n, Int_t *a)
Simple recursive algorithm to find the permutations of n natural numbers, not necessarily all distinc...
Definition TMath.cxx:2560
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754.
Definition TMath.h:913
Double_t Floor(Double_t x)
Rounds x downward, returning the largest integral value that is not greater than x.
Definition TMath.h:691
Double_t PoissonI(Double_t x, Double_t par)
Computes the Discrete Poisson distribution function for (x,par).
Definition TMath.cxx:615
Double_t CauchyDist(Double_t x, Double_t t=0, Double_t s=1)
Computes the density of Cauchy distribution at point x by default, standard Cauchy distribution is us...
Definition TMath.cxx:2178
Double_t ATan(Double_t)
Returns the principal value of the arc tangent of x, expressed in radians.
Definition TMath.h:651
Double_t StruveL1(Double_t x)
Modified Struve functions of order 0.
Definition TMath.cxx:1970
constexpr Double_t Gn()
Standard acceleration of gravity in .
Definition TMath.h:177
Double_t ASinH(Double_t)
Returns the area hyperbolic sine of x.
Definition TMath.cxx:67
Double_t LaplaceDistI(Double_t x, Double_t alpha=0, Double_t beta=1)
Computes the cumulative distribution function (lower tail integral) of Laplace distribution at point ...
Definition TMath.cxx:2383
ULong_t Hash(const void *txt, Int_t ntxt)
Calculates hash index from any char string.
Definition TMath.cxx:1408
constexpr Double_t QeUncertainty()
Elementary charge uncertainty.
Definition TMath.h:349
constexpr Double_t K()
Boltzmann's constant in : .
Definition TMath.h:250
Double_t BreitWigner(Double_t x, Double_t mean=0, Double_t gamma=1)
Calculates a Breit Wigner function with mean and gamma.
Definition TMath.cxx:442
constexpr Double_t Sqrt2()
Definition TMath.h:89
constexpr Double_t KUncertainty()
Boltzmann's constant uncertainty.
Definition TMath.h:264
constexpr Double_t Hbarcgs()
Definition TMath.h:221
Double_t Landau(Double_t x, Double_t mpv=0, Double_t sigma=1, Bool_t norm=kFALSE)
The LANDAU function.
Definition TMath.cxx:492
Double_t Voigt(Double_t x, Double_t sigma, Double_t lg, Int_t r=4)
Computation of Voigt function (normalised).
Definition TMath.cxx:898
Double_t Student(Double_t T, Double_t ndf)
Computes density function for Student's t- distribution.
Definition TMath.cxx:2625
constexpr Double_t CUncertainty()
Speed of light uncertainty.
Definition TMath.h:131
constexpr Double_t Qe()
Elementary charge in .
Definition TMath.h:342
Double_t Ceil(Double_t x)
Rounds x upward, returning the smallest integral value that is not less than x.
Definition TMath.h:679
constexpr Double_t PiOver2()
Definition TMath.h:54
constexpr Double_t HCcgs()
Definition TMath.h:243
T MinElement(Long64_t n, const T *a)
Returns minimum of array a of length n.
Definition TMath.h:971
Double_t BetaDistI(Double_t x, Double_t p, Double_t q)
Computes the cumulative distribution function of the Beta distribution, i.e.
Definition TMath.cxx:2090
T NormCross(const T v1[3], const T v2[3], T out[3])
Calculates the Normalized Cross Product of two vectors.
Definition TMath.h:962
Int_t Finite(Double_t x)
Check if it is finite with a mask in order to be consistent in presence of fast math.
Definition TMath.h:781
Double_t TanH(Double_t)
Returns the hyperbolic tangent of x.
Definition TMath.h:629
Int_t FloorNint(Double_t x)
Returns the nearest integer of TMath::Floor(x).
Definition TMath.h:697
Double_t ACosH(Double_t)
Returns the nonnegative area hyperbolic cosine of x.
Definition TMath.cxx:81
Double_t BesselK0(Double_t x)
Modified Bessel function I_0(x)
Definition TMath.cxx:1460
Double_t BesselY0(Double_t x)
Bessel function J1(x) for any real x.
Definition TMath.cxx:1705
Double_t ATan2(Double_t y, Double_t x)
Returns the principal value of the arc tangent of y/x, expressed in radians.
Definition TMath.h:657
constexpr Double_t RUncertainty()
Universal gas constant uncertainty.
Definition TMath.h:312
Double_t BetaCf(Double_t x, Double_t a, Double_t b)
Continued fraction evaluation by modified Lentz's method used in calculation of incomplete Beta funct...
Definition TMath.cxx:2020
Long64_t LocMax(Long64_t n, const T *a)
Returns index of array with the maximum element.
Definition TMath.h:1099
Double_t ErfInverse(Double_t x)
Returns the inverse error function.
Definition TMath.cxx:208
Double_t LaplaceDist(Double_t x, Double_t alpha=0, Double_t beta=1)
Computes the probability density function of Laplace distribution at point x, with location parameter...
Definition TMath.cxx:2367
constexpr Double_t E()
Base of natural log: .
Definition TMath.h:96
constexpr Double_t GnUncertainty()
Standard acceleration of gravity uncertainty.
Definition TMath.h:184
constexpr Double_t Hcgs()
Definition TMath.h:198
constexpr Double_t HUncertainty()
Planck's constant uncertainty.
Definition TMath.h:205
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Definition TMath.h:767
constexpr Double_t DegToRad()
Conversion from degree to radian: .
Definition TMath.h:82
Double_t Erfc(Double_t x)
Computes the complementary error function erfc(x).
Definition TMath.cxx:199
Double_t VavilovI(Double_t x, Double_t kappa, Double_t beta2)
Returns the value of the Vavilov cumulative distribution function (lower tail integral of the probabi...
Definition TMath.cxx:2817
constexpr Double_t Sigma()
Stefan-Boltzmann constant in : .
Definition TMath.h:273
Double_t Beta(Double_t p, Double_t q)
Calculates Beta-function Gamma(p)*Gamma(q)/Gamma(p+q).
Definition TMath.cxx:2011
constexpr Double_t Kcgs()
Definition TMath.h:257
Double_t Sq(Double_t x)
Returns x*x.
Definition TMath.h:667
Double_t Poisson(Double_t x, Double_t par)
Computes the Poisson distribution function for (x,par).
Definition TMath.cxx:587
constexpr Double_t H()
Planck's constant in : .
Definition TMath.h:191
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:732
Int_t CeilNint(Double_t x)
Returns the nearest integer of TMath::Ceil(x).
Definition TMath.h:685
Double_t Ldexp(Double_t x, Int_t exp)
Returns the result of multiplying x (the significant) by 2 raised to the power of exp (the exponent).
Definition TMath.h:726
Double_t BesselJ0(Double_t x)
Modified Bessel function K_1(x)
Definition TMath.cxx:1634
constexpr Double_t LogE()
Base-10 log of e (to convert ln to log)
Definition TMath.h:110
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition TMath.cxx:353
constexpr Double_t MWair()
Molecular weight of dry air 1976 US Standard Atmosphere in or
Definition TMath.h:320
constexpr Double_t Gcgs()
Definition TMath.h:146
Double_t StruveL0(Double_t x)
Struve functions of order 1.
Definition TMath.cxx:1923
Double_t NormQuantile(Double_t p)
Computes quantiles for standard normal distribution N(0, 1) at probability p.
Definition TMath.cxx:2459
constexpr Double_t Ln10()
Natural log of 10 (to convert log to ln)
Definition TMath.h:103
Double_t Hypot(Double_t x, Double_t y)
Returns sqrt(x*x + y*y)
Definition TMath.cxx:59
constexpr Double_t EulerGamma()
Euler-Mascheroni Constant.
Definition TMath.h:335
constexpr Double_t PiOver4()
Definition TMath.h:61
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:605
void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted=kTRUE, Int_t *index=nullptr, Int_t type=7)
Computes sample quantiles, corresponding to the given probabilities.
Definition TMath.cxx:1207
constexpr Double_t Pi()
Definition TMath.h:40
Double_t StruveH0(Double_t x)
Bessel function Y1(x) for positive x.
Definition TMath.cxx:1777
constexpr Double_t R()
Universal gas constant ( ) in
Definition TMath.h:305
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition TMath.cxx:509
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Comparing floating points.
Definition TMath.h:429
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Comparing floating points.
Definition TMath.h:421
Double_t Mean(Long64_t n, const T *a, const Double_t *w=nullptr)
Returns the weighted mean of an array a with length n.
Definition TMath.h:1176
Double_t KolmogorovProb(Double_t z)
Calculates the Kolmogorov distribution function,.
Definition TMath.cxx:679
constexpr Double_t InvPi()
Definition TMath.h:68
Bool_t RootsCubic(const Double_t coef[4], Double_t &a, Double_t &b, Double_t &c)
Calculates roots of polynomial of 3rd order a*x^3 + b*x^2 + c*x + d, where.
Definition TMath.cxx:1107
Double_t ChisquareQuantile(Double_t p, Double_t ndf)
Evaluate the quantiles of the chi-squared probability distribution function.
Definition TMath.cxx:2196
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:599
Double_t FDist(Double_t F, Double_t N, Double_t M)
Computes the density function of F-distribution (probability function, integral of density,...
Definition TMath.cxx:2280
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN).
Definition TMath.h:921
Double_t BreitWignerRelativistic(Double_t x, Double_t median=0, Double_t gamma=1)
Calculates a Relativistic Breit Wigner function with median and gamma.
Definition TMath.cxx:452
T * Cross(const T v1[3], const T v2[3], T out[3])
Calculates the Cross Product of two vectors: out = [v1 x v2].
Definition TMath.h:1284
void BubbleLow(Int_t Narr, Double_t *arr1, Int_t *arr2)
Opposite ordering of the array arr2[] to that of BubbleHigh.
Definition TMath.cxx:1353
Double_t BesselK(Int_t n, Double_t x)
Integer order modified Bessel function I_n(x)
Definition TMath.cxx:1561
constexpr Double_t Na()
Avogadro constant (Avogadro's Number) in .
Definition TMath.h:287
T MaxElement(Long64_t n, const T *a)
Returns maximum of array a of length n.
Definition TMath.h:979
Double_t BesselJ1(Double_t x)
Bessel function J0(x) for any real x.
Definition TMath.cxx:1669
Double_t BetaIncomplete(Double_t x, Double_t a, Double_t b)
Calculates the incomplete Beta-function.
Definition TMath.cxx:2103
constexpr Double_t Rgair()
Dry Air Gas Constant (R / MWair) in
Definition TMath.h:328
constexpr Double_t Hbar()
in :
Definition TMath.h:214
Double_t StruveH1(Double_t x)
Struve functions of order 0.
Definition TMath.cxx:1846
Double_t Freq(Double_t x)
Computation of the normal frequency function freq(x).
Definition TMath.cxx:270
Double_t Tan(Double_t)
Returns the tangent of an angle of x radians.
Definition TMath.h:611
Double_t LandauI(Double_t x)
Returns the cumulative (lower tail integral) of the Landau distribution function at point x.
Definition TMath.cxx:2847
Double_t StdDev(Long64_t n, const T *a, const Double_t *w=nullptr)
Definition TMath.h:534
Double_t ATanH(Double_t)
Returns the area hyperbolic tangent of x.
Definition TMath.cxx:95
constexpr Double_t RadToDeg()
Conversion from radian to degree: .
Definition TMath.h:75
Double_t BesselI0(Double_t x)
Integer order modified Bessel function K_n(x)
Definition TMath.cxx:1426
Double_t Log10(Double_t x)
Returns the common (base-10) logarithm of x.
Definition TMath.h:773
Double_t StudentI(Double_t T, Double_t ndf)
Calculates the cumulative distribution function of Student's t-distribution second parameter stands f...
Definition TMath.cxx:2648
Double_t StudentQuantile(Double_t p, Double_t ndf, Bool_t lower_tail=kTRUE)
Computes quantiles of the Student's t-distribution 1st argument is the probability,...
Definition TMath.cxx:2676
Double_t BesselY1(Double_t x)
Bessel function Y0(x) for positive x.
Definition TMath.cxx:1739
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124
Double_t GammaDist(Double_t x, Double_t gamma, Double_t mu=0, Double_t beta=1)
Computes the density function of Gamma distribution at point x.
Definition TMath.cxx:2350
constexpr Double_t GhbarC()
in
Definition TMath.h:161
constexpr Double_t HC()
in
Definition TMath.h:236
constexpr Double_t TwoPi()
Definition TMath.h:47
constexpr Double_t HbarUncertainty()
uncertainty.
Definition TMath.h:228
Double_t Infinity()
Returns an infinity as defined by the IEEE standard.
Definition TMath.h:928
Double_t ErfcInverse(Double_t x)
Returns the inverse of the complementary error function.
Definition TMath.cxx:242
Double_t SinH(Double_t)
Returns the hyperbolic sine of `x.
Definition TMath.h:617
static T Min()
Returns maximum representation for type T.
Definition TMath.h:936
static T Epsilon()
Returns minimum double representation.
Definition TMath.h:952
static T Max()
Returns minimum double representation.
Definition TMath.h:944
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2339