ROOT   Reference Guide
RooMath.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
15 *****************************************************************************/
16
17// -- CLASS DESCRIPTION [MISC] --
18// RooMath is a singleton class implementing various mathematical
20
21
22#include <RooMath.h>
23
25
26#include <algorithm>
27#include <cmath>
28#include <complex>
29#include <iostream>
30
32{
34}
35
37{
39}
40
41std::complex<double> RooMath::erfc(const std::complex<double> z)
42{
43 double re = -z.real() * z.real() + z.imag() * z.imag();
44 double im = -2. * z.real() * z.imag();
46 return (z.real() >= 0.) ? (std::complex<double>(re, im) * faddeeva(std::complex<double>(-z.imag(), z.real())))
47 : (2. - std::complex<double>(re, im) * faddeeva(std::complex<double>(z.imag(), -z.real())));
48}
49
50std::complex<double> RooMath::erfc_fast(const std::complex<double> z)
51{
52 double re = -z.real() * z.real() + z.imag() * z.imag();
53 double im = -2. * z.real() * z.imag();
55 return (z.real() >= 0.)
56 ? (std::complex<double>(re, im) * faddeeva_fast(std::complex<double>(-z.imag(), z.real())))
57 : (2. - std::complex<double>(re, im) * faddeeva_fast(std::complex<double>(z.imag(), -z.real())));
58}
59
60std::complex<double> RooMath::erf(const std::complex<double> z)
61{
62 double re = -z.real() * z.real() + z.imag() * z.imag();
63 double im = -2. * z.real() * z.imag();
65 return (z.real() >= 0.) ? (1. - std::complex<double>(re, im) * faddeeva(std::complex<double>(-z.imag(), z.real())))
66 : (std::complex<double>(re, im) * faddeeva(std::complex<double>(z.imag(), -z.real())) - 1.);
67}
68
69std::complex<double> RooMath::erf_fast(const std::complex<double> z)
70{
71 double re = -z.real() * z.real() + z.imag() * z.imag();
72 double im = -2. * z.real() * z.imag();
74 return (z.real() >= 0.)
75 ? (1. - std::complex<double>(re, im) * faddeeva_fast(std::complex<double>(-z.imag(), z.real())))
76 : (std::complex<double>(re, im) * faddeeva_fast(std::complex<double>(z.imag(), -z.real())) - 1.);
77}
78
79double RooMath::interpolate(double ya[], Int_t n, double x)
80{
81 // Interpolate array 'ya' with 'n' elements for 'x' (between 0 and 'n'-1)
82
83 // Int to Double conversion is faster via array lookup than type conversion!
84 static double itod[20] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
85 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0};
86 int i, m, ns = 1;
87 double den, dif, dift /*,ho,hp,w*/, y, dy;
88 double c[20], d[20];
89
90 dif = fabs(x);
91 for (i = 1; i <= n; i++) {
92 if ((dift = fabs(x - itod[i - 1])) < dif) {
93 ns = i;
94 dif = dift;
95 }
96 c[i] = ya[i - 1];
97 d[i] = ya[i - 1];
98 }
99
100 y = ya[--ns];
101 for (m = 1; m < n; m++) {
102 for (i = 1; i <= n - m; i++) {
103 den = (c[i + 1] - d[i]) / itod[m];
104 d[i] = (x - itod[i + m - 1]) * den;
105 c[i] = (x - itod[i - 1]) * den;
106 }
107 dy = (2 * ns) < (n - m) ? c[ns + 1] : d[ns--];
108 y += dy;
109 }
110 return y;
111}
112
113double RooMath::interpolate(double xa[], double ya[], Int_t n, double x)
114{
115 // Interpolate array 'ya' with 'n' elements for 'xa'
116
117 // Int to Double conversion is faster via array lookup than type conversion!
118 // static double itod[20] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
119 // 10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0} ;
120 int i, m, ns = 1;
121 double den, dif, dift, ho, hp, w, y, dy;
122 double c[20], d[20];
123
124 dif = fabs(x - xa[0]);
125 for (i = 1; i <= n; i++) {
126 if ((dift = fabs(x - xa[i - 1])) < dif) {
127 ns = i;
128 dif = dift;
129 }
130 c[i] = ya[i - 1];
131 d[i] = ya[i - 1];
132 }
133
134 y = ya[--ns];
135 for (m = 1; m < n; m++) {
136 for (i = 1; i <= n - m; i++) {
137 ho = xa[i - 1] - x;
138 hp = xa[i - 1 + m] - x;
139 w = c[i + 1] - d[i];
140 den = ho - hp;
141 if (den == 0.) {
142 std::cerr << "In " << __func__ << "(), " << __FILE__ << ":" << __LINE__
143 << ": Zero distance between points not allowed." << std::endl;
144 return 0;
145 }
146 den = w / den;
147 d[i] = hp * den;
148 c[i] = ho * den;
149 }
150 dy = (2 * ns) < (n - m) ? c[ns + 1] : d[ns--];
151 y += dy;
152 }
153 return y;
154}
#define d(i)
Definition: RSha256.hxx:102
#define c(i)
Definition: RSha256.hxx:101
static std::complex< double > erfc(const std::complex< double > z)
complex erfc function
Definition: RooMath.cxx:41
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:60
static std::complex< double > faddeeva(std::complex< double > z)
evaluate Faddeeva function for complex argument
Definition: RooMath.cxx:31
static double interpolate(double yArr[], Int_t nOrder, double x)
Definition: RooMath.cxx:79
static std::complex< double > faddeeva_fast(std::complex< double > z)
evaluate Faddeeva function for complex argument (fast version)
Definition: RooMath.cxx:36
static std::complex< double > erfc_fast(const std::complex< double > z)
complex erfc function (fast version)
Definition: RooMath.cxx:50
static std::complex< double > erf_fast(const std::complex< double > z)
complex erf function (fast version)
Definition: RooMath.cxx:69
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
static constexpr double ns
__roodevice__ __roohost__ std::complex< double > faddeeva(std::complex< double > z)
__roodevice__ __roohost__ std::complex< double > faddeeva_fast(std::complex< double > z)
__roodevice__ static __roohost__ void cexp(double &re, double &im)