diff --git a/src/external/libGapIntegrals/GapIntegrals.pdf b/src/external/libGapIntegrals/GapIntegrals.pdf new file mode 100644 index 00000000..86112612 Binary files /dev/null and b/src/external/libGapIntegrals/GapIntegrals.pdf differ diff --git a/src/external/libGapIntegrals/INSTALL b/src/external/libGapIntegrals/INSTALL index ad927704..e8163714 100644 --- a/src/external/libGapIntegrals/INSTALL +++ b/src/external/libGapIntegrals/INSTALL @@ -1,6 +1,7 @@ #--------------------------------------------------------------------- # INSTALL -# Bastian M. Wojek, 2009/09/07 +# Bastian M. Wojek +# $Id$ #--------------------------------------------------------------------- Installation of the musrfit-plugin "libGapIntegrals.so" @@ -22,9 +23,15 @@ Installation of the musrfit-plugin "libGapIntegrals.so" * The plugin classes can now be accessed inside THEORY blocks of musrfit msr-files via: - userFcn libGapIntegrals.so TGapSWave 4 5 - userFcn libGapIntegrals.so TGapDWave 4 5 - userFcn libGapIntegrals.so TGapAnSWave 4 5 6 + userFcn libGapIntegrals.so TGapSWave 1 2 (Tc Delta0) + userFcn libGapIntegrals.so TGapDWave 1 2 (Tc Delta0) + userFcn libGapIntegrals.so TGapAnSWave 1 2 3 (Tc Delta0 a) + userFcn libGapIntegrals.so TGapNonMonDWave1 1 2 3 (Tc Delta0 a) + userFcn libGapIntegrals.so TGapNonMonDWave2 1 2 3 (Tc Delta0 a) + userFcn libGapIntegrals.so TGapPowerLaw 1 2 (Tc n) + +* For an introductory discussion on which functions are actually calculated, please refer to + GapIntegrals.pdf #--------------------------------------------------------------------- # this is the end ... diff --git a/src/external/libGapIntegrals/Makefile.libGapIntegrals b/src/external/libGapIntegrals/Makefile.libGapIntegrals index cebbf628..ff85507b 100644 --- a/src/external/libGapIntegrals/Makefile.libGapIntegrals +++ b/src/external/libGapIntegrals/Makefile.libGapIntegrals @@ -1,41 +1,126 @@ #--------------------------------------------------- -# get compilation flags from root-config +# Makefile.libGapIntegrals +# +# Author: Bastian M. Wojek +# e-mail: bastian.wojek@psi.ch +# +# $Id$ +# +#--------------------------------------------------- + +#--------------------------------------------------- +# get compilation and library flags from root-config ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags) -# ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs) +ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs) +ROOTGLIBS = $(shell $(ROOTSYS)/bin/root-config --glibs) #--------------------------------------------------- +# depending on the architecture, choose the compiler, +# linker, and the flags to use +# -OS = LINUX +ARCH = $(shell $(ROOTSYS)/bin/root-config --arch) + +ifeq ($(ARCH),linux) +OS = LINUX +endif +ifeq ($(ARCH),linuxx8664gcc) +OS = LINUX +endif +ifeq ($(ARCH),win32gcc) +OS = WIN32GCC +endif +ifeq ($(ARCH),macosx) +OS = DARWIN +endif + +# -- Linux +ifeq ($(OS),LINUX) CXX = g++ CXXFLAGS = -O3 -Wall -Wno-trigraphs -fPIC -MUSRFITINCLUDE = ../../include -LOCALINCLUDE = . -ROOTINCLUDE = $(ROOTSYS)/include -INCLUDES = -I$(LOCALINCLUDE) -I$(MUSRFITINCLUDE) -I$(ROOTINCLUDE) +PMUSRPATH = ../../include +MNPATH = $(ROOTSYS)/include +GSLPATH = /usr/include/gsl +CUBAPATH = /usr/local/include +INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(GSLPATH) -I$(CUBAPATH) LD = g++ LDFLAGS = -O -SOFLAGS = -shared +SOFLAGS = -shared +SHLIB = libGapIntegrals.so +endif + +# -- Windows/Cygwin +ifeq ($(OS),WIN32GCC) +CXX = g++ +CXXFLAGS = -O3 -Wall -Wno-trigraphs -D_DLL +PMUSRPATH = ../../include +MNPATH = $(ROOTSYS)/include +GSLPATH = /usr/include/gsl +CUBAPATH = /usr/local/include +INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(GSLPATH) -I$(CUBAPATH) +LD = g++ +LDFLAGS = -O -Wl,--enable-auto-import -Wl,--enable-runtime-pseudo-reloc +SOFLAGS = -shared -Wl,--export-all-symbols +SHLIB = libGapIntegrals.dll +endif + +# -- MacOSX/Darwin +ifeq ($(OS),DARWIN) +CXX = g++ +CXXFLAGS = -O3 -Wall -Wno-trigraphs -fPIC +PMUSRPATH = ../../include +MNPATH = $(ROOTSYS)/include +FINKPATH = /sw/include +GSLPATH = $(FINKPATH)/gsl +CUBAPATH = $(FINKPATH) +INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(GSLPATH) -I$(CUBAPATH) +LD = g++ +LDFLAGS = -O -Xlinker -bind_at_load +SOFLAGS = -dynamiclib -flat_namespace -undefined suppress -Wl,-x +SHLIB = libGapIntegrals.dylib +endif # the output from the root-config script: CXXFLAGS += $(ROOTCFLAGS) -LDFLAGS += -L/usr/local/lib -lcuba -lm +LDFLAGS += + +# the ROOT libraries (G = graphic) +LIBS = $(ROOTLIBS) +GLIBS = $(ROOTGLIBS) + +# GSL lib +GSLLIB = -lgslcblas -lgsl +# Cuba lib +CUBALIB = -L/usr/local/lib -lcuba -lm + +ifeq ($(OS),WIN32GCC) +# GSL lib +GSLLIB = -L/usr/lib -lgslcblas -lgsl +# Cuba lib +CUBALIB = -L/usr/local/lib -lcuba -lm +endif + +ifeq ($(OS),DARWIN) +# GSL lib +GSLLIB = -L/sw/lib -lgslcblas -lgsl +# Cuba lib +CUBALIB = -L/sw/lib -lcuba -lm +endif # some definitions: headers (used to generate *Dict* stuff), sources, objects,... OBJS = OBJS += TIntegrator.o OBJS += TGapIntegrals.o TGapIntegralsDict.o -SHLIB = libGapIntegrals.so +# make the shared libs: -# make the shared lib: -# all: $(SHLIB) $(SHLIB): $(OBJS) @echo "---> Building shared library $(SHLIB) ..." /bin/rm -f $(SHLIB) - $(LD) $(SOFLAGS) -o $(SHLIB) $(OBJS) $(LDFLAGS) + $(LD) $(SOFLAGS) -o $(SHLIB) $(OBJS) $(LDFLAGS) $(LIBS) $(GSLLIB) $(CUBALIB) @echo "done" # clean up: remove all object file (and core files) @@ -53,7 +138,7 @@ $(OBJS): %.o: %.cpp TGapIntegralsDict.cpp: ./TIntegrator.h ./TGapIntegrals.h ./TGapIntegralsLinkDef.h @echo "Generating dictionary $@..." - rootcint -f $@ -c -p -I$(MUSRFITINCLUDE) $^ + rootcint -f $@ -c -p -I$(PMUSRPATH) $^ install: all @echo "Installing shared lib: $(SHLIB)" @@ -61,3 +146,14 @@ ifeq ($(OS),LINUX) cp -pv $(SHLIB) $(ROOTSYS)/lib cp -pv TIntegrator.h TGapIntegrals.h $(ROOTSYS)/include endif +ifeq ($(OS),WIN32GCC) + cp -pv $(SHLIB) $(ROOTSYS)/bin + ln -sf $(ROOTSYS)/bin/$(SHLIB) $(ROOTSYS)/lib/$(SHLIB) + cp -pv TIntegrator.h TGapIntegrals.h $(ROOTSYS)/include +endif +ifeq ($(OS),DARWIN) + cp -pv $(SHLIB) $(ROOTSYS)/lib + cp -pv TIntegrator.h TGapIntegrals.h $(ROOTSYS)/include +endif + +cleaninstall: clean install diff --git a/src/external/libGapIntegrals/PMusrCanvas.cpp.BMWpatch b/src/external/libGapIntegrals/PMusrCanvas.cpp.BMWpatch new file mode 100644 index 00000000..4ffdbad3 --- /dev/null +++ b/src/external/libGapIntegrals/PMusrCanvas.cpp.BMWpatch @@ -0,0 +1,101 @@ +--- PMusrCanvas.cpp.orig 2009-09-18 18:58:50.000000000 +0200 ++++ PMusrCanvas.cpp 2009-09-16 21:51:45.000000000 +0200 +@@ -5,7 +5,7 @@ + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +- $Id$ ++ $Id: PMusrCanvas.cpp 4106 2009-08-31 11:19:03Z nemu $ + + ***************************************************************************/ + +@@ -3291,6 +3291,36 @@ + return; + } + ++// BMW: Writing an extra file with the complete theory curve only until the issue of writing only single theory points has been resolved ++ ++ // generate output filename ++ ++ // in order to handle names with "." correctly this slightly odd data-filename generation ++ tokens = fMsrHandler->GetFileName().Tokenize("."); ++ str = ""; ++ for (int i=0; iGetEntries()-1; i++) { ++ ostr = dynamic_cast(tokens->At(i)); ++ str += ostr->GetString() + TString("."); ++ } ++ str += "dataTheory"; ++ ++ if (tokens) { ++ delete tokens; ++ tokens = 0; ++ } ++ ++ // open file ++ ofstream foutTheory; ++ ++ // open data-file ++ foutTheory.open(str.Data(), iostream::out); ++ if (!fout.is_open()) { ++ cout << endl << ">> PMusrCanvas::SaveDataAscii: **ERROR** couldn't open file " << str.Data() << " for writing." << endl; ++ return; ++ } ++ ++// endofBMW ++ + // extract data + Double_t time, xval, yval; + Int_t xminBin; +@@ -3688,11 +3718,39 @@ + fout << "Data" << fNonMusrData.size()-1 << ", eData" << fNonMusrData.size()-1 << ", Theo" << fNonMusrData.size()-1; + fout << endl; + // write data ++ ++// BMW: Writing an extra file with the complete theory curve only until the issue of writing only single theory points has been resolved ++ ++ // get current x-range ++ xminBin = fNonMusrData[0].theory->GetXaxis()->GetFirst(); // first bin of the zoomed range ++ xmaxBin = fNonMusrData[0].theory->GetXaxis()->GetLast(); // last bin of the zoomed range ++ xmin = fNonMusrData[0].theory->GetXaxis()->GetBinCenter(xminBin); ++ xmax = fNonMusrData[0].theory->GetXaxis()->GetBinCenter(xmaxBin); ++ ++ // get data ++ for (int i=0; iGetN(); i++) { ++ fNonMusrData[0].theory->GetPoint(i,xval,yval); // get values ++ if ((xval < xmin) || (xval > xmax)) ++ continue; ++ foutTheory << xval << " "; ++ for (unsigned int j=0; jGetPoint(i,xval,yval); // get values ++ foutTheory << yval << " "; ++ } ++ // write last data set ++ fNonMusrData[fNonMusrData.size()-1].theory->GetPoint(i,xval,yval); // get values ++ foutTheory << yval << " "; ++ foutTheory << endl; ++ } ++ ++// endofBMW ++ + // get current x-range + xminBin = fNonMusrData[0].data->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fNonMusrData[0].data->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fNonMusrData[0].data->GetXaxis()->GetBinCenter(xminBin); + xmax = fNonMusrData[0].data->GetXaxis()->GetBinCenter(xmaxBin); ++ + // get data + for (int i=0; iGetN(); i++) { + fNonMusrData[0].data->GetPoint(i,xval,yval); // get values +@@ -3739,6 +3797,12 @@ + // close file + fout.close(); + ++// BMW: Writing an extra file with the complete theory curve only until the issue of writing only single theory points has been resolved ++ ++ foutTheory.close(); ++ ++// endofBMW ++ + cout << endl << ">> Data windows saved in ascii format ..." << endl; + } + diff --git a/src/external/libGapIntegrals/TGapIntegrals.cpp b/src/external/libGapIntegrals/TGapIntegrals.cpp index 3555cec2..3016607a 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.cpp +++ b/src/external/libGapIntegrals/TGapIntegrals.cpp @@ -5,7 +5,7 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2009/09/06 + $Id$ ***************************************************************************/ @@ -44,6 +44,9 @@ using namespace std; ClassImp(TGapSWave) ClassImp(TGapDWave) ClassImp(TGapAnSWave) +ClassImp(TGapNonMonDWave1) +ClassImp(TGapNonMonDWave2) +ClassImp(TGapPowerLaw) // s wave gap integral @@ -86,6 +89,32 @@ TGapAnSWave::TGapAnSWave() { fPar.clear(); } +TGapNonMonDWave1::TGapNonMonDWave1() { + TNonMonDWave1GapIntegralCuhre *gapint = new TNonMonDWave1GapIntegralCuhre(); + fGapIntegral = gapint; + gapint = 0; + delete gapint; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + +TGapNonMonDWave2::TGapNonMonDWave2() { + TNonMonDWave2GapIntegralCuhre *gapint = new TNonMonDWave2GapIntegralCuhre(); + fGapIntegral = gapint; + gapint = 0; + delete gapint; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + TGapSWave::~TGapSWave() { delete fGapIntegral; fGapIntegral = 0; @@ -119,6 +148,28 @@ TGapAnSWave::~TGapAnSWave() { fPar.clear(); } +TGapNonMonDWave1::~TGapNonMonDWave1() { + delete fGapIntegral; + fGapIntegral = 0; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + +TGapNonMonDWave2::~TGapNonMonDWave2() { + delete fGapIntegral; + fGapIntegral = 0; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + double TGapSWave::operator()(double t, const vector &par) const { assert(par.size() == 2); // two parameters: Tc, Delta(0) @@ -165,11 +216,11 @@ double TGapSWave::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))); fGapIntegral->SetParameters(intPar); - ds = 1.0+2.0*fGapIntegral->IntegrateFunc(0.0, 2.0*(t+intPar[1])); + ds = 1.0-1.0/intPar[0]*fGapIntegral->IntegrateFunc(0.0, 2.0*(t+intPar[1])); intPar.clear(); @@ -327,3 +378,154 @@ double TGapAnSWave::operator()(double t, const vector &par) const { return fIntegralValues[vectorIndex]; } + +double TGapNonMonDWave1::operator()(double t, const vector &par) const { + + assert(par.size() == 3); // two parameters: Tc, Delta(0), a + + 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: 2 k_B T, Delta(T), a, E_c, phi_c + 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+intPar[1])); // upper limit of energy-integration: cutoff energy + intPar.push_back(TMath::PiOver2()); // upper limit of phi-integration + + fGapIntegral->SetParameters(intPar); + + ds = 1.0-intPar[3]/intPar[0]*fGapIntegral->IntegrateFunc(); + + intPar.clear(); + + if (newTemp) + fIntegralValues.push_back(ds); + else + fIntegralValues[vectorIndex] = ds; + + fCalcNeeded[vectorIndex] = false; + } + + return fIntegralValues[vectorIndex]; + +} + +double TGapNonMonDWave2::operator()(double t, const vector &par) const { + + assert(par.size() == 3); // two parameters: Tc, Delta(0), a + + 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: 2 k_B T, Delta(T), a, E_c, phi_c + 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+intPar[1])); // upper limit of energy-integration: cutoff energy + intPar.push_back(TMath::PiOver2()); // upper limit of phi-integration + + fGapIntegral->SetParameters(intPar); + + ds = 1.0-intPar[3]/intPar[0]*fGapIntegral->IntegrateFunc(); + + intPar.clear(); + + if (newTemp) + fIntegralValues.push_back(ds); + else + fIntegralValues[vectorIndex] = ds; + + fCalcNeeded[vectorIndex] = false; + } + + return fIntegralValues[vectorIndex]; + +} + +double TGapPowerLaw::operator()(double t, const vector &par) const { + + assert(par.size() == 2); // two parameters: Tc, N + + if(t<=0.0) + return 1.0; + else if (t >= par[0]) + return 0.0; + + return 1.0 - pow(t/par[0], par[1]); + +} diff --git a/src/external/libGapIntegrals/TGapIntegrals.h b/src/external/libGapIntegrals/TGapIntegrals.h index 4efa847f..3d2ec2a5 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.h +++ b/src/external/libGapIntegrals/TGapIntegrals.h @@ -5,7 +5,7 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2009/09/06 + $Id$ ***************************************************************************/ @@ -100,4 +100,59 @@ private: ClassDef(TGapAnSWave,1) }; +class TGapNonMonDWave1 : public PUserFcnBase { + +public: + TGapNonMonDWave1(); + ~TGapNonMonDWave1(); + + double operator()(double, const vector&) const; + +private: + TNonMonDWave1GapIntegralCuhre *fGapIntegral; + mutable vector fTemp; + mutable vector::const_iterator fTempIter; + mutable vector fIntegralValues; + mutable vector fCalcNeeded; + + mutable vector fPar; + + ClassDef(TGapNonMonDWave1,1) +}; + +class TGapNonMonDWave2 : public PUserFcnBase { + +public: + TGapNonMonDWave2(); + ~TGapNonMonDWave2(); + + double operator()(double, const vector&) const; + +private: + TNonMonDWave2GapIntegralCuhre *fGapIntegral; + mutable vector fTemp; + mutable vector::const_iterator fTempIter; + mutable vector fIntegralValues; + mutable vector fCalcNeeded; + + mutable vector fPar; + + ClassDef(TGapNonMonDWave2,1) +}; + + +class TGapPowerLaw : public PUserFcnBase { + +public: + TGapPowerLaw() {} + ~TGapPowerLaw() {} + + double operator()(double, const vector&) const; + +private: + + ClassDef(TGapPowerLaw,1) +}; + + #endif //_TGapIntegrals_H_ diff --git a/src/external/libGapIntegrals/TGapIntegralsLinkDef.h b/src/external/libGapIntegrals/TGapIntegralsLinkDef.h index c112003b..ddb102a4 100644 --- a/src/external/libGapIntegrals/TGapIntegralsLinkDef.h +++ b/src/external/libGapIntegrals/TGapIntegralsLinkDef.h @@ -5,7 +5,7 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2009/09/06 + $Id$ ***************************************************************************/ @@ -40,6 +40,9 @@ #pragma link C++ class TGapSWave+; #pragma link C++ class TGapDWave+; #pragma link C++ class TGapAnSWave+; +#pragma link C++ class TGapNonMonDWave1+; +#pragma link C++ class TGapNonMonDWave2+; +#pragma link C++ class TGapPowerLaw+; #endif //__CINT__ // root dictionary stuff -------------------------------------------------- diff --git a/src/external/libGapIntegrals/TIntegrator.cpp b/src/external/libGapIntegrals/TIntegrator.cpp index c17700ec..129bcaa2 100644 --- a/src/external/libGapIntegrals/TIntegrator.cpp +++ b/src/external/libGapIntegrals/TIntegrator.cpp @@ -5,7 +5,7 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2009/09/07 + $Id$ ***************************************************************************/ @@ -75,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 (1000); + const unsigned int MINEVAL (100); const unsigned int MAXEVAL (1000000); const unsigned int KEY (13); @@ -173,3 +173,70 @@ void TAnSWaveGapIntegralSuave::Integrand(const int *ndim, const double x[], 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 TNonMonDWave1GapIntegralCuhre::fPar; + +double TNonMonDWave1GapIntegralCuhre::IntegrateFunc() +{ + const unsigned int NCOMP(1); + const double EPSREL (1e-4); + const double EPSABS (1e-6); + const unsigned int VERBOSE (0); + const unsigned int LAST (4); + const unsigned int MINEVAL (100); + const unsigned int MAXEVAL (1000000); + + const unsigned int KEY (13); + + int nregions, neval, fail; + double integral[NCOMP], error[NCOMP], prob[NCOMP]; + + Cuhre(fNDim, NCOMP, Integrand, + EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, + KEY, + &nregions, &neval, &fail, integral, error, prob); + + return integral[0]; +} + +void TNonMonDWave1GapIntegralCuhre::Integrand(const int *ndim, const double x[], + const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic} +{ + double deltasq(TMath::Power(fPar[1]*(fPar[2]*TMath::Cos(2.0*x[1]*fPar[4])+(1.0-fPar[2])*TMath::Cos(6.0*x[1]*fPar[4])),2.0)); + f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0); + return; +} + +std::vector TNonMonDWave2GapIntegralCuhre::fPar; + +double TNonMonDWave2GapIntegralCuhre::IntegrateFunc() +{ + const unsigned int NCOMP(1); + const double EPSREL (1e-4); + const double EPSABS (1e-6); + const unsigned int VERBOSE (0); + const unsigned int LAST (4); + const unsigned int MINEVAL (100); + const unsigned int MAXEVAL (1000000); + + const unsigned int KEY (13); + + int nregions, neval, fail; + double integral[NCOMP], error[NCOMP], prob[NCOMP]; + + Cuhre(fNDim, NCOMP, Integrand, + EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, + KEY, + &nregions, &neval, &fail, integral, error, prob); + + return integral[0]; +} + +void TNonMonDWave2GapIntegralCuhre::Integrand(const int *ndim, const double x[], + const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic} +{ + double deltasq(4.0*fPar[2]/27.0*TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]*fPar[4]), 2.0) \ + / TMath::Power(1.0 + fPar[2]*TMath::Cos(2.0*x[1]*fPar[4])*TMath::Cos(2.0*x[1]*fPar[4]), 3.0)); + f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0); + return; +} diff --git a/src/external/libGapIntegrals/TIntegrator.h b/src/external/libGapIntegrals/TIntegrator.h index 1fadf779..0a80899b 100644 --- a/src/external/libGapIntegrals/TIntegrator.h +++ b/src/external/libGapIntegrals/TIntegrator.h @@ -5,7 +5,7 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2009/09/07 + $Id$ ***************************************************************************/ @@ -184,6 +184,34 @@ class TAnSWaveGapIntegralSuave { }; +class TNonMonDWave1GapIntegralCuhre { + public: + TNonMonDWave1GapIntegralCuhre() : fNDim(2) {} + ~TNonMonDWave1GapIntegralCuhre() { 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; + +}; + +class TNonMonDWave2GapIntegralCuhre { + public: + TNonMonDWave2GapIntegralCuhre() : fNDim(2) {} + ~TNonMonDWave2GapIntegralCuhre() { 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 { @@ -230,7 +258,6 @@ 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)); } - // To be integrated: Bessel function times Exponential class TIntBesselJ0Exp : public TIntegrator { @@ -259,7 +286,7 @@ inline double TIntSinGss::FuncAtX(double x) const return TMath::Sin(TMath::TwoPi()*fPar[0]*x) * TMath::Exp(-0.5*fPar[1]*fPar[1]*x*x); } -// To be integrated: df/dE * E / sqrt(E^2+ Delta^2) +// To be integrated: df/dE * E / sqrt(E^2 - Delta^2) class TGapIntegral : public TIntegrator { public: @@ -271,9 +298,7 @@ class TGapIntegral : public TIntegrator { inline double TGapIntegral::FuncAtX(double e) const { - double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K -// return -e/(2.0*twokt*TMath::CosH(e/twokt)*TMath::CosH(e/twokt)*TMath::Sqrt(e*e-fPar[1]*fPar[1])); - return -1.0/(2.0*twokt*TMath::CosH(TMath::Sqrt(e*e+fPar[1]*fPar[1])/twokt)*TMath::CosH(TMath::Sqrt(e*e+fPar[1]*fPar[1])/twokt)); + return 1.0/(TMath::Power(TMath::CosH(TMath::Sqrt(e*e+fPar[1]*fPar[1])/fPar[0]),2.0)); } #endif //_TIntegrator_H_