LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

BMWlibs: add two p-wave gap funtions which can be used to evaluate the superfluid density.

Browse Source
This commit is contained in:
suter_a 2020-11-04 15:03:58 +01:00
parent 6ef53c6b6a
commit 32cf3221d9
7 changed files with 836 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

148
src/external/BMWtools/BMWIntegrator.cpp vendored
View File

@ -35,8 +35,127 @@
#define SEED 0
#define STATEFILE NULL
//-----------------------------------------------------------------------------
std::vector<double> TPointPWaveGapIntegralCuhre::fPar;
//-----------------------------------------------------------------------------
/**
* <p>Integrate the function using the Cuhre interface
*
* <p><b>return:</b>
* - value of the integral
*/
double TPointPWaveGapIntegralCuhre::IntegrateFunc()
{
const unsigned int NCOMP(1);
const unsigned int NVEC(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, USERDATA, NVEC,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY, STATEFILE, SPIN,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value for the use with Cuhre---actual implementation of the function
*
* <p><b>return:</b>
* - 0
*
* \param ndim number of dimensions of the integral (2 here)
* \param x point where the function should be evaluated
* \param ncomp number of components of the integrand (1 here)
* \param f function value
* \param userdata additional user parameters (required by the interface, NULL here)
*/
int TPointPWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
const int *ncomp, double f[], void *userdata) // x = {E, theta}, fPar = {twokBT, Delta(T), Ec, thetac}
{
double sinTheta = TMath::Sin(x[1]*fPar[3]);
double deltasq(TMath::Power(fPar[1]*sinTheta,2.0));
f[0] = sinTheta/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0);
return 0;
}
//-----------------------------------------------------------------------------
std::vector<double> TLinePWaveGapIntegralCuhre::fPar;
//-----------------------------------------------------------------------------
/**
* <p>Integrate the function using the Cuhre interface
*
* <p><b>return:</b>
* - value of the integral
*/
double TLinePWaveGapIntegralCuhre::IntegrateFunc()
{
const unsigned int NCOMP(1);
const unsigned int NVEC(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, USERDATA, NVEC,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY, STATEFILE, SPIN,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value for the use with Cuhre---actual implementation of the function
*
* <p><b>return:</b>
* - 0
*
* \param ndim number of dimensions of the integral (2 here)
* \param x point where the function should be evaluated
* \param ncomp number of components of the integrand (1 here)
* \param f function value
* \param userdata additional user parameters (required by the interface, NULL here)
*/
int TLinePWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
const int *ncomp, double f[], void *userdata) // x = {E, theta}, fPar = {twokBT, Delta(T), Ec, thetac}
{
double sinTheta = TMath::Sin(x[1]*fPar[3]);
double cosTheta = TMath::Cos(x[1]*fPar[3]);
double deltasq(TMath::Power(fPar[1]*cosTheta,2.0));
f[0] = sinTheta/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0);
return 0;
}
//-----------------------------------------------------------------------------
std::vector<double> TDWaveGapIntegralCuhre::fPar;
//-----------------------------------------------------------------------------
/**
* <p>Integrate the function using the Cuhre interface
*
@ -67,6 +186,7 @@ double TDWaveGapIntegralCuhre::IntegrateFunc()
return integral[0];
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value for the use with Cuhre---actual implementation of the function
*
@ -87,8 +207,11 @@ int TDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
return 0;
}
//-----------------------------------------------------------------------------
std::vector<double> TCosSqDWaveGapIntegralCuhre::fPar;
//-----------------------------------------------------------------------------
/**
* <p>Integrate the function using the Cuhre interface
*
@ -119,6 +242,7 @@ double TCosSqDWaveGapIntegralCuhre::IntegrateFunc()
return integral[0];
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value for the use with Cuhre---actual implementation of the function
*
@ -139,8 +263,11 @@ int TCosSqDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
return 0;
}
//-----------------------------------------------------------------------------
std::vector<double> TSinSqDWaveGapIntegralCuhre::fPar;
//-----------------------------------------------------------------------------
/**
* <p>Integrate the function using the Cuhre interface
*
@ -171,6 +298,7 @@ double TSinSqDWaveGapIntegralCuhre::IntegrateFunc()
return integral[0];
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value for the use with Cuhre---actual implementation of the function
*
@ -191,8 +319,11 @@ int TSinSqDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
return 0;
}
//-----------------------------------------------------------------------------
std::vector<double> TAnSWaveGapIntegralCuhre::fPar;
//-----------------------------------------------------------------------------
/**
* <p>Integrate the function using the Cuhre interface
*
@ -223,6 +354,7 @@ double TAnSWaveGapIntegralCuhre::IntegrateFunc()
return integral[0];
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value for the use with Cuhre---actual implementation of the function
*
@ -243,8 +375,11 @@ int TAnSWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
return 0;
}
//-----------------------------------------------------------------------------
std::vector<double> TAnSWaveGapIntegralDivonne::fPar;
//-----------------------------------------------------------------------------
/**
* <p>Integrate the function using the Divonne interface
*
@ -283,6 +418,7 @@ double TAnSWaveGapIntegralDivonne::IntegrateFunc()
return integral[0];
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value for the use with Divonne---actual implementation of the function
*
@ -303,8 +439,11 @@ int TAnSWaveGapIntegralDivonne::Integrand(const int *ndim, const double x[],
return 0;
}
//-----------------------------------------------------------------------------
std::vector<double> TAnSWaveGapIntegralSuave::fPar;
//-----------------------------------------------------------------------------
/**
* <p>Integrate the function using the Suave interface
*
@ -337,6 +476,7 @@ double TAnSWaveGapIntegralSuave::IntegrateFunc()
return integral[0];
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value for the use with Suave---actual implementation of the function
*
@ -357,8 +497,11 @@ int TAnSWaveGapIntegralSuave::Integrand(const int *ndim, const double x[],
return 0;
}
//-----------------------------------------------------------------------------
std::vector<double> TNonMonDWave1GapIntegralCuhre::fPar;
//-----------------------------------------------------------------------------
/**
* <p>Integrate the function using the Cuhre interface
*
@ -389,6 +532,7 @@ double TNonMonDWave1GapIntegralCuhre::IntegrateFunc()
return integral[0];
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value for the use with Cuhre---actual implementation of the function
*
@ -409,8 +553,11 @@ int TNonMonDWave1GapIntegralCuhre::Integrand(const int *ndim, const double x[],
return 0;
}
//-----------------------------------------------------------------------------
std::vector<double> TNonMonDWave2GapIntegralCuhre::fPar;
//-----------------------------------------------------------------------------
/**
* <p>Integrate the function using the Cuhre interface
*
@ -441,6 +588,7 @@ double TNonMonDWave2GapIntegralCuhre::IntegrateFunc()
return integral[0];
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value for the use with Cuhre---actual implementation of the function
*

149
src/external/BMWtools/BMWIntegrator.h vendored
View File

@ -36,6 +36,7 @@
#include <cmath>
#include <vector>
//-----------------------------------------------------------------------------
/**
* <p>Alternative base class for 1D integrations using the GNU Scientific Library integrator.
* The difference to the other class is that here the integration parameters have to be supplied directly to the IntegrateFunc method.
@ -54,6 +55,7 @@ class T2Integrator {
static double FuncAtXgsl(double, void *);
};
//-----------------------------------------------------------------------------
/**
* <p>Method for passing the integrand function value to the integrator.
*
@ -69,6 +71,7 @@ inline double T2Integrator::FuncAtXgsl(double x, void *ptrPair)
return pairOfPointers->first->FuncAtX(x, *(pairOfPointers->second));
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the integral of the function between the given boundaries
*
@ -93,20 +96,7 @@ inline double T2Integrator::IntegrateFunc(double x1, double x2, const std::vecto
return value;
}
//-----------------------------------------------------------------------------
/**
* <p>Base class for 1D integrations using the GNU Scientific Library integrator.
* The function which should be integrated has to be implemented in a derived class.
@ -129,6 +119,7 @@ class TIntegrator {
mutable double (*fFunc)(double, void *); ///< pointer to the integrand function
};
//-----------------------------------------------------------------------------
/**
* <p>Constructor of the base class for 1D integrations
* Allocation of memory for an integration using the adaptive 31 point Gauss-Kronrod rule
@ -137,6 +128,7 @@ inline TIntegrator::TIntegrator() : fFunc(0) {
fIntegrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS31);
}
//-----------------------------------------------------------------------------
/**
* <p>Destructor of the base class for 1D integrations
* Clean up.
@ -148,6 +140,7 @@ inline TIntegrator::~TIntegrator(){
fFunc=0;
}
//-----------------------------------------------------------------------------
/**
* <p>Method for passing the integrand function value to the integrator.
*
@ -162,6 +155,7 @@ inline double TIntegrator::FuncAtXgsl(double x, void *obj)
return ((TIntegrator*)obj)->FuncAtX(x);
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the integral of the function between the given boundaries
*
@ -177,6 +171,7 @@ inline double TIntegrator::IntegrateFunc(double x1, double x2)
return fIntegrator->Integral(fFunc, (this), x1, x2);
}
//-----------------------------------------------------------------------------
/**
* <p>Base class for multidimensional Monte-Carlo integrations using the GNU Scientific Library integrator.
* The function which should be integrated has to be implemented in a derived class.
@ -199,6 +194,7 @@ class TMCIntegrator {
mutable double (*fFunc)(double *, size_t, void *); ///< pointer to the integrand function
};
//-----------------------------------------------------------------------------
/**
* <p>Constructor of the base class for multidimensional Monte-Carlo integrations
* Allocation of memory for an integration using the MISER algorithm of Press and Farrar
@ -207,6 +203,7 @@ inline TMCIntegrator::TMCIntegrator() : fFunc(0) {
fMCIntegrator = new ROOT::Math::GSLMCIntegrator(ROOT::Math::MCIntegration::kMISER, 1.E-6, 1.E-4, 500000);
}
//-----------------------------------------------------------------------------
/**
* <p>Destructor of the base class for 1D integrations
* Clean up.
@ -218,6 +215,7 @@ inline TMCIntegrator::~TMCIntegrator(){
fFunc=0;
}
//-----------------------------------------------------------------------------
/**
* <p>Method for passing the integrand function value to the integrator.
*
@ -233,6 +231,7 @@ inline double TMCIntegrator::FuncAtXgsl(double *x, size_t dim, void *obj)
return ((TMCIntegrator*)obj)->FuncAtX(x);
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the integral of the function between the given boundaries
*
@ -249,6 +248,45 @@ inline double TMCIntegrator::IntegrateFunc(size_t dim, double *x1, double *x2)
return fMCIntegrator->Integral(fFunc, dim, x1, x2, (this));
}
//-----------------------------------------------------------------------------
/**
* <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and a point p symmetry of the superconducting order parameter.
* The integration uses the Cuhre algorithm of the Cuba library.
*/
class TPointPWaveGapIntegralCuhre {
public:
TPointPWaveGapIntegralCuhre() : fNDim(2) {}
~TPointPWaveGapIntegralCuhre() { fPar.clear(); }
void SetParameters(const std::vector<double> &par) { fPar=par; }
static int Integrand(const int*, const double[], const int*, double[], void*);
double IntegrateFunc();
protected:
static std::vector<double> fPar; ///< parameters of the integrand
unsigned int fNDim; ///< dimension of the integral
};
//-----------------------------------------------------------------------------
/**
* <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and a line p symmetry of the superconducting order parameter.
* The integration uses the Cuhre algorithm of the Cuba library.
*/
class TLinePWaveGapIntegralCuhre {
public:
TLinePWaveGapIntegralCuhre() : fNDim(2) {}
~TLinePWaveGapIntegralCuhre() { fPar.clear(); }
void SetParameters(const std::vector<double> &par) { fPar=par; }
static int Integrand(const int*, const double[], const int*, double[], void*);
double IntegrateFunc();
protected:
static std::vector<double> fPar; ///< parameters of the integrand
unsigned int fNDim; ///< dimension of the integral
};
//-----------------------------------------------------------------------------
/**
* <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and a d_{x^2-y^2} symmetry of the superconducting order parameter.
@ -267,6 +305,7 @@ class TDWaveGapIntegralCuhre {
unsigned int fNDim; ///< dimension of the integral
};
//-----------------------------------------------------------------------------
/**
* <p>Two-dimensional integrator class for the efficient calculation of the superfluid density along the a-axis
* within the semi-classical model assuming a cylindrical Fermi surface and a mixed d_{x^2-y^2} + s symmetry of the
@ -286,6 +325,7 @@ class TCosSqDWaveGapIntegralCuhre {
unsigned int fNDim; ///< dimension of the integral
};
//-----------------------------------------------------------------------------
/**
* <p>Two-dimensional integrator class for the efficient calculation of the superfluid density along the b-axis
* within the semi-classical model assuming a cylindrical Fermi surface and a mixed d_{x^2-y^2} + s symmetry of the
@ -305,6 +345,7 @@ class TSinSqDWaveGapIntegralCuhre {
unsigned int fNDim; ///< dimension of the integral
};
//-----------------------------------------------------------------------------
/**
* <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and an "anisotropic s-wave" symmetry of the superconducting order parameter.
@ -323,6 +364,7 @@ class TAnSWaveGapIntegralCuhre {
unsigned int fNDim; ///< dimension of the integral
};
//-----------------------------------------------------------------------------
/**
* <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and an "anisotropic s-wave" symmetry of the superconducting order parameter.
@ -341,6 +383,7 @@ class TAnSWaveGapIntegralDivonne {
unsigned int fNDim; ///< dimension of the integral
};
//-----------------------------------------------------------------------------
/**
* <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and an "anisotropic s-wave" symmetry of the superconducting order parameter.
@ -359,6 +402,7 @@ class TAnSWaveGapIntegralSuave {
unsigned int fNDim; ///< dimension of the integral
};
//-----------------------------------------------------------------------------
/**
* <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and an "non-monotonic d-wave" symmetry of the superconducting order parameter.
@ -377,6 +421,7 @@ class TNonMonDWave1GapIntegralCuhre {
unsigned int fNDim; ///< dimension of the integral
};
//-----------------------------------------------------------------------------
/**
* <p>Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and an "non-monotonic d-wave" symmetry of the superconducting order parameter.
@ -395,6 +440,7 @@ class TNonMonDWave2GapIntegralCuhre {
unsigned int fNDim; ///< dimension of the integral
};
//-----------------------------------------------------------------------------
/**
* <p>Test class for the 2D MC integration
* Integral: x*y dx dy
@ -406,6 +452,7 @@ class T2DTest : public TMCIntegrator {
double FuncAtX(double *) const;
};
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value---actual implementation of the function x*y
*
@ -419,6 +466,65 @@ inline double T2DTest::FuncAtX(double *x) const
return x[0]*x[1];
}
//-----------------------------------------------------------------------------
/**
* <p>Class for the 2D Monte-Carlo integration for the calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and a point p symmetry of the superconducting order parameter.
* The integration uses the GSL integration routines.
*/
class TPointPWaveGapIntegral : public TMCIntegrator {
public:
TPointPWaveGapIntegral() {}
~TPointPWaveGapIntegral() {}
double FuncAtX(double *) const;
};
//-----------------------------------------------------------------------------
/**
* <p>Class for the 2D Monte-Carlo integration for the calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and a line p symmetry of the superconducting order parameter.
* The integration uses the GSL integration routines.
*/
class TLinePWaveGapIntegral : public TMCIntegrator {
public:
TLinePWaveGapIntegral() {}
~TLinePWaveGapIntegral() {}
double FuncAtX(double *) const;
};
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value---actual implementation of the function
*
* <p><b>return:</b>
* - function value
*
* \param x point where the function should be evaluated
*/
inline double TPointPWaveGapIntegral::FuncAtX(double *x) const // x = {E, theta}, fPar = {T, Delta(T)}
{
double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K
double deltasq(TMath::Power(fPar[1]*TMath::Sin(x[1]),2.0));
return -TMath::Sin(x[1])/(4.0*twokt*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt)*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt));
}
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value---actual implementation of the function
*
* <p><b>return:</b>
* - function value
*
* \param x point where the function should be evaluated
*/
inline double TLinePWaveGapIntegral::FuncAtX(double *x) const // x = {E, theta}, fPar = {T, Delta(T)}
{
double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K
double deltasq(TMath::Power(fPar[1]*TMath::Cos(x[1]),2.0));
return -TMath::Sin(x[1])/(4.0*twokt*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt)*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt));
}
//-----------------------------------------------------------------------------
/**
* <p>Class for the 2D Monte-Carlo integration for the calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and a d_{x^2-y^2} symmetry of the superconducting order parameter.
@ -431,6 +537,7 @@ class TDWaveGapIntegral : public TMCIntegrator {
double FuncAtX(double *) const;
};
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value---actual implementation of the function
*
@ -446,6 +553,7 @@ inline double TDWaveGapIntegral::FuncAtX(double *x) const // x = {E, phi}, fPar
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));
}
//-----------------------------------------------------------------------------
/**
* <p>Class for the 2D Monte-Carlo integration for the calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and an "anisotropic s-wave" symmetry of the superconducting order parameter.
@ -458,6 +566,7 @@ class TAnSWaveGapIntegral : public TMCIntegrator {
double FuncAtX(double *) const;
};
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value---actual implementation of the function
*
@ -473,6 +582,7 @@ inline double TAnSWaveGapIntegral::FuncAtX(double *x) const // x = {E, phi}, fPa
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));
}
//-----------------------------------------------------------------------------
/**
* <p>Class for the 1D integration of j0(a*x)*exp(-b*x)
* The integration uses the GSL integration routines.
@ -486,6 +596,7 @@ class TIntBesselJ0Exp : public T2Integrator {
double FuncAtX(double, const std::vector<double>&) const;
};
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value---actual implementation of the function j0(a*x)*exp(-b*x)
*
@ -506,6 +617,7 @@ inline double TIntBesselJ0Exp::FuncAtX(double x, const std::vector<double> &par)
return j0 * TMath::Exp(-par[1]*x);
}
//-----------------------------------------------------------------------------
/**
* <p>Class for the 1D integration of sin(a*x)*exp(-b*x*x)
* The integration uses the GSL integration routines.
@ -519,6 +631,7 @@ class TIntSinGss : public T2Integrator {
double FuncAtX(double, const std::vector<double>&) const;
};
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value---actual implementation of the function sin(a*x)*exp(-b*x*x)
*
@ -532,6 +645,7 @@ inline double TIntSinGss::FuncAtX(double x, const std::vector<double> &par) cons
return TMath::Sin(TMath::TwoPi()*par[0]*x) * TMath::Exp(-0.5*par[1]*par[1]*x*x);
}
//-----------------------------------------------------------------------------
/**
* <p>Class for the 1D integration of the "DeRenzi Spin Glass Interpolation Integrand"
* See Eq. (5) of R. De Renzi and S. Fanesi, Physica B 289-290, 209-212 (2000).
@ -547,6 +661,7 @@ class TIntSGInterpolation : public T2Integrator {
double FuncAtX(double, const std::vector<double>&) const;
};
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value---actual implementation of the function
*
@ -563,6 +678,7 @@ inline double TIntSGInterpolation::FuncAtX(double x, const std::vector<double> &
return (wt*TMath::Cos(wt)-TMath::Sin(wt))/(wt*wt)*TMath::Exp(-TMath::Power(expo,par[3]))/TMath::Power(expo,(1.0-par[3]));
}
//-----------------------------------------------------------------------------
/**
* <p>Class for the 1D integration for the calculation of the superfluid density within the semi-classical model
* assuming a cylindrical Fermi surface and an isotropic s-wave symmetry of the superconducting order parameter.
@ -575,6 +691,7 @@ class TGapIntegral : public TIntegrator {
double FuncAtX(double) const; // variable: E
};
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value---actual implementation of the function df/dE * E / sqrt(E^2 - Delta^2)
*
@ -588,6 +705,7 @@ 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));
}
//-----------------------------------------------------------------------------
/**
* <p>Class for the 1D integration for the calculation of the uniaxial static Gauss-Kubo-Toyabe function
* The integration uses the GSL integration routines.
@ -601,6 +719,7 @@ class TFirstUniaxialGssKTIntegral : public T2Integrator {
virtual double FuncAtX(double, const std::vector<double>&) const; // variable: x
};
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value---actual implementation of the integrand in Eq. (7) of Solt's article
*
@ -618,6 +737,7 @@ inline double TFirstUniaxialGssKTIntegral::FuncAtX(double x, const std::vector<d
return (1.0 - x*x)*(p - SsqTsq)/TMath::Power(p, 2.5)*TMath::Exp(-0.5*SsqTsq/p);
}
//-----------------------------------------------------------------------------
/**
* <p>Class for the 1D integration for the calculation of the uniaxial static Gauss-Kubo-Toyabe function
* The integration uses the GSL integration routines.
@ -631,6 +751,7 @@ class TSecondUniaxialGssKTIntegral : public T2Integrator {
virtual double FuncAtX(double, const std::vector<double>&) const; // variable: x
};
//-----------------------------------------------------------------------------
/**
* <p>Calculate the function value---actual implementation of the integrand in Eq. (7) of Solt's article
*

BIN
src/external/libGapIntegrals/GapIntegrals.pdf vendored
View File

Binary file not shown.

48
src/external/libGapIntegrals/GapIntegrals.tex vendored
View File

@ -74,15 +74,16 @@ E-Mail: & \verb?andreas.suter@psi.ch? &&
\section*{\musrfithead plug-in for the calculation of the temperature dependence of $\bm{1/\lambda^2}$ for various gap symmetries}%
This memo is intended to give a short summary of the background on which the \gapint plug-in for \musrfit \cite{musrfit} has been developed. The aim of this implementation is the efficient calculation of integrals of the form
\begin{equation}\label{int}
\begin{equation}\label{eq:int_phi}
I(T) = 1 + \frac{1}{\pi}\int_0^{2\pi}\int_{\Delta(\varphi,T)}^{\infty}\left(\frac{\partial f}{\partial E}\right) \frac{E}{\sqrt{E^2-\Delta^2(\varphi,T)}}\mathrm{d}E\mathrm{d}\varphi\,,
\end{equation}
where $f = (1+\exp(E/k_{\mathrm B}T))^{-1}$, like they appear e.g. in the theoretical temperature dependence of $1/\lambda^2$~\cite{Manzano}.
In order not to do too many unnecessary function calls during the final numerical evaluation we simplify the integral (\ref{int}) as far as possible analytically. The derivative of $f$ is given by
For gap symmetries which involve not only a $E$- and $\varphi$-dependence but also a $\theta$-dependence, see the special section towards the end of the memo.
In order not to do too many unnecessary function calls during the final numerical evaluation we simplify the integral (\ref{eq:int_phi}) as far as possible analytically. The derivative of $f$ is given by
\begin{equation}\label{derivative}
\frac{\partial f}{\partial E} = -\frac{1}{k_{\mathrm B}T}\frac{\exp(E/k_{\mathrm B}T)}{\left(1+\exp(E/k_{\mathrm B}T)\right)^2} = -\frac{1}{4k_{\mathrm B}T} \frac{1}{\cosh^2\left(E/2k_{\mathrm B}T\right)}.
\end{equation}
Using (\ref{derivative}) and doing the substitution $E'^2 = E^2-\Delta^2(\varphi,T)$, equation (\ref{int}) can be written as
Using (\ref{derivative}) and doing the substitution $E'^2 = E^2-\Delta^2(\varphi,T)$, equation (\ref{eq:int_phi}) can be written as
\begin{equation}
\begin{split}
I(T) & = 1 - \frac{1}{4\pi k_{\mathrm B}T}\int_0^{2\pi}\int_{\Delta(\varphi,T)}^{\infty}\frac{1}{\cosh^2\left(E/2k_{\mathrm B}T\right)}\frac{E}{\sqrt{E^2-\Delta^2(\varphi,T)}}\mathrm{d}E\mathrm{d}\varphi \\
@ -167,6 +168,20 @@ The \gapint plug-in calculates $\tilde{I}(T)$ for the following $\Delta(\varphi)
Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $a~(1)$, $[c_0~(1),~a_{\rm G}~(1)]$.
If $c_0$ and $a_{\rm G}$ are provided, the temperature dependence according to Eq.(\ref{eq:gapT_Prozorov}) will be used,
otherwise Eq.(\ref{eq:gapT_Manzano}) will be utilized.
\item[\textit{p}-wave (point) \cite{Pang2015}:]
\begin{equation}
\Delta(\theta, T) = \Delta(T) \sin(\theta)
\end{equation}
\musrfit theory line: \verb?userFcn libGapIntegrals TGapPointPWave 1 2 [3 4]?\\[1.5ex]
Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $[c_0~(1),~a_{\rm G}~(1)]$. If $c_0$ and $a_{\rm G}$ are provided,
the temperature dependence according to Eq.(\ref{eq:gapT_Prozorov}) will be used, otherwise Eq.(\ref{eq:gapT_Manzano}) will be utilized.
\item[\textit{p}-wave (line) \cite{Ozaki1986}:]
\begin{equation}
\Delta(\theta, T) = \Delta(T) \cos(\theta)
\end{equation}
\musrfit theory line: \verb?userFcn libGapIntegrals TGapLinePWave 1 2 [3 4]?\\[1.5ex]
Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $[c_0~(1),~a_{\rm G}~(1)]$. If $c_0$ and $a_{\rm G}$ are provided,
the temperature dependence according to Eq.(\ref{eq:gapT_Prozorov}) will be used, otherwise Eq.(\ref{eq:gapT_Manzano}) will be utilized.
\end{description}
\noindent It is also possible to calculate a power law temperature dependence (in the two fluid approximation $n=4$) and the dirty \textit{s}-wave expression.
@ -193,6 +208,31 @@ Obviously for this no integration is needed.
within the semi-classical model assuming a cylindrical Fermi surface and a mixed $d_{x^2-y^2} + s$ symmetry of the superconducting order parameter
(effectively: $d_{x^2-y^2}$ with shifted nodes and \textit{a}-\textit{b}-anisotropy)) see the source code.
\subsection*{Gap Integrals for $\bm{\theta}$-, and $\bm{(\theta, \varphi)}$-dependent Gaps}%
\noindent For the most general case for which the gap-symmetry is $(E,\varphi,\theta)$-dependent, the integral to be calculate takes the form
\begin{equation}\label{eq:int_phi_theta}
I = 1 + \frac{1}{2\pi} \int_0^\pi \sin(\theta)\, \mathrm{d}\theta \int_0^{2\pi} \mathrm{d}\varphi \int_{\Delta(\varphi,\theta,T)}^\infty \mathrm{d}E \left\{ \left(\frac{\partial f}{\partial E}\right) \frac{E}{\sqrt{E^2 - \Delta^2(\varphi, \theta, T)}} \right\}
\end{equation}
\subsubsection*{Gap Integrals for $\bm{\theta}$-dpendent Gaps}%
\noindent In case the gap-symmetry is only depending on $(E,\theta)$ the $\varphi$-integral can be carried out and Eq.(\ref{eq:int_phi_theta}) simplifies to
\begin{equation}\label{eq:int_theta}
I(T) = 1 + \int_0^\pi \sin(\theta)\, \mathrm{d}\theta \int_{\Delta(\varphi,\theta,T)}^\infty \mathrm{d}E \left\{ \left(\frac{\partial f}{\partial E}\right) \frac{E}{\sqrt{E^2 - \Delta^2(\varphi, \theta, T)}} \right\}
\end{equation}
\noindent Following the same simplification steps as for Eq.(\ref{eq:int_phi}) we find
\begin{equation}\label{eq:int_tilde_theta}
\tilde{I}(T) = 1 - \frac{\pi E_c}{4 k_{\rm B}T}\, \int_0^1 \mathrm{d}x \sin(\pi x) \int_0^1 \mathrm{d}F\, \frac{1}{\cosh^2\left(\sqrt{(E_c\cdot F)^2 + \Delta^2(\pi\cdot x, T)}/(2 k_{\rm B} T)\right)}
\end{equation}
\noindent where $x = \theta / \pi$, and $F = E^\prime / E_c$.
\subsection*{License}
The \gapint library 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 \cite{GPL}; either version 2 of the License, or (at your option) any later version.
@ -208,6 +248,8 @@ The \gapint library is free software; you can redistribute it and/or modify it u
\bibitem{Matsui} H.~Matsui~\textit{et al.}, \textit{Direct Observation of a Nonmonotonic $d_{x^2-y^2}$-Wave Superconducting Gap in the Electron-Doped High-T$_{\mathrm c}$ Superconductor Pr$_{0.89}$LaCe$_{0.11}$CuO$_4$}, Phys.~Rev.~Lett.~\textbf{95}~(2005)~017003
\bibitem{Eremin} I.~Eremin, E.~Tsoncheva, and A.V.~Chubukov, \textit{Signature of the nonmonotonic $d$-wave gap in electron-doped cuprates}, Phys.~Rev.~B~\textbf{77}~(2008)~024508
\bibitem{AnisotropicSWave} ??
\bibitem{Pang2015} G.M.~Pang, \emph{et al.}, Phys.~Rev.~B~\textbf{91}~(2015)~220502(R), and references in there.
\bibitem{Ozaki1986} M.~Ozaki, \emph{et al.}, Prog.~Theor.~Phys.~\textbf{75}~(1986)~442.
\bibitem{Tinkham} M.~Tinkham, \textit{Introduction to Superconductivity} $2^{\rm nd}$ ed. (Dover Publications, New York, 2004).
\bibitem{GPL} http://www.gnu.org/licenses/old-licenses/gpl-2.0.html

369
src/external/libGapIntegrals/TGapIntegrals.cpp vendored
View File

@ -36,6 +36,8 @@
#define TWOPI 6.28318530717958647692
ClassImp(TGapSWave)
ClassImp(TGapPointPWave)
ClassImp(TGapLinePWave)
ClassImp(TGapDWave)
ClassImp(TGapCosSqDWave)
ClassImp(TGapSinSqDWave)
@ -46,6 +48,8 @@ ClassImp(TGapPowerLaw)
ClassImp(TGapDirtySWave)
ClassImp(TLambdaSWave)
ClassImp(TLambdaPointPWave)
ClassImp(TLambdaLinePWave)
ClassImp(TLambdaDWave)
ClassImp(TLambdaAnSWave)
ClassImp(TLambdaNonMonDWave1)
@ -53,6 +57,8 @@ ClassImp(TLambdaNonMonDWave2)
ClassImp(TLambdaPowerLaw)
ClassImp(TLambdaInvSWave)
ClassImp(TLambdaInvPointPWave)
ClassImp(TLambdaInvLinePWave)
ClassImp(TLambdaInvDWave)
ClassImp(TLambdaInvAnSWave)
ClassImp(TLambdaInvNonMonDWave1)
@ -77,6 +83,38 @@ TGapSWave::TGapSWave() {
fPar.clear();
}
//--------------------------------------------------------------------
/**
* <p> point p wave gap integral
*/
TGapPointPWave::TGapPointPWave() {
TPointPWaveGapIntegralCuhre *gapint = new TPointPWaveGapIntegralCuhre();
fGapIntegral = gapint;
gapint = nullptr;
fTemp.clear();
fTempIter = fTemp.end();
fIntegralValues.clear();
fCalcNeeded.clear();
fPar.clear();
}
//--------------------------------------------------------------------
/**
* <p> line p wave gap integral
*/
TGapLinePWave::TGapLinePWave() {
TLinePWaveGapIntegralCuhre *gapint = new TLinePWaveGapIntegralCuhre();
fGapIntegral = gapint;
gapint = nullptr;
fTemp.clear();
fTempIter = fTemp.end();
fIntegralValues.clear();
fCalcNeeded.clear();
fPar.clear();
}
//--------------------------------------------------------------------
/**
* <p>
@ -181,6 +219,22 @@ TLambdaSWave::TLambdaSWave() {
fLambdaInvSq = new TGapSWave();
}
//--------------------------------------------------------------------
/**
* <p>
*/
TLambdaPointPWave::TLambdaPointPWave() {
fLambdaInvSq = new TGapPointPWave();
}
//--------------------------------------------------------------------
/**
* <p>
*/
TLambdaLinePWave::TLambdaLinePWave() {
fLambdaInvSq = new TGapLinePWave();
}
//--------------------------------------------------------------------
/**
* <p>
@ -221,6 +275,22 @@ TLambdaInvSWave::TLambdaInvSWave() {
fLambdaInvSq = new TGapSWave();
}
//--------------------------------------------------------------------
/**
* <p>
*/
TLambdaInvPointPWave::TLambdaInvPointPWave() {
fLambdaInvSq = new TGapPointPWave();
}
//--------------------------------------------------------------------
/**
* <p>
*/
TLambdaInvLinePWave::TLambdaInvLinePWave() {
fLambdaInvSq = new TGapLinePWave();
}
//--------------------------------------------------------------------
/**
* <p>
@ -268,6 +338,36 @@ TGapSWave::~TGapSWave() {
fPar.clear();
}
//--------------------------------------------------------------------
/**
* <p>
*/
TGapPointPWave::~TGapPointPWave() {
delete fGapIntegral;
fGapIntegral = nullptr;
fTemp.clear();
fTempIter = fTemp.end();
fIntegralValues.clear();
fCalcNeeded.clear();
fPar.clear();
}
//--------------------------------------------------------------------
/**
* <p>
*/
TGapLinePWave::~TGapLinePWave() {
delete fGapIntegral;
fGapIntegral = nullptr;
fTemp.clear();
fTempIter = fTemp.end();
fIntegralValues.clear();
fCalcNeeded.clear();
fPar.clear();
}
//--------------------------------------------------------------------
/**
* <p>
@ -367,6 +467,24 @@ TLambdaSWave::~TLambdaSWave() {
fLambdaInvSq = nullptr;
}
//--------------------------------------------------------------------
/**
* <p>
*/
TLambdaPointPWave::~TLambdaPointPWave() {
delete fLambdaInvSq;
fLambdaInvSq = nullptr;
}
//--------------------------------------------------------------------
/**
* <p>
*/
TLambdaLinePWave::~TLambdaLinePWave() {
delete fLambdaInvSq;
fLambdaInvSq = nullptr;
}
//--------------------------------------------------------------------
/**
* <p>
@ -412,6 +530,24 @@ TLambdaInvSWave::~TLambdaInvSWave() {
fLambdaInvSq = nullptr;
}
//--------------------------------------------------------------------
/**
* <p>
*/
TLambdaInvPointPWave::~TLambdaInvPointPWave() {
delete fLambdaInvSq;
fLambdaInvSq = nullptr;
}
//--------------------------------------------------------------------
/**
* <p>
*/
TLambdaInvLinePWave::~TLambdaInvLinePWave() {
delete fLambdaInvSq;
fLambdaInvSq = nullptr;
}
//--------------------------------------------------------------------
/**
* <p>
@ -528,6 +664,162 @@ double TGapSWave::operator()(double t, const std::vector<double> &par) const {
}
//--------------------------------------------------------------------
/**
* <p>prepare the needed parameters for the integration carried out in TPointPWaveGapIntegralCuhre.
* For details see also the Memo GapIntegrals.pdf, , especially Eq.(19) and (20).
*/
double TGapPointPWave::operator()(double t, const std::vector<double> &par) const {
// parameters: [0] Tc (K), [1] Delta0 (meV), [[2] c0 (1), [3] aG (1)]
assert((par.size() == 2) || (par.size() == 4)); // 2 parameters: see A.~Carrington and F.~Manzano, Physica~C~\textbf{385}~(2003)~205
// 4 parameters: see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67
// and Erratum Supercond. Sci. Technol. 21 (2008) 082003
// c0 in the original context is c0 = (pi kB Tc) / Delta0
if (t <= 0.0)
return 1.0;
else if (t >= par[0])
return 0.0;
bool integralParChanged(false);
if (fPar.empty()) { // first time calling this routine
fPar = par;
integralParChanged = true;
} else { // check if Tc or Delta0 have changed
for (unsigned int i(0); i<par.size(); i++) {
if (par[i] != fPar[i]) {
fPar[i] = par[i];
integralParChanged = true;
}
}
}
bool newTemp(false);
unsigned int vectorIndex;
if (integralParChanged) {
fCalcNeeded.clear();
fCalcNeeded.resize(fTemp.size(), true);
}
fTempIter = find(fTemp.begin(), fTemp.end(), t);
if(fTempIter == fTemp.end()) {
fTemp.push_back(t);
vectorIndex = fTemp.size() - 1;
fCalcNeeded.push_back(true);
newTemp = true;
} else {
vectorIndex = fTempIter - fTemp.begin();
}
if (fCalcNeeded[vectorIndex]) {
double ds;
std::vector<double> intPar; // parameters for the integral, T & Delta(T)
intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K
if (par.size() == 2) { // Carrington/Manzano
intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
} else { // Prozorov/Giannetta
intPar.push_back(par[1]*tanh(par[2]*sqrt(par[3]*(par[0]/t-1.0)))); // Delta0*tanh(c0*sqrt(aG*(Tc/T-1)))
}
intPar.push_back(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy
intPar.push_back(TMath::Pi()); // upper limit of theta-integration
fGapIntegral->SetParameters(intPar);
ds = 1.0-(intPar[2]*PI)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc();
intPar.clear();
if (newTemp)
fIntegralValues.push_back(ds);
else
fIntegralValues[vectorIndex] = ds;
fCalcNeeded[vectorIndex] = false;
}
return fIntegralValues[vectorIndex];
}
//--------------------------------------------------------------------
/**
* <p>prepare the needed parameters for the integration carried out in TLinePWaveGapIntegralCuhre.
* For details see also the Memo GapIntegrals.pdf, , especially Eq.(19) and (20).
*/
double TGapLinePWave::operator()(double t, const std::vector<double> &par) const {
// parameters: [0] Tc (K), [1] Delta0 (meV), [[2] c0 (1), [3] aG (1)]
assert((par.size() == 2) || (par.size() == 4)); // 2 parameters: see A.~Carrington and F.~Manzano, Physica~C~\textbf{385}~(2003)~205
// 4 parameters: see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67
// and Erratum Supercond. Sci. Technol. 21 (2008) 082003
// c0 in the original context is c0 = (pi kB Tc) / Delta0
if (t <= 0.0)
return 1.0;
else if (t >= par[0])
return 0.0;
bool integralParChanged(false);
if (fPar.empty()) { // first time calling this routine
fPar = par;
integralParChanged = true;
} else { // check if Tc or Delta0 have changed
for (unsigned int i(0); i<par.size(); i++) {
if (par[i] != fPar[i]) {
fPar[i] = par[i];
integralParChanged = true;
}
}
}
bool newTemp(false);
unsigned int vectorIndex;
if (integralParChanged) {
fCalcNeeded.clear();
fCalcNeeded.resize(fTemp.size(), true);
}
fTempIter = find(fTemp.begin(), fTemp.end(), t);
if(fTempIter == fTemp.end()) {
fTemp.push_back(t);
vectorIndex = fTemp.size() - 1;
fCalcNeeded.push_back(true);
newTemp = true;
} else {
vectorIndex = fTempIter - fTemp.begin();
}
if (fCalcNeeded[vectorIndex]) {
double ds;
std::vector<double> intPar; // parameters for the integral, T & Delta(T)
intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K
if (par.size() == 2) { // Carrington/Manzano
intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
} else { // Prozorov/Giannetta
intPar.push_back(par[1]*tanh(par[2]*sqrt(par[3]*(par[0]/t-1.0)))); // Delta0*tanh(c0*sqrt(aG*(Tc/T-1)))
}
intPar.push_back(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy
intPar.push_back(TMath::Pi()); // upper limit of theta-integration
fGapIntegral->SetParameters(intPar);
ds = 1.0-(intPar[2]*PI)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc();
intPar.clear();
if (newTemp)
fIntegralValues.push_back(ds);
else
fIntegralValues[vectorIndex] = ds;
fCalcNeeded[vectorIndex] = false;
}
return fIntegralValues[vectorIndex];
}
//--------------------------------------------------------------------
/**
* <p>prepare the needed parameters for the integration carried out in TDWaveGapIntegralCuhre.
@ -609,7 +901,6 @@ double TGapDWave::operator()(double t, const std::vector<double> &par) const {
}
return fIntegralValues[vectorIndex];
}
//--------------------------------------------------------------------
@ -1103,7 +1394,40 @@ double TLambdaSWave::operator()(double t, const std::vector<double> &par) const
return 1.0;
return 1.0/sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
/**
* <p>
*/
double TLambdaPointPWave::operator()(double t, const std::vector<double> &par) const
{
assert(par.size() == 2); // two parameters: Tc, Delta0
if (t >= par[0])
return -1.0;
if (t <= 0.0)
return 1.0;
return 1.0/sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
/**
* <p>
*/
double TLambdaLinePWave::operator()(double t, const std::vector<double> &par) const
{
assert(par.size() == 2); // two parameters: Tc, Delta0
if (t >= par[0])
return -1.0;
if (t <= 0.0)
return 1.0;
return 1.0/sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
@ -1121,7 +1445,6 @@ double TLambdaDWave::operator()(double t, const std::vector<double> &par) const
return 1.0;
return 1.0/sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
@ -1139,7 +1462,6 @@ double TLambdaAnSWave::operator()(double t, const std::vector<double> &par) cons
return 1.0;
return 1.0/sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
@ -1157,7 +1479,6 @@ double TLambdaNonMonDWave1::operator()(double t, const std::vector<double> &par)
return 1.0;
return 1.0/sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
@ -1175,7 +1496,6 @@ double TLambdaNonMonDWave2::operator()(double t, const std::vector<double> &par)
return 1.0;
return 1.0/sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
@ -1192,7 +1512,6 @@ double TLambdaPowerLaw::operator()(double t, const std::vector<double> &par) con
return -1.0;
return 1.0/sqrt(1.0 - pow(t/par[0], par[1]));
}
@ -1211,7 +1530,40 @@ double TLambdaInvSWave::operator()(double t, const std::vector<double> &par) con
return 1.0;
return sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
/**
* <p>
*/
double TLambdaInvPointPWave::operator()(double t, const std::vector<double> &par) const
{
assert(par.size() == 2); // two parameters: Tc, Delta0
if (t >= par[0])
return 0.0;
if (t <= 0.0)
return 1.0;
return sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
/**
* <p>
*/
double TLambdaInvLinePWave::operator()(double t, const std::vector<double> &par) const
{
assert(par.size() == 2); // two parameters: Tc, Delta0
if (t >= par[0])
return 0.0;
if (t <= 0.0)
return 1.0;
return sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
@ -1229,7 +1581,6 @@ double TLambdaInvDWave::operator()(double t, const std::vector<double> &par) con
return 1.0;
return sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
@ -1247,7 +1598,6 @@ double TLambdaInvAnSWave::operator()(double t, const std::vector<double> &par) c
return 1.0;
return sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
@ -1265,7 +1615,6 @@ double TLambdaInvNonMonDWave1::operator()(double t, const std::vector<double> &p
return 1.0;
return sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
@ -1283,7 +1632,6 @@ double TLambdaInvNonMonDWave2::operator()(double t, const std::vector<double> &p
return 1.0;
return sqrt((*fLambdaInvSq)(t, par));
}
//--------------------------------------------------------------------
@ -1300,7 +1648,6 @@ double TLambdaInvPowerLaw::operator()(double t, const std::vector<double> &par)
return 0.0;
return sqrt(1.0 - pow(t/par[0], par[1]));
}
//--------------------------------------------------------------------

144
src/external/libGapIntegrals/TGapIntegrals.h vendored
View File

@ -63,6 +63,62 @@ private:
ClassDef(TGapSWave,1)
};
//--------------------------------------------------------------------
/**
* <p>
*/
class TGapPointPWave : public PUserFcnBase {
public:
TGapPointPWave();
virtual ~TGapPointPWave();
virtual Bool_t NeedGlobalPart() const { return false; }
virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
virtual Bool_t GlobalPartIsValid() const { return true; }
double operator()(double, const std::vector<double>&) const;
private:
TPointPWaveGapIntegralCuhre *fGapIntegral;
mutable std::vector<double> fTemp;
mutable std::vector<double>::const_iterator fTempIter;
mutable std::vector<double> fIntegralValues;
mutable std::vector<bool> fCalcNeeded;
mutable std::vector<double> fPar;
ClassDef(TGapPointPWave,1)
};
//--------------------------------------------------------------------
/**
* <p>
*/
class TGapLinePWave : public PUserFcnBase {
public:
TGapLinePWave();
virtual ~TGapLinePWave();
virtual Bool_t NeedGlobalPart() const { return false; }
virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
virtual Bool_t GlobalPartIsValid() const { return true; }
double operator()(double, const std::vector<double>&) const;
private:
TLinePWaveGapIntegralCuhre *fGapIntegral;
mutable std::vector<double> fTemp;
mutable std::vector<double>::const_iterator fTempIter;
mutable std::vector<double> fIntegralValues;
mutable std::vector<bool> fCalcNeeded;
mutable std::vector<double> fPar;
ClassDef(TGapLinePWave,1)
};
//--------------------------------------------------------------------
/**
* <p>
@ -297,6 +353,50 @@ private:
ClassDef(TLambdaSWave,1)
};
//--------------------------------------------------------------------
/**
* <p>
*/
class TLambdaPointPWave : public PUserFcnBase {
public:
TLambdaPointPWave();
virtual ~TLambdaPointPWave();
virtual Bool_t NeedGlobalPart() const { return false; }
virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
virtual Bool_t GlobalPartIsValid() const { return true; }
double operator()(double, const std::vector<double>&) const;
private:
TGapPointPWave *fLambdaInvSq;
ClassDef(TLambdaPointPWave,1)
};
//--------------------------------------------------------------------
/**
* <p>
*/
class TLambdaLinePWave : public PUserFcnBase {
public:
TLambdaLinePWave();
virtual ~TLambdaLinePWave();
virtual Bool_t NeedGlobalPart() const { return false; }
virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
virtual Bool_t GlobalPartIsValid() const { return true; }
double operator()(double, const std::vector<double>&) const;
private:
TGapLinePWave *fLambdaInvSq;
ClassDef(TLambdaLinePWave,1)
};
//--------------------------------------------------------------------
/**
* <p>
@ -428,6 +528,50 @@ private:
ClassDef(TLambdaInvSWave,1)
};
//--------------------------------------------------------------------
/**
* <p>
*/
class TLambdaInvPointPWave : public PUserFcnBase {
public:
TLambdaInvPointPWave();
virtual ~TLambdaInvPointPWave();
virtual Bool_t NeedGlobalPart() const { return false; }
virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
virtual Bool_t GlobalPartIsValid() const { return true; }
double operator()(double, const std::vector<double>&) const;
private:
TGapPointPWave *fLambdaInvSq;
ClassDef(TLambdaInvPointPWave,1)
};
//--------------------------------------------------------------------
/**
* <p>
*/
class TLambdaInvLinePWave : public PUserFcnBase {
public:
TLambdaInvLinePWave();
virtual ~TLambdaInvLinePWave();
virtual Bool_t NeedGlobalPart() const { return false; }
virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { }
virtual Bool_t GlobalPartIsValid() const { return true; }
double operator()(double, const std::vector<double>&) const;
private:
TGapLinePWave *fLambdaInvSq;
ClassDef(TLambdaInvLinePWave,1)
};
//--------------------------------------------------------------------
/**
* <p>

6
src/external/libGapIntegrals/TGapIntegralsLinkDef.h vendored
View File

@ -35,6 +35,8 @@
#pragma link off all functions;
#pragma link C++ class TGapSWave+;
#pragma link C++ class TGapPointPWave+;
#pragma link C++ class TGapLinePWave+;
#pragma link C++ class TGapDWave+;
#pragma link C++ class TGapCosSqDWave+;
#pragma link C++ class TGapSinSqDWave+;
@ -44,12 +46,16 @@
#pragma link C++ class TGapPowerLaw+;
#pragma link C++ class TGapDirtySWave+;
#pragma link C++ class TLambdaSWave+;
#pragma link C++ class TLambdaPointPWave+;
#pragma link C++ class TLambdaLinePWave+;
#pragma link C++ class TLambdaDWave+;
#pragma link C++ class TLambdaAnSWave+;
#pragma link C++ class TLambdaNonMonDWave1+;
#pragma link C++ class TLambdaNonMonDWave2+;
#pragma link C++ class TLambdaPowerLaw+;
#pragma link C++ class TLambdaInvSWave+;
#pragma link C++ class TLambdaInvPointPWave+;
#pragma link C++ class TLambdaInvLinePWave+;
#pragma link C++ class TLambdaInvDWave+;
#pragma link C++ class TLambdaInvAnSWave+;
#pragma link C++ class TLambdaInvNonMonDWave1+;