diff --git a/src/external/libBNMR/TBNMR.cpp b/src/external/libBNMR/TBNMR.cpp
index 8365b7fd..b49c741a 100644
--- a/src/external/libBNMR/TBNMR.cpp
+++ b/src/external/libBNMR/TBNMR.cpp
@@ -30,9 +30,6 @@
  ***************************************************************************/
 
 #include "TBNMR.h"
-#include "TF1.h"
-#include "Math/WrappedTF1.h"
-#include "Math/GaussIntegrator.h"
 
 #define tau_Li 1210
 #define gamma_Li 6.3018 // In units kHz/mT
@@ -65,6 +62,11 @@ double ExpRlx::operator()(double x, const vector<double> &par) const {
   return y;
 }
 
+
+//initialize Integrators
+TF1 SExpRlx::sexp1=TF1("sexp", "exp(-([0]-x)/[3])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20000.0);
+TF1 SExpRlx::sexp2=TF1("sexp", "exp(-([3]-x)/[4])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20000.0);
+
 double SExpRlx::operator()(double x, const vector<double> &par) const {
   assert(par.size()==3); // make sure the number of parameters handed to the function is correct
   
@@ -72,23 +74,14 @@ double SExpRlx::operator()(double x, const vector<double> &par) const {
   // par[1] is the relaxation rate
   // par[2] is the exponent
 
-  double y;
-
-
-  if ( x >= 0 && x <= par[0] ) { 
-    TF1 sexp("sexp", "exp(-([0]-x)/[3])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20000.0);
-    sexp.SetParameters(x, par[1], par[2],tau_Li);
-    sexp.SetNpx(1000);
-    y=sexp.Integral(0.0,x)/(1-exp(-x/tau_Li))/tau_Li;
+  if ( x >= 0 && x <= par[0] ) {
+    sexp1.SetParameters(x, par[1], par[2],tau_Li);
+    return sexp1.Integral(0.0,x)/(1-exp(-x/tau_Li))/tau_Li;
   } else if ( x > par[0] ) {
-    TF1 sexp("sexp", "exp(-([3]-x)/[4])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20000.0);
-    sexp.SetParameters(x, par[1], par[2], par[0],tau_Li);
-    sexp.SetNpx(1000);
-    y=sexp.Integral(0.0,par[0])/(1-exp(-par[0]/tau_Li))/tau_Li;
-  }  else {
-    y = 0;
+    sexp2.SetParameters(x, par[1], par[2], par[0],tau_Li);
+    return sexp2.Integral(0.0,par[0])/(1-exp(-par[0]/tau_Li))/tau_Li;
   }
-  return y;
+  return 0;
 }
 
 double MLRes::operator()(double x, const vector<double> &par) const {
diff --git a/src/external/libBNMR/TBNMR.h b/src/external/libBNMR/TBNMR.h
index 81c827a1..79dbe22d 100644
--- a/src/external/libBNMR/TBNMR.h
+++ b/src/external/libBNMR/TBNMR.h
@@ -31,6 +31,9 @@
 
 
 #include "PUserFcnBase.h"
+#include "TF1.h"
+#include "Math/WrappedTF1.h"
+#include "Math/GaussIntegrator.h"
 #include <cassert>
 #include <cmath>
 #include <vector>
@@ -63,7 +66,7 @@ class SExpRlx : public PUserFcnBase {
 
 public:
   // default constructor and destructor
-  SExpRlx(){}
+  SExpRlx(){sexp1.SetNpx(1000); sexp2.SetNpx(1000);}
   ~SExpRlx(){}
 
   Bool_t NeedGlobalPart() const { return false; }
@@ -72,6 +75,9 @@ public:
 
   // function operator
   double operator()(double, const vector<double>&) const;
+private:
+  static TF1 sexp1;
+  static TF1 sexp2;
 
   // definition of the class for the ROOT-dictionary
   ClassDef(SExpRlx,1)
@@ -95,4 +101,4 @@ public:
   ClassDef(MLRes,1)
 };
 
-#endif //LIBBNMRH
\ No newline at end of file
+#endif //LIBBNMRH