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

Resolve conflict from Jonas' commits

Browse Source
This commit is contained in:
salman
2018-08-18 14:17:41 +02:00
parent 3cdf3d356e
commit ef653c11b2
2 changed files with 26 additions and 93 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -10,8 +10,8 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2010 by Zaher Salman *
* zaher.salman@psi.ch *
* Copyright (C) 2010 by Zaher Salman *
* zaher.salman@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
@ -20,7 +20,7 @@
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
@ -30,29 +30,15 @@
***************************************************************************/
#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
#define PI 3.14159265358979323846
#define TWOPI 6.28318530717958647692
ClassImp(TBNMR) // for the ROOT-dictionary
ClassImp(ExpRlx)
ClassImp(ExpRlx) // for the ROOT-dictionary
ClassImp(SExpRlx)
ClassImp(MLRes)
double TBNMR::operator()(double x, const vector<double> &par) const {
assert(par.size()==1); // make sure the number of parameters handed to the function is correct
double arg(par[0]*x);
if(!arg)
return 1.0;
return sin(arg)/arg;
}
double ExpRlx::operator()(double x, const vector<double> &par) const {
assert(par.size()==2); // make sure the number of parameters handed to the function is correct
@ -75,55 +61,25 @@ 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
// par[0] time of beam off
// par[1] is the relaxation rate
// par[2] is the exponent
double tau_p;
double y;
tau_p = (tau_Li/(1.+par[1]*tau_Li));
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;
}
double MLRes::operator()(double x, const vector<double> &par) const {
assert(par.size()==3); // make sure the number of parameters handed to the function is correct
// par[0] time of beam off
// 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, 10000.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;
} else if ( x > par[0] ) {
TF1 sexp("sexp", "exp(-([3]-x)/[4])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 10000.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;
}
return y;
return 0;
}

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

@ -31,29 +31,18 @@
#include "PUserFcnBase.h"
#include "TF1.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
#include <cassert>
#include <cmath>
#include <vector>
#ifndef LIBBNMRH
#define LIBBNMRH
using namespace std;
class TBNMR : public PUserFcnBase {
public:
// default constructor and destructor
TBNMR(){}
~TBNMR(){}
Bool_t NeedGlobalPart() const { return false; }
void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
Bool_t GlobalPartIsValid() const { return true; }
// function operator
double operator()(double, const vector<double>&) const;
// definition of the class for the ROOT-dictionary
ClassDef(TBNMR,1)
};
class ExpRlx : public PUserFcnBase {
@ -77,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; }
@ -86,25 +75,13 @@ 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)
};
class MLRes : public PUserFcnBase {
public:
// default constructor and destructor
MLRes(){}
~MLRes(){}
Bool_t NeedGlobalPart() const { return false; }
void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) { }
Bool_t GlobalPartIsValid() const { return true; }
// function operator
double operator()(double, const vector<double>&) const;
// definition of the class for the ROOT-dictionary
ClassDef(MLRes,1)
};
#endif //LIBBNMRH