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

Speed up of SExpRlx by a factor of ~5 by using static TF1s. Might not be thread save.

Browse Source
This commit is contained in:
Jonas A. Krieger 2018-08-15 14:51:51 +02:00
parent e73ded0078
commit 32f8535926
2 changed files with 19 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
src/external/libBNMR/TBNMR.cpp vendored
View File

@ -30,9 +30,6 @@
***************************************************************************/ ***************************************************************************/
#include "TBNMR.h" #include "TBNMR.h"
#include "TF1.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
#define tau_Li 1210 #define tau_Li 1210
#define gamma_Li 6.3018 // In units kHz/mT #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; 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 { 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 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[1] is the relaxation rate
// par[2] is the exponent // par[2] is the exponent
double y;
if ( x >= 0 && x <= par[0] ) { if ( x >= 0 && x <= par[0] ) {
TF1 sexp("sexp", "exp(-([0]-x)/[3])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20000.0); sexp1.SetParameters(x, par[1], par[2],tau_Li);
sexp.SetParameters(x, par[1], par[2],tau_Li); return sexp1.Integral(0.0,x)/(1-exp(-x/tau_Li))/tau_Li;
sexp.SetNpx(1000);
y=sexp.Integral(0.0,x)/(1-exp(-x/tau_Li))/tau_Li;
} else if ( x > par[0] ) { } else if ( x > par[0] ) {
TF1 sexp("sexp", "exp(-([3]-x)/[4])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20000.0); sexp2.SetParameters(x, par[1], par[2], par[0],tau_Li);
sexp.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;
sexp.SetNpx(1000);
y=sexp.Integral(0.0,par[0])/(1-exp(-par[0]/tau_Li))/tau_Li;
} else {
y = 0;
} }
return y; return 0;
} }
double MLRes::operator()(double x, const vector<double> &par) const { double MLRes::operator()(double x, const vector<double> &par) const {

8
src/external/libBNMR/TBNMR.h vendored
View File

@ -31,6 +31,9 @@
#include "PUserFcnBase.h" #include "PUserFcnBase.h"
#include "TF1.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
#include <cassert> #include <cassert>
#include <cmath> #include <cmath>
#include <vector> #include <vector>
@ -63,7 +66,7 @@ class SExpRlx : public PUserFcnBase {
public: public:
// default constructor and destructor // default constructor and destructor
SExpRlx(){} SExpRlx(){sexp1.SetNpx(1000); sexp2.SetNpx(1000);}
~SExpRlx(){} ~SExpRlx(){}
Bool_t NeedGlobalPart() const { return false; } Bool_t NeedGlobalPart() const { return false; }
@ -72,6 +75,9 @@ public:
// function operator // function operator
double operator()(double, const vector<double>&) const; double operator()(double, const vector<double>&) const;
private:
static TF1 sexp1;
static TF1 sexp2;
// definition of the class for the ROOT-dictionary // definition of the class for the ROOT-dictionary
ClassDef(SExpRlx,1) ClassDef(SExpRlx,1)