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

Beta-NMR helicity subtraction seems to work in the new class. Still needs debugging and tests. Some of the changes from Jonas were reverted since SExp does not work otherwise.

Browse Source
This commit is contained in:
salman 2018-08-20 17:51:39 +02:00
parent 52a5bc25b1
commit e9ec126d9e
5 changed files with 131 additions and 93 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
src/classes/PRunAsymmetryBNMR.cpp
View File

@ -539,11 +539,11 @@ Bool_t PRunAsymmetryBNMR::PrepareData()
return false;
}
// keep the time resolution in (us)
// keep the time resolution in (ms)
// possibility to rescale for betaNMR
fTimeResolution = runData->GetTimeResolution()/1.0e3;
cout.precision(10);
cout << endl << ">> PRunAsymmetryBNMR::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ns)" << endl;
cout << endl << ">> PRunAsymmetryBNMR::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ms)" << endl;
// get all the proper t0's and addt0's for the current RUN block
if (!GetProperT0(runData, globalBlock, forwardHistoNo, backwardHistoNo)) {
@ -602,6 +602,8 @@ Bool_t PRunAsymmetryBNMR::PrepareData()
// set forward/backward histo data of the first group
fForwardp.resize(forward[0].size());
fBackwardp.resize(backward[0].size());
fForwardm.resize(forward[0].size());
fBackwardm.resize(backward[0].size());
for (UInt_t i=0; i<fForwardp.size(); i++) {
fForwardp[i] = forward[0][i];
fBackwardp[i] = backward[0][i];
@ -609,26 +611,6 @@ Bool_t PRunAsymmetryBNMR::PrepareData()
fBackwardm[i] = backward[1][i];
}
// // group histograms, add all the remaining forward histograms of the group
// for (UInt_t i=1; i<forwardHistoNo.size(); i++) { // loop over the groupings
// for (UInt_t j=0; j<runData->GetDataBin(forwardHistoNo[i])->size(); j++) { // loop over the bin indices
// // make sure that the index stays within proper range
// if ((j+fT0s[2*i]-fT0s[0] >= 0) && (j+fT0s[2*i]-fT0s[0] < runData->GetDataBin(forwardHistoNo[i])->size())) {
// fForward[j] += forward[i][j+(Int_t)fT0s[2*i]-(Int_t)fT0s[0]];
// }
// }
// }
// // group histograms, add all the remaining backward histograms of the group
// for (UInt_t i=1; i<backwardHistoNo.size(); i++) { // loop over the groupings
// for (UInt_t j=0; j<runData->GetDataBin(backwardHistoNo[i])->size(); j++) { // loop over the bin indices
// // make sure that the index stays within proper range
// if ((j+fT0s[2*i+1]-fT0s[1] >= 0) && (j+fT0s[2*i+1]-fT0s[1] < runData->GetDataBin(backwardHistoNo[i])->size())) {
// fBackward[j] += backward[i][j+(Int_t)fT0s[2*i+1]-(Int_t)fT0s[1]];
// }
// }
// }
// subtract background from histogramms ------------------------------------------
if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given
if (fRunInfo->GetBkgRange(0) >= 0) { // background range given
@ -706,22 +688,40 @@ Bool_t PRunAsymmetryBNMR::PrepareData()
Bool_t PRunAsymmetryBNMR::SubtractFixBkg()
{
Double_t dval;
for (UInt_t i=0; i<fForward.size(); i++) {
// Order in RunInfo structure Fp, Fm, Bp, Bm
for (UInt_t i=0; i<fForwardp.size(); i++) {
// keep the error, and make sure that the bin is NOT empty
if (fForward[i] != 0.0)
dval = TMath::Sqrt(fForward[i]);
if (fForwardp[i] != 0.0)
dval = TMath::Sqrt(fForwardp[i]);
else
dval = 1.0;
fForwardErr.push_back(dval);
fForward[i] -= fRunInfo->GetBkgFix(0);
fForwardpErr.push_back(dval);
fForwardp[i] -= fRunInfo->GetBkgFix(0);
// keep the error, and make sure that the bin is NOT empty
if (fBackward[i] != 0.0)
dval = TMath::Sqrt(fBackward[i]);
if (fForwardm[i] != 0.0)
dval = TMath::Sqrt(fForwardm[i]);
else
dval = 1.0;
fBackwardErr.push_back(dval);
fBackward[i] -= fRunInfo->GetBkgFix(1);
fForwardmErr.push_back(dval);
fForwardm[i] -= fRunInfo->GetBkgFix(1);
// keep the error, and make sure that the bin is NOT empty
if (fBackwardp[i] != 0.0)
dval = TMath::Sqrt(fBackwardp[i]);
else
dval = 1.0;
fBackwardpErr.push_back(dval);
fBackwardp[i] -= fRunInfo->GetBkgFix(2);
// keep the error, and make sure that the bin is NOT empty
if (fBackwardm[i] != 0.0)
dval = TMath::Sqrt(fBackwardm[i]);
else
dval = 1.0;
fBackwardmErr.push_back(dval);
fBackwardm[i] -= fRunInfo->GetBkgFix(3);
}
return true;
@ -807,8 +807,8 @@ Bool_t PRunAsymmetryBNMR::SubtractEstimatedBkg()
// calculate background
Double_t bkgp[2] = {0.0, 0.0};
Double_t errBkgp[2] = {0.0, 0.0};
Double_t bkgn[2] = {0.0, 0.0};
Double_t errBkgn[2] = {0.0, 0.0};
Double_t bkgm[2] = {0.0, 0.0};
Double_t errBkgm[2] = {0.0, 0.0};
// forward
for (UInt_t i=start[0]; i<=end[0]; i++) {
@ -898,7 +898,7 @@ Bool_t PRunAsymmetryBNMR::PrepareFitData()
PRunData forwardpPacked;
PRunData backwardpPacked;
PRunData forwardmPacked;
PRunData backwarmpPacked;
PRunData backwardmPacked;
Double_t valuep = 0.0;
Double_t errorp = 0.0;
Double_t valuem = 0.0;
@ -1150,6 +1150,8 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2
Double_t errorp = 0.0;
Double_t valuem = 0.0;
Double_t errorm = 0.0;
Double_t value = 0.0;
Double_t error = 0.0;
// forward
for (Int_t i=start[0]; i<end[0]; i++) {

42
src/external/libBNMR/ExpRlx-MUD.msr vendored
View File

@ -4,52 +4,36 @@
FITPARAMETER
###############################################################
# No Name Value Err Min Max
1 Alpha 1.11587 0.00038 none
2 Asy 0.0570 0.0011 none
1 Alpha 1 0 none
2 Asy 0.05615 0.00073 none
3 T 1 0 none
4 Rlx 1.015 0.023 none
5 One 1 0 none
6 Bet -1.052 0.027 none
4 Rlx 1.046 0.024 none
###############################################################
THEORY
###############################################################
asymmetry fun1
asymmetry 2
userFcn .libs/libBNMR.so ExpRlx 3 4
###############################################################
FUNCTIONS
#FUNCTIONS
###############################################################
fun1 = 0.5 * map1 * map2
#fun1 = 0.5 * map1 * map2
###############################################################
RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format)
fittype 5 (beta-NMR fit)
alpha 1
forward 3
backward 4
forward 3 5
backward 4 6
data 11 800 11 800
#backgr.fix 0
background 1 9 1 9 # estimated bkg: 8.0000 / 9.6250
t0 10.0 10.0
map 2 5 0 0 0 0 0 0 0 0
t0 10.0 10.0 10.0 10.0
map 0 0 0 0 0 0 0 0 0 0
fit 0.5 8
packing 5
RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format)
fittype 5 (beta-NMR fit)
alpha 1
forward 5
backward 6
data 11 800 11 800
#backgr.fix 0
background 1 9 1 9 # estimated bkg: 11.6250 / 15.6250
t0 10.0 10.0
map 2 6 0 0 0 0 0 0 0 0
fit 0.5 8
packing 5
###############################################################
COMMANDS
MINIMIZE
@ -58,7 +42,7 @@ SAVE
###############################################################
PLOT 5 (beta-NMR asymmetry plot)
runs 1 2
runs 1
use_fit_ranges
view_packing 10
@ -72,5 +56,5 @@ plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_RE
phase 8
#range FRQMIN FRQMAX
###############################################################
STATISTIC --- 2018-08-18 23:43:35
chisq = 399.5, NDF = 290, chisq/NDF = 1.377736
STATISTIC --- 2018-08-20 17:31:40
chisq = 309.1, NDF = 145, chisq/NDF = 2.131395

61
src/external/libBNMR/SExpRlx-MUD.msr vendored Normal file
View File

@ -0,0 +1,61 @@
# Run Numbers: 1111
###############################################################
FITPARAMETER
###############################################################
# No Name Value Err Min Max
1 Alpha 1 0 none
2 Asy 0.108 0.011 none 0 0.2
3 T 1 0 none
4 Rlx 3.63 0.89 none 0 15000
5 Beta 0.452 0.039 none 0.3 2
###############################################################
THEORY
###############################################################
asymmetry 2
userFcn .libs/libBNMR.so SExpRlx 3 4 5
###############################################################
#FUNCTIONS
###############################################################
#fun1 = 0.5 * map1 * map2
###############################################################
RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format)
fittype 5 (beta-NMR fit)
alpha 1
forward 3 5
backward 4 6
data 11 800 11 800
#backgr.fix 0
background 1 9 1 9 # estimated bkg: 8.0000 / 9.6250
t0 10.0 10.0 10.0 10.0
map 0 0 0 0 0 0 0 0 0 0
fit 0.5 8
packing 5
###############################################################
COMMANDS
MINIMIZE
HESSE
SAVE
###############################################################
PLOT 5 (beta-NMR asymmetry plot)
runs 1
use_fit_ranges
view_packing 10
###############################################################
FOURIER
units MHz # units either 'Gauss', 'Tesla', 'MHz', or 'Mc/s'
fourier_power 12
apodization STRONG # NONE, WEAK, MEDIUM, STRONG
plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL
phase 8
#range FRQMIN FRQMAX
###############################################################
STATISTIC --- 2018-08-20 17:48:41
chisq = 177.2, NDF = 144, chisq/NDF = 1.230711

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

@ -31,53 +31,43 @@
#include "TBNMR.h"
#define tau_Li 1210
#define gamma_Li 6.3018 // In units kHz/mT
#define PI 3.14159265358979323846
#define TWOPI 6.28318530717958647692
ClassImp(ExpRlx) // for the ROOT-dictionary
ClassImp(SExpRlx)
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
// par[0] time of beam off
// par[1] is the relaxation rate
// par[0] time of beam on (pulse length) in seconds
// par[1] is the relaxation rate in 1/s
double tau_p;
double y;
tau_p = (tau_Li/(1.+par[1]*tau_Li));
// x should be in seconds, otherwise it should be rescaled here
if ( x <= par[0] && x >= 0) {
y=(tau_p/tau_Li)*(1-exp(-x/tau_p))/(1-exp(-x/tau_Li));
return (tau_p/tau_Li)*(1-exp(-x/tau_p))/(1-exp(-x/tau_Li));
} else if ( x > par[0] ){
y=(tau_p/tau_Li)*(1-exp(-par[0]/tau_p))/(1-exp(-par[0]/tau_Li))*exp(-par[1]*(x-par[0]));
} else {
y = 0;
return (tau_p/tau_Li)*(1-exp(-par[0]/tau_p))/(1-exp(-par[0]/tau_Li))*exp(-par[1]*(x-par[0]));
}
return y;
return 0;
}
//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[0] beam of beam on (pulse length) in seconds
// par[1] is the relaxation rate in 1/s
// par[2] is the exponent
// x should be in seconds, otherwise it should be rescaled here
if ( x >= 0 && x <= par[0] ) {
TF1 sexp1("sexp1", "exp(-([0]-x)/[3])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20.0);
sexp1.SetParameters(x, par[1], par[2],tau_Li);
sexp1.SetNpx(1000);
return sexp1.Integral(0.0,x)/(1-exp(-x/tau_Li))/tau_Li;
} else if ( x > par[0] ) {
TF1 sexp2("sexp2", "exp(-([3]-x)/[4])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20.0);
sexp2.SetParameters(x, par[1], par[2], par[0],tau_Li);
sexp2.SetNpx(1000);
return sexp2.Integral(0.0,par[0])/(1-exp(-par[0]/tau_Li))/tau_Li;
}
return 0;

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

@ -41,6 +41,10 @@
#ifndef LIBBNMRH
#define LIBBNMRH
#define tau_Li 1.210 // In seconds
#define PI 3.14159265358979323846
#define TWOPI 6.28318530717958647692
using namespace std;
@ -66,7 +70,7 @@ class SExpRlx : public PUserFcnBase {
public:
// default constructor and destructor
SExpRlx(){sexp1.SetNpx(1000); sexp2.SetNpx(1000);}
SExpRlx(){}
~SExpRlx(){}
Bool_t NeedGlobalPart() const { return false; }
@ -75,9 +79,6 @@ 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)