From ae5e2761565d720e28c9b14c41cda205eada5874 Mon Sep 17 00:00:00 2001 From: "Bastian M. Wojek" Date: Wed, 9 Sep 2009 16:12:01 +0000 Subject: [PATCH] Removed a redundant check in msr2data --- .../libGapIntegrals/TGapIntegrals.cpp | 18 ++--- src/external/libGapIntegrals/TGapIntegrals.h | 2 +- src/external/libGapIntegrals/TIntegrator.cpp | 65 +++++++++++++------ src/external/libGapIntegrals/TIntegrator.h | 14 ++++ src/msr2data.cpp | 4 -- 5 files changed, 70 insertions(+), 33 deletions(-) diff --git a/src/external/libGapIntegrals/TGapIntegrals.cpp b/src/external/libGapIntegrals/TGapIntegrals.cpp index 6eaf85de..e3c5a5cb 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.cpp +++ b/src/external/libGapIntegrals/TGapIntegrals.cpp @@ -74,7 +74,7 @@ TGapDWave::TGapDWave() { } TGapAnSWave::TGapAnSWave() { - TAnSWaveGapIntegralDivonne *gapint = new TAnSWaveGapIntegralDivonne(); + TAnSWaveGapIntegralSuave *gapint = new TAnSWaveGapIntegralSuave(); fGapIntegral = gapint; gapint = 0; delete gapint; @@ -230,16 +230,17 @@ double TGapDWave::operator()(double t, const vector &par) const { double ds; vector intPar; // parameters for the integral, T & Delta(T) - intPar.push_back(t); + intPar.push_back(2.0*0.08617384436*t); // 2 kB T, kB in meV/K intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51))); - + intPar.push_back(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy + intPar.push_back(TMath::PiOver2()); // upper limit of phi-integration // double xl[] = {0.0, 0.0}; // lower bound E, phi -// double xu[] = {2.0*(t+intPar[1]), 0.5*PI}; // upper bound E, phi +// double xu[] = {4.0*(t+intPar[1]), 0.5*PI}; // upper bound E, phi fGapIntegral->SetParameters(intPar); // ds = 1.0+4.0/PI*fGapIntegral->IntegrateFunc(2, xl, xu); - ds = 1.0+4.0/PI*fGapIntegral->IntegrateFunc(); + ds = 1.0-intPar[2]/intPar[0]*fGapIntegral->IntegrateFunc(); intPar.clear(); @@ -300,17 +301,18 @@ double TGapAnSWave::operator()(double t, const vector &par) const { double ds; vector intPar; // parameters for the integral, T & Delta(T) - intPar.push_back(t); + intPar.push_back(2.0*0.08617384436*t); // 2 kB T, kB in meV/K intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51))); intPar.push_back(par[2]); - + intPar.push_back(4.0*(t+(1.0+par[2])*intPar[1])); // upper limit of energy-integration: cutoff energy + intPar.push_back(TMath::PiOver2()); // upper limit of phi-integration // double xl[] = {0.0, 0.0}; // lower bound E, phi // double xu[] = {4.0*(t+intPar[1]), 0.5*PI}; // upper bound E, phi fGapIntegral->SetParameters(intPar); // ds = 1.0+4.0/PI*fGapIntegral->IntegrateFunc(2, xl, xu); - ds = 1.0+4.0/PI*fGapIntegral->IntegrateFunc(); + ds = 1.0-intPar[3]/intPar[0]*fGapIntegral->IntegrateFunc(); intPar.clear(); diff --git a/src/external/libGapIntegrals/TGapIntegrals.h b/src/external/libGapIntegrals/TGapIntegrals.h index 897a330b..d49093ca 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.h +++ b/src/external/libGapIntegrals/TGapIntegrals.h @@ -89,7 +89,7 @@ public: double operator()(double, const vector&) const; private: - TAnSWaveGapIntegralDivonne *fGapIntegral; + TAnSWaveGapIntegralSuave *fGapIntegral; mutable vector fTemp; mutable vector::const_iterator fTempIter; mutable vector fIntegralValues; diff --git a/src/external/libGapIntegrals/TIntegrator.cpp b/src/external/libGapIntegrals/TIntegrator.cpp index 4a22f36a..40364bf1 100644 --- a/src/external/libGapIntegrals/TIntegrator.cpp +++ b/src/external/libGapIntegrals/TIntegrator.cpp @@ -59,13 +59,10 @@ double TDWaveGapIntegralCuhre::IntegrateFunc() } void TDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[], - const int *ncomp, double f[]) // x = {E, phi}, fPar = {T, Delta(T)} + const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T), Ec, phic} { - double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K - double Ec(4.0*(fPar[0]+fPar[1])); // upper limit of energy-integration: cutoff energy - double phic(TMath::PiOver2()); // upper limit of phi-integration - double deltasq(TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]*phic),2.0)); - f[0] = -phic*Ec/(2.0*twokt*TMath::CosH(TMath::Sqrt(x[0]*x[0]*Ec*Ec+deltasq)/twokt)*TMath::CosH(TMath::Sqrt(x[0]*x[0]*Ec*Ec+deltasq)/twokt)); + double deltasq(TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]*fPar[3]),2.0)); + f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0); return; } @@ -78,7 +75,7 @@ double TAnSWaveGapIntegralCuhre::IntegrateFunc() const double EPSABS (1e-6); const unsigned int VERBOSE (0); const unsigned int LAST (4); - const unsigned int MINEVAL (0); + const unsigned int MINEVAL (10000); const unsigned int MAXEVAL (1000000); const unsigned int KEY (13); @@ -95,13 +92,10 @@ double TAnSWaveGapIntegralCuhre::IntegrateFunc() } void TAnSWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[], - const int *ncomp, double f[]) // x = {E, phi}, fPar = {T, Delta(T),a} + const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic} { - double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K - double Ec(4.0*(fPar[0]+(1.0+fPar[2])*fPar[1])); // upper limit of energy-integration: cutoff energy - double phic(TMath::PiOver2()); // upper limit of phi-integration - double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*phic)),2.0)); - f[0] = -phic*Ec/(2.0*twokt*TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*Ec*Ec+deltasq)/twokt),2.0)); + double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0)); + f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0); return; } @@ -113,7 +107,7 @@ double TAnSWaveGapIntegralDivonne::IntegrateFunc() const double EPSREL (1e-4); const double EPSABS (1e-6); const unsigned int VERBOSE (0); - const unsigned int MINEVAL (0); + const unsigned int MINEVAL (10000); const unsigned int MAXEVAL (1000000); const unsigned int KEY1 (47); const unsigned int KEY2 (1); @@ -139,12 +133,43 @@ double TAnSWaveGapIntegralDivonne::IntegrateFunc() } void TAnSWaveGapIntegralDivonne::Integrand(const int *ndim, const double x[], - const int *ncomp, double f[]) // x = {E, phi}, fPar = {T, Delta(T),a} + const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic} { - double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K - double Ec(4.0*(fPar[0]+(1.0+fPar[2])*fPar[1])); // upper limit of energy-integration: cutoff energy - double phic(TMath::PiOver2()); // upper limit of phi-integration - double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*phic)),2.0)); - f[0] = -phic*Ec/(2.0*twokt*TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*Ec*Ec+deltasq)/twokt),2.0)); + double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0)); + f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0); + return; +} + +std::vector TAnSWaveGapIntegralSuave::fPar; + +double TAnSWaveGapIntegralSuave::IntegrateFunc() +{ + const unsigned int NCOMP(1); + const double EPSREL (1e-4); + const double EPSABS (1e-6); + const unsigned int VERBOSE (0); + const unsigned int LAST (4); + const unsigned int MINEVAL (10000); + const unsigned int MAXEVAL (1000000); + + const unsigned int NNEW (1000); + const double FLATNESS (25.); + + int nregions, neval, fail; + double integral[NCOMP], error[NCOMP], prob[NCOMP]; + + Suave(fNDim, NCOMP, Integrand, + EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, + NNEW, FLATNESS, + &nregions, &neval, &fail, integral, error, prob); + + return integral[0]; +} + +void TAnSWaveGapIntegralSuave::Integrand(const int *ndim, const double x[], + const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic} +{ + double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0)); + f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0); return; } diff --git a/src/external/libGapIntegrals/TIntegrator.h b/src/external/libGapIntegrals/TIntegrator.h index 0d72ea07..1fadf779 100644 --- a/src/external/libGapIntegrals/TIntegrator.h +++ b/src/external/libGapIntegrals/TIntegrator.h @@ -170,6 +170,20 @@ class TAnSWaveGapIntegralDivonne { }; +class TAnSWaveGapIntegralSuave { + public: + TAnSWaveGapIntegralSuave() : fNDim(2) {} + ~TAnSWaveGapIntegralSuave() { fPar.clear(); } + void SetParameters(const std::vector &par) { fPar=par; } + static void Integrand(const int*, const double[], const int*, double[]); + double IntegrateFunc(); + + protected: + static vector fPar; + unsigned int fNDim; + +}; + // To be integrated: x*y dx dy class T2DTest : public TMCIntegrator { diff --git a/src/msr2data.cpp b/src/msr2data.cpp index 326064f1..a4f0cddf 100644 --- a/src/msr2data.cpp +++ b/src/msr2data.cpp @@ -178,8 +178,6 @@ int msr2data_doFitting(vector &arg, bool &chainfit) temp = -1; chainfit = false; iter = arg.erase(iter); - if (iter == arg.end()) - break; } else if (!iter->substr(0,4).compare("fit-")) { if (temp) { @@ -194,8 +192,6 @@ int msr2data_doFitting(vector &arg, bool &chainfit) iss.str(s); iss >> temp; iter = arg.erase(iter); - if (iter == arg.end()) - break; } else { iter++; }