diff --git a/src/external/BMWtools/BMWIntegrator.cpp b/src/external/BMWtools/BMWIntegrator.cpp index 96ce4154..fb8ead3a 100644 --- a/src/external/BMWtools/BMWIntegrator.cpp +++ b/src/external/BMWtools/BMWIntegrator.cpp @@ -35,8 +35,127 @@ #define SEED 0 #define STATEFILE NULL + +//----------------------------------------------------------------------------- + +std::vector TPointPWaveGapIntegralCuhre::fPar; + +//----------------------------------------------------------------------------- +/** + *

Integrate the function using the Cuhre interface + * + *

return: + * - 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]; +} + +//----------------------------------------------------------------------------- +/** + *

Calculate the function value for the use with Cuhre---actual implementation of the function + * + *

return: + * - 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 TLinePWaveGapIntegralCuhre::fPar; + +//----------------------------------------------------------------------------- +/** + *

Integrate the function using the Cuhre interface + * + *

return: + * - 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]; +} + +//----------------------------------------------------------------------------- +/** + *

Calculate the function value for the use with Cuhre---actual implementation of the function + * + *

return: + * - 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 TDWaveGapIntegralCuhre::fPar; +//----------------------------------------------------------------------------- /** *

Integrate the function using the Cuhre interface * @@ -67,6 +186,7 @@ double TDWaveGapIntegralCuhre::IntegrateFunc() return integral[0]; } +//----------------------------------------------------------------------------- /** *

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 TCosSqDWaveGapIntegralCuhre::fPar; +//----------------------------------------------------------------------------- /** *

Integrate the function using the Cuhre interface * @@ -119,6 +242,7 @@ double TCosSqDWaveGapIntegralCuhre::IntegrateFunc() return integral[0]; } +//----------------------------------------------------------------------------- /** *

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 TSinSqDWaveGapIntegralCuhre::fPar; +//----------------------------------------------------------------------------- /** *

Integrate the function using the Cuhre interface * @@ -171,6 +298,7 @@ double TSinSqDWaveGapIntegralCuhre::IntegrateFunc() return integral[0]; } +//----------------------------------------------------------------------------- /** *

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 TAnSWaveGapIntegralCuhre::fPar; +//----------------------------------------------------------------------------- /** *

Integrate the function using the Cuhre interface * @@ -223,6 +354,7 @@ double TAnSWaveGapIntegralCuhre::IntegrateFunc() return integral[0]; } +//----------------------------------------------------------------------------- /** *

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 TAnSWaveGapIntegralDivonne::fPar; +//----------------------------------------------------------------------------- /** *

Integrate the function using the Divonne interface * @@ -283,6 +418,7 @@ double TAnSWaveGapIntegralDivonne::IntegrateFunc() return integral[0]; } +//----------------------------------------------------------------------------- /** *

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 TAnSWaveGapIntegralSuave::fPar; +//----------------------------------------------------------------------------- /** *

Integrate the function using the Suave interface * @@ -337,6 +476,7 @@ double TAnSWaveGapIntegralSuave::IntegrateFunc() return integral[0]; } +//----------------------------------------------------------------------------- /** *

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 TNonMonDWave1GapIntegralCuhre::fPar; +//----------------------------------------------------------------------------- /** *

Integrate the function using the Cuhre interface * @@ -389,6 +532,7 @@ double TNonMonDWave1GapIntegralCuhre::IntegrateFunc() return integral[0]; } +//----------------------------------------------------------------------------- /** *

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 TNonMonDWave2GapIntegralCuhre::fPar; +//----------------------------------------------------------------------------- /** *

Integrate the function using the Cuhre interface * @@ -441,6 +588,7 @@ double TNonMonDWave2GapIntegralCuhre::IntegrateFunc() return integral[0]; } +//----------------------------------------------------------------------------- /** *

Calculate the function value for the use with Cuhre---actual implementation of the function * diff --git a/src/external/BMWtools/BMWIntegrator.h b/src/external/BMWtools/BMWIntegrator.h index a901004a..5ecb7406 100644 --- a/src/external/BMWtools/BMWIntegrator.h +++ b/src/external/BMWtools/BMWIntegrator.h @@ -36,6 +36,7 @@ #include #include +//----------------------------------------------------------------------------- /** *

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 *); }; +//----------------------------------------------------------------------------- /** *

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)); } +//----------------------------------------------------------------------------- /** *

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; } - - - - - - - - - - - - - - +//----------------------------------------------------------------------------- /** *

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 }; +//----------------------------------------------------------------------------- /** *

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); } +//----------------------------------------------------------------------------- /** *

