diff --git a/src/external/libGbGLF/Makefile.am b/src/external/libGbGLF/Makefile.am index 2c77bbd1..cd04bf0f 100644 --- a/src/external/libGbGLF/Makefile.am +++ b/src/external/libGbGLF/Makefile.am @@ -18,7 +18,7 @@ dict_cpp_sources = \ include_HEADERS = $(h_sources) noinst_HEADERS = $(h_linkdef) $(dict_h_sources) -AM_CPPFLAGS = -I$(top_srcdir)/src/include $(PMUSR_CFLAGS) $(FFTW3_CFLAGS) $(GSL_CFLAGS) -I$(ROOTINCDIR) +AM_CPPFLAGS = -I$(top_srcdir)/src/include $(PMUSR_CFLAGS) -I$(ROOTINCDIR) AM_CXXFLAGS = $(LOCAL_LIB_CXXFLAGS) BUILT_SOURCES = $(dict_cpp_sources) $(dict_h_sources) @@ -31,7 +31,7 @@ CLEANFILES = *Dict.cpp *Dict.h *~ core lib_LTLIBRARIES = libGbGLF.la libGbGLF_la_SOURCES = $(h_sources) $(cpp_sources) $(dict_h_sources) $(dict_cpp_sources) -libGbGLF_la_LIBADD = $(USERFCN_LIBS) $(GSL_LIBS) +libGbGLF_la_LIBADD = $(USERFCN_LIBS) libGbGLF_la_LDFLAGS = -version-info $(PLUGIN_LIBRARY_VERSION) -release $(PLUGIN_RELEASE) $(AM_LDFLAGS) ## For the moment do not build pkgconfig files for musrfit plug-ins... diff --git a/src/external/libGbGLF/PGbGLF.cpp b/src/external/libGbGLF/PGbGLF.cpp index edfafba5..a3e07fde 100644 --- a/src/external/libGbGLF/PGbGLF.cpp +++ b/src/external/libGbGLF/PGbGLF.cpp @@ -43,45 +43,6 @@ using namespace std; ClassImp(PGbGLF) -//-------------------------------------------------------------------------- -// Constructor -//-------------------------------------------------------------------------- -/** - * - */ -PGbGLF::PGbGLF() -{ - fWs = 0; - fWs = gsl_integration_workspace_alloc(65536); - if (fWs == 0) { - cerr << "debug> PGbGLF::PGbGLF(): **ERROR** couldn't allocate fWs ..." << endl; - } - - fIntTab = 0; - - fGslFunParam.wExt = 0.0; - fGslFunParam.s0 = 0.0; - fGslFunParam.s1 = 0.0; - - fPrevParam[0] = 0.0; - fPrevParam[1] = 0.0; - fPrevParam[2] = 0.0; -} - -//-------------------------------------------------------------------------- -// Destructor -//-------------------------------------------------------------------------- -/** - * - */ -PGbGLF::~PGbGLF() -{ - if (fWs) - gsl_integration_workspace_free(fWs); - if (fIntTab) - gsl_integration_qawo_table_free(fIntTab); -} - //-------------------------------------------------------------------------- // operator() //-------------------------------------------------------------------------- @@ -93,59 +54,50 @@ PGbGLF::~PGbGLF() */ Double_t PGbGLF::operator()(Double_t t, const std::vector &par) const { - // param: [0] Bext (G), [1] Delta0 (1/us), [2] DeltaGbG (1/us) + // param: [0] Bext (G), [1] Delta0 (1/us), [2] Rb (1) assert(par.size()==3); if (t<0.0) return 1.0; - if (t>fTmax) - t = fTmax; Double_t dval = 0.0; Double_t wExt = TWO_PI * GAMMA_BAR_MUON * par[0]; Double_t s0 = par[1]; - Double_t s1 = par[2]; + Double_t s1 = s0*par[2]; // sigma0 * Rb Double_t wExt2 = wExt*wExt; Double_t s02 = s0*s0; Double_t s12 = s1*s1; Double_t t2 = t*t; - Double_t result = 0.0; - Double_t err = 0.0; + if (par[0] < 3.0) { // take ZF solution for Bext <= 3G. This is good enough + Double_t aa = 1.0 + s12*t2; + dval = 0.333333333333333333 + 0.666666666666667 / pow(aa, 1.5)*(1.0-s02*t2/aa)*exp(-0.5*s02*t2/aa); + } else { + // P_z^LF (GbG, 1st part) + Double_t aa = 1.0+t2*s12; + dval = 1.0 - 2.0*(s02+s12)/wExt2 + 2.0*(s02+s12*aa)/(wExt2*pow(aa,2.5))*cos(wExt*t)*exp(-0.5*s02*t2/aa); - // if integration table is not already setup, do it - if (fIntTab == 0) { - fIntTab = gsl_integration_qawo_table_alloc(wExt, fTmax, GSL_INTEG_SINE, 16); - if (fIntTab == 0) { - cerr << ">> PGbGLF::operator(): **ERROR** couldn't invoke qawo table." << endl; - return 0.0; - } - fFun.function = &pz_GbG_2; - fGslFunParam.wExt = wExt; - fGslFunParam.s0 = s0; - fGslFunParam.s1 = s1; - fFun.params = &fGslFunParam; + // P_z^LF (GbG, 2nd part) + Double_t dt = t; + Int_t n=1; + Double_t sumT = dt * pz_GbG_2(t, par) * 0.5; + Double_t sumM = 0.0; + Double_t tt = 0.0; + do { + sumM = 0.0; + for (Int_t i=0; i 1.0e-5) && (n < 8192)); + + dval += sumT; } - // if parameter set has changed, adopted fFun parameters - if (fPrevParam[1] != par[1]) { - fPrevParam[0] = par[0]; - fPrevParam[1] = par[1]; - fPrevParam[2] = par[2]; - fGslFunParam.wExt = wExt; - fGslFunParam.s0 = s0; - fGslFunParam.s1 = s1; - fFun.params = &fGslFunParam; - } - - // P_z^LF (GbG, 1st part) - Double_t aa = 1.0+t2*s12; - dval = 1.0 - 2.0*(s02+s12)/wExt2 + 2.0*(s02+s12*aa)/(wExt2*pow(aa,2.5))*cos(wExt*t)*exp(-0.5*s02*t2/aa); - - // P_z^LF (GbG, 2nd part) - gsl_integration_qawo_table_set_length(fIntTab, t); - gsl_integration_qawo(&fFun, 0.0, 1.0e-9, 1.0e-9, 1024, fWs, fIntTab, &result, &err); - dval += result; return dval; } @@ -154,20 +106,20 @@ Double_t PGbGLF::operator()(Double_t t, const std::vector &par) const // pz_GbG_2 (private) //-------------------------------------------------------------------------- /** - *

