Logo ROOT  
Reference Guide
 
Loading...
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 *
10 * and Stanford University. All rights reserved. *
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
30std::complex<double> RooMath::faddeeva(std::complex<double> z)
31{
33}
34
35std::complex<double> RooMath::faddeeva_fast(std::complex<double> z)
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, m, ns = 1;
86 double den, dif, dift /*,ho,hp,w*/, y, dy;
87 double c[20], d[20];
88
89 dif = std::abs(x);
90 for (i = 1; i <= n; i++) {
91 if ((dift = std::abs(x - itod[i - 1])) < dif) {
92 ns = i;
93 dif = dift;
94 }
95 c[i] = ya[i - 1];
96 d[i] = ya[i - 1];
97 }
98
99 y = ya[--ns];
100 for (m = 1; m < n; m++) {
101 for (i = 1; i <= n - m; i++) {
102 den = (c[i + 1] - d[i]) / itod[m];
103 d[i] = (x - itod[i + m - 1]) * den;
104 c[i] = (x - itod[i - 1]) * den;
105 }
106 dy = (2 * ns) < (n - m) ? c[ns + 1] : d[ns--];
107 y += dy;
108 }
109 return y;
110}
111
112double RooMath::interpolate(double xa[], double ya[], Int_t n, double x)
113{
114 // Interpolate array 'ya' with 'n' elements for 'xa'
115
116 // Int to Double conversion is faster via array lookup than type conversion!
117 // static double itod[20] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
118 // 10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0} ;
119 int i, m, ns = 1;
120 double den, dif, dift, ho, hp, w, y, dy;
121 double c[20], d[20];
122
123 dif = std::abs(x - xa[0]);
124 for (i = 1; i <= n; i++) {
125 if ((dift = std::abs(x - xa[i - 1])) < dif) {
126 ns = i;
127 dif = dift;
128 }
129 c[i] = ya[i - 1];
130 d[i] = ya[i - 1];
131 }
132
133 y = ya[--ns];
134 for (m = 1; m < n; m++) {
135 for (i = 1; i <= n - m; i++) {
136 ho = xa[i - 1] - x;
137 hp = xa[i - 1 + m] - x;
138 w = c[i + 1] - d[i];
139 den = ho - hp;
140 if (den == 0.) {
141 std::cerr << "In " << __func__ << "(), " << __FILE__ << ":" << __LINE__
142 << ": Zero distance between points not allowed." << std::endl;
143 return 0;
144 }
145 den = w / den;
146 d[i] = hp * den;
147 c[i] = ho * den;
148 }
149 dy = (2 * ns) < (n - m) ? c[ns + 1] : d[ns--];
150 y += dy;
151 }
152 return y;
153}
#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