Destructor of the base class for 1D integrations * Clean up. @@ -148,6 +140,7 @@ inline TIntegrator::~TIntegrator(){ fFunc=0; } +//----------------------------------------------------------------------------- /** *

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); } +//----------------------------------------------------------------------------- /** *

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); } +//----------------------------------------------------------------------------- /** *

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 }; +//----------------------------------------------------------------------------- /** *

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); } +//----------------------------------------------------------------------------- /** *

Destructor of the base class for 1D integrations * Clean up. @@ -218,6 +215,7 @@ inline TMCIntegrator::~TMCIntegrator(){ fFunc=0; } +//----------------------------------------------------------------------------- /** *

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); } +//----------------------------------------------------------------------------- /** *

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)); } +//----------------------------------------------------------------------------- +/** + *

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 &par) { fPar=par; } + static int Integrand(const int*, const double[], const int*, double[], void*); + double IntegrateFunc(); + + protected: + static std::vector fPar; ///< parameters of the integrand + unsigned int fNDim; ///< dimension of the integral +}; + +//----------------------------------------------------------------------------- +/** + *

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 &par) { fPar=par; } + static int Integrand(const int*, const double[], const int*, double[], void*); + double IntegrateFunc(); + + protected: + static std::vector fPar; ///< parameters of the integrand + unsigned int fNDim; ///< dimension of the integral +}; + +//----------------------------------------------------------------------------- /** *

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 }; +//----------------------------------------------------------------------------- /** *

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 }; +//----------------------------------------------------------------------------- /** *

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 }; +//----------------------------------------------------------------------------- /** *

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 }; +//----------------------------------------------------------------------------- /** *

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 }; +//----------------------------------------------------------------------------- /** *

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 }; +//----------------------------------------------------------------------------- /** *

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 }; +//----------------------------------------------------------------------------- /** *

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 }; +//----------------------------------------------------------------------------- /** *

Test class for the 2D MC integration * Integral: x*y dx dy @@ -406,6 +452,7 @@ class T2DTest : public TMCIntegrator { double FuncAtX(double *) const; }; +//----------------------------------------------------------------------------- /** *

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]; } +//----------------------------------------------------------------------------- +/** + *

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; +}; + +//----------------------------------------------------------------------------- +/** + *

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; +}; + +//----------------------------------------------------------------------------- +/** + *

Calculate the function value---actual implementation of the function + * + *

return: + * - 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)); +} + +//----------------------------------------------------------------------------- +/** + *

Calculate the function value---actual implementation of the function + * + *

return: + * - 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)); +} + +//----------------------------------------------------------------------------- /** *

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; }; +//----------------------------------------------------------------------------- /** *

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)); } +//----------------------------------------------------------------------------- /** *

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; }; +//----------------------------------------------------------------------------- /** *

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)); } +//----------------------------------------------------------------------------- /** *

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&) const; }; +//----------------------------------------------------------------------------- /** *

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 &par) return j0 * TMath::Exp(-par[1]*x); } +//----------------------------------------------------------------------------- /** *

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&) const; }; +//----------------------------------------------------------------------------- /** *

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 &par) cons return TMath::Sin(TMath::TwoPi()*par[0]*x) * TMath::Exp(-0.5*par[1]*par[1]*x*x); } +//----------------------------------------------------------------------------- /** *

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&) const; }; +//----------------------------------------------------------------------------- /** *

Calculate the function value---actual implementation of the function * @@ -563,6 +678,7 @@ inline double TIntSGInterpolation::FuncAtX(double x, const std::vector & return (wt*TMath::Cos(wt)-TMath::Sin(wt))/(wt*wt)*TMath::Exp(-TMath::Power(expo,par[3]))/TMath::Power(expo,(1.0-par[3])); } +//----------------------------------------------------------------------------- /** *

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 }; +//----------------------------------------------------------------------------- /** *

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)); } +//----------------------------------------------------------------------------- /** *

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&) const; // variable: x }; +//----------------------------------------------------------------------------- /** *

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::vectorClass 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&) const; // variable: x }; +//----------------------------------------------------------------------------- /** *

