ROOT   Reference Guide
Searching...
No Matches
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 *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17// -- CLASS DESCRIPTION [MISC] --
18// RooMath is a singleton class implementing various mathematical
19// functions not found in TMath, mostly involving complex algebra
20
21#include <RooMath.h>
22
24
25#include <algorithm>
26#include <cmath>
27#include <complex>
28#include <iostream>
29
31{
33}
34
36{
38}
39
40std::complex<double> RooMath::erfc(const std::complex<double> z)
41{
42 double re = -z.real() * z.real() + z.imag() * z.imag();
43 double im = -2. * z.real() * z.imag();
45 return (z.real() >= 0.) ? (std::complex<double>(re, im) * faddeeva(std::complex<double>(-z.imag(), z.real())))
46 : (2. - std::complex<double>(re, im) * faddeeva(std::complex<double>(z.imag(), -z.real())));
47}
48
49std::complex<double> RooMath::erfc_fast(const std::complex<double> z)
50{
51 double re = -z.real() * z.real() + z.imag() * z.imag();
52 double im = -2. * z.real() * z.imag();
54 return (z.real() >= 0.)
55 ? (std::complex<double>(re, im) * faddeeva_fast(std::complex<double>(-z.imag(), z.real())))
56 : (2. - std::complex<double>(re, im) * faddeeva_fast(std::complex<double>(z.imag(), -z.real())));
57}
58
59std::complex<double> RooMath::erf(const std::complex<double> z)
60{
61 double re = -z.real() * z.real() + z.imag() * z.imag();
62 double im = -2. * z.real() * z.imag();
64 return (z.real() >= 0.) ? (1. - std::complex<double>(re, im) * faddeeva(std::complex<double>(-z.imag(), z.real())))
65 : (std::complex<double>(re, im) * faddeeva(std::complex<double>(z.imag(), -z.real())) - 1.);
66}
67
68std::complex<double> RooMath::erf_fast(const std::complex<double> z)
69{
70 double re = -z.real() * z.real() + z.imag() * z.imag();
71 double im = -2. * z.real() * z.imag();
73 return (z.real() >= 0.)
74 ? (1. - std::complex<double>(re, im) * faddeeva_fast(std::complex<double>(-z.imag(), z.real())))
75 : (std::complex<double>(re, im) * faddeeva_fast(std::complex<double>(z.imag(), -z.real())) - 1.);
76}
77
78double RooMath::interpolate(double ya[], Int_t n, double x)
79{
80 // Interpolate array 'ya' with 'n' elements for 'x' (between 0 and 'n'-1)
81
82 // Int to Double conversion is faster via array lookup than type conversion!
83 static double itod[20] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
84 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0};
85 int i;
86 int m;
87 int ns = 1;
88 double den;
89 double dif;
90 double dift /*,ho,hp,w*/;
91 double y;
92 double dy;
93 double c[20];
94 double d[20];
95
96 dif = std::abs(x);
97 for (i = 1; i <= n; i++) {
98 if ((dift = std::abs(x - itod[i - 1])) < dif) {
99 ns = i;
100 dif = dift;
101 }
102 c[i] = ya[i - 1];
103 d[i] = ya[i - 1];
104 }
105
106 y = ya[--ns];
107 for (m = 1; m < n; m++) {
108 for (i = 1; i <= n - m; i++) {
109 den = (c[i + 1] - d[i]) / itod[m];
110 d[i] = (x - itod[i + m - 1]) * den;
111 c[i] = (x - itod[i - 1]) * den;
112 }
113 dy = (2 * ns) < (n - m) ? c[ns + 1] : d[ns--];
114 y += dy;
115 }
116 return y;
117}
118
119double RooMath::interpolate(double xa[], double ya[], Int_t n, double x)
120{
121 // Interpolate array 'ya' with 'n' elements for 'xa'
122
123 // Int to Double conversion is faster via array lookup than type conversion!
124 // static double itod[20] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
125 // 10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0} ;
126 int i;
127 int m;
128 int ns = 1;
129 double den;
130 double dif;
131 double dift;
132 double ho;
133 double hp;
134 double w;
135 double y;
136 double dy;
137 double c[20];
138 double d[20];
139
140 dif = std::abs(x - xa[0]);
141 for (i = 1; i <= n; i++) {
142 if ((dift = std::abs(x - xa[i - 1])) < dif) {
143 ns = i;
144 dif = dift;
145 }
146 c[i] = ya[i - 1];
147 d[i] = ya[i - 1];
148 }
149
150 y = ya[--ns];
151 for (m = 1; m < n; m++) {
152 for (i = 1; i <= n - m; i++) {
153 ho = xa[i - 1] - x;
154 hp = xa[i - 1 + m] - x;
155 w = c[i + 1] - d[i];
156 den = ho - hp;
157 if (den == 0.) {
158 std::cerr << "In " << __func__ << "(), " << __FILE__ << ":" << __LINE__
159 << ": Zero distance between points not allowed." << std::endl;
160 return 0;
161 }
162 den = w / den;
163 d[i] = hp * den;
164 c[i] = ho * den;
165 }
166 dy = (2 * ns) < (n - m) ? c[ns + 1] : d[ns--];
167 y += dy;
168 }
169 return y;
170}
#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:40
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition RooMath.cxx:59
static std::complex< double > faddeeva(std::complex< double > z)
evaluate Faddeeva function for complex argument
Definition RooMath.cxx:30
static double interpolate(double yArr[], Int_t nOrder, double x)
Definition RooMath.cxx:78
static std::complex< double > faddeeva_fast(std::complex< double > z)
evaluate Faddeeva function for complex argument (fast version)
Definition RooMath.cxx:35
static std::complex< double > erfc_fast(const std::complex< double > z)
complex erfc function (fast version)
Definition RooMath.cxx:49
static std::complex< double > erf_fast(const std::complex< double > z)
complex erf function (fast version)
Definition RooMath.cxx:68
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
__roodevice__ static __roohost__ void cexp(double &re, double &im)
__roodevice__ __roohost__ STD::complex< double > faddeeva(STD::complex< double > z)
__roodevice__ __roohost__ STD::complex< double > faddeeva_fast(STD::complex< double > z)
TMarker m
Definition textangle.C:8