As always forgot some files...
This commit is contained in:
Bastian M. Wojek
parent
ccaba04d51
commit
68092aad6b
242
src/external/BMWIntegrator/BMWIntegrator.cpp
vendored
Normal file
242
src/external/BMWIntegrator/BMWIntegrator.cpp
vendored
Normal file
@ -0,0 +1,242 @@
|
|||||||
|
/***************************************************************************
|
||||||
|
|
||||||
|
BMWIntegrator.cpp
|
||||||
|
|
||||||
|
Author: Bastian M. Wojek
|
||||||
|
e-mail: bastian.wojek@psi.ch
|
||||||
|
|
||||||
|
$Id$
|
||||||
|
|
||||||
|
***************************************************************************/
|
||||||
|
|
||||||
|
/***************************************************************************
|
||||||
|
* Copyright (C) 2009 by Bastian M. Wojek *
|
||||||
|
* bastian.wojek@psi.ch *
|
||||||
|
* *
|
||||||
|
* This program is free software; you can redistribute it and/or modify *
|
||||||
|
* it under the terms of the GNU General Public License as published by *
|
||||||
|
* the Free Software Foundation; either version 2 of the License, or *
|
||||||
|
* (at your option) any later version. *
|
||||||
|
* *
|
||||||
|
* This program is distributed in the hope that it will be useful, *
|
||||||
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||||
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||||
|
* GNU General Public License for more details. *
|
||||||
|
* *
|
||||||
|
* You should have received a copy of the GNU General Public License *
|
||||||
|
* along with this program; if not, write to the *
|
||||||
|
* Free Software Foundation, Inc., *
|
||||||
|
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
|
||||||
|
***************************************************************************/
|
||||||
|
|
||||||
|
# include "BMWIntegrator.h"
|
||||||
|
|
||||||
|
# include "cuba.h"
|
||||||
|
|
||||||
|
std::vector<double> TDWaveGapIntegralCuhre::fPar;
|
||||||
|
|
||||||
|
double TDWaveGapIntegralCuhre::IntegrateFunc()
|
||||||
|
{
|
||||||
|
const unsigned int NCOMP(1);
|
||||||
|
const double EPSREL (1e-4);
|
||||||
|
const double EPSABS (1e-6);
|
||||||
|
const unsigned int VERBOSE (0);
|
||||||
|
const unsigned int LAST (4);
|
||||||
|
const unsigned int MINEVAL (0);
|
||||||
|
const unsigned int MAXEVAL (50000);
|
||||||
|
|
||||||
|
const unsigned int KEY (13);
|
||||||
|
|
||||||
|
int nregions, neval, fail;
|
||||||
|
double integral[NCOMP], error[NCOMP], prob[NCOMP];
|
||||||
|
|
||||||
|
Cuhre(fNDim, NCOMP, Integrand,
|
||||||
|
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
|
||||||
|
KEY,
|
||||||
|
&nregions, &neval, &fail, integral, error, prob);
|
||||||
|
|
||||||
|
return integral[0];
|
||||||
|
}
|
||||||
|
|
||||||
|
void TDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
|
||||||
|
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T), Ec, phic}
|
||||||
|
{
|
||||||
|
double deltasq(TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]*fPar[3]),2.0));
|
||||||
|
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<double> TAnSWaveGapIntegralCuhre::fPar;
|
||||||
|
|
||||||
|
double TAnSWaveGapIntegralCuhre::IntegrateFunc()
|
||||||
|
{
|
||||||
|
const unsigned int NCOMP(1);
|
||||||
|
const double EPSREL (1e-4);
|
||||||
|
const double EPSABS (1e-6);
|
||||||
|
const unsigned int VERBOSE (0);
|
||||||
|
const unsigned int LAST (4);
|
||||||
|
const unsigned int MINEVAL (100);
|
||||||
|
const unsigned int MAXEVAL (1000000);
|
||||||
|
|
||||||
|
const unsigned int KEY (13);
|
||||||
|
|
||||||
|
int nregions, neval, fail;
|
||||||
|
double integral[NCOMP], error[NCOMP], prob[NCOMP];
|
||||||
|
|
||||||
|
Cuhre(fNDim, NCOMP, Integrand,
|
||||||
|
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
|
||||||
|
KEY,
|
||||||
|
&nregions, &neval, &fail, integral, error, prob);
|
||||||
|
|
||||||
|
return integral[0];
|
||||||
|
}
|
||||||
|
|
||||||
|
void TAnSWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
|
||||||
|
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
|
||||||
|
{
|
||||||
|
double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0));
|
||||||
|
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<double> TAnSWaveGapIntegralDivonne::fPar;
|
||||||
|
|
||||||
|
double TAnSWaveGapIntegralDivonne::IntegrateFunc()
|
||||||
|
{
|
||||||
|
const unsigned int NCOMP(1);
|
||||||
|
const double EPSREL (1e-4);
|
||||||
|
const double EPSABS (1e-6);
|
||||||
|
const unsigned int VERBOSE (0);
|
||||||
|
const unsigned int MINEVAL (1000);
|
||||||
|
const unsigned int MAXEVAL (1000000);
|
||||||
|
const unsigned int KEY1 (47);
|
||||||
|
const unsigned int KEY2 (1);
|
||||||
|
const unsigned int KEY3 (1);
|
||||||
|
const unsigned int MAXPASS (5);
|
||||||
|
const double BORDER (0.);
|
||||||
|
const double MAXCHISQ (10.);
|
||||||
|
const double MINDEVIATION (.25);
|
||||||
|
const unsigned int NGIVEN (0);
|
||||||
|
const unsigned int LDXGIVEN (fNDim);
|
||||||
|
const unsigned int NEXTRA (0);
|
||||||
|
|
||||||
|
int nregions, neval, fail;
|
||||||
|
double integral[NCOMP], error[NCOMP], prob[NCOMP];
|
||||||
|
|
||||||
|
Divonne(fNDim, NCOMP, Integrand,
|
||||||
|
EPSREL, EPSABS, VERBOSE, MINEVAL, MAXEVAL,
|
||||||
|
KEY1, KEY2, KEY3, MAXPASS, BORDER, MAXCHISQ, MINDEVIATION,
|
||||||
|
NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL,
|
||||||
|
&nregions, &neval, &fail, integral, error, prob);
|
||||||
|
|
||||||
|
return integral[0];
|
||||||
|
}
|
||||||
|
|
||||||
|
void TAnSWaveGapIntegralDivonne::Integrand(const int *ndim, const double x[],
|
||||||
|
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
|
||||||
|
{
|
||||||
|
double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0));
|
||||||
|
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<double> TAnSWaveGapIntegralSuave::fPar;
|
||||||
|
|
||||||
|
double TAnSWaveGapIntegralSuave::IntegrateFunc()
|
||||||
|
{
|
||||||
|
const unsigned int NCOMP(1);
|
||||||
|
const double EPSREL (1e-4);
|
||||||
|
const double EPSABS (1e-6);
|
||||||
|
const unsigned int VERBOSE (0);
|
||||||
|
const unsigned int LAST (4);
|
||||||
|
const unsigned int MINEVAL (1000);
|
||||||
|
const unsigned int MAXEVAL (1000000);
|
||||||
|
|
||||||
|
const unsigned int NNEW (1000);
|
||||||
|
const double FLATNESS (25.);
|
||||||
|
|
||||||
|
int nregions, neval, fail;
|
||||||
|
double integral[NCOMP], error[NCOMP], prob[NCOMP];
|
||||||
|
|
||||||
|
Suave(fNDim, NCOMP, Integrand,
|
||||||
|
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
|
||||||
|
NNEW, FLATNESS,
|
||||||
|
&nregions, &neval, &fail, integral, error, prob);
|
||||||
|
|
||||||
|
return integral[0];
|
||||||
|
}
|
||||||
|
|
||||||
|
void TAnSWaveGapIntegralSuave::Integrand(const int *ndim, const double x[],
|
||||||
|
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
|
||||||
|
{
|
||||||
|
double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0));
|
||||||
|
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<double> TNonMonDWave1GapIntegralCuhre::fPar;
|
||||||
|
|
||||||
|
double TNonMonDWave1GapIntegralCuhre::IntegrateFunc()
|
||||||
|
{
|
||||||
|
const unsigned int NCOMP(1);
|
||||||
|
const double EPSREL (1e-4);
|
||||||
|
const double EPSABS (1e-6);
|
||||||
|
const unsigned int VERBOSE (0);
|
||||||
|
const unsigned int LAST (4);
|
||||||
|
const unsigned int MINEVAL (100);
|
||||||
|
const unsigned int MAXEVAL (1000000);
|
||||||
|
|
||||||
|
const unsigned int KEY (13);
|
||||||
|
|
||||||
|
int nregions, neval, fail;
|
||||||
|
double integral[NCOMP], error[NCOMP], prob[NCOMP];
|
||||||
|
|
||||||
|
Cuhre(fNDim, NCOMP, Integrand,
|
||||||
|
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
|
||||||
|
KEY,
|
||||||
|
&nregions, &neval, &fail, integral, error, prob);
|
||||||
|
|
||||||
|
return integral[0];
|
||||||
|
}
|
||||||
|
|
||||||
|
void TNonMonDWave1GapIntegralCuhre::Integrand(const int *ndim, const double x[],
|
||||||
|
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
|
||||||
|
{
|
||||||
|
double deltasq(TMath::Power(fPar[1]*(fPar[2]*TMath::Cos(2.0*x[1]*fPar[4])+(1.0-fPar[2])*TMath::Cos(6.0*x[1]*fPar[4])),2.0));
|
||||||
|
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<double> TNonMonDWave2GapIntegralCuhre::fPar;
|
||||||
|
|
||||||
|
double TNonMonDWave2GapIntegralCuhre::IntegrateFunc()
|
||||||
|
{
|
||||||
|
const unsigned int NCOMP(1);
|
||||||
|
const double EPSREL (1e-4);
|
||||||
|
const double EPSABS (1e-6);
|
||||||
|
const unsigned int VERBOSE (0);
|
||||||
|
const unsigned int LAST (4);
|
||||||
|
const unsigned int MINEVAL (100);
|
||||||
|
const unsigned int MAXEVAL (1000000);
|
||||||
|
|
||||||
|
const unsigned int KEY (13);
|
||||||
|
|
||||||
|
int nregions, neval, fail;
|
||||||
|
double integral[NCOMP], error[NCOMP], prob[NCOMP];
|
||||||
|
|
||||||
|
Cuhre(fNDim, NCOMP, Integrand,
|
||||||
|
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
|
||||||
|
KEY,
|
||||||
|
&nregions, &neval, &fail, integral, error, prob);
|
||||||
|
|
||||||
|
return integral[0];
|
||||||
|
}
|
||||||
|
|
||||||
|
void TNonMonDWave2GapIntegralCuhre::Integrand(const int *ndim, const double x[],
|
||||||
|
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
|
||||||
|
{
|
||||||
|
double deltasq(4.0*fPar[2]/27.0*TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]*fPar[4]), 2.0) \
|
||||||
|
/ TMath::Power(1.0 + fPar[2]*TMath::Cos(2.0*x[1]*fPar[4])*TMath::Cos(2.0*x[1]*fPar[4]), 3.0));
|
||||||
|
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
|
||||||
|
return;
|
||||||
|
}
|
||||||
329
src/external/BMWIntegrator/BMWIntegrator.h
vendored
Normal file
329
src/external/BMWIntegrator/BMWIntegrator.h
vendored
Normal file
@ -0,0 +1,329 @@
|
|||||||
|
/***************************************************************************
|
||||||
|
|
||||||
|
BMWIntegrator.h
|
||||||
|
|
||||||
|
Author: Bastian M. Wojek
|
||||||
|
e-mail: bastian.wojek@psi.ch
|
||||||
|
|
||||||
|
$Id:
|
||||||
|
|
||||||
|
***************************************************************************/
|
||||||
|
|
||||||
|
/***************************************************************************
|
||||||
|
* Copyright (C) 2009 by Bastian M. Wojek *
|
||||||
|
* bastian.wojek@psi.ch *
|
||||||
|
* *
|
||||||
|
* This program is free software; you can redistribute it and/or modify *
|
||||||
|
* it under the terms of the GNU General Public License as published by *
|
||||||
|
* the Free Software Foundation; either version 2 of the License, or *
|
||||||
|
* (at your option) any later version. *
|
||||||
|
* *
|
||||||
|
* This program is distributed in the hope that it will be useful, *
|
||||||
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||||
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||||
|
* GNU General Public License for more details. *
|
||||||
|
* *
|
||||||
|
* You should have received a copy of the GNU General Public License *
|
||||||
|
* along with this program; if not, write to the *
|
||||||
|
* Free Software Foundation, Inc., *
|
||||||
|
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
|
||||||
|
***************************************************************************/
|
||||||
|
|
||||||
|
#ifndef _BMWIntegrator_H_
|
||||||
|
#define _BMWIntegrator_H_
|
||||||
|
|
||||||
|
#include "Math/GSLIntegrator.h"
|
||||||
|
#include "Math/GSLMCIntegrator.h"
|
||||||
|
#include "TMath.h"
|
||||||
|
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
using namespace std;
|
||||||
|
|
||||||
|
// 1D Integrator base class - the function to be integrated have to be implemented in a derived class
|
||||||
|
|
||||||
|
class TIntegrator {
|
||||||
|
public:
|
||||||
|
TIntegrator();
|
||||||
|
virtual ~TIntegrator();
|
||||||
|
void SetParameters(const std::vector<double> &par) const { fPar=par; }
|
||||||
|
virtual double FuncAtX(double) const = 0;
|
||||||
|
double IntegrateFunc(double, double);
|
||||||
|
|
||||||
|
protected:
|
||||||
|
mutable vector<double> fPar;
|
||||||
|
|
||||||
|
private:
|
||||||
|
static double FuncAtXgsl(double, void *);
|
||||||
|
ROOT::Math::GSLIntegrator *fIntegrator;
|
||||||
|
mutable double (*fFunc)(double, void *);
|
||||||
|
};
|
||||||
|
|
||||||
|
inline TIntegrator::TIntegrator() : fFunc(0) {
|
||||||
|
ROOT::Math::GSLIntegrator *integrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS31);
|
||||||
|
fIntegrator = integrator;
|
||||||
|
integrator = 0;
|
||||||
|
delete integrator;
|
||||||
|
}
|
||||||
|
|
||||||
|
inline TIntegrator::~TIntegrator(){
|
||||||
|
delete fIntegrator;
|
||||||
|
fIntegrator=0;
|
||||||
|
fFunc=0;
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double TIntegrator::FuncAtXgsl(double x, void *obj)
|
||||||
|
{
|
||||||
|
return ((TIntegrator*)obj)->FuncAtX(x);
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double TIntegrator::IntegrateFunc(double x1, double x2)
|
||||||
|
{
|
||||||
|
fFunc = &TIntegrator::FuncAtXgsl;
|
||||||
|
return fIntegrator->Integral(fFunc, (this), x1, x2);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Multi dimensional GSL Monte Carlo Integrations
|
||||||
|
|
||||||
|
class TMCIntegrator {
|
||||||
|
public:
|
||||||
|
TMCIntegrator();
|
||||||
|
virtual ~TMCIntegrator();
|
||||||
|
void SetParameters(const std::vector<double> &par) const { fPar=par; }
|
||||||
|
virtual double FuncAtX(double *) const = 0;
|
||||||
|
double IntegrateFunc(unsigned int, double *, double *);
|
||||||
|
|
||||||
|
protected:
|
||||||
|
mutable vector<double> fPar;
|
||||||
|
|
||||||
|
private:
|
||||||
|
static double FuncAtXgsl(double *, unsigned int, void *);
|
||||||
|
ROOT::Math::GSLMCIntegrator *fMCIntegrator;
|
||||||
|
mutable double (*fFunc)(double *, unsigned int, void *);
|
||||||
|
};
|
||||||
|
|
||||||
|
inline TMCIntegrator::TMCIntegrator() : fFunc(0) {
|
||||||
|
ROOT::Math::GSLMCIntegrator *integrator = new ROOT::Math::GSLMCIntegrator(ROOT::Math::MCIntegration::kMISER, 1.E-6, 1.E-4, 500000);
|
||||||
|
fMCIntegrator = integrator;
|
||||||
|
integrator = 0;
|
||||||
|
delete integrator;
|
||||||
|
}
|
||||||
|
|
||||||
|
inline TMCIntegrator::~TMCIntegrator(){
|
||||||
|
delete fMCIntegrator;
|
||||||
|
fMCIntegrator=0;
|
||||||
|
fFunc=0;
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double TMCIntegrator::FuncAtXgsl(double *x, unsigned int dim, void *obj)
|
||||||
|
{
|
||||||
|
return ((TMCIntegrator*)obj)->FuncAtX(x);
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double TMCIntegrator::IntegrateFunc(unsigned int dim, double *x1, double *x2)
|
||||||
|
{
|
||||||
|
fFunc = &TMCIntegrator::FuncAtXgsl;
|
||||||
|
return fMCIntegrator->Integral(fFunc, dim, x1, x2, (this));
|
||||||
|
}
|
||||||
|
|
||||||
|
// Multidimensional Integrator class for a d-wave gap integral using the Cuhre algorithm
|
||||||
|
|
||||||
|
class TDWaveGapIntegralCuhre {
|
||||||
|
public:
|
||||||
|
TDWaveGapIntegralCuhre() : fNDim(2) {}
|
||||||
|
~TDWaveGapIntegralCuhre() { fPar.clear(); }
|
||||||
|
void SetParameters(const std::vector<double> &par) { fPar=par; }
|
||||||
|
static void Integrand(const int*, const double[], const int*, double[]);
|
||||||
|
double IntegrateFunc();
|
||||||
|
|
||||||
|
protected:
|
||||||
|
static vector<double> fPar;
|
||||||
|
unsigned int fNDim;
|
||||||
|
|
||||||
|
};
|
||||||
|
|
||||||
|
class TAnSWaveGapIntegralCuhre {
|
||||||
|
public:
|
||||||
|
TAnSWaveGapIntegralCuhre() : fNDim(2) {}
|
||||||
|
~TAnSWaveGapIntegralCuhre() { fPar.clear(); }
|
||||||
|
void SetParameters(const std::vector<double> &par) { fPar=par; }
|
||||||
|
static void Integrand(const int*, const double[], const int*, double[]);
|
||||||
|
double IntegrateFunc();
|
||||||
|
|
||||||
|
protected:
|
||||||
|
static vector<double> fPar;
|
||||||
|
unsigned int fNDim;
|
||||||
|
|
||||||
|
};
|
||||||
|
|
||||||
|
class TAnSWaveGapIntegralDivonne {
|
||||||
|
public:
|
||||||
|
TAnSWaveGapIntegralDivonne() : fNDim(2) {}
|
||||||
|
~TAnSWaveGapIntegralDivonne() { fPar.clear(); }
|
||||||
|
void SetParameters(const std::vector<double> &par) { fPar=par; }
|
||||||
|
static void Integrand(const int*, const double[], const int*, double[]);
|
||||||
|
double IntegrateFunc();
|
||||||
|
|
||||||
|
protected:
|
||||||
|
static vector<double> fPar;
|
||||||
|
unsigned int fNDim;
|
||||||
|
|
||||||
|
};
|
||||||
|
|
||||||
|
class TAnSWaveGapIntegralSuave {
|
||||||
|
public:
|
||||||
|
TAnSWaveGapIntegralSuave() : fNDim(2) {}
|
||||||
|
~TAnSWaveGapIntegralSuave() { fPar.clear(); }
|
||||||
|
void SetParameters(const std::vector<double> &par) { fPar=par; }
|
||||||
|
static void Integrand(const int*, const double[], const int*, double[]);
|
||||||
|
double IntegrateFunc();
|
||||||
|
|
||||||
|
protected:
|
||||||
|
static vector<double> fPar;
|
||||||
|
unsigned int fNDim;
|
||||||
|
|
||||||
|
};
|
||||||
|
|
||||||
|
class TNonMonDWave1GapIntegralCuhre {
|
||||||
|
public:
|
||||||
|
TNonMonDWave1GapIntegralCuhre() : fNDim(2) {}
|
||||||
|
~TNonMonDWave1GapIntegralCuhre() { fPar.clear(); }
|
||||||
|
void SetParameters(const std::vector<double> &par) { fPar=par; }
|
||||||
|
static void Integrand(const int*, const double[], const int*, double[]);
|
||||||
|
double IntegrateFunc();
|
||||||
|
|
||||||
|
protected:
|
||||||
|
static vector<double> fPar;
|
||||||
|
unsigned int fNDim;
|
||||||
|
|
||||||
|
};
|
||||||
|
|
||||||
|
class TNonMonDWave2GapIntegralCuhre {
|
||||||
|
public:
|
||||||
|
TNonMonDWave2GapIntegralCuhre() : fNDim(2) {}
|
||||||
|
~TNonMonDWave2GapIntegralCuhre() { fPar.clear(); }
|
||||||
|
void SetParameters(const std::vector<double> &par) { fPar=par; }
|
||||||
|
static void Integrand(const int*, const double[], const int*, double[]);
|
||||||
|
double IntegrateFunc();
|
||||||
|
|
||||||
|
protected:
|
||||||
|
static vector<double> fPar;
|
||||||
|
unsigned int fNDim;
|
||||||
|
|
||||||
|
};
|
||||||
|
|
||||||
|
// To be integrated: x*y dx dy
|
||||||
|
|
||||||
|
class T2DTest : public TMCIntegrator {
|
||||||
|
public:
|
||||||
|
T2DTest() {}
|
||||||
|
~T2DTest() {}
|
||||||
|
double FuncAtX(double *) const;
|
||||||
|
};
|
||||||
|
|
||||||
|
inline double T2DTest::FuncAtX(double *x) const
|
||||||
|
{
|
||||||
|
return x[0]*x[1];
|
||||||
|
}
|
||||||
|
|
||||||
|
// To be integrated: d wave gap integral
|
||||||
|
|
||||||
|
class TDWaveGapIntegral : public TMCIntegrator {
|
||||||
|
public:
|
||||||
|
TDWaveGapIntegral() {}
|
||||||
|
~TDWaveGapIntegral() {}
|
||||||
|
double FuncAtX(double *) const;
|
||||||
|
};
|
||||||
|
|
||||||
|
inline double TDWaveGapIntegral::FuncAtX(double *x) const // x = {E, phi}, fPar = {T, Delta(T)}
|
||||||
|
{
|
||||||
|
double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K
|
||||||
|
double deltasq(TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]),2.0));
|
||||||
|
return -1.0/(2.0*twokt*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt)*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt));
|
||||||
|
}
|
||||||
|
|
||||||
|
// To be integrated: anisotropic s wave gap integral
|
||||||
|
|
||||||
|
class TAnSWaveGapIntegral : public TMCIntegrator {
|
||||||
|
public:
|
||||||
|
TAnSWaveGapIntegral() {}
|
||||||
|
~TAnSWaveGapIntegral() {}
|
||||||
|
double FuncAtX(double *) const;
|
||||||
|
};
|
||||||
|
|
||||||
|
inline double TAnSWaveGapIntegral::FuncAtX(double *x) const // x = {E, phi}, fPar = {T, Delta(T), a}
|
||||||
|
{
|
||||||
|
double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K
|
||||||
|
double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1])),2.0));
|
||||||
|
return -1.0/(2.0*twokt*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt)*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt));
|
||||||
|
}
|
||||||
|
|
||||||
|
// To be integrated: Bessel function times Exponential
|
||||||
|
|
||||||
|
class TIntBesselJ0Exp : public TIntegrator {
|
||||||
|
public:
|
||||||
|
TIntBesselJ0Exp() {}
|
||||||
|
~TIntBesselJ0Exp() {}
|
||||||
|
double FuncAtX(double) const;
|
||||||
|
};
|
||||||
|
|
||||||
|
inline double TIntBesselJ0Exp::FuncAtX(double x) const
|
||||||
|
{
|
||||||
|
double w0t(TMath::TwoPi()*fPar[0]*x);
|
||||||
|
double j0;
|
||||||
|
if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
|
||||||
|
j0 = 1.0;
|
||||||
|
} else {
|
||||||
|
j0 = TMath::Sin(w0t)/w0t;
|
||||||
|
}
|
||||||
|
return j0 * TMath::Exp(-fPar[1]*x);
|
||||||
|
}
|
||||||
|
|
||||||
|
// To be integrated: Sine times Gaussian
|
||||||
|
|
||||||
|
class TIntSinGss : public TIntegrator {
|
||||||
|
public:
|
||||||
|
TIntSinGss() {}
|
||||||
|
~TIntSinGss() {}
|
||||||
|
double FuncAtX(double) const;
|
||||||
|
};
|
||||||
|
|
||||||
|
inline double TIntSinGss::FuncAtX(double x) const
|
||||||
|
{
|
||||||
|
return TMath::Sin(TMath::TwoPi()*fPar[0]*x) * TMath::Exp(-0.5*fPar[1]*fPar[1]*x*x);
|
||||||
|
}
|
||||||
|
|
||||||
|
// To be integrated: DeRenzi Spin Glass Interpolation Integrand
|
||||||
|
|
||||||
|
class TIntSGInterpolation : public TIntegrator {
|
||||||
|
public:
|
||||||
|
TIntSGInterpolation() {}
|
||||||
|
~TIntSGInterpolation() {}
|
||||||
|
double FuncAtX(double) const;
|
||||||
|
};
|
||||||
|
|
||||||
|
inline double TIntSGInterpolation::FuncAtX(double x) const
|
||||||
|
{
|
||||||
|
// Parameters: nu_L [MHz], a [1/us], lambda [1/us], beta [1], t [us]
|
||||||
|
double wt(TMath::TwoPi()*fPar[0]*x);
|
||||||
|
double expo(0.5*fPar[1]*fPar[1]*x*x/fPar[3]+fPar[2]*fPar[4]);
|
||||||
|
return (wt*TMath::Cos(wt)-TMath::Sin(wt))/(wt*wt)*TMath::Exp(-TMath::Power(expo,fPar[3]))/TMath::Power(expo,(1.0-fPar[3]));
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
// To be integrated: df/dE * E / sqrt(E^2 - Delta^2)
|
||||||
|
|
||||||
|
class TGapIntegral : public TIntegrator {
|
||||||
|
public:
|
||||||
|
TGapIntegral() {}
|
||||||
|
~TGapIntegral() {}
|
||||||
|
double FuncAtX(double) const; // variable: E
|
||||||
|
|
||||||
|
};
|
||||||
|
|
||||||
|
inline double TGapIntegral::FuncAtX(double e) const
|
||||||
|
{
|
||||||
|
return 1.0/(TMath::Power(TMath::CosH(TMath::Sqrt(e*e+fPar[1]*fPar[1])/fPar[0]),2.0));
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif //_TIntegrator_H_
|
||||||
110
src/external/libLFRelaxation/TIntegrator.h
vendored
110
src/external/libLFRelaxation/TIntegrator.h
vendored
@ -1,110 +0,0 @@
|
|||||||
/***************************************************************************
|
|
||||||
|
|
||||||
TIntegrator.h
|
|
||||||
|
|
||||||
Author: Bastian M. Wojek
|
|
||||||
e-mail: bastian.wojek@psi.ch
|
|
||||||
|
|
||||||
2008/12/24
|
|
||||||
|
|
||||||
***************************************************************************/
|
|
||||||
|
|
||||||
#ifndef _TIntegrator_H_
|
|
||||||
#define _TIntegrator_H_
|
|
||||||
|
|
||||||
#include "Math/GSLIntegrator.h"
|
|
||||||
#include "TMath.h"
|
|
||||||
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
using namespace std;
|
|
||||||
|
|
||||||
// Integrator base class - the function to be integrated have to be implemented in a derived class
|
|
||||||
|
|
||||||
class TIntegrator {
|
|
||||||
public:
|
|
||||||
TIntegrator();
|
|
||||||
virtual ~TIntegrator();
|
|
||||||
void SetParameters(const std::vector<double> &par) const { fPar=par; }
|
|
||||||
virtual double FuncAtX(double) const = 0;
|
|
||||||
double IntegrateFunc(double, double);
|
|
||||||
|
|
||||||
protected:
|
|
||||||
mutable vector<double> fPar;
|
|
||||||
|
|
||||||
private:
|
|
||||||
static double FuncAtXgsl(double, void *);
|
|
||||||
ROOT::Math::GSLIntegrator *fIntegrator;
|
|
||||||
mutable double (*fFunc)(double, void *);
|
|
||||||
};
|
|
||||||
|
|
||||||
inline TIntegrator::TIntegrator() : fFunc(0) {
|
|
||||||
ROOT::Math::GSLIntegrator *integrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS31);
|
|
||||||
fIntegrator = integrator;
|
|
||||||
integrator = 0;
|
|
||||||
delete integrator;
|
|
||||||
}
|
|
||||||
|
|
||||||
inline TIntegrator::~TIntegrator(){
|
|
||||||
delete fIntegrator;
|
|
||||||
fIntegrator=0;
|
|
||||||
fFunc=0;
|
|
||||||
}
|
|
||||||
|
|
||||||
inline double TIntegrator::FuncAtXgsl(double x, void *obj)
|
|
||||||
{
|
|
||||||
return ((TIntegrator*)obj)->FuncAtX(x);
|
|
||||||
}
|
|
||||||
|
|
||||||
inline double TIntegrator::IntegrateFunc(double x1, double x2)
|
|
||||||
{
|
|
||||||
fFunc = &TIntegrator::FuncAtXgsl;
|
|
||||||
return fIntegrator->Integral(fFunc, (this), x1, x2);
|
|
||||||
}
|
|
||||||
|
|
||||||
// To be integrated: Bessel function times Exponential
|
|
||||||
|
|
||||||
class TIntBesselJ0Exp : public TIntegrator {
|
|
||||||
public:
|
|
||||||
TIntBesselJ0Exp() {}
|
|
||||||
~TIntBesselJ0Exp() {}
|
|
||||||
double FuncAtX(double) const;
|
|
||||||
};
|
|
||||||
|
|
||||||
inline double TIntBesselJ0Exp::FuncAtX(double x) const
|
|
||||||
{
|
|
||||||
return TMath::BesselJ0(TMath::TwoPi()*fPar[0]*x) * TMath::Exp(-fPar[1]*x);
|
|
||||||
}
|
|
||||||
|
|
||||||
// To be integrated: Sine times Gaussian
|
|
||||||
|
|
||||||
class TIntSinGss : public TIntegrator {
|
|
||||||
public:
|
|
||||||
TIntSinGss() {}
|
|
||||||
~TIntSinGss() {}
|
|
||||||
double FuncAtX(double) const;
|
|
||||||
};
|
|
||||||
|
|
||||||
inline double TIntSinGss::FuncAtX(double x) const
|
|
||||||
{
|
|
||||||
return TMath::Sin(TMath::TwoPi()*fPar[0]*x) * TMath::Exp(-0.5*fPar[1]*fPar[1]*x*x);
|
|
||||||
}
|
|
||||||
|
|
||||||
// To be integrated: df/dE * E / sqrt(E^2+ Delta^2)
|
|
||||||
|
|
||||||
class TGapIntegral : public TIntegrator {
|
|
||||||
public:
|
|
||||||
TGapIntegral() {}
|
|
||||||
~TGapIntegral() {}
|
|
||||||
double FuncAtX(double) const; // parameter: E
|
|
||||||
|
|
||||||
};
|
|
||||||
|
|
||||||
inline double TGapIntegral::FuncAtX(double e) const
|
|
||||||
{
|
|
||||||
double kt(0.08617384436*fPar[0]); // kB in meV/K
|
|
||||||
double expekt(TMath::Exp(e/kt));
|
|
||||||
return -expekt*e/(kt*(1.0+expekt)*TMath::Sqrt(e*e+fPar[1]*fPar[1]));
|
|
||||||
}
|
|
||||||
|
|
||||||
#endif //_TIntegrator_H_
|
|
||||||
Loading…
x
Reference in New Issue
Block a user