Calculate the function value---actual implementation of the integrand in Eq. (7) of Solt's article * diff --git a/src/external/libGapIntegrals/GapIntegrals.pdf b/src/external/libGapIntegrals/GapIntegrals.pdf index b03edf6c..ee36ca2c 100644 Binary files a/src/external/libGapIntegrals/GapIntegrals.pdf and b/src/external/libGapIntegrals/GapIntegrals.pdf differ diff --git a/src/external/libGapIntegrals/GapIntegrals.tex b/src/external/libGapIntegrals/GapIntegrals.tex index 75812971..29c1180d 100644 --- a/src/external/libGapIntegrals/GapIntegrals.tex +++ b/src/external/libGapIntegrals/GapIntegrals.tex @@ -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 diff --git a/src/external/libGapIntegrals/TGapIntegrals.cpp b/src/external/libGapIntegrals/TGapIntegrals.cpp index 702f1589..b90ee836 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.cpp +++ b/src/external/libGapIntegrals/TGapIntegrals.cpp @@ -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(); } +//-------------------------------------------------------------------- +/** + *

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(); +} + +//-------------------------------------------------------------------- +/** + *

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(); +} + //-------------------------------------------------------------------- /** *

@@ -181,6 +219,22 @@ TLambdaSWave::TLambdaSWave() { fLambdaInvSq = new TGapSWave(); } +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaPointPWave::TLambdaPointPWave() { + fLambdaInvSq = new TGapPointPWave(); +} + +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaLinePWave::TLambdaLinePWave() { + fLambdaInvSq = new TGapLinePWave(); +} + //-------------------------------------------------------------------- /** *

@@ -221,6 +275,22 @@ TLambdaInvSWave::TLambdaInvSWave() { fLambdaInvSq = new TGapSWave(); } +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaInvPointPWave::TLambdaInvPointPWave() { + fLambdaInvSq = new TGapPointPWave(); +} + +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaInvLinePWave::TLambdaInvLinePWave() { + fLambdaInvSq = new TGapLinePWave(); +} + //-------------------------------------------------------------------- /** *

@@ -268,6 +338,36 @@ TGapSWave::~TGapSWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ +TGapPointPWave::~TGapPointPWave() { + delete fGapIntegral; + fGapIntegral = nullptr; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + +//-------------------------------------------------------------------- +/** + *

+ */ +TGapLinePWave::~TGapLinePWave() { + delete fGapIntegral; + fGapIntegral = nullptr; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + //-------------------------------------------------------------------- /** *

@@ -367,6 +467,24 @@ TLambdaSWave::~TLambdaSWave() { fLambdaInvSq = nullptr; } +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaPointPWave::~TLambdaPointPWave() { + delete fLambdaInvSq; + fLambdaInvSq = nullptr; +} + +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaLinePWave::~TLambdaLinePWave() { + delete fLambdaInvSq; + fLambdaInvSq = nullptr; +} + //-------------------------------------------------------------------- /** *

@@ -412,6 +530,24 @@ TLambdaInvSWave::~TLambdaInvSWave() { fLambdaInvSq = nullptr; } +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaInvPointPWave::~TLambdaInvPointPWave() { + delete fLambdaInvSq; + fLambdaInvSq = nullptr; +} + +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaInvLinePWave::~TLambdaInvLinePWave() { + delete fLambdaInvSq; + fLambdaInvSq = nullptr; +} + //-------------------------------------------------------------------- /** *

@@ -528,6 +664,162 @@ double TGapSWave::operator()(double t, const std::vector &par) const { } +//-------------------------------------------------------------------- +/** + *

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 &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 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]; +} + +//-------------------------------------------------------------------- +/** + *

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 &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 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]; +} + //-------------------------------------------------------------------- /** *

prepare the needed parameters for the integration carried out in TDWaveGapIntegralCuhre. @@ -609,7 +901,6 @@ double TGapDWave::operator()(double t, const std::vector &par) const { } return fIntegralValues[vectorIndex]; - } //-------------------------------------------------------------------- @@ -1103,7 +1394,40 @@ double TLambdaSWave::operator()(double t, const std::vector &par) const return 1.0; return 1.0/sqrt((*fLambdaInvSq)(t, par)); +} +//-------------------------------------------------------------------- +/** + *

+ */ +double TLambdaPointPWave::operator()(double t, const std::vector &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)); +} + +//-------------------------------------------------------------------- +/** + *

