- Removed the extensions of the external libraries in the test files.

ROOT is anyway checking multiple extensions for dynamic libraries, therefore leaving them out yields platform-independent msr files.

- Minor changes in libTFitPofB

- Added a user function for Uemura's ZF/LF dynamical spin-glass relaxation function
  see, e.g. Y.J. Uemura et al., Phys. Rev. B 31, 546-563 (1985)
  or Y.J. Uemura, Hyperfine Interact. 8, 739 (1981)

  This yields an overall Lorentzian field distribution with motional narrowing.
  However, the present implementation is at best some design study:
  The Lorentzian distribution of width "a" is modeled by Gaussian distributions of widths "sigma" from 0 to infinity.
  Some "resonable cut-offs" would be: 0.1*a < sigma < 10000*a
  Due to finite memory and computing power, in the present version these cut-offs are: 0.1*a < sigma < 40*a
  This yields depolarization functions overall similar to those in Uemura's articles.
  However, due to the low cut-off the first derivative of the depolarization function is zero in the limit t->0
  (what should not be the case for a true Lorentzan).
  Furthermore, the calculation is rather slow and the resulting functions should only be regarded as crude approximations.
  Therefore, at the moment this is far from being of practical use in analyzing experimental data.

configure.ac

@@ -8,7 +8,9 @@ dnl -----------------------------------------------
 dnl Check if pkg-config is installed
 dnl -----------------------------------------------
 
-m4_ifdef([PKG_CHECK_MODULES],[],AC_MSG_ERROR([Please install pkg-config before configuring musrfit!]))
+PKG_PROG_PKG_CONFIG([0.10])
+#The above macro still needs to be tested under Cygwin without pkg-config---but it is promising to be a better solution than the macro below.
+#m4_ifdef([PKG_CHECK_MODULES],[],AC_MSG_ERROR([Please install pkg-config before configuring musrfit!]))
 
 dnl -----------------------------------------------
 dnl Package names and version numbers
@@ -501,6 +503,9 @@ if test "${BUILD_BMW_LIBS}" = "1"; then
   LIBS="$LIBS $FFTW3_LIBS"
   AC_SEARCH_LIBS([fftw_init_threads], [fftw3_threads], [FFTW3_LIBS="$FFTW3_LIBS -lfftw3_threads -lpthread"
                                       AC_DEFINE([HAVE_LIBFFTW3_THREADS], [1], [Define to 1 if fftw3_threads are available])], [], [-lpthread])
+  # Check for fftw3f library. If it is not available the BMWlibs will not be built!
+  AC_SEARCH_LIBS([fftwf_malloc], [fftw3f], [FFTW3_LIBS="$FFTW3_LIBS -lfftw3f"], [BUILD_BMW_LIBS=0
+                                      AC_MSG_WARN([The float version of FFTW3 is not available. The BMWlibs will not be built!])], [])
   CFLAGS="$SAVED_CFLAGS"
   LIBS="$SAVED_LIBS"
 

src/external/BMWIntegrator/BMWIntegrator.cpp

@@ -5,7 +5,7 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  $Id
+  $Id$
 
 ***************************************************************************/
 

src/external/BMWIntegrator/BMWIntegrator.h

@@ -5,7 +5,7 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  $Id
+  $Id$
 
 ***************************************************************************/
 

src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp

