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

As always forgot some more files...

Browse Source
This commit is contained in:
Bastian M. Wojek 2009-10-17 19:20:17 +00:00
parent 68092aad6b
commit 0a5d35944d
5 changed files with 15 additions and 557 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
src/external/libGapIntegrals/Makefile.libGapIntegrals vendored
View File

@ -116,33 +116,37 @@ endif
# some definitions: headers (used to generate *Dict* stuff), sources, objects,...
OBJS =
OBJS += TIntegrator.o
OBJS += TGapIntegrals.o TGapIntegralsDict.o
EXTOBJS = BMWIntegrator.o
# make the shared libs:
all: $(SHLIB)
$(SHLIB): $(OBJS)
$(SHLIB): $(EXTOBJS) $(OBJS)
@echo "---> Building shared library $(SHLIB) ..."
/bin/rm -f $(SHLIB)
$(LD) $(SOFLAGS) $(LDFLAGS) -o $(SHLIB) $(OBJS) $(LIBS) $(PMUSRLIB) $(GSLLIB) $(CUBALIB)
$(LD) $(SOFLAGS) $(LDFLAGS) -o $(SHLIB) $(EXTOBJS) $(OBJS) $(LIBS) $(PMUSRLIB) $(GSLLIB) $(CUBALIB)
@echo "done"
# clean up: remove all object file (and core files)
# semicolon needed to tell make there is no source
# for this target!
#
clean:; @rm -f $(OBJS) $(SHLIB) *Dict* core* *~
@echo "---> removing $(OBJS) $(SHLIB)"
clean:; @rm -f $(OBJS) $(EXTOBJS) $(SHLIB) *Dict* core* *~
@echo "---> removing $(OBJS) $(EXTOBJS) $(SHLIB)"
#
BMWIntegrator.o: ../BMWIntegrator/BMWIntegrator.cpp
$(CXX) $(INCLUDES) $(CXXFLAGS) -c $<
$(OBJS): %.o: %.cpp
$(CXX) $(INCLUDES) $(CXXFLAGS) -c $<
# Generate the ROOT CINT dictionary
TGapIntegralsDict.cpp: ./TIntegrator.h ./TGapIntegrals.h ./TGapIntegralsLinkDef.h
TGapIntegralsDict.cpp: ../BMWIntegrator/BMWIntegrator.h ./TGapIntegrals.h ./TGapIntegralsLinkDef.h
@echo "Generating dictionary $@..."
rootcint -f $@ -c -p $^
@ -150,16 +154,16 @@ install: all
@echo "Installing shared lib: $(SHLIB)"
ifeq ($(OS),LINUX)
cp -pv $(SHLIB) $(ROOTSYS)/lib
cp -pv TIntegrator.h TGapIntegrals.h $(ROOTSYS)/include
cp -pv ../BMWIntegrator/BMWIntegrator.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
cp -pv ../BMWIntegrator/BMWIntegrator.h TGapIntegrals.h $(ROOTSYS)/include
endif
ifeq ($(OS),DARWIN)
cp -pv $(SHLIB) $(ROOTSYS)/lib
cp -pv TIntegrator.h TGapIntegrals.h $(ROOTSYS)/include
cp -pv ../BMWIntegrator/BMWIntegrator.h TGapIntegrals.h $(ROOTSYS)/include
endif
cleaninstall: clean install

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

@ -34,7 +34,7 @@
using namespace std;
#include "TIntegrator.h"
#include "../BMWIntegrator/BMWIntegrator.h"
#include "TGapIntegrals.h"
#define PI 3.14159265358979323846

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