+ */ +double TLambdaLinePWave::operator()(double t, const std::vector &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 &par) const return 1.0; return 1.0/sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1139,7 +1462,6 @@ double TLambdaAnSWave::operator()(double t, const std::vector &par) cons return 1.0; return 1.0/sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1157,7 +1479,6 @@ double TLambdaNonMonDWave1::operator()(double t, const std::vector &par) return 1.0; return 1.0/sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1175,7 +1496,6 @@ double TLambdaNonMonDWave2::operator()(double t, const std::vector &par) return 1.0; return 1.0/sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1192,7 +1512,6 @@ double TLambdaPowerLaw::operator()(double t, const std::vector &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 &par) con return 1.0; return sqrt((*fLambdaInvSq)(t, par)); +} +//-------------------------------------------------------------------- +/** + *

+ */ +double TLambdaInvPointPWave::operator()(double t, const std::vector &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)); +} + +//-------------------------------------------------------------------- +/** + *

+ */ +double TLambdaInvLinePWave::operator()(double t, const std::vector &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 &par) con return 1.0; return sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1247,7 +1598,6 @@ double TLambdaInvAnSWave::operator()(double t, const std::vector &par) c return 1.0; return sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1265,7 +1615,6 @@ double TLambdaInvNonMonDWave1::operator()(double t, const std::vector &p return 1.0; return sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1283,7 +1632,6 @@ double TLambdaInvNonMonDWave2::operator()(double t, const std::vector &p return 1.0; return sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1300,7 +1648,6 @@ double TLambdaInvPowerLaw::operator()(double t, const std::vector &par) return 0.0; return sqrt(1.0 - pow(t/par[0], par[1])); - } //-------------------------------------------------------------------- diff --git a/src/external/libGapIntegrals/TGapIntegrals.h b/src/external/libGapIntegrals/TGapIntegrals.h index a0d344ee..07c498be 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.h +++ b/src/external/libGapIntegrals/TGapIntegrals.h @@ -63,6 +63,62 @@ private: ClassDef(TGapSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ +class TGapPointPWave : public PUserFcnBase { + +public: + TGapPointPWave(); + virtual ~TGapPointPWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TPointPWaveGapIntegralCuhre *fGapIntegral; + mutable std::vector fTemp; + mutable std::vector::const_iterator fTempIter; + mutable std::vector fIntegralValues; + mutable std::vector fCalcNeeded; + + mutable std::vector fPar; + + ClassDef(TGapPointPWave,1) +}; + +//-------------------------------------------------------------------- +/** + *

+ */ +class TGapLinePWave : public PUserFcnBase { + +public: + TGapLinePWave(); + virtual ~TGapLinePWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TLinePWaveGapIntegralCuhre *fGapIntegral; + mutable std::vector fTemp; + mutable std::vector::const_iterator fTempIter; + mutable std::vector fIntegralValues; + mutable std::vector fCalcNeeded; + + mutable std::vector fPar; + + ClassDef(TGapLinePWave,1) +}; + //-------------------------------------------------------------------- /** *

@@ -297,6 +353,50 @@ private: ClassDef(TLambdaSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ +class TLambdaPointPWave : public PUserFcnBase { + +public: + TLambdaPointPWave(); + virtual ~TLambdaPointPWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TGapPointPWave *fLambdaInvSq; + + ClassDef(TLambdaPointPWave,1) +}; + +//-------------------------------------------------------------------- +/** + *

+ */ +class TLambdaLinePWave : public PUserFcnBase { + +public: + TLambdaLinePWave(); + virtual ~TLambdaLinePWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TGapLinePWave *fLambdaInvSq; + + ClassDef(TLambdaLinePWave,1) +}; + //-------------------------------------------------------------------- /** *

@@ -428,6 +528,50 @@ private: ClassDef(TLambdaInvSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ +class TLambdaInvPointPWave : public PUserFcnBase { + +public: + TLambdaInvPointPWave(); + virtual ~TLambdaInvPointPWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TGapPointPWave *fLambdaInvSq; + + ClassDef(TLambdaInvPointPWave,1) +}; + +//-------------------------------------------------------------------- +/** + *

+ */ +class TLambdaInvLinePWave : public PUserFcnBase { + +public: + TLambdaInvLinePWave(); + virtual ~TLambdaInvLinePWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TGapLinePWave *fLambdaInvSq; + + ClassDef(TLambdaInvLinePWave,1) +}; + //-------------------------------------------------------------------- /** *

diff --git a/src/external/libGapIntegrals/TGapIntegralsLinkDef.h b/src/external/libGapIntegrals/TGapIntegralsLinkDef.h index 1deab362..03f17895 100644 --- a/src/external/libGapIntegrals/TGapIntegralsLinkDef.h +++ b/src/external/libGapIntegrals/TGapIntegralsLinkDef.h @@ -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+;