@@ -1994,6 +1994,7 @@ void TBulkAnisotropicTriVortexLondonFieldCalc::CalculateGrid() const {
   const int NFFT(fSteps);
   const int NFFT_2(fSteps/2);
   const int NFFTsq(fSteps*fSteps);
+  const int NFFTsq_2((fSteps/2 + 1) * fSteps);
 
   // fill the field Fourier components in the matrix
 
@@ -2015,70 +2016,102 @@ void TBulkAnisotropicTriVortexLondonFieldCalc::CalculateGrid() const {
   double kk, ll;
   int k, l, lNFFT_2;
 
-  for (l = 0; l < NFFT_2; l += 2) {
+  // zero first everything since the r2c FFT changes the input, too
-    lNFFT_2 = l*(NFFT_2 + 1);
+  #ifdef HAVE_GOMP
-    ll = static_cast<double>(l*l);
+  #pragma omp parallel for default(shared) private(k) schedule(dynamic)
-    for (k = 0; k < NFFT_2; k += 2) {
+  #endif
-      kk = static_cast<double>(k*k);
+  for (k = 0; k < NFFTsq_2; ++k) {
-      fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+    fFFTin[k][0] = 0.0;
-      fFFTin[lNFFT_2 + k][1] = 0.0;
+    fFFTin[k][1] = 0.0;
-      fFFTin[lNFFT_2 + k + 1][0] = 0.0;
-      fFFTin[lNFFT_2 + k + 1][1] = 0.0;
-    }
-    k = NFFT_2;
-    kk = static_cast<double>(k*k);
-    fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
-    fFFTin[lNFFT_2 + k][1] = 0.0;
   }
 
+  #ifdef HAVE_GOMP
+  #pragma omp parallel private(k, l, lNFFT_2, kk, ll) num_threads(4)
+  {
+    switch(omp_get_thread_num()) {
+      case 0:
+  #endif
+        for (l = 0; l < NFFT_2; l += 2) {
+          lNFFT_2 = l*(NFFT_2 + 1);
+          ll = static_cast<double>(l*l);
+          for (k = 0; k < NFFT_2; k += 2) {
+            kk = static_cast<double>(k*k);
+            fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+//            fFFTin[lNFFT_2 + k][1] = 0.0;
+//            fFFTin[lNFFT_2 + k + 1][0] = 0.0;
+//            fFFTin[lNFFT_2 + k + 1][1] = 0.0;
+          }
 
-  for (l = NFFT_2; l < NFFT; l += 2) {
+          k = NFFT_2;
-    lNFFT_2 = l*(NFFT_2 + 1);
+          kk = static_cast<double>(k*k);
-    ll = static_cast<double>((NFFT-l)*(NFFT-l));
+          fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
-    for (k = 0; k < NFFT_2; k += 2) {
+//          fFFTin[lNFFT_2 + k][1] = 0.0;
-      kk = static_cast<double>(k*k);
+
-      fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+        }
-      fFFTin[lNFFT_2 + k][1] = 0.0;
+  #ifdef HAVE_GOMP
-      fFFTin[lNFFT_2 + k + 1][0] = 0.0;
+      break;
-      fFFTin[lNFFT_2 + k + 1][1] = 0.0;
+    case 1:
+  #endif
+        for (l = NFFT_2; l < NFFT; l += 2) {
+          lNFFT_2 = l*(NFFT_2 + 1);
+          ll = static_cast<double>((NFFT-l)*(NFFT-l));
+          for (k = 0; k < NFFT_2; k += 2) {
+            kk = static_cast<double>(k*k);
+            fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+//            fFFTin[lNFFT_2 + k][1] = 0.0;
+//            fFFTin[lNFFT_2 + k + 1][0] = 0.0;
+//            fFFTin[lNFFT_2 + k + 1][1] = 0.0;
+          }
+
+          k = NFFT_2;
+          kk = static_cast<double>(k*k);
+          fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+//          fFFTin[lNFFT_2 + k][1] = 0.0;
+        }
+  #ifdef HAVE_GOMP
+      break;
+    case 2:
+  #endif
+        // intermediate rows
+        for (l = 1; l < NFFT_2; l += 2) {
+          lNFFT_2 = l*(NFFT_2 + 1);
+          ll = static_cast<double>(l*l);
+          for (k = 0; k < NFFT_2; k += 2) {
+            kk = static_cast<double>((k + 1)*(k + 1));
+//            fFFTin[lNFFT_2 + k][0] = 0.0;
+//            fFFTin[lNFFT_2 + k][1] = 0.0;
+            fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+//            fFFTin[lNFFT_2 + k + 1][1] = 0.0;
+          }
+//          k = NFFT_2;
+//          fFFTin[lNFFT_2 + k][0] = 0.0;
+//          fFFTin[lNFFT_2 + k][1] = 0.0;
+        }
+  #ifdef HAVE_GOMP
+      break;
+    case 3:
+  #endif
+        for (l = NFFT_2 + 1; l < NFFT; l += 2) {
+          lNFFT_2 = l*(NFFT_2 + 1);
+          ll = static_cast<double>((NFFT-l)*(NFFT-l));
+          for (k = 0; k < NFFT_2; k += 2) {
+            kk = static_cast<double>((k+1)*(k+1));
+//            fFFTin[lNFFT_2 + k][0] = 0.0;
+//            fFFTin[lNFFT_2 + k][1] = 0.0;
+            fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
+//           fFFTin[lNFFT_2 + k + 1][1] = 0.0;
+          }
+          k = NFFT_2;
+//          fFFTin[lNFFT_2 + k][0] = 0.0;
+//          fFFTin[lNFFT_2 + k][1] = 0.0;
+        }
+  #ifdef HAVE_GOMP
+      break;
+    default:
+      break;
     }
-    k = NFFT_2;
-    kk = static_cast<double>(k*k);
-    fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
-    fFFTin[lNFFT_2 + k][1] = 0.0;
-  }
-
-  // intermediate rows
-
-  for (l = 1; l < NFFT_2; l += 2) {
-    lNFFT_2 = l*(NFFT_2 + 1);
-    ll = static_cast<double>(l*l);
-    for (k = 0; k < NFFT_2; k += 2) {
-      kk = static_cast<double>((k + 1)*(k + 1));
-      fFFTin[lNFFT_2 + k][0] = 0.0;
-      fFFTin[lNFFT_2 + k][1] = 0.0;
-      fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
-      fFFTin[lNFFT_2 + k + 1][1] = 0.0;
-    }
-    k = NFFT_2;
-    fFFTin[lNFFT_2 + k][0] = 0.0;
-    fFFTin[lNFFT_2 + k][1] = 0.0;
-  }
-
-  for (l = NFFT_2 + 1; l < NFFT; l += 2) {
-    lNFFT_2 = l*(NFFT_2 + 1);
-    ll = static_cast<double>((NFFT-l)*(NFFT-l));
-    for (k = 0; k < NFFT_2; k += 2) {
-      kk = static_cast<double>((k+1)*(k+1));
-      fFFTin[lNFFT_2 + k][0] = 0.0;
-      fFFTin[lNFFT_2 + k][1] = 0.0;
-      fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk);
-      fFFTin[lNFFT_2 + k + 1][1] = 0.0;
-    }
-    k = NFFT_2;
-    fFFTin[lNFFT_2 + k][0] = 0.0;
-    fFFTin[lNFFT_2 + k][1] = 0.0;
   }
+  #endif
 
   // Do the Fourier transform to get B(x,y)
 

src/external/TFitPofB-lib/classes/TPofBCalc.cpp

@@ -62,7 +62,10 @@ TPofBCalc::TPofBCalc(const vector<double> &para) : fBmin(0.0), fBmax(0.0), fDT(p
 
   int i;
 
-  for (i = 0; i < static_cast<int>(fPBSize); i++) {
+  #ifdef HAVE_GOMP
+  #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+  #endif
+  for (i = 0; i < static_cast<int>(fPBSize); ++i) {
     fB[i] = static_cast<double>(i)*fDB;
     fPB[i] = 0.0;
   }
@@ -78,7 +81,12 @@ TPofBCalc::TPofBCalc(const vector<double>& b, const vector<double>& pb, double d
   fB = new double[fPBSize];
   fPB = new double[fPBSize];
 
-  for (unsigned int i(0); i < fPBSize; i++) {
+  int i;
+
+  #ifdef HAVE_GOMP
+  #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+  #endif
+  for (i = 0; i < static_cast<int>(fPBSize); ++i) {
     fB[i] = b[i];
     fPB[i] = pb[i];
   }
@@ -86,22 +94,22 @@ TPofBCalc::TPofBCalc(const vector<double>& b, const vector<double>& pb, double d
   vector<double>::const_iterator iter, iterB;
   iterB = b.begin();
 
-  for(iter = pb.begin(); iter != pb.end(); iter++){
+  for(iter = pb.begin(); iter != pb.end(); ++iter){
     if(*iter != 0.0) {
       fBmin = *iterB;
 //      cout << fBmin << endl;
       break;
     }
-    iterB++;
+    ++iterB;
   }
 
-  for( ; iter != b.end(); iter++){
+  for( ; iter != b.end(); ++iter){
     if(*iter == 0.0) {
       fBmax = *(iterB-1);
 //      cout << fBmax << endl;
       break;
     }
-    iterB++;
+    ++iterB;
   }
 
   fDT = dt; // needed if a convolution should be done
@@ -117,12 +125,35 @@ void TPofBCalc::UnsetPBExists() {
 #ifdef HAVE_GOMP
 #pragma omp parallel for default(shared) private(i) schedule(dynamic)
 #endif
-  for (i = 0; i < static_cast<int>(fPBSize); i++) {
+  for (i = 0; i < static_cast<int>(fPBSize); ++i) {
     fPB[i] = 0.0;
   }
   fPBExists = false;
 }
 
+void TPofBCalc::Normalize(unsigned int minFilledIndex = 0, unsigned int maxFilledIndex = 0) const {
+
+    if (!maxFilledIndex)
+      maxFilledIndex = fPBSize - 1;
+
+    int i;
+
+    double pBsum(0.0);
+
+    #ifdef HAVE_GOMP
+    #pragma omp parallel for default(shared) private(i) schedule(dynamic) reduction(+:pBsum)
+    #endif
+    for (i = minFilledIndex; i <= static_cast<int>(maxFilledIndex); ++i)
+      pBsum += fPB[i];
+    pBsum *= fDB;
+
+    #ifdef HAVE_GOMP
+    #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+    #endif
+    for (i = minFilledIndex; i <= static_cast<int>(maxFilledIndex); ++i)
+      fPB[i] /= pBsum;
+}
+
 void TPofBCalc::Calculate(const string &type, const vector<double> &para) {
 
   if (type == "skg"){ // skewed Gaussian
@@ -133,29 +164,30 @@ void TPofBCalc::Calculate(const string &type, const vector<double> &para) {
     int a3(static_cast<int>(floor(fBmax/fDB)));
     int a4(static_cast<int>(ceil(fBmax/fDB)));
 
-    unsigned int BmaxIndex((a3 < a4) ? a4 : (a4 + 1));
+    int BmaxIndex((a3 < a4) ? a4 : (a4 + 1));
-    unsigned int B0Index(static_cast<int>(ceil(para[2]/(gBar*fDB))));
+    int B0Index(static_cast<int>(ceil(para[2]/(gBar*fDB))));
 
     double expominus(para[3]*para[3]/(2.0*pi*pi*gBar*gBar));
     double expoplus(para[4]*para[4]/(2.0*pi*pi*gBar*gBar));
     double B0(para[2]/(gBar));
 
-    for (unsigned int i(0); i < B0Index; i++) {
+    int i;
+
+    #ifdef HAVE_GOMP
+    #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+    #endif
+    for (i = 0; i < B0Index; ++i) {
       fPB[i] = exp(-(fB[i]-B0)*(fB[i]-B0)/expominus);
     }
 
-    for (unsigned int i(B0Index); i <= BmaxIndex; i++) {
+    #ifdef HAVE_GOMP
+    #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+    #endif
+    for (i = B0Index; i <= BmaxIndex; ++i) {
       fPB[i] = exp(-(fB[i]-B0)*(fB[i]-B0)/expoplus);
     }
 
-    // normalize p(B)
+    Normalize(0, BmaxIndex);
-
-    double pBsum = 0.0;
-    for (unsigned int i(0); i <= BmaxIndex; i++)
-      pBsum += fPB[i];
-    pBsum *= fDB;
-    for (unsigned int i(0); i <= BmaxIndex; i++)
-      fPB[i] /= pBsum;
 
   }
 
@@ -192,12 +224,12 @@ void TPofBCalc::Calculate(const TBofZCalcInverse *BofZ, const TTrimSPData *dataT
 
   // calculate p(B) from the inverse of B(z)
 
-  for (i = firstZerosEnd; i <= lastZerosStart; i++) {
+  for (i = firstZerosEnd; i <= lastZerosStart; ++i) {
 
     vector< pair<double, double> > inv;
     inv = BofZ->GetInverseAndDerivative(fB[i]);
 
-    for (unsigned int j(0); j < inv.size(); j++) {
+    for (unsigned int j(0); j < inv.size(); ++j) {
       fPB[i] += dataTrimSP->GetNofZ(inv[j].first, para[2])*fabs(inv[j].second);
     }
 //    if (fPB[i])
@@ -206,12 +238,7 @@ void TPofBCalc::Calculate(const TBofZCalcInverse *BofZ, const TTrimSPData *dataT
 
   // normalize p(B)
 
-  double pBsum = 0.0;
+  Normalize(firstZerosEnd, lastZerosStart);
-  for (i = firstZerosEnd; i<=lastZerosStart; i++)
-    pBsum += fPB[i];
-  pBsum *= fDB;
-  for (i = firstZerosEnd; i<=lastZerosStart; i++)
-    fPB[i] /= pBsum;
 
   if(para.size() == 6 && para[5] != 0.0)
     AddBackground(para[3], para[4], para[5]);
@@ -366,13 +393,7 @@ void TPofBCalc::Calculate(const TBofZCalc *BofZ, const TTrimSPData *dataTrimSP,
   bofzBZ = 0;
 
   // normalize p(B)
-
+  Normalize(firstZerosEnd, lastZerosStart);
-  double pBsum = 0.0;
-  for (unsigned int i(firstZerosEnd); i<=lastZerosStart; i++)
-    pBsum += fPB[i];
-  pBsum *= fDB;
-  for (unsigned int i(firstZerosEnd); i<=lastZerosStart; i++)
-    fPB[i] /= pBsum;
 
   fPBExists = true;
   return;
@@ -391,8 +412,8 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
   fBmin = vortexLattice->GetBmin();
   fBmax = vortexLattice->GetBmax();
 
-  int a1(static_cast<int>(floor(fBmin/fDB)));
+//  int a1(static_cast<int>(floor(fBmin/fDB)));
-  int a2(static_cast<int>(ceil(fBmin/fDB)));
+//  int a2(static_cast<int>(ceil(fBmin/fDB)));
   int a3(static_cast<int>(floor(fBmax/fDB)));
   int a4(static_cast<int>(ceil(fBmax/fDB)));
 
@@ -459,9 +480,10 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
                + (numberOfSteps_2 - j)*(numberOfSteps_2 - j))/static_cast<double>(numberOfStepsSq);
           fPB[fill_index] += 1.0/(1.0+sigmaSq*Rsq1) + 1.0/(1.0+sigmaSq*Rsq2) + 1.0/(1.0+sigmaSq*Rsq3) \
                            + 1.0/(1.0+sigmaSq*Rsq4) + 1.0/(1.0+sigmaSq*Rsq5) + 1.0/(1.0+sigmaSq*Rsq6);
-
+/*
-//          of << 1.0/(1.0+sigmaSq*Rsq1) + 1.0/(1.0+sigmaSq*Rsq2) + 1.0/(1.0+sigmaSq*Rsq3) \
+          of << 1.0/(1.0+sigmaSq*Rsq1) + 1.0/(1.0+sigmaSq*Rsq2) + 1.0/(1.0+sigmaSq*Rsq3) \
-//                           + 1.0/(1.0+sigmaSq*Rsq4) + 1.0/(1.0+sigmaSq*Rsq5) + 1.0/(1.0+sigmaSq*Rsq6) << " ";
+                           + 1.0/(1.0+sigmaSq*Rsq4) + 1.0/(1.0+sigmaSq*Rsq5) + 1.0/(1.0+sigmaSq*Rsq6) << " ";
+*/
         }
       }
 //      of << endl;
@@ -488,10 +510,10 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
         field = vortexFields[i + numberOfSteps*j] \
               + para[5]*(exp(Rsq1*one_xiSq) + exp(Rsq2*one_xiSq) + exp(Rsq3*one_xiSq) \
                         +exp(Rsq4*one_xiSq) + exp(Rsq5*one_xiSq) + exp(Rsq6*one_xiSq));
-
+/*
-//          of << para[5]*(exp(Rsq1*one_xiSq) - exp(Rsq2*one_xiSq) + exp(Rsq3*one_xiSq) \
+          of << para[5]*(exp(Rsq1*one_xiSq) - exp(Rsq2*one_xiSq) + exp(Rsq3*one_xiSq) \
-//                        -exp(Rsq4*one_xiSq) + exp(Rsq5*one_xiSq) - exp(Rsq6*one_xiSq)) << " ";
+                        -exp(Rsq4*one_xiSq) + exp(Rsq5*one_xiSq) - exp(Rsq6*one_xiSq)) << " ";
-
+*/
         fill_index = static_cast<unsigned int>(ceil(fabs((field/fDB))));
         if (fill_index < fPBSize) {
           fPB[fill_index] += 1.0;
@@ -503,30 +525,68 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
     }
 //    of.close();
   } else {
-    for (unsigned int j(0); j < numberOfSteps_2; ++j) {
+    int i,j;
-      for (unsigned int i(0); i < numberOfSteps_2; ++i) {
+    #ifdef HAVE_GOMP
+    // cannot use a reduction clause here (like e.g. in Normalize()), since pBvec[] is not a scalar variable
+    // therefore, we need to work on it a bit more
+    int n(omp_get_num_procs()), tid, offset;
+    vector< vector<unsigned int> > pBvec(n, vector<unsigned int>(fPBSize, 0));
+
+    int indexStep(static_cast<int>(floor(static_cast<float>(numberOfSteps_2)/static_cast<float>(n))));
+
+    #pragma omp parallel private(tid, i, j, offset, fill_index) num_threads(n)
+    {
+      tid = omp_get_thread_num();
+      offset = tid*indexStep;
+
+      if (tid == n-1) {
+        for (j = offset; j < static_cast<int>(numberOfSteps_2); ++j) {
+          for (i = 0; i < static_cast<int>(numberOfSteps_2); ++i) {
+            fill_index = static_cast<unsigned int>(ceil(fabs((vortexFields[i + numberOfSteps*j]/fDB))));
+            if (fill_index < fPBSize) {
+              pBvec[tid][fill_index] += 1;
+            }
+          }
+        }
+      } else {
+        for (j = 0; j < indexStep; ++j) {
+          for (i = 0; i < static_cast<int>(numberOfSteps_2); ++i) {
+            fill_index = static_cast<unsigned int>(ceil(fabs((vortexFields[offset + i + numberOfSteps*j]/fDB))));
+            if (fill_index < fPBSize) {
+              pBvec[tid][fill_index] += 1;
+            }
+          }
+        }
+      }
+    }
+
+    for (j = 0; j < n; ++j) {
+      #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+      for (i = 0; i < static_cast<int>(fPBSize); ++i) {
+        fPB[i] += static_cast<double>(pBvec[j][i]);
+      }
+      pBvec[j].clear();
+    }
+    pBvec.clear();
+
+    #else
+
+    for (j = 0; j < static_cast<int>(numberOfSteps_2); ++j) {
+      for (i = 0; i < static_cast<int>(numberOfSteps_2); ++i) {
         fill_index = static_cast<unsigned int>(ceil(fabs((vortexFields[i + numberOfSteps*j]/fDB))));
         if (fill_index < fPBSize) {
           fPB[fill_index] += 1.0;
         }
       }
     }
+
+    #endif
   }
 
   vortexFields = 0;
 
   // normalize P(B)
-  double sum(0.0);
+  Normalize();
-  for (unsigned int i(0); i < fPBSize; ++i)
-    sum += fPB[i];
-  sum *= fDB;
-  int i;
-#ifdef HAVE_GOMP
-#pragma omp parallel for default(shared) private(i) schedule(dynamic)
-#endif
-  for (i = 0; i < static_cast<int>(fPBSize); ++i)
-    fPB[i] /= sum;
-// end pragma omp parallel
 
   if(para.size() == 5)
     AddBackground(para[2], para[3], para[4]);
@@ -542,24 +602,40 @@ void TPofBCalc::AddBackground(double B, double s, double w) {
   if(!s || w<0.0 || w>1.0 || B<0.0)
     return;
 
+  int i;
   double BsSq(s*s/(gBar*gBar*4.0*pi*pi));
 
   // calculate Gaussian background
   double bg[fPBSize];
-  for(unsigned int i(0); i < fPBSize; i++) {
+  #ifdef HAVE_GOMP
+  #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+  #endif
+  for(i = 0; i < static_cast<int>(fPBSize); ++i) {
     bg[i] = exp(-(fB[i]-B)*(fB[i]-B)/(2.0*BsSq));
   }
 
   // normalize background
+
   double bgsum(0.0);
-  for (unsigned int i(0); i < fPBSize; i++)
+  #ifdef HAVE_GOMP
+  #pragma omp parallel for default(shared) private(i) schedule(dynamic) reduction(+:bgsum)
+  #endif
+  for (i = 0; i < static_cast<int>(fPBSize); ++i)
     bgsum += bg[i];
+
   bgsum *= fDB;
-  for (unsigned int i(0); i < fPBSize; i++)
+
+  #ifdef HAVE_GOMP
+  #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+  #endif
+  for (i = 0; i < static_cast<int>(fPBSize); ++i)
     bg[i] /= bgsum;
 
   // add background to P(B)
-  for (unsigned int i(0); i < fPBSize; i++)
+  #ifdef HAVE_GOMP
+  #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+  #endif
+  for (i = 0; i < static_cast<int>(fPBSize); ++i)
     fPB[i] = (1.0 - w)*fPB[i] + w*bg[i];
 
 //   // check if normalization is still valid
@@ -582,7 +658,7 @@ void TPofBCalc::ConvolveGss(double w) {
   fftw_plan FFTplanToFieldDomain;
   fftw_complex *FFTout;
 
-  TBin = 1.0/(gBar*double(NFFT-1)*fDB);
+  TBin = 1.0/(gBar*static_cast<double>(NFFT-1)*fDB);
 
   FFTout = new fftw_complex[NFFT/2 + 1]; //(fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (NFFT/2+1));
 
@@ -597,7 +673,11 @@ void TPofBCalc::ConvolveGss(double w) {
   double GssInTimeDomain;
   double expo(-2.0*PI*PI*gBar*gBar*w*w*TBin*TBin);
 
-  for (unsigned int i(0); i < NFFT/2+1; i++) {
+  int i;
+  #ifdef HAVE_GOMP
+  #pragma omp parallel for default(shared) private(GssInTimeDomain, i) schedule(dynamic)
+  #endif
+  for (i = 0; i < static_cast<int>(NFFT/2+1); ++i) {
     GssInTimeDomain = exp(expo*static_cast<double>(i*i));
     FFTout[i][0] *= GssInTimeDomain;
     FFTout[i][1] *= GssInTimeDomain;
@@ -618,24 +698,57 @@ void TPofBCalc::ConvolveGss(double w) {
 //  fftw_cleanup();
 
   // normalize p(B)
-
+  Normalize();
-  double pBsum = 0.0;
-  for (unsigned int i(0); i < NFFT; i++)
-    pBsum += fPB[i];
-  pBsum *= fDB;
-  for (unsigned int i(0); i < NFFT; i++)
-    fPB[i] /= pBsum;
 
   return;
 }
 
 double TPofBCalc::GetFirstMoment() const {
 
+  int i;
+
   double pBsum(0.0);
-  for (unsigned int i(0); i < fPBSize; i++)
+
+  #ifdef HAVE_GOMP
+  #pragma omp parallel for default(shared) private(i) schedule(dynamic) reduction(+:pBsum)
+  #endif
+  for (i = 0; i < static_cast<int>(fPBSize); ++i)
     pBsum += fB[i]*fPB[i];
+
   pBsum *= fDB;
 
   return pBsum;
 
 }
+
+double TPofBCalc::GetCentralMoment(unsigned int n) const {
+
+  double firstMoment(GetFirstMoment());
+  double diff;
+
+  int i;
+
+  double pBsum(0.0);
+
+  #ifdef HAVE_GOMP
+  #pragma omp parallel for default(shared) private(i, diff) schedule(dynamic) reduction(+:pBsum)
+  #endif
+  for (i = 0; i < static_cast<int>(fPBSize); ++i) {
+    diff = fB[i]-firstMoment;
+    pBsum += pow(diff, static_cast<double>(n))*fPB[i];
+  }
+
+  pBsum *= fDB;
+
+  return pBsum;
+
+}
+
+double TPofBCalc::GetSkewnessAlpha() const {
+
+  double M2(GetCentralMoment(2));
+  double M3(GetCentralMoment(3));
+
+  return M3 > 0.0 ? pow(M3, 1.0/3.0)/pow(M2, 0.5) : -pow(-M3, 1.0/3.0)/pow(M2, 0.5);
+
+}

src/external/TFitPofB-lib/classes/TVortex.cpp

@@ -45,6 +45,10 @@ ClassImp(TBulkTriVortexAGL)
 ClassImp(TBulkTriVortexNGL)
 ClassImp(TBulkAnisotropicTriVortexLondonGlobal)
 ClassImp(TBulkAnisotropicTriVortexLondon)
+ClassImp(TBulkAnisotropicTriVortexMLGlobal)
+ClassImp(TBulkAnisotropicTriVortexML)
+ClassImp(TBulkAnisotropicTriVortexAGLGlobal)
+ClassImp(TBulkAnisotropicTriVortexAGL)
 
 //------------------
 // Destructor of the TBulkTriVortexLondon class -- cleaning up

src/external/TFitPofB-lib/include/TPofBCalc.h

@@ -68,7 +68,10 @@ public:
   void ConvolveGss(double);
   void AddBackground(double, double, double);
   double GetFirstMoment() const;
+  double GetCentralMoment(unsigned int) const;
+  double GetSkewnessAlpha() const;
   void UnsetPBExists();
+  void Normalize(unsigned int, unsigned int) const;
 
 private:
   double *fB; ///< array of discrete points in the field domain

src/external/libGapIntegrals/TGapIntegrals.h

@@ -44,7 +44,7 @@ class TGapSWave : public PUserFcnBase {
 
 public:
   TGapSWave();
-  ~TGapSWave();
+  virtual ~TGapSWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -68,7 +68,7 @@ class TGapDWave : public PUserFcnBase {
 
 public:
   TGapDWave();
-  ~TGapDWave();
+  virtual ~TGapDWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -92,7 +92,7 @@ class TGapCosSqDWave : public PUserFcnBase {
 
 public:
   TGapCosSqDWave();
-  ~TGapCosSqDWave();
+  virtual ~TGapCosSqDWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -116,7 +116,7 @@ class TGapSinSqDWave : public PUserFcnBase {
 
 public:
   TGapSinSqDWave();
-  ~TGapSinSqDWave();
+  virtual ~TGapSinSqDWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -141,7 +141,7 @@ class TGapAnSWave : public PUserFcnBase {
 
 public:
   TGapAnSWave();
-  ~TGapAnSWave();
+  virtual ~TGapAnSWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -165,7 +165,7 @@ class TGapNonMonDWave1 : public PUserFcnBase {
 
 public:
   TGapNonMonDWave1();
-  ~TGapNonMonDWave1();
+  virtual ~TGapNonMonDWave1();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -189,7 +189,7 @@ class TGapNonMonDWave2 : public PUserFcnBase {
 
 public:
   TGapNonMonDWave2();
-  ~TGapNonMonDWave2();
+  virtual ~TGapNonMonDWave2();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -214,7 +214,7 @@ class TGapPowerLaw : public PUserFcnBase {
 
 public:
   TGapPowerLaw() {}
-  ~TGapPowerLaw() {}
+  virtual ~TGapPowerLaw() {}
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -231,7 +231,7 @@ class TGapDirtySWave : public PUserFcnBase {
 
 public:
   TGapDirtySWave() {}
-  ~TGapDirtySWave() {}
+  virtual ~TGapDirtySWave() {}
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -249,7 +249,7 @@ class TLambdaSWave : public PUserFcnBase {
 
 public:
   TLambdaSWave();
-  ~TLambdaSWave();
+  virtual ~TLambdaSWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -267,7 +267,7 @@ class TLambdaDWave : public PUserFcnBase {
 
 public:
   TLambdaDWave();
-  ~TLambdaDWave();
+  virtual ~TLambdaDWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -285,7 +285,7 @@ class TLambdaAnSWave : public PUserFcnBase {
 
 public:
   TLambdaAnSWave();
-  ~TLambdaAnSWave();
+  virtual ~TLambdaAnSWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -303,7 +303,7 @@ class TLambdaNonMonDWave1 : public PUserFcnBase {
 
 public:
   TLambdaNonMonDWave1();
-  ~TLambdaNonMonDWave1();
+  virtual ~TLambdaNonMonDWave1();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -321,7 +321,7 @@ class TLambdaNonMonDWave2 : public PUserFcnBase {
 
 public:
   TLambdaNonMonDWave2();
-  ~TLambdaNonMonDWave2();
+  virtual ~TLambdaNonMonDWave2();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -340,7 +340,7 @@ class TLambdaPowerLaw : public PUserFcnBase {
 
 public:
   TLambdaPowerLaw() {}
-  ~TLambdaPowerLaw() {}
+  virtual ~TLambdaPowerLaw() {}
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -357,7 +357,7 @@ class TLambdaInvSWave : public PUserFcnBase {
 
 public:
   TLambdaInvSWave();
-  ~TLambdaInvSWave();
+  virtual ~TLambdaInvSWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -375,7 +375,7 @@ class TLambdaInvDWave : public PUserFcnBase {
 
 public:
   TLambdaInvDWave();
-  ~TLambdaInvDWave();
+  virtual ~TLambdaInvDWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -393,7 +393,7 @@ class TLambdaInvAnSWave : public PUserFcnBase {
 
 public:
   TLambdaInvAnSWave();
-  ~TLambdaInvAnSWave();
+  virtual ~TLambdaInvAnSWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -411,7 +411,7 @@ class TLambdaInvNonMonDWave1 : public PUserFcnBase {
 
 public:
   TLambdaInvNonMonDWave1();
-  ~TLambdaInvNonMonDWave1();
+  virtual ~TLambdaInvNonMonDWave1();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -429,7 +429,7 @@ class TLambdaInvNonMonDWave2 : public PUserFcnBase {
 
 public:
   TLambdaInvNonMonDWave2();
-  ~TLambdaInvNonMonDWave2();
+  virtual ~TLambdaInvNonMonDWave2();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -448,7 +448,7 @@ class TLambdaInvPowerLaw : public PUserFcnBase {
 
 public:
   TLambdaInvPowerLaw() {}
-  ~TLambdaInvPowerLaw() {}
+  virtual ~TLambdaInvPowerLaw() {}
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
@@ -465,7 +465,7 @@ class TFilmMagnetizationDWave : public PUserFcnBase {
 
 public:
   TFilmMagnetizationDWave();
-  ~TFilmMagnetizationDWave();
+  virtual ~TFilmMagnetizationDWave();
 
   virtual Bool_t NeedGlobalPart() const { return false; }
   virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }

src/external/libLFRelaxation/README

@@ -3,17 +3,18 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/12/05
+  2011/03/12
 
 ***************************************************************************/
 
 Implementation of a userFcn-interface to Gaussian and Lorentzian static and dynamic LF relaxation functions.
-At the moment this is localized to l_wojek@pc5405, because an absolute path had to be set. Of course this can be easily 
+At the moment this is a simple alternative implementation to the functions provided by musrfit itself.
-changed in the code if needed.
+Mostly, this effort should be regarded as a design study which is not really indended for production use.
 
-The functions are then called from within musrfit as:
+The functions are called from within musrfit as:
 
-userFcn  libLFRelaxation.so TLFStatGssKT 1 2       (frequency rate)
+userFcn  libLFRelaxation TLFStatGssKT 1 2       (frequency rate)
-userFcn  libLFRelaxation.so TLFStatLorKT 1 2       (frequency rate)
+userFcn  libLFRelaxation TLFStatLorKT 1 2       (frequency rate)
-userFcn  libLFRelaxation.so TLFDynGssKT 1 2 3      (frequency rate fluct.rate)
+userFcn  libLFRelaxation TLFDynGssKT  1 2 3     (frequency rate fluct.rate)
-userFcn  libLFRelaxation.so TLFDynLorKT 1 2 3      (frequency rate fluct.rate)
+userFcn  libLFRelaxation TLFDynLorKT  1 2 3     (frequency rate fluct.rate)
+userFcn  libLFRelaxation TLFDynSG     1 2 3     (frequency rate fluct.rate)

src/external/libLFRelaxation/TLFRelaxation.cpp

@@ -41,7 +41,7 @@ using namespace std;
 
 #define PI 3.14159265358979323846
 #define TWOPI 6.28318530717958647692
-
+#define Nsteps 100
 
 ClassImp(TLFStatGssKT)
 ClassImp(TLFStatLorKT)
@@ -140,7 +140,7 @@ TLFDynGssKT::TLFDynGssKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("Words
   if (wordsOfWisdomR == NULL) {
     cout << "TLFDynGssKT::TLFDynGssKT: Couldn't open wisdom file ..." << endl;
   } else {
-    wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR);
+    wisdomLoaded = fftwf_import_wisdom_from_file(wordsOfWisdomR);
     fclose(wordsOfWisdomR);
   }
 
@@ -150,11 +150,11 @@ TLFDynGssKT::TLFDynGssKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("Words
   // END of WisdomLoading
 
   // allocating memory for the FFtransform pairs and create the FFT plans
-
+ 
-  fFFTtime = (double *)malloc(sizeof(double) * fNSteps);
+  fFFTtime = (float *)malloc(sizeof(float) * fNSteps);
-  fFFTfreq = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (fNSteps/2+1));
+  fFFTfreq = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * (fNSteps/2+1));
-  fFFTplanFORW = fftw_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_ESTIMATE);
+  fFFTplanFORW = fftwf_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_ESTIMATE);
-  fFFTplanBACK = fftw_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_ESTIMATE);
+  fFFTplanBACK = fftwf_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_ESTIMATE);
 }
 
 TLFDynGssKT::~TLFDynGssKT() {
@@ -164,22 +164,22 @@ TLFDynGssKT::~TLFDynGssKT() {
   if (wordsOfWisdomW == NULL) {
     cout << "TLFDynGssKT::~TLFDynGssKT: Could not open file ... No wisdom is exported..." << endl;
   } else {
-    fftw_export_wisdom_to_file(wordsOfWisdomW);
+    fftwf_export_wisdom_to_file(wordsOfWisdomW);
     fclose(wordsOfWisdomW);
   }
   // END of Wisdom Export
 
   // clean up
-  fftw_destroy_plan(fFFTplanFORW);
+  fftwf_destroy_plan(fFFTplanFORW);
-  fftw_destroy_plan(fFFTplanBACK);
+  fftwf_destroy_plan(fFFTplanBACK);
   free(fFFTtime);
-  fftw_free(fFFTfreq);
+  fftwf_free(fFFTfreq);
   cout << "TLFDynGssKT::~TLFDynGssKT(): " << fCounter << " full FFT cycles needed..." << endl;
 }
 
 double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 3); // three parameters nuL=gbar*B,sigma,fluct.rate nu
+  assert(par.size() == 3 || par.size() == 4); // three parameters nuL=gbar*B,sigma,fluct.rate nu (fourth parameter: t[us] for SG-function)
 
   if(t<0.0)
     return 1.0;
@@ -192,7 +192,7 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
     fFirstCall=false;
   }
 
-  for (unsigned int i(0); i<par.size(); i++) {
+  for (unsigned int i(0); i<3; i++) {
     if( fPar[i]-par[i] ) {
       fPar[i] = par[i];
       fCalcNeeded=true;
@@ -263,7 +263,7 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
 */
     // Transform to frequency domain
 
-    fftw_execute(fFFTplanFORW);
+    fftwf_execute(fFFTplanFORW);
 
     // calculate F(s)
 
@@ -279,7 +279,7 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
 
     // Transform back to time domain
 
-    fftw_execute(fFFTplanBACK);
+    fftwf_execute(fFFTplanBACK);
 
 //    for (unsigned int i(0); i<fNSteps; i++) {
 //      fFFTtime[i]=(fDw*TMath::Exp(fC*i*fDt)/TMath::Pi()*fFFTtime[i]);
@@ -296,6 +296,57 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
   return fDw*exp(fC*t)/PI*fFFTtime[int(t/fDt)];
 }
 
+// LF Dynamic Spin Glass Relaxation
+
+TLFDynSG::TLFDynSG() {
+  for(unsigned int i(0); i<Nsteps; ++i)
+    fLFDynGssIntegral.push_back(new TLFDynGssKT());
+  fPar.resize(3, 0.0);
+}
+
+TLFDynSG::~TLFDynSG() {
+  for(unsigned int i(0); i<Nsteps; ++i)
+    delete fLFDynGssIntegral[i];
+  fLFDynGssIntegral.clear();
+  fPar.clear();
+}
+
+double TLFDynSG::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 3); // three parameters nuL=gbar*B, a ,fluct.rate nu
+
+  if(t<0.0)
+    return 1.0;
+
+  if(t>20.0)
+    return 0.0;
+
+  // Parameters for the integration
+  for (unsigned int i(0); i<3; ++i)
+    fPar[i] = par[i];
+
+  double sigma_step((40.0-0.1)*par[1]/static_cast<double>(Nsteps));
+  // at the moment, integrate between 0.1a and 40a (a "real upper limit" would be 10000a)
+  double integral(0.0);
+  double a_sigmaSq(0.0);
+  double rho(0.0);
+  double rhoNormalizer(0.0); // is [1/]sqrt(2./PI) when the integration is done between 0 and infinity
+
+  for (unsigned int i(0); i<Nsteps; ++i) {
+    fPar[1] = 0.1*par[1] + (static_cast<double>(i) + 0.5)*sigma_step;
+    a_sigmaSq = (par[1]/(fPar[1]*fPar[1]));
+    rho = a_sigmaSq*exp(-0.5*par[1]*a_sigmaSq);
+    rhoNormalizer += rho;
+    integral += (*(fLFDynGssIntegral[i]))(t, fPar)*rho;
+  }
+  integral /= rhoNormalizer;
+
+  //cout << "value at " << t << ": " << integral << endl;
+
+  return integral;
+
+}
+
 // LF Dynamic Lorentz KT
 
 TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("WordsOfWisdom.dat"), fNSteps(524288), fDt(0.000040), fCounter(0), fL1(0.0), fL2(0.0) {
@@ -311,7 +362,7 @@ TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("Words
   if (wordsOfWisdomR == NULL) {
     cout << "TLFDynLorKT::TLFDynLorKT: Couldn't open wisdom file ..." << endl;
   } else {
-    wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR);
+    wisdomLoaded = fftwf_import_wisdom_from_file(wordsOfWisdomR);
     fclose(wordsOfWisdomR);
   }
 
@@ -322,10 +373,10 @@ TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("Words
 
   // allocating memory for the FFtransform pairs and create the FFT plans
 
-  fFFTtime = (double *)malloc(sizeof(double) * fNSteps);
+  fFFTtime = (float *)malloc(sizeof(float) * fNSteps);
-  fFFTfreq = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (fNSteps/2+1));
+  fFFTfreq = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * (fNSteps/2+1));
-  fFFTplanFORW = fftw_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_ESTIMATE);
+  fFFTplanFORW = fftwf_plan_dft_r2c_1d(fNSteps, fFFTtime, fFFTfreq, FFTW_ESTIMATE);
-  fFFTplanBACK = fftw_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_ESTIMATE);
+  fFFTplanBACK = fftwf_plan_dft_c2r_1d(fNSteps, fFFTfreq, fFFTtime, FFTW_ESTIMATE);
 }
 
 TLFDynLorKT::~TLFDynLorKT() {
@@ -335,16 +386,16 @@ TLFDynLorKT::~TLFDynLorKT() {
   if (wordsOfWisdomW == NULL) {
     cout << "TLFDynLorKT::~TLFDynLorKT: Could not open file ... No wisdom is exported..." << endl;
   } else {
-    fftw_export_wisdom_to_file(wordsOfWisdomW);
+    fftwf_export_wisdom_to_file(wordsOfWisdomW);
     fclose(wordsOfWisdomW);
   }
   // END of Wisdom Export
 
   // clean up
-  fftw_destroy_plan(fFFTplanFORW);
+  fftwf_destroy_plan(fFFTplanFORW);
-  fftw_destroy_plan(fFFTplanBACK);
+  fftwf_destroy_plan(fFFTplanBACK);
   free(fFFTtime);
-  fftw_free(fFFTfreq);
+  fftwf_free(fFFTfreq);
   cout << "TLFDynLorKT::~TLFDynLorKT(): " << fCounter << " full FFT cyles needed..." << endl;
 }
 
@@ -457,7 +508,7 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
 
     // Transform to frequency domain
 
-    fftw_execute(fFFTplanFORW);
+    fftwf_execute(fFFTplanFORW);
 
     // calculate F(s)
 
@@ -474,7 +525,7 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
 
     // Transform back to time domain
 
-    fftw_execute(fFFTplanBACK);
+    fftwf_execute(fFFTplanBACK);
 
 //    for (unsigned int i(0); i<fNSteps; i++) {
 //      fFFTtime[i]=(fDw*TMath::Exp(fC*i*fDt)/TMath::Pi()*fFFTtime[i]);

src/external/libLFRelaxation/TLFRelaxation.h

@@ -107,10 +107,10 @@ private:
   double fDt;
   double fDw;
   double fC;
-  fftw_plan fFFTplanFORW;
+  fftwf_plan fFFTplanFORW;
-  fftw_plan fFFTplanBACK;
+  fftwf_plan fFFTplanBACK;
-  double *fFFTtime;
+  float *fFFTtime;
-  fftw_complex *fFFTfreq;
+  fftwf_complex *fFFTfreq;
   mutable unsigned int fCounter;
 
   ClassDef(TLFDynGssKT,1)
@@ -137,10 +137,10 @@ private:
   double fDt;
   double fDw;
   double fC;
-  fftw_plan fFFTplanFORW;
+  fftwf_plan fFFTplanFORW;
-  fftw_plan fFFTplanBACK;
+  fftwf_plan fFFTplanBACK;
-  double *fFFTtime;
+  float *fFFTtime;
-  fftw_complex *fFFTfreq;
+  fftwf_complex *fFFTfreq;
   mutable unsigned int fCounter;
   mutable double fL1;
   mutable double fL2;
@@ -148,6 +148,25 @@ private:
   ClassDef(TLFDynLorKT,1)
 };
 
+class TLFDynSG : public PUserFcnBase {
+
+public:
+  TLFDynSG();
+  ~TLFDynSG();
+
+  Bool_t NeedGlobalPart() const { return false; }
+  void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
+  Bool_t GlobalPartIsValid() const { return true; }
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  mutable vector<double> fPar;
+  vector<TLFDynGssKT*> fLFDynGssIntegral;
+
+  ClassDef(TLFDynSG,1)
+};
+
 class TLFSGInterpolation : public PUserFcnBase {
 
 public:

src/external/libLFRelaxation/TLFRelaxationLinkDef.h

@@ -40,6 +40,7 @@
 #pragma link C++ class TLFStatLorKT+;
 #pragma link C++ class TLFDynGssKT+;
 #pragma link C++ class TLFDynLorKT+;
+#pragma link C++ class TLFDynSG+;
 #pragma link C++ class TLFSGInterpolation+;
 
 #endif //__CINT__

src/tests/CheckMusrfitInstall/ASlibs/test-nonlocal.msr

@@ -24,7 +24,7 @@ FITPARAMETER
 ###############################################################
 THEORY
 asymmetry      3
-userFcn  libPNL_PippardFitter.so PNL_PippardFitter 4 5 6 7 8 9 10 fun2 11
+userFcn  libPNL_PippardFitter PNL_PippardFitter 4 5 6 7 8 9 10 fun2 11
 simpleGss     12          (rate)
 
 ###############################################################

src/tests/CheckMusrfitInstall/BMWlibs/test-libCalcMeanFieldsLEM-DB.msr

@@ -10,7 +10,7 @@ FITPARAMETER
          
 ###############################################################
 THEORY
-userFcn  libCalcMeanFieldsLEM.so TMeanFieldsForScSingleLayer 1 2 3 4 5 5 5
+userFcn  libCalcMeanFieldsLEM TMeanFieldsForScSingleLayer 1 2 3 4 5 5 5
 
 ###############################################################
 #FUNCTIONS

src/tests/CheckMusrfitInstall/BMWlibs/test-libGapIntegrals-ASCII.msr

@@ -9,7 +9,7 @@ FITPARAMETER
 ###############################################################
 THEORY
 asymmetry      1
-userFcn  libGapIntegrals.so TGapDWave 2 3
+userFcn  libGapIntegrals TGapDWave 2 3
 
 ###############################################################
 #FUNCTIONS

src/tests/CheckMusrfitInstall/BMWlibs/test-libLFRelaxation-MUD.msr

@@ -17,7 +17,7 @@ asymmetry   fun1
 dynExpKTLF  fun3     4     5    (frequency damping hopping-rate)
 +
 asymmetry   fun2
-userFcn  libLFRelaxation.so TLFDynLorKT fun3 4 5
+userFcn  libLFRelaxation TLFDynLorKT fun3 4 5
 
 ###############################################################
 FUNCTIONS

src/tests/CheckMusrfitInstall/BMWlibs/test-libTFitPofB_Abrikosov-PSI-BIN.msr

@@ -25,7 +25,7 @@ FITPARAMETER
 THEORY
 asymmetry   fun3
 simpleGss      6          (rate)
-userFcn  libTFitPofB.so TBulkTriVortexML fun2 2 3 4
+userFcn  libTFitPofB TBulkTriVortexML fun2 2 3 4
 +
 asymmetry   fun4
 simpleGss      8          (rate)

src/tests/CheckMusrfitInstall/BMWlibs/test-libTFitPofB_Meissner-ROOT-NPP.msr

@@ -26,7 +26,7 @@ FITPARAMETER
 THEORY
 asymmetry   fun1
 simpleGss   map3          (rate)
-userFcn  libTFitPofB.so TLondon1DHS fun2 map6 1 2 3
+userFcn  libTFitPofB TLondon1DHS fun2 map6 1 2 3
 
 ###############################################################
 FUNCTIONS