Added Solt's ZF depolarization function to the libZFRelaxation

While doing that I realized that all LF functions within libLFRelaxation cannot be used at the moment
with a parallelized musrfit installation. Therefore, the built-in functions should be used for now.
One day this should be fixed ...

src/external/BMWtools/BMWIntegrator.h

@@ -40,6 +40,91 @@
 
 using namespace std;
 
+/**
+ * <p>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.
+ *    Therefore, integrals with different parameters can be calculated in parallel (at least I hope so).
+ *    The function which should be integrated has to be implemented in a derived class.
+ *    Note: The purpose of this is to offer an easy-to-use interface---not the most efficient integration routine.
+ */
+class T2Integrator {
+  public:
+    T2Integrator();
+    virtual ~T2Integrator();
+    virtual double FuncAtX(double, const std::vector<double> &par) const = 0;
+    virtual double IntegrateFunc(double, double, const std::vector<double> &par);
+
+  private:
+    static double FuncAtXgsl(double, void *);
+    ROOT::Math::GSLIntegrator *fIntegrator; ///< pointer to the GSL integrator
+};
+
+/**
+ * <p>Constructor of the alternative base class for 1D integrations
+ *    Allocation of memory for an integration using the adaptive 31 point Gauss-Kronrod rule
+ */
+inline T2Integrator::T2Integrator() {
+  fIntegrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS31);
+}
+
+/**
+ * <p>Destructor of the alternative base class for 1D integrations
+ *    Clean up.
+ */
+inline T2Integrator::~T2Integrator(){
+  delete fIntegrator;
+  fIntegrator=0;
+}
+
+/**
+ * <p>Method for passing the integrand function value to the integrator.
+ *
+ * <p><b>return:</b>
+ * - function value of the integrand
+ *
+ * \param x point at which the function value is calculated
+ * \param ptrPair pointers to the function object and the parameters for the integration
+ */
+inline double T2Integrator::FuncAtXgsl(double x, void *ptrPair)
+{
+  pair<T2Integrator*, const vector<double>*> *pairOfPointers = static_cast<pair<T2Integrator*, const vector<double>*>*>(ptrPair);
+  return pairOfPointers->first->FuncAtX(x, *(pairOfPointers->second));
+}
+
+/**
+ * <p>Calculate the integral of the function between the given boundaries
+ *
+ * <p><b>return:</b>
+ * - value of the integral
+ *
+ * \param x1 lower boundary
+ * \param x2 upper boundary
+ * \param par additional parameters for the integration
+ */
+inline double T2Integrator::IntegrateFunc(double x1, double x2, const std::vector<double> &par)
+{
+  pair<T2Integrator*, const vector<double>*> ptrPair;
+  ptrPair.first = (this);
+  ptrPair.second = &par;
+  return fIntegrator->Integral(&T2Integrator::FuncAtXgsl, static_cast<void*>(&ptrPair), x1, x2);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 /**
  * <p>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.
@@ -515,4 +600,60 @@ 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));
 }
 
+/**
+ * <p>Class for the 1D integration for the calculation of the uniaxial static Gauss-Kubo-Toyabe function
+ *    The integration uses the GSL integration routines.
+ */
+class TFirstUniaxialGssKTIntegral : public T2Integrator {
+  public:
+    TFirstUniaxialGssKTIntegral() {}
+    ~TFirstUniaxialGssKTIntegral() {}
+    double FuncAtX(double, const vector<double>&) const; // variable: x
+};
+
+/**
+ * <p>Calculate the function value---actual implementation of the integrand in Eq. (7) of Solt's article
+ *
+ * <p><b>return:</b>
+ * - function value
+ *
+ * \param x point where the function should be evaluated
+ */
+inline double TFirstUniaxialGssKTIntegral::FuncAtX(double x, const vector<double> &par) const
+{
+  double eps(par[0]*par[0]/(par[1]*par[1]) - 1.0);
+  double p(1.0 + eps*x*x);
+  double SsqTsq(par[0]*par[0]*par[2]*par[2]);
+
+  return (1.0 - x*x)*(p - SsqTsq)/TMath::Power(p, 2.5)*TMath::Exp(-0.5*SsqTsq/p);
+}
+
+/**
+ * <p>Class for the 1D integration for the calculation of the uniaxial static Gauss-Kubo-Toyabe function
+ *    The integration uses the GSL integration routines.
+ */
+class TSecondUniaxialGssKTIntegral : public T2Integrator {
+  public:
+    TSecondUniaxialGssKTIntegral() {}
+    ~TSecondUniaxialGssKTIntegral() {}
+    double FuncAtX(double, const vector<double>&) const; // variable: x
+};
+
+/**
+ * <p>Calculate the function value---actual implementation of the integrand in Eq. (7) of Solt's article
+ *
+ * <p><b>return:</b>
+ * - function value
+ *
+ * \param x point where the function should be evaluated
+ */
+inline double TSecondUniaxialGssKTIntegral::FuncAtX(double x, const vector<double> &par) const
+{
+  double eps(par[0]*par[0]/(par[1]*par[1]) - 1.0);
+  double p(1.0 + eps*x*x);
+  double SsqTsq(par[0]*par[0]*par[2]*par[2]);
+
+  return (p - SsqTsq)/TMath::Power(p, 2.5)*TMath::Exp(-0.5*SsqTsq/p);
+}
+
 #endif //_TIntegrator_H_

src/external/libLFRelaxation/TLFRelaxation.h

@@ -38,14 +38,6 @@
 
 using namespace std;
 
-//#include "gsl/gsl_math.h"
-//#include "gsl/gsl_sf_exp.h"
-//#include "gsl/gsl_sf_log.h"
-//#include "gsl/gsl_sf_trig.h"
-//#include "gsl/gsl_sf_bessel.h"
-
-
-//#include "TMath.h"
 #include "PUserFcnBase.h"
 #include "fftw3.h"
 #include "BMWIntegrator.h"

src/external/libZFRelaxation/Makefile.am

@@ -19,7 +19,7 @@ dict_cpp_sources = \
 include_HEADERS = $(h_sources)
 noinst_HEADERS = $(h_linkdef) $(dict_h_sources)
 
-INCLUDES = -I$(top_srcdir)/src/include $(PMUSR_CFLAGS) $(ROOT_CFLAGS)
+INCLUDES = -I$(top_srcdir)/src/include $(BMWTOOLS_CFLAGS) $(PMUSR_CFLAGS) $(ROOT_CFLAGS)
 AM_CXXFLAGS = $(LOCAL_LIB_CXXFLAGS)
 
 BUILT_SOURCES = $(dict_cpp_sources) $(dict_h_sources)
@@ -32,7 +32,7 @@ CLEANFILES = *Dict.cpp *Dict.h *~ core
 lib_LTLIBRARIES = libZFRelaxation.la
 
 libZFRelaxation_la_SOURCES = $(h_sources) $(cpp_sources) $(dict_h_sources) $(dict_cpp_sources)
-libZFRelaxation_la_LIBADD = $(USERFCN_LIBS) $(ROOT_LIBS)
+libZFRelaxation_la_LIBADD = $(BMWTOOLS_LIBS) $(USERFCN_LIBS) $(GSL_LIBS) $(ROOT_LIBS)
 libZFRelaxation_la_LDFLAGS = -version-info $(PLUGIN_LIBRARY_VERSION) -release $(PLUGIN_RELEASE) $(AM_LDFLAGS)
 
 ## For the moment do not build pkgconfig files for musrfit plug-ins...

src/external/libZFRelaxation/ZFRelaxation.cpp

@@ -71,3 +71,86 @@ double ZFMagExp::operator()(double t, const vector<double> &par) const {
   double w(TWO_PI * par[1]);
   return par[0] * (cos(w*t) - par[2]/w*sin(w*t)) * exp(-par[2]*t) + (1.0-par[0]) * exp(-par[3]*t);
 }
+
+//--------------------------------------------------------------------------
+// UniaxialStatGssKT::UniaxialStatGssKT()
+//--------------------------------------------------------------------------
+/**
+ * <p>Constructor
+ */
+UniaxialStatGssKT::UniaxialStatGssKT() {
+  fIntFirst = new TFirstUniaxialGssKTIntegral();
+  fIntSecond = new TSecondUniaxialGssKTIntegral();
+}
+
+//--------------------------------------------------------------------------
+// UniaxialStatGssKT::~UniaxialStatGssKT()
+//--------------------------------------------------------------------------
+/**
+ * <p>Destructor
+ */
+UniaxialStatGssKT::~UniaxialStatGssKT() {
+  delete fIntFirst; fIntFirst = 0;
+  delete fIntSecond; fIntSecond = 0;
+}
+
+//--------------------------------------------------------------------------
+// UniaxialStatGssKT::operator()
+//--------------------------------------------------------------------------
+/**
+ * <p>Method returning the function value at a given time for a given set of parameters.
+ *
+ * \param t time \htmlonly (&mu;s) \endhtmlonly \latexonly ($\mu\mathrm{s}$) \endlatexonly
+ * \param par parameters [\htmlonly &sigma;<sub>1</sub> (&mu;s<sup>-1</sup>), &sigma;<sub>2</sub> (&mu;s<sup>-1</sup>), &Theta; (&deg;)\endhtmlonly \latexonly $\sigma_{1}~(\mu\mathrm{s}^{-1})$, $\sigma_{2}~(\mu\mathrm{s}^{-1})$, $\Theta~(^{\circ})$ \endlatexonly]
+ */
+double UniaxialStatGssKT::operator()(double t, const vector<double> &par) const {
+  assert(par.size() == 3);
+
+  if (t < 0.0)
+    return 1.0;
+
+  // set the parameters for the integrations
+  vector<double> intParam(3);
+  intParam[0] = par[0];
+  intParam[1] = par[1];
+  intParam[2] = t;
+
+  if (((fabs(par[1]) < 1.0e-5) && (fabs(par[0]) > 1.0e-2)) || (fabs(par[0]/par[1]) > 1000.0)) {
+    cerr << endl;
+    cerr << ">> UniaxialStatGssKT: WARNING Ratio sigma1/sigma2 unreasonably large! Set it internally to 1000. Please check your parameters!";
+    cerr << endl;
+    intParam[1] = 1.0e-3*intParam[0];
+  } else if (fabs(par[1]) < 1.0e-10) {
+    cerr << endl;
+    cerr << ">> UniaxialStatGssKT: WARNING sigma2 too small! Set it internally to 1.0E-10. Please check your parameters!";
+    cerr << endl;
+    intParam[1] = 1.0e-10;
+  }
+
+  double eps(intParam[0]*intParam[0]/(intParam[1]*intParam[1]) - 1.0);
+
+  // check if the anisotropy is so small that the normal statGssKT can be used
+  if (fabs(eps) < 1.0e-6) {
+    double sSqtSq(intParam[0]*intParam[0]*t*t);
+    return 0.333333333333333 + 0.666666666666667*(1.0 - sSqtSq)*exp(-0.5*sSqtSq);
+  }
+
+  double sinSqTh(pow(sin(DEG_TO_RAD*par[2]), 2.0)), cosSqTh(1.5*pow(cos(DEG_TO_RAD*par[2]), 2.0) - 0.5);
+
+  // check if the limit sigma2 >> sigma1 applies
+  if (eps < -0.999999) {
+    return 0.5*sinSqTh + cosSqTh;
+  }
+
+  // otherwise calculate the full depolarization function
+  double sqrtOnePlusEps(sqrt(1.0 + eps));
+  double avg;
+  if (eps > 0.0) {
+    avg = 1.0 - sqrtOnePlusEps/eps*(sqrtOnePlusEps - log(sqrt(eps) + sqrtOnePlusEps)/sqrt(eps));
+  } else {
+    avg = 1.0 - sqrtOnePlusEps/eps*(sqrtOnePlusEps - asin(sqrt(-eps))/sqrt(-eps));
+  }
+
+  return sinSqTh*(0.5 + sqrtOnePlusEps*fIntSecond->IntegrateFunc(0.0, 1.0, intParam)) + cosSqTh*(avg + sqrtOnePlusEps*fIntFirst->IntegrateFunc(0.0, 1.0, intParam));
+
+}

src/external/libZFRelaxation/ZFRelaxation.h

@@ -32,13 +32,14 @@
 #define _ZFRelaxation_H_
 
 #include "PUserFcnBase.h"
+#include "BMWIntegrator.h"
 
 #include <vector>
 using namespace std;
 
 //-----------------------------------------------------------------------------------------------------------------
 /**
- * <p>User function for the muon spin depolarization resulting from static Gaussian broadened randomly oriented internal fields
+ * <p>User function for the muon-spin depolarization resulting from static Gaussian broadened randomly oriented internal fields
  * See also: E.I. Kornilov and V.Yu. Pomjakushin
  * \htmlonly Phys. Lett. A <b>153</b>, 364&#150;367 (1991), doi:<a href="http://dx.doi.org/10.1016/0375-9601(91)90959-C">10.1016/0375-9601(91)90959-C</a>
  * \endhtmlonly
@@ -63,7 +64,7 @@ public:
 
 //-----------------------------------------------------------------------------------------------------------------
 /**
- * <p>User function for the muon spin depolarization resulting from static Lorentzian broadened randomly oriented internal fields
+ * <p>User function for the muon-spin depolarization resulting from static Lorentzian broadened randomly oriented internal fields
  * See also: M. I. Larkin, Y. Fudamoto, I. M. Gat, A. Kinkhabwala, K. M. Kojima, G. M. Luke, J. Merrin, B. Nachumi, Y. J. Uemura, M. Azuma, T. Saito, and M. Takano
  * \htmlonly Physica B <b>289&#150;290</b>, 153&#150;156 (2000), doi:<a href="http://dx.doi.org/10.1016/S0921-4526(00)00337-9">10.1016/S0921-4526(00)00337-9</a>
  * \endhtmlonly
@@ -86,4 +87,32 @@ public:
   ClassDef(ZFMagExp,1)
 };
 
+//-----------------------------------------------------------------------------------------------------------------
+/**
+ * <p>User function for the muon-spin depolarization resulting from static Gaussian distributed fields with uniaxial anisotropy
+ * See also: G. Solt
+ * \htmlonly Hyperfine Interactions <b>96</b>, 167&#150;175 (1995), doi:<a href="http://dx.doi.org/10.1007/BF02066280">10.1007/BF02066280</a>
+ * \endhtmlonly
+ * \latexonly Hyperfine Interactions \textbf{96}, 167--175 (1995), \texttt{http://dx.doi.org/10.1007/BF02066280}
+ * \endlatexonly
+ */
+class UniaxialStatGssKT : public PUserFcnBase {
+
+public:
+  UniaxialStatGssKT();
+  virtual ~UniaxialStatGssKT();
+
+  virtual Bool_t NeedGlobalPart() const { return false; }
+  virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
+  virtual Bool_t GlobalPartIsValid() const { return true; }
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  TFirstUniaxialGssKTIntegral *fIntFirst;
+  TSecondUniaxialGssKTIntegral *fIntSecond;
+
+  ClassDef(UniaxialStatGssKT,1)
+};
+
 #endif //_ZFRelaxation_H_

src/external/libZFRelaxation/ZFRelaxationLinkDef.h

@@ -37,6 +37,7 @@
 
 #pragma link C++ class ZFMagGss+;
 #pragma link C++ class ZFMagExp+;
+#pragma link C++ class UniaxialStatGssKT+;
 
 #endif //__CINT__
 // root dictionary stuff --------------------------------------------------