Revision 716fb731a3c8c9b583235feddac21006994d9626 authored by Kyle Cranmer on 03 August 2011, 22:09:36 UTC, committed by Kyle Cranmer on 03 August 2011, 22:09:36 UTC
git-svn-id: http://root.cern.ch/svn/root/branches/v5-30-00-patches@40446 27541ba8-7e3a-0410-8455-c3a389f83636
1 parent 0db9ccb
RooMath.cxx
/*****************************************************************************
* Project: RooFit *
* Package: RooFitCore *
* @(#)root/roofitcore:$Id$
* Authors: *
* WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
* DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
* *
* Copyright (c) 2000-2005, Regents of the University of California *
* and Stanford University. All rights reserved. *
* *
* Redistribution and use in source and binary forms, *
* with or without modification, are permitted according to the terms *
* listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
*****************************************************************************/
// -- CLASS DESCRIPTION [MISC] --
// RooMath is a singleton class implementing various mathematical
// functions not found in TMath, mostly involving complex algebra
//
//
#include "RooFit.h"
#include "RooMath.h"
#include "RooMath.h"
#include "TMath.h"
#include <math.h>
#include "Riostream.h"
#include "RooMsgService.h"
#include "RooSentinel.h"
#include "TSystem.h"
ClassImp(RooMath)
;
RooComplex RooMath::FastComplexErrFunc(const RooComplex& z)
{
return ITPComplexErrFunc(z,z.im()>0?3:4) ;
}
Double_t RooMath::FastComplexErrFuncRe(const RooComplex& z)
{
return ITPComplexErrFuncRe(z,z.im()>0?3:4) ;
}
Double_t RooMath::FastComplexErrFuncIm(const RooComplex& z)
{
return ITPComplexErrFuncIm(z,z.im()>0?3:4) ;
}
void RooMath::cacheCERF(Bool_t flag)
{
_cacheTable = flag ;
}
RooComplex RooMath::ComplexErrFunc(Double_t re, Double_t im) {
// Return CERNlib complex error function for Z(re,im)
return ComplexErrFunc(RooComplex(re,im));
}
RooComplex RooMath::ComplexErrFunc(const RooComplex& Z) {
// Return CERNlib complex error function
//
// This code is translated from the fortran version in the CERN mathlib.
// (see ftp://asisftp.cern.ch/cernlib/share/pro/src/mathlib/gen/c/cwerf64.F)
RooComplex ZH,S,T,V;
static RooComplex R[38];
static const Double_t Z1= 1, HF= Z1/2, Z10= 10;
static const Double_t C1= 74/Z10, C2= 83/Z10, C3= Z10/32, C4 = 16/Z10;
static const Double_t C = 1.12837916709551257, P = pow(2*C4,33);
static const RooComplex zero(0);
Double_t X(Z.re()),Y(Z.im()), XA(fabs(X)), YA(fabs(Y));
int N;
if((YA < C1) && (XA < C2)) {
ZH= RooComplex(YA+C4,XA);
R[37]=zero;
N= 36;
while(N > 0) {
T=ZH+R[N+1].conj()*N;
R[N--]=(T*HF)/T.abs2();
}
Double_t XL=P;
S=zero;
N= 33;
while(N > 0) {
XL=C3*XL;
S=R[N--]*(S+XL);
}
V=S*C;
}
else {
ZH=RooComplex(YA,XA);
R[1]=zero;
N= 9;
while(N > 0) {
T=ZH+R[1].conj()*N;
R[1]=(T*HF)/T.abs2();
N--;
}
V=R[1]*C;
}
if(YA==0) V=RooComplex(exp(-(XA*XA)),V.im());
if(Y < 0) {
RooComplex tmp(XA,YA);
tmp= -tmp*tmp;
V=tmp.exp()*2-V;
if(X > 0) V= V.conj();
}
else {
if(X < 0) V= V.conj();
}
return V;
}
void RooMath::initFastCERF(Int_t reBins, Double_t reMin, Double_t reMax, Int_t imBins, Double_t imMin, Double_t imMax)
{
// Allocate and initialize lookup table for interpolated complex error function
// for given grid parameters
// Store grid dimensions and related parameters
_reBins = reBins ;
_imBins = imBins ;
_reMin = reMin ;
_reMax = reMax ;
_reRange = _reMax-_reMin ;
_reStep = _reRange/_reBins ;
_imMin = imMin ;
_imMax = imMax ;
_imRange = _imMax-_imMin ;
_imStep = _imRange/_imBins ;
oocxcoutD((TObject*)0,Eval) << endl
<< "RooMath::initFastCERF: Allocating Complex Error Function lookup table" << endl
<< " Re: " << _reBins << " bins in range (" << _reMin << "," << _reMax << ")" << endl
<< " Im: " << _imBins << " bins in range (" << _imMin << "," << _imMax << ")" << endl
<< " Allocation size : " << _reBins*_imBins * 2 * sizeof(Double_t) / 1024 << " kB" << endl ;
// Allocate storage matrix for Im(cerf) and Re(cerf) and fill it using ComplexErrFunc()
RooSentinel::activate() ;
Int_t imIdx,reIdx ;
_reCerfArray = new pDouble_t[_imBins] ;
_imCerfArray = new pDouble_t[_imBins] ;
for (imIdx=0 ; imIdx<_imBins ; imIdx++) {
_reCerfArray[imIdx] = new Double_t[_reBins] ;
_imCerfArray[imIdx] = new Double_t[_reBins] ;
}
Bool_t cacheLoaded(kFALSE) ;
if (!_cacheTable || !(cacheLoaded=loadCache())) {
oocxcoutD((TObject*)0,Eval) << endl
<< " Filling table: |..................................................|\r"
<< " Filling table: |" ;
// Allocate storage matrix for Im(cerf) and Re(cerf) and fill it using ComplexErrFunc()
for (imIdx=0 ; imIdx<_imBins ; imIdx++) {
if (imIdx % (_imBins/50) ==0) {
ooccxcoutD((TObject*)0,Eval) << ">" ; cout.flush() ;
}
for (reIdx=0 ; reIdx<_reBins ; reIdx++) {
RooComplex val=ComplexErrFunc(_reMin+reIdx*_reStep,_imMin+imIdx*_imStep) ;
_reCerfArray[imIdx][reIdx] = val.re();
_imCerfArray[imIdx][reIdx] = val.im() ;
}
}
ooccoutI((TObject*)0,Eval) << endl ;
}
if (_cacheTable && !cacheLoaded) storeCache() ;
}
void RooMath::cleanup()
{
if (_reCerfArray) {
for (Int_t imIdx=0 ; imIdx<_imBins ; imIdx++) {
delete[] _reCerfArray[imIdx] ;
delete[] _imCerfArray[imIdx] ;
}
delete[] _reCerfArray ;
delete[] _imCerfArray ;
_reCerfArray = 0 ;
_imCerfArray = 0 ;
}
}
RooComplex RooMath::ITPComplexErrFunc(const RooComplex& z, Int_t nOrder)
{
// Return complex error function interpolated from lookup tabel created
// by initFastCERF(). Interpolation is performed in Im and Re plane
// to specified order.
// Initialize grid
if (!_reCerfArray) initFastCERF() ;
// Located nearest grid point
Double_t imPrime = (z.im()-_imMin) / _imStep ;
Double_t rePrime = (z.re()-_reMin) / _reStep ;
// Calculate corners of nOrder X nOrder grid box
Int_t imIdxLo = Int_t(imPrime - 1.0*nOrder/2 + 0.5) ;
Int_t imIdxHi = imIdxLo+nOrder-1 ;
Int_t reIdxLo = Int_t(rePrime - 1.0*nOrder/2 + 0.5) ;
Int_t reIdxHi = reIdxLo+nOrder-1 ;
// Check if the box is fully contained in the grid
if (imIdxLo<0 || imIdxHi>_imBins-1 || reIdxLo<0 || reIdxHi>_reBins-1) {
return ComplexErrFunc(z) ;
}
// Allocate temporary array space for interpolation
Int_t imIdx /*, reIdx*/ ;
Double_t imYR[10] ;
Double_t imYI[10] ;
// Loop over imaginary grid points
for (imIdx=imIdxLo ; imIdx<=imIdxHi ; imIdx++) {
// Interpolate real array and store as array point for imaginary interpolation
imYR[imIdx-imIdxLo] = interpolate(&_reCerfArray[imIdx][reIdxLo],nOrder,rePrime-reIdxLo) ;
imYI[imIdx-imIdxLo] = interpolate(&_imCerfArray[imIdx][reIdxLo],nOrder,rePrime-reIdxLo) ;
}
// Interpolate imaginary arrays and construct complex return value
Double_t re = interpolate(imYR,nOrder,imPrime-imIdxLo) ;
Double_t im = interpolate(imYI,nOrder,imPrime-imIdxLo) ;
return RooComplex(re,im) ;
}
Double_t RooMath::ITPComplexErrFuncRe(const RooComplex& z, Int_t nOrder)
{
// Return real component of complex error function interpolated from
// lookup table created by initFastCERF(). Interpolation is performed in
// Im and Re plane to specified order. This functions is noticably faster
// than ITPComplexErrrFunc().re() because only the real lookup table
// is interpolated
// Initialize grid
if (!_reCerfArray) initFastCERF() ;
// Located nearest grid point
Double_t imPrime = (z.im()-_imMin) / _imStep ;
Double_t rePrime = (z.re()-_reMin) / _reStep ;
// Calculate corners of nOrder X nOrder grid box
Int_t imIdxLo = Int_t(imPrime - 1.0*nOrder/2 + 0.5) ;
Int_t imIdxHi = imIdxLo+nOrder-1 ;
Int_t reIdxLo = Int_t(rePrime - 1.0*nOrder/2 + 0.5) ;
Int_t reIdxHi = reIdxLo+nOrder-1 ;
// Check if the box is fully contained in the grid
if (imIdxLo<0 || imIdxHi>_imBins-1 || reIdxLo<0 || reIdxHi>_reBins-1) {
//cout << "RooMath::ITPComplexErrFuncRe: (" << z.re() << "," << z.im() << ") outside interpolation grid" << endl ;
return ComplexErrFunc(z).re() ;
}
if (nOrder==1) return _reCerfArray[imIdxLo][reIdxLo] ;
Int_t imIdx ;
Double_t imYR[10] ;
// Allocate temporary array space for interpolation
for (imIdx=imIdxLo ; imIdx<=imIdxHi ; imIdx++) {
// Interpolate real array and store as array point for imaginary interpolation
imYR[imIdx-imIdxLo] = interpolate(&_reCerfArray[imIdx][reIdxLo],nOrder,rePrime-reIdxLo) ;
}
// Interpolate imaginary arrays and construct complex return value
return interpolate(imYR,nOrder,imPrime-imIdxLo) ;
}
Double_t RooMath::ITPComplexErrFuncIm(const RooComplex& z, Int_t nOrder)
{
// Return real component of complex error function interpolated from
// lookup table created by initFastCERF(). Interpolation is performed in
// Im and Re plane to specified order. This functions is noticably faster
// than ITPComplexErrrFunc().im() because only the imaginary lookup table
// is interpolated
// Initialize grid
if (!_reCerfArray) initFastCERF() ;
// Located nearest grid point
Double_t imPrime = (z.im()-_imMin) / _imStep ;
Double_t rePrime = (z.re()-_reMin) / _reStep ;
// Calculate corners of nOrder X nOrder grid box
Int_t imIdxLo = Int_t(imPrime - 1.0*nOrder/2 + 0.5) ;
Int_t imIdxHi = imIdxLo+nOrder-1 ;
Int_t reIdxLo = Int_t(rePrime - 1.0*nOrder/2 + 0.5) ;
Int_t reIdxHi = reIdxLo+nOrder-1 ;
// Check if the box is fully contained in the grid
if (imIdxLo<0 || imIdxHi>_imBins-1 || reIdxLo<0 || reIdxHi>_reBins-1) {
return ComplexErrFunc(z).im() ;
}
// Allocate temporary array space for interpolation
Int_t imIdx ;
Double_t imYI[10] ;
// Loop over imaginary grid points
for (imIdx=imIdxLo ; imIdx<=imIdxHi ; imIdx++) {
// Interpolate real array and store as array point for imaginary interpolation
imYI[imIdx-imIdxLo] = interpolate(&_imCerfArray[imIdx][reIdxLo],nOrder,rePrime-reIdxLo) ;
}
// Interpolate imaginary arrays and construct complex return value
return interpolate(imYI,nOrder,imPrime-imIdxLo) ;
}
Double_t RooMath::interpolate(Double_t ya[], Int_t n, Double_t x)
{
// Interpolate array 'ya' with 'n' elements for 'x' (between 0 and 'n'-1)
// Int to Double conversion is faster via array lookup than type conversion!
static Double_t itod[20] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0} ;
int i,m,ns=1 ;
Double_t den,dif,dift/*,ho,hp,w*/,y,dy ;
Double_t c[20], d[20] ;
dif = fabs(x) ;
for(i=1 ; i<=n ; i++) {
if ((dift=fabs(x-itod[i-1]))<dif) {
ns=i ;
dif=dift ;
}
c[i] = ya[i-1] ;
d[i] = ya[i-1] ;
}
y=ya[--ns] ;
for(m=1 ; m<n; m++) {
for(i=1 ; i<=n-m ; i++) {
den=(c[i+1] - d[i])/itod[m] ;
d[i]=(x-itod[i+m-1])*den ;
c[i]=(x-itod[i-1])*den ;
}
dy = (2*ns)<(n-m) ? c[ns+1] : d[ns--] ;
y += dy ;
}
return y ;
}
Double_t RooMath::interpolate(Double_t xa[], Double_t ya[], Int_t n, Double_t x)
{
// Interpolate array 'ya' with 'n' elements for 'xa'
// Int to Double conversion is faster via array lookup than type conversion!
// static Double_t itod[20] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
// 10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0} ;
int i,m,ns=1 ;
Double_t den,dif,dift,ho,hp,w,y,dy ;
Double_t c[20], d[20] ;
dif = fabs(x-xa[0]) ;
for(i=1 ; i<=n ; i++) {
if ((dift=fabs(x-xa[i-1]))<dif) {
ns=i ;
dif=dift ;
}
c[i] = ya[i-1] ;
d[i] = ya[i-1] ;
}
y=ya[--ns] ;
for(m=1 ; m<n; m++) {
for(i=1 ; i<=n-m ; i++) {
ho=xa[i-1]-x ;
hp=xa[i-1+m]-x ;
w=c[i+1]-d[i] ;
den=ho-hp ;
if (den==0.) {
oocoutE((TObject*)0,Eval) << "RooMath::interpolate ERROR: zero distance between points not allowed" << endl ;
return 0 ;
}
den = w/den ;
d[i]=hp*den ;
c[i]=ho*den;
}
dy = (2*ns)<(n-m) ? c[ns+1] : d[ns--] ;
y += dy ;
}
return y ;
}
Bool_t RooMath::loadCache()
{
// Load the complex error function lookup table from the cache file
const char* fName = cacheFileName() ;
// Open cache file
ifstream ifs(fName) ;
// Return immediately if file doesn't exist
if (ifs.fail()) return kFALSE ;
oocxcoutD((TObject*)0,Eval) << endl << " Filling table from cache file " << fName << endl ;
// Load data in memory arrays
Bool_t ok(kTRUE) ;
Int_t i ;
for (i=0 ; i<_imBins ; i++) {
ifs.read((char*)_imCerfArray[i],_reBins*sizeof(Double_t)) ;
if (ifs.fail()) ok=kFALSE ;
ifs.read((char*)_reCerfArray[i],_reBins*sizeof(Double_t)) ;
if (ifs.fail()) ok=kFALSE ;
}
// Issue error message on partial read failure
if (!ok) {
oocoutE((TObject*)0,Eval) << "RooMath::loadCERFCache: error reading file " << cacheFileName() << endl ;
}
return ok ;
}
void RooMath::storeCache()
{
// Store the complex error function lookup table in the cache file
ofstream ofs(cacheFileName()) ;
oocxcoutI((TObject*)0,Eval) << endl << " Writing table to cache file " << cacheFileName() << endl ;
Int_t i ;
for (i=0 ; i<_imBins ; i++) {
ofs.write((char*)_imCerfArray[i],_reBins*sizeof(Double_t)) ;
ofs.write((char*)_reCerfArray[i],_reBins*sizeof(Double_t)) ;
}
}
const char* RooMath::cacheFileName()
{
// Construct and return the name of the complex error function cache file
static char fileName[1024] ;
snprintf(fileName,1024,"/tmp/RooMath_CERFcache_R%04d_I%04d_%d.dat",_reBins,_imBins,gSystem->GetUid()) ;
return fileName ;
}
Double_t RooMath::erf(Double_t x)
{
return TMath::Erf(x) ;
}
Double_t RooMath::erfc(Double_t x)
{
return TMath::Erfc(x) ;
}
// Instantiation of static members
Int_t RooMath::_reBins(0) ;
Int_t RooMath::_imBins(0) ;
Double_t RooMath::_reMin(0) ;
Double_t RooMath::_reMax(0) ;
Double_t RooMath::_reRange(0) ;
Double_t RooMath::_reStep(0) ;
Double_t RooMath::_imMin(0) ;
Double_t RooMath::_imMax(0) ;
Double_t RooMath::_imRange(0) ;
Double_t RooMath::_imStep(0) ;
pDouble_t* RooMath::_reCerfArray = 0;
pDouble_t* RooMath::_imCerfArray = 0;
Bool_t RooMath::_cacheTable(kTRUE) ;
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...