@ -38,7 +38,7 @@
using namespace std;
#include "PUserFcnBase.h"
#include "TIntegrator.h"
#include "../BMWIntegrator/BMWIntegrator.h"
class TGapSWave : public PUserFcnBase {

242
src/external/libGapIntegrals/TIntegrator.cpp vendored
View File

@ -1,242 +0,0 @@
/***************************************************************************
TIntegrator.cpp
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2009 by Bastian M. Wojek *
* bastian.wojek@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
# include "TIntegrator.h"
# include "cuba.h"
std::vector<double> TDWaveGapIntegralCuhre::fPar;
double TDWaveGapIntegralCuhre::IntegrateFunc()
{
const unsigned int NCOMP(1);
const double EPSREL (1e-4);
const double EPSABS (1e-6);
const unsigned int VERBOSE (0);
const unsigned int LAST (4);
const unsigned int MINEVAL (0);
const unsigned int MAXEVAL (50000);
const unsigned int KEY (13);
int nregions, neval, fail;
double integral[NCOMP], error[NCOMP], prob[NCOMP];
Cuhre(fNDim, NCOMP, Integrand,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
}
void TDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T), Ec, phic}
{
double deltasq(TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]*fPar[3]),2.0));
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0);
return;
}
std::vector<double> TAnSWaveGapIntegralCuhre::fPar;
double TAnSWaveGapIntegralCuhre::IntegrateFunc()
{
const unsigned int NCOMP(1);
const double EPSREL (1e-4);
const double EPSABS (1e-6);
const unsigned int VERBOSE (0);
const unsigned int LAST (4);
const unsigned int MINEVAL (100);
const unsigned int MAXEVAL (1000000);
const unsigned int KEY (13);
int nregions, neval, fail;
double integral[NCOMP], error[NCOMP], prob[NCOMP];
Cuhre(fNDim, NCOMP, Integrand,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
}
void TAnSWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
{
double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0));
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
return;
}
std::vector<double> TAnSWaveGapIntegralDivonne::fPar;
double TAnSWaveGapIntegralDivonne::IntegrateFunc()
{
const unsigned int NCOMP(1);
const double EPSREL (1e-4);
const double EPSABS (1e-6);
const unsigned int VERBOSE (0);
const unsigned int MINEVAL (1000);
const unsigned int MAXEVAL (1000000);
const unsigned int KEY1 (47);
const unsigned int KEY2 (1);
const unsigned int KEY3 (1);
const unsigned int MAXPASS (5);
const double BORDER (0.);
const double MAXCHISQ (10.);
const double MINDEVIATION (.25);
const unsigned int NGIVEN (0);
const unsigned int LDXGIVEN (fNDim);
const unsigned int NEXTRA (0);
int nregions, neval, fail;
double integral[NCOMP], error[NCOMP], prob[NCOMP];
Divonne(fNDim, NCOMP, Integrand,
EPSREL, EPSABS, VERBOSE, MINEVAL, MAXEVAL,
KEY1, KEY2, KEY3, MAXPASS, BORDER, MAXCHISQ, MINDEVIATION,
NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
}
void TAnSWaveGapIntegralDivonne::Integrand(const int *ndim, const double x[],
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
{
double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0));
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
return;
}
std::vector<double> TAnSWaveGapIntegralSuave::fPar;
double TAnSWaveGapIntegralSuave::IntegrateFunc()
{
const unsigned int NCOMP(1);
const double EPSREL (1e-4);
const double EPSABS (1e-6);
const unsigned int VERBOSE (0);
const unsigned int LAST (4);
const unsigned int MINEVAL (1000);
const unsigned int MAXEVAL (1000000);
const unsigned int NNEW (1000);
const double FLATNESS (25.);
int nregions, neval, fail;
double integral[NCOMP], error[NCOMP], prob[NCOMP];
Suave(fNDim, NCOMP, Integrand,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
NNEW, FLATNESS,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
}
void TAnSWaveGapIntegralSuave::Integrand(const int *ndim, const double x[],
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
{
double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0));
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
return;
}
std::vector<double> TNonMonDWave1GapIntegralCuhre::fPar;
double TNonMonDWave1GapIntegralCuhre::IntegrateFunc()
{
const unsigned int NCOMP(1);
const double EPSREL (1e-4);
const double EPSABS (1e-6);
const unsigned int VERBOSE (0);
const unsigned int LAST (4);
const unsigned int MINEVAL (100);
const unsigned int MAXEVAL (1000000);
const unsigned int KEY (13);
int nregions, neval, fail;
double integral[NCOMP], error[NCOMP], prob[NCOMP];
Cuhre(fNDim, NCOMP, Integrand,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
}
void TNonMonDWave1GapIntegralCuhre::Integrand(const int *ndim, const double x[],
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
{
double deltasq(TMath::Power(fPar[1]*(fPar[2]*TMath::Cos(2.0*x[1]*fPar[4])+(1.0-fPar[2])*TMath::Cos(6.0*x[1]*fPar[4])),2.0));
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
return;
}
std::vector<double> TNonMonDWave2GapIntegralCuhre::fPar;
double TNonMonDWave2GapIntegralCuhre::IntegrateFunc()
{
const unsigned int NCOMP(1);
const double EPSREL (1e-4);
const double EPSABS (1e-6);
const unsigned int VERBOSE (0);
const unsigned int LAST (4);
const unsigned int MINEVAL (100);
const unsigned int MAXEVAL (1000000);
const unsigned int KEY (13);
int nregions, neval, fail;
double integral[NCOMP], error[NCOMP], prob[NCOMP];
Cuhre(fNDim, NCOMP, Integrand,
EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
KEY,
&nregions, &neval, &fail, integral, error, prob);
return integral[0];
}
void TNonMonDWave2GapIntegralCuhre::Integrand(const int *ndim, const double x[],
const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
{
double deltasq(4.0*fPar[2]/27.0*TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]*fPar[4]), 2.0) \
/ TMath::Power(1.0 + fPar[2]*TMath::Cos(2.0*x[1]*fPar[4])*TMath::Cos(2.0*x[1]*fPar[4]), 3.0));
f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
return;
}

304
src/external/libGapIntegrals/TIntegrator.h vendored
View File

@ -1,304 +0,0 @@
/***************************************************************************
TIntegrator.h
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2009 by Bastian M. Wojek *
* bastian.wojek@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef _TIntegrator_H_
#define _TIntegrator_H_
#include "Math/GSLIntegrator.h"
#include "Math/GSLMCIntegrator.h"
#include "TMath.h"
#include <vector>
using namespace std;
// 1D Integrator base class - the function to be integrated have to be implemented in a derived class
class TIntegrator {
public:
TIntegrator();
virtual ~TIntegrator();
void SetParameters(const std::vector<double> &par) const { fPar=par; }
virtual double FuncAtX(double) const = 0;
double IntegrateFunc(double, double);
protected:
mutable vector<double> fPar;
private:
static double FuncAtXgsl(double, void *);
ROOT::Math::GSLIntegrator *fIntegrator;
mutable double (*fFunc)(double, void *);
};
inline TIntegrator::TIntegrator() : fFunc(0) {
ROOT::Math::GSLIntegrator *integrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS31);
fIntegrator = integrator;
integrator = 0;
delete integrator;
}
inline TIntegrator::~TIntegrator(){
delete fIntegrator;
fIntegrator=0;
fFunc=0;
}
inline double TIntegrator::FuncAtXgsl(double x, void *obj)
{
return ((TIntegrator*)obj)->FuncAtX(x);
}
inline double TIntegrator::IntegrateFunc(double x1, double x2)
{
fFunc = &TIntegrator::FuncAtXgsl;
return fIntegrator->Integral(fFunc, (this), x1, x2);
}
// Multi dimensional GSL Monte Carlo Integrations
class TMCIntegrator {
public:
TMCIntegrator();
virtual ~TMCIntegrator();
void SetParameters(const std::vector<double> &par) const { fPar=par; }
virtual double FuncAtX(double *) const = 0;
double IntegrateFunc(unsigned int, double *, double *);
protected:
mutable vector<double> fPar;
private:
static double FuncAtXgsl(double *, unsigned int, void *);
ROOT::Math::GSLMCIntegrator *fMCIntegrator;
mutable double (*fFunc)(double *, unsigned int, void *);
};
inline TMCIntegrator::TMCIntegrator() : fFunc(0) {
ROOT::Math::GSLMCIntegrator *integrator = new ROOT::Math::GSLMCIntegrator(ROOT::Math::MCIntegration::kMISER, 1.E-6, 1.E-4, 500000);
fMCIntegrator = integrator;
integrator = 0;
delete integrator;
}
inline TMCIntegrator::~TMCIntegrator(){
delete fMCIntegrator;
fMCIntegrator=0;
fFunc=0;
}
inline double TMCIntegrator::FuncAtXgsl(double *x, unsigned int dim, void *obj)
{
return ((TMCIntegrator*)obj)->FuncAtX(x);
}
inline double TMCIntegrator::IntegrateFunc(unsigned int dim, double *x1, double *x2)
{
fFunc = &TMCIntegrator::FuncAtXgsl;
return fMCIntegrator->Integral(fFunc, dim, x1, x2, (this));
}
// Multidimensional Integrator class for a d-wave gap integral using the Cuhre algorithm
class TDWaveGapIntegralCuhre {
public:
TDWaveGapIntegralCuhre() : fNDim(2) {}
~TDWaveGapIntegralCuhre() { fPar.clear(); }
void SetParameters(const std::vector<double> &par) { fPar=par; }
static void Integrand(const int*, const double[], const int*, double[]);
double IntegrateFunc();
protected:
static vector<double> fPar;
unsigned int fNDim;
};
class TAnSWaveGapIntegralCuhre {
public:
TAnSWaveGapIntegralCuhre() : fNDim(2) {}
~TAnSWaveGapIntegralCuhre() { fPar.clear(); }
void SetParameters(const std::vector<double> &par) { fPar=par; }
static void Integrand(const int*, const double[], const int*, double[]);
double IntegrateFunc();
protected:
static vector<double> fPar;
unsigned int fNDim;
};
class TAnSWaveGapIntegralDivonne {
public:
TAnSWaveGapIntegralDivonne() : fNDim(2) {}
~TAnSWaveGapIntegralDivonne() { fPar.clear(); }
void SetParameters(const std::vector<double> &par) { fPar=par; }
static void Integrand(const int*, const double[], const int*, double[]);
double IntegrateFunc();
protected:
static vector<double> fPar;
unsigned int fNDim;
};
class TAnSWaveGapIntegralSuave {
public:
TAnSWaveGapIntegralSuave() : fNDim(2) {}
~TAnSWaveGapIntegralSuave() { fPar.clear(); }
void SetParameters(const std::vector<double> &par) { fPar=par; }
static void Integrand(const int*, const double[], const int*, double[]);
double IntegrateFunc();
protected:
static vector<double> fPar;
unsigned int fNDim;
};
class TNonMonDWave1GapIntegralCuhre {
public:
TNonMonDWave1GapIntegralCuhre() : fNDim(2) {}
~TNonMonDWave1GapIntegralCuhre() { fPar.clear(); }
void SetParameters(const std::vector<double> &par) { fPar=par; }
static void Integrand(const int*, const double[], const int*, double[]);
double IntegrateFunc();
protected:
static vector<double> fPar;
unsigned int fNDim;
};
class TNonMonDWave2GapIntegralCuhre {
public:
TNonMonDWave2GapIntegralCuhre() : fNDim(2) {}
~TNonMonDWave2GapIntegralCuhre() { fPar.clear(); }
void SetParameters(const std::vector<double> &par) { fPar=par; }
static void Integrand(const int*, const double[], const int*, double[]);
double IntegrateFunc();
protected:
static vector<double> fPar;
unsigned int fNDim;
};
// To be integrated: x*y dx dy
class T2DTest : public TMCIntegrator {
public:
T2DTest() {}
~T2DTest() {}
double FuncAtX(double *) const;
};
inline double T2DTest::FuncAtX(double *x) const
{
return x[0]*x[1];
}
// To be integrated: d wave gap integral
class TDWaveGapIntegral : public TMCIntegrator {
public:
TDWaveGapIntegral() {}
~TDWaveGapIntegral() {}
double FuncAtX(double *) const;
};
inline double TDWaveGapIntegral::FuncAtX(double *x) const // x = {E, phi}, fPar = {T, Delta(T)}
{
double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K
double deltasq(TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]),2.0));
return -1.0/(2.0*twokt*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt)*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt));
}
// To be integrated: anisotropic s wave gap integral
class TAnSWaveGapIntegral : public TMCIntegrator {
public:
TAnSWaveGapIntegral() {}
~TAnSWaveGapIntegral() {}
double FuncAtX(double *) const;
};
inline double TAnSWaveGapIntegral::FuncAtX(double *x) const // x = {E, phi}, fPar = {T, Delta(T), a}
{
double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K
double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1])),2.0));
return -1.0/(2.0*twokt*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt)*TMath::CosH(TMath::Sqrt(x[0]*x[0]+deltasq)/twokt));
}
// To be integrated: Bessel function times Exponential
class TIntBesselJ0Exp : public TIntegrator {
public:
TIntBesselJ0Exp() {}
~TIntBesselJ0Exp() {}
double FuncAtX(double) const;
};
inline double TIntBesselJ0Exp::FuncAtX(double x) const
{
return TMath::BesselJ0(TMath::TwoPi()*fPar[0]*x) * TMath::Exp(-fPar[1]*x);
}
// To be integrated: Sine times Gaussian
class TIntSinGss : public TIntegrator {
public:
TIntSinGss() {}
~TIntSinGss() {}
double FuncAtX(double) const;
};
inline double TIntSinGss::FuncAtX(double x) const
{
return TMath::Sin(TMath::TwoPi()*fPar[0]*x) * TMath::Exp(-0.5*fPar[1]*fPar[1]*x*x);
}
// To be integrated: df/dE * E / sqrt(E^2 - Delta^2)
class TGapIntegral : public TIntegrator {
public:
TGapIntegral() {}
~TGapIntegral() {}
double FuncAtX(double) const; // variable: E
};
inline double TGapIntegral::FuncAtX(double e) const
{
return 1.0/(TMath::Power(TMath::CosH(TMath::Sqrt(e*e+fPar[1]*fPar[1])/fPar[0]),2.0));
}
#endif //_TIntegrator_H_