Integrand of the non-analytic part of the GbG-LF polarization + *

Integrand of the non-analytic part * - * \param x - * \param param + * \param t time */ - -double pz_GbG_2(double x, void *param) +Double_t PGbGLF::pz_GbG_2(Double_t t, const std::vector &par) const { - gslFunParam *p = (gslFunParam*) param; + Double_t wExt = TWO_PI * GAMMA_BAR_MUON * par[0]; + Double_t s0 = par[1]; + Double_t s1 = s0*par[2]; // sigma0 * Rb - double s02 = p->s0 * p->s0; - double s12 = p->s1 * p->s1; - double x2 = x*x; - double aa = 1.0+x2*s12; + Double_t s02 = s0*s0; + Double_t s12 = s1*s1; + Double_t t2 = t*t; + Double_t aa = 1.0+t2*s12; - return 2.0*(s02*s02+3.0*s12*s12*aa*aa+6.0*s02*s12*aa)/(pow(p->wExt, 3.0)*pow(aa, 4.5))*exp(-0.5*s02*x2/aa); + return 2.0*(s02*s02+3.0*s12*s12*aa*aa+6.0*s02*s12*aa)/(pow(wExt,3.0)*pow(aa,4.5))*exp(-0.5*s02*t2/aa)*sin(wExt*t); } diff --git a/src/external/libGbGLF/PGbGLF.h b/src/external/libGbGLF/PGbGLF.h index 5a07f748..aa01f369 100644 --- a/src/external/libGbGLF/PGbGLF.h +++ b/src/external/libGbGLF/PGbGLF.h @@ -33,19 +33,8 @@ #include using namespace std; -#include - #include "PUserFcnBase.h" -typedef struct { - double wExt; - double s0; - double s1; -} gslFunParam; - - -double pz_GbG_2(double x, void *param); - //-------------------------------------------------------------------------------------------- /** *

Interface class for the user function. @@ -53,8 +42,8 @@ double pz_GbG_2(double x, void *param); class PGbGLF : public PUserFcnBase { public: - PGbGLF(); - virtual ~PGbGLF(); + PGbGLF() {} + virtual ~PGbGLF() {} virtual Bool_t NeedGlobalPart() const { return false; } virtual void SetGlobalPart(vector &globalPart, UInt_t idx) {} @@ -63,14 +52,9 @@ class PGbGLF : public PUserFcnBase virtual Double_t operator()(Double_t t, const std::vector ¶m) const; private: - static const Double_t fTmax = 12.0; - mutable Double_t fPrevParam[3]; - mutable gslFunParam fGslFunParam; - gsl_integration_workspace *fWs; - mutable gsl_function fFun; - mutable gsl_integration_qawo_table *fIntTab; + Double_t pz_GbG_2(Double_t t, const std::vector &par) const; - ClassDef(PGbGLF, 1) + ClassDef(PGbGLF, 1) }; #endif // _PGBGLF_H_