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.

src/classes/PRunAsymmetryBNMR.cpp

@@ -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)
+    if (fForwardp[i] != 0.0)
-      dval = TMath::Sqrt(fForward[i]);
+      dval = TMath::Sqrt(fForwardp[i]);
     else
       dval = 1.0;
-    fForwardErr.push_back(dval);
+    fForwardpErr.push_back(dval);
-    fForward[i] -= fRunInfo->GetBkgFix(0);
+    fForwardp[i] -= fRunInfo->GetBkgFix(0);
 
     // keep the error, and make sure that the bin is NOT empty
-    if (fBackward[i] != 0.0)
+    if (fForwardm[i] != 0.0)
-      dval = TMath::Sqrt(fBackward[i]);
+      dval = TMath::Sqrt(fForwardm[i]);
     else
       dval = 1.0;
-    fBackwardErr.push_back(dval);
+    fForwardmErr.push_back(dval);
-    fBackward[i] -= fRunInfo->GetBkgFix(1);
+    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 bkgm[2]    = {0.0, 0.0};
-  Double_t errBkgn[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++) {

src/external/libBNMR/ExpRlx-MUD.msr

@@ -4,52 +4,36 @@
 FITPARAMETER
 ###############################################################
 #     No     Name       Value    Err     Min  Max                  
-        1 Alpha       1.11587   0.00038     none        
+        1 Alpha       1         0           none        
-        2 Asy         0.0570    0.0011      none        
+        2 Asy         0.05615   0.00073     none        
         3 T           1         0           none        
-        4 Rlx         1.015     0.023       none        
+        4 Rlx         1.046     0.024       none        
-        5 One         1         0           none        
-        6 Bet         -1.052    0.027       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 
+forward         3 5 
-backward        4 
+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    
+t0              10.0    10.0    10.0    10.0    
-map             2    5    0    0    0    0    0    0    0    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
+STATISTIC --- 2018-08-20 17:31:40
-  chisq = 399.5, NDF = 290, chisq/NDF = 1.377736
+  chisq = 309.1, NDF = 145, chisq/NDF = 2.131395

src/external/libBNMR/SExpRlx-MUD.msr

@@ -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

src/external/libBNMR/TBNMR.cpp

@@ -31,53 +31,48 @@
 
 #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 std::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[0] time of beam on (pulse length) in seconds
-  // par[1] is the relaxation rate
+  // 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]));
+    return (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 0;
-  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 std::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[0] beam of beam on (pulse length) in seconds 
-  // par[1] is the relaxation rate
+  // 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;

src/external/libBNMR/TBNMR.h

@@ -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; }