diff --git a/CMakeLists.txt b/CMakeLists.txt
index aef7326a..3875fe16 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -5,7 +5,7 @@ if (CMAKE_VERSION GREATER_EQUAL 3.12)
   cmake_policy(SET CMP0075 NEW)
 endif (CMAKE_VERSION GREATER_EQUAL 3.12)
 
-project(musrfit VERSION 1.6.0 LANGUAGES C CXX)
+project(musrfit VERSION 1.6.1 LANGUAGES C CXX)
 
 #--- musrfit specific options -------------------------------------------------
 option(nexus "build optional NeXus support. Needed for ISIS" OFF)
diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp
index 5bedae20..fd4f70fb 100644
--- a/src/classes/PMsrHandler.cpp
+++ b/src/classes/PMsrHandler.cpp
@@ -819,15 +819,24 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
         } else if (sstr.BeginsWith("alpha ")) {
           fout.width(16);
           fout << std::left << "alpha";
-          fout << fRuns[runNo].GetAlphaParamNo() << std::endl;
+          // check of alpha is given as a function
+          if (fRuns[runNo].GetAlphaParamNo() >= MSR_PARAM_FUN_OFFSET)
+            fout << "fun" << fRuns[runNo].GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+          else
+            fout << fRuns[runNo].GetAlphaParamNo() << std::endl;
+          fout << std::endl;
         } else if (sstr.BeginsWith("beta ")) {
           fout.width(16);
           fout << std::left << "beta";
-          fout << fRuns[runNo].GetBetaParamNo()  << std::endl;
+          if (fRuns[runNo].GetBetaParamNo() >= MSR_PARAM_FUN_OFFSET)
+            fout << "fun" << fRuns[runNo].GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+          else
+            fout << fRuns[runNo].GetBetaParamNo()  << std::endl;
+          fout << std::endl;
         } else if (sstr.BeginsWith("norm")) {
           fout.width(16);
           fout << std::left << "norm";
-          // check if norm is give as a function
+          // check if norm is given as a function
           if (fRuns[runNo].GetNormParamNo() >= MSR_PARAM_FUN_OFFSET)
             fout << "fun" << fRuns[runNo].GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
           else
@@ -1939,14 +1948,24 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, std::map<UInt_t, TString
     if (fRuns[i].GetAlphaParamNo() != -1) {
       fout.width(16);
       fout << std::left << "alpha";
-      fout << fRuns[i].GetAlphaParamNo() << std::endl;
+      // check if alpha is give as a function
+      if (fRuns[i].GetAlphaParamNo() >= MSR_PARAM_FUN_OFFSET)
+        fout << "fun" << fRuns[i].GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+      else
+        fout << fRuns[i].GetAlphaParamNo();
+      fout << std::endl;
     }
 
     // beta
     if (fRuns[i].GetBetaParamNo() != -1) {
       fout.width(16);
       fout << std::left << "beta";
-      fout << fRuns[i].GetBetaParamNo()  << std::endl;
+      // check if beta is give as a function
+      if (fRuns[i].GetBetaParamNo() >= MSR_PARAM_FUN_OFFSET)
+        fout << "fun" << fRuns[i].GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+      else
+        fout << fRuns[i].GetBetaParamNo();
+      fout << std::endl;
     }
 
     // norm
@@ -3379,6 +3398,12 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines)
             param.SetAlphaParamNo(ival);
           else
             error = true;
+        } else if (str.Contains("fun")) {
+          Int_t no;
+          if (FilterNumber(str, "fun", MSR_PARAM_FUN_OFFSET, no))
+            param.SetAlphaParamNo(no);
+          else
+            error = true;
         } else {
           error = true;
         }
@@ -3399,6 +3424,12 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines)
           ival = str.Atoi();
           if (ival > 0)
             param.SetBetaParamNo(ival);
+          else
+            error = true;          
+        } else if (str.Contains("fun")) {
+          Int_t no;
+          if (FilterNumber(str, "fun", MSR_PARAM_FUN_OFFSET, no))
+            param.SetBetaParamNo(no);
           else
             error = true;
         } else {
@@ -4582,9 +4613,12 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
 
     // handle a single PLOT block
     while ((iter1 != iter2) && !error) {
+      TString line = iter1->fLine;
+      if (line.First('#') != -1) // remove trailing comment before proceed
+        line.Resize(line.First('#'));
 
-      if (iter1->fLine.Contains("PLOT")) { // handle plot header
-        tokens = iter1->fLine.Tokenize(" \t");
+      if (line.Contains("PLOT")) { // handle plot header
+        tokens = line.Tokenize(" \t");
         if (!tokens) {
           std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo;
           std::cerr << std::endl << std::endl;
@@ -4605,9 +4639,9 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
           delete tokens;
           tokens = nullptr;
         }
-      } else if (iter1->fLine.Contains("lifetimecorrection", TString::kIgnoreCase)) {
+      } else if (line.Contains("lifetimecorrection", TString::kIgnoreCase)) {
         param.fLifeTimeCorrection = true;
-      } else if (iter1->fLine.Contains("runs", TString::kIgnoreCase)) { // handle plot runs
+      } else if (line.Contains("runs", TString::kIgnoreCase)) { // handle plot runs
         TComplex run;
         PStringNumberList *rl;
         std::string errorMsg;
@@ -4623,7 +4657,7 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
           case MSR_PLOT_ASYM_RRF:
           case MSR_PLOT_NON_MUSR:
           case MSR_PLOT_MU_MINUS:
-            rl = new PStringNumberList(iter1->fLine.Data());
+            rl = new PStringNumberList(line.Data());
             if (!rl->Parse(errorMsg, true)) {
               std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo;
               std::cerr << std::endl << ">>   Error Message: " << errorMsg;
@@ -4643,14 +4677,14 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
             error = true;
             break;
         }
-      } else if (iter1->fLine.Contains("range ", TString::kIgnoreCase)) { // handle plot range
+      } else if (line.Contains("range ", TString::kIgnoreCase)) { // handle plot range
         // remove previous entries
         param.fTmin.clear();
         param.fTmax.clear();
         param.fYmin.clear();
         param.fYmax.clear();
 
-        tokens = iter1->fLine.Tokenize(" \t");
+        tokens = line.Tokenize(" \t");
         if (!tokens) {
           std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo;
           std::cerr << std::endl << std::endl;
@@ -4700,14 +4734,14 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
           delete tokens;
           tokens = nullptr;
         }
-      } else if (iter1->fLine.Contains("sub_ranges", TString::kIgnoreCase)) {
+      } else if (line.Contains("sub_ranges", TString::kIgnoreCase)) {
         // remove previous entries
         param.fTmin.clear();
         param.fTmax.clear();
         param.fYmin.clear();
         param.fYmax.clear();
 
-        tokens = iter1->fLine.Tokenize(" \t");
+        tokens = line.Tokenize(" \t");
         if (!tokens) {
           std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo;
           std::cerr << std::endl << std::endl;
@@ -4762,11 +4796,11 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
           delete tokens;
           tokens = nullptr;
         }
-      } else if (iter1->fLine.Contains("use_fit_ranges", TString::kIgnoreCase)) {
+      } else if (line.Contains("use_fit_ranges", TString::kIgnoreCase)) {
         param.fUseFitRanges = true;
         // check if y-ranges are given
 
-        tokens = iter1->fLine.Tokenize(" \t");
+        tokens = line.Tokenize(" \t");
         if (!tokens) {
           std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo;
           std::cerr << std::endl << std::endl;
@@ -5207,12 +5241,20 @@ UInt_t PMsrHandler::GetNoOfFitParameters(UInt_t idx)
     paramVector.push_back(fRuns[idx].GetBkgFitParamNo());
 
   // get alpha parameter if present (asymmetry fit)
-  if (fRuns[idx].GetAlphaParamNo() != -1)
-    paramVector.push_back(fRuns[idx].GetAlphaParamNo());
+  if (fRuns[idx].GetAlphaParamNo() != -1) {
+    if (fRuns[idx].GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) // parameter
+      paramVector.push_back(fRuns[idx].GetAlphaParamNo());
+    else // function
+      funVector.push_back(fRuns[idx].GetAlphaParamNo() - MSR_PARAM_FUN_OFFSET);
+  }
 
   // get beta parameter if present (asymmetry fit)
-  if (fRuns[idx].GetBetaParamNo() != -1)
-    paramVector.push_back(fRuns[idx].GetBetaParamNo());
+  if (fRuns[idx].GetBetaParamNo() != -1) {
+    if (fRuns[idx].GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) // parameter
+      paramVector.push_back(fRuns[idx].GetBetaParamNo());
+    else // function
+      funVector.push_back(fRuns[idx].GetBetaParamNo() - MSR_PARAM_FUN_OFFSET);
+  }
 
   // go through the theory block and collect parameters
   // possible entries: number -> parameter, fun<number> -> function, map<number> -> maps
diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp
index e1586867..d4fa8e91 100644
--- a/src/classes/PRunAsymmetry.cpp
+++ b/src/classes/PRunAsymmetry.cpp
@@ -111,7 +111,16 @@ PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, UIn
     return;
   }
   // check if alpha parameter is within proper bounds
-  if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > static_cast<Int_t>(param->size()))) {
+  if (fRunInfo->GetAlphaParamNo() >= MSR_PARAM_FUN_OFFSET) { // alpha seems to be a function
+    if ((fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET < 0) ||
+        (fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET > msrInfo->GetNoOfFuncs())) {
+      std::cerr << std::endl << ">> PRunAsymmetry::PRunAsymmetry(): **ERROR** alpha parameter is a function with no = " << fRunInfo->GetAlphaParamNo();
+      std::cerr << std::endl << ">> This is out of bound, since there are only " << msrInfo->GetNoOfFuncs() << " functions.";
+      std::cerr << std::endl;
+      fValid = false;
+      return;
+    }
+  } else if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > static_cast<Int_t>(param->size()))) {
     std::cerr << std::endl << ">> PRunAsymmetry::PRunAsymmetry(): **ERROR** alpha parameter no = " << fRunInfo->GetAlphaParamNo();
     std::cerr << std::endl << ">> This is out of bound, since there are only " << param->size() << " parameters.";
     std::cerr << std::endl;
@@ -120,10 +129,11 @@ PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, UIn
   }
   // check if alpha is fixed
   Bool_t alphaFixedToOne = false;
-  if (((*param)[fRunInfo->GetAlphaParamNo()-1].fStep == 0.0) &&
-      ((*param)[fRunInfo->GetAlphaParamNo()-1].fValue == 1.0))
-    alphaFixedToOne = true;
-
+  if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is parameter
+    if (((*param)[fRunInfo->GetAlphaParamNo()-1].fStep == 0.0) &&
+        ((*param)[fRunInfo->GetAlphaParamNo()-1].fValue == 1.0))
+      alphaFixedToOne = true;
+  }
   // check if beta is given
   Bool_t betaFixedToOne = false;
   if (fRunInfo->GetBetaParamNo() == -1) { // no beta given hence assuming beta == 1
@@ -196,6 +206,58 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
   Double_t time(1.0);
   Int_t i;
 
+  // determine alpha/beta
+  switch (fAlphaBetaTag) {
+    case 1: // alpha == 1, beta == 1
+      a = 1.0;
+      b = 1.0;
+      break;
+    case 2: // alpha != 1, beta == 1
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        a = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      b = 1.0;
+      break;
+    case 3: // alpha == 1, beta != 1
+      a = 1.0;
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        b = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      break;
+    case 4: // alpha != 1, beta != 1
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        a = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        b = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      break;
+    default:
+      a = 1.0;
+      b = 1.0;
+      break;
+  }
+
   // Calculate the theory function once to ensure one function evaluation for the current set of parameters.
   // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
   // for a given set of parameters---which should be done outside of the parallelized loop.
@@ -206,34 +268,12 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
   Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
   if (chunk < 10)
     chunk = 10;
-  #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq)
+  #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,f) schedule(dynamic,chunk) reduction(+:chisq)
   #endif
   for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
     time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
-    switch (fAlphaBetaTag) {
-      case 1: // alpha == 1, beta == 1
-        asymFcnValue = fTheory->Func(time, par, fFuncValues);
-        break;
-      case 2: // alpha != 1, beta == 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
-        f = fTheory->Func(time, par, fFuncValues);
-        asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0));
-        break;
-      case 3: // alpha == 1, beta != 1
-        b = par[fRunInfo->GetBetaParamNo()-1];
-        f = fTheory->Func(time, par, fFuncValues);
-        asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0));
-        break;
-      case 4: // alpha != 1, beta != 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
-        b = par[fRunInfo->GetBetaParamNo()-1];
-        f = fTheory->Func(time, par, fFuncValues);
-        asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
-        break;
-      default:
-        asymFcnValue = 0.0;
-        break;
-    }
+    f = fTheory->Func(time, par, fFuncValues);
+    asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
     diff = fData.GetValue()->at(i) - asymFcnValue;
     chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
   }
@@ -429,18 +469,46 @@ void PRunAsymmetry::CalcTheory()
         asymFcnValue = fTheory->Func(time, par, fFuncValues);
         break;
       case 2: // alpha != 1, beta == 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
+        if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+          a = par[fRunInfo->GetAlphaParamNo()-1];
+        } else { // alpha is function
+          // get function number
+          UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0));
         break;
       case 3: // alpha == 1, beta != 1
-        b = par[fRunInfo->GetBetaParamNo()-1];
+        if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+          b = par[fRunInfo->GetBetaParamNo()-1];
+        } else { // beta is a function
+          // get function number
+          UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0));
         break;
       case 4: // alpha != 1, beta != 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
-        b = par[fRunInfo->GetBetaParamNo()-1];
+        if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+          a = par[fRunInfo->GetAlphaParamNo()-1];
+        } else { // alpha is function
+          // get function number
+          UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
+        if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+          b = par[fRunInfo->GetBetaParamNo()-1];
+        } else { // beta is a function
+          // get function number
+          UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
         break;
@@ -1138,16 +1206,44 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2])
       beta  = 1.0;
       break;
     case 2: // alpha != 1, beta == 1
-      alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       beta  = 1.0;
       break;
     case 3: // alpha == 1, beta != 1
       alpha = 1.0;
-      beta  = par[fRunInfo->GetBetaParamNo()-1];
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        beta = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       break;
     case 4: // alpha != 1, beta != 1
-      alpha = par[fRunInfo->GetAlphaParamNo()-1];
-      beta  = par[fRunInfo->GetBetaParamNo()-1];
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        beta = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       break;
     default:
       break;
@@ -1349,16 +1445,44 @@ Bool_t PRunAsymmetry::PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]
       beta  = 1.0;
       break;
     case 2: // alpha != 1, beta == 1
-      alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       beta  = 1.0;
       break;
     case 3: // alpha == 1, beta != 1
       alpha = 1.0;
-      beta  = par[fRunInfo->GetBetaParamNo()-1];
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        beta = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       break;
     case 4: // alpha != 1, beta != 1
-      alpha = par[fRunInfo->GetAlphaParamNo()-1];
-      beta  = par[fRunInfo->GetBetaParamNo()-1];
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        beta = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       break;
     default:
       break;
diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp
index 9b95c074..66b7626c 100644
--- a/src/classes/PRunAsymmetryBNMR.cpp
+++ b/src/classes/PRunAsymmetryBNMR.cpp
@@ -206,56 +206,90 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector<Double_t>& par)
   Double_t time(1.0),alphaest;
   Int_t i;
 
+  // determine alpha/beta
+  alphaest = fRunInfo->GetEstimatedAlpha();
+  switch (fAlphaBetaTag) {
+    case 1: // alpha == 1, beta == 1
+      a = 1.0;
+      b = 1.0;
+      break;
+    case 2: // alpha != 1, beta == 1
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        a = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      b = 1.0;
+      break;
+    case 3: // alpha == 1, beta != 1
+      a = 1.0;
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        b = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      break;
+    case 4: // alpha != 1, beta != 1
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        a = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        b = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      break;
+    case 5: // alpha ?? , beta == 1
+      a = alphaest;
+      b = 1.0;
+      break;
+    case 6: // alpha ??, beta != 1
+      a = alphaest;
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        b = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      break;
+    default:
+      a = 1.0;
+      b = 1.0;
+      break;
+  }
+
   // Calculate the theory function once to ensure one function evaluation for the current set of parameters.
   // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
   // for a given set of parameters---which should be done outside of the parallelized loop.
   // For all other functions it means a tiny and acceptable overhead.
   asymFcnValue = fTheory->Func(time, par, fFuncValues);
-  alphaest = fRunInfo->GetEstimatedAlpha();
 
   #ifdef HAVE_GOMP
   Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
   if (chunk < 10)
     chunk = 10;
-  #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq)
+  #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,f) schedule(dynamic,chunk) reduction(+:chisq)
   #endif
   for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
     time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
-    switch (fAlphaBetaTag) {
-      case 1: // alpha == 1, beta == 1
-        asymFcnValue = fTheory->Func(time, par, fFuncValues);
-        break;
-      case 2: // alpha != 1, beta == 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
-        f = fTheory->Func(time, par, fFuncValues)/2.0;
-        asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0));
-        break;
-      case 3: // alpha == 1, beta != 1
-        b = par[fRunInfo->GetBetaParamNo()-1];
-        f = fTheory->Func(time, par, fFuncValues)/2.0;
-        asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0))-f*(b+1.0)/(2.0+f*(b-1.0));
-        break;
-      case 4: // alpha != 1, beta != 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
-        b = par[fRunInfo->GetBetaParamNo()-1];
-        f = fTheory->Func(time, par, fFuncValues)/2.0;
-        asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0));
-        break;
-      case 5: // alpha ?? , beta == 1
-        a = alphaest;
-        f = fTheory->Func(time, par, fFuncValues)/2.0;
-        asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0));
-        break;
-      case 6: // alpha ??, beta != 1
-        a = alphaest;
-        b = par[fRunInfo->GetBetaParamNo()-1];
-        f = fTheory->Func(time, par, fFuncValues)/2.0;
-        asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0));
-        break;
-      default:
-        asymFcnValue = 0.0;
-        break;
-    }
+    f = fTheory->Func(time, par, fFuncValues)/2.0;
+    asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0));
     diff = fData.GetValue()->at(i) - asymFcnValue;
     chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
   }
@@ -455,18 +489,46 @@ void PRunAsymmetryBNMR::CalcTheory()
         asymFcnValue = fTheory->Func(time, par, fFuncValues);
         break;
       case 2: // alpha != 1, beta == 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
+        if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+          a = par[fRunInfo->GetAlphaParamNo()-1];
+        } else { // alpha is function
+          // get function number
+          UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0));
         break;
       case 3: // alpha == 1, beta != 1
-        b = par[fRunInfo->GetBetaParamNo()-1];
+        if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+          b = par[fRunInfo->GetBetaParamNo()-1];
+        } else { // beta is a function
+          // get function number
+          UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0))-f*(b+1.0)/(2.0+f*(b-1.0));
         break;
       case 4: // alpha != 1, beta != 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
-        b = par[fRunInfo->GetBetaParamNo()-1];
+        if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+          a = par[fRunInfo->GetAlphaParamNo()-1];
+        } else { // alpha is function
+          // get function number
+          UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
+        if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+          b = par[fRunInfo->GetBetaParamNo()-1];
+        } else { // beta is a function
+          // get function number
+          UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0));
         break;
@@ -477,7 +539,14 @@ void PRunAsymmetryBNMR::CalcTheory()
         break;
       case 6: // alpha ??, beta != 1
         a = alphaest;
-        b = par[fRunInfo->GetBetaParamNo()-1];
+        if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+          b = par[fRunInfo->GetBetaParamNo()-1];
+        } else { // beta is a function
+          // get function number
+          UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0));
         break;
@@ -1287,16 +1356,44 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2
       beta  = 1.0;
       break;
     case 2: // alpha != 1, beta == 1
-      alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       beta  = 1.0;
       break;
     case 3: // alpha == 1, beta != 1
       alpha = 1.0;
-      beta  = par[fRunInfo->GetBetaParamNo()-1];
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        beta = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       break;
     case 4: // alpha != 1, beta != 1
-      alpha = par[fRunInfo->GetAlphaParamNo()-1];
-      beta  = par[fRunInfo->GetBetaParamNo()-1];
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        beta = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       break;
     case 5: // alpha ?? , beta == 1
       // use estimated value
@@ -1304,7 +1401,14 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2
       break;
     case 6: // alpha ??, beta != 1
       // use estimated value
-      beta  = par[fRunInfo->GetBetaParamNo()-1];
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        beta = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       break;
     default:
       break;
diff --git a/src/classes/PRunAsymmetryRRF.cpp b/src/classes/PRunAsymmetryRRF.cpp
index 8294c28e..4c15952b 100644
--- a/src/classes/PRunAsymmetryRRF.cpp
+++ b/src/classes/PRunAsymmetryRRF.cpp
@@ -188,6 +188,58 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
   Double_t time(1.0);
   Int_t i;
 
+  // determine alpha/beta
+  switch (fAlphaBetaTag) {
+    case 1: // alpha == 1, beta == 1
+      a = 1.0;
+      b = 1.0;
+      break;
+    case 2: // alpha != 1, beta == 1
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        a = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      b = 1.0;
+      break;
+    case 3: // alpha == 1, beta != 1
+      a = 1.0;
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        b = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      break;
+    case 4: // alpha != 1, beta != 1
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        a = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        b = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      break;
+    default:
+      a = 1.0;
+      b = 1.0;
+      break;
+  }
+
   // Calculate the theory function once to ensure one function evaluation for the current set of parameters.
   // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
   // for a given set of parameters---which should be done outside of the parallelized loop.
@@ -198,34 +250,12 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
   Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
   if (chunk < 10)
     chunk = 10;
-  #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq)
+  #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,f) schedule(dynamic,chunk) reduction(+:chisq)
   #endif
   for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
     time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
-    switch (fAlphaBetaTag) {
-      case 1: // alpha == 1, beta == 1
-        asymFcnValue = fTheory->Func(time, par, fFuncValues);
-        break;
-      case 2: // alpha != 1, beta == 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
-        f = fTheory->Func(time, par, fFuncValues);
-        asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0));
-        break;
-      case 3: // alpha == 1, beta != 1
-        b = par[fRunInfo->GetBetaParamNo()-1];
-        f = fTheory->Func(time, par, fFuncValues);
-        asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0));
-        break;
-      case 4: // alpha != 1, beta != 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
-        b = par[fRunInfo->GetBetaParamNo()-1];
-        f = fTheory->Func(time, par, fFuncValues);
-        asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
-        break;
-      default:
-        asymFcnValue = 0.0;
-        break;
-    }
+    f = fTheory->Func(time, par, fFuncValues);
+    asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
     diff = fData.GetValue()->at(i) - asymFcnValue;
     chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
   }
@@ -421,18 +451,46 @@ void PRunAsymmetryRRF::CalcTheory()
         asymFcnValue = fTheory->Func(time, par, fFuncValues);
         break;
       case 2: // alpha != 1, beta == 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
+        if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+          a = par[fRunInfo->GetAlphaParamNo()-1];
+        } else { // alpha is function
+          // get function number
+          UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0));
         break;
       case 3: // alpha == 1, beta != 1
-        b = par[fRunInfo->GetBetaParamNo()-1];
+        if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+          b = par[fRunInfo->GetBetaParamNo()-1];
+        } else { // beta is a function
+          // get function number
+          UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0));
         break;
       case 4: // alpha != 1, beta != 1
-        a = par[fRunInfo->GetAlphaParamNo()-1];
-        b = par[fRunInfo->GetBetaParamNo()-1];
+        if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+          a = par[fRunInfo->GetAlphaParamNo()-1];
+        } else { // alpha is function
+          // get function number
+          UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
+        if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+          b = par[fRunInfo->GetBetaParamNo()-1];
+        } else { // beta is a function
+          // get function number
+          UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+          // evaluate function
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
         break;
@@ -1034,16 +1092,44 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]
       beta  = 1.0;
       break;
     case 2: // alpha != 1, beta == 1
-      alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       beta  = 1.0;
       break;
     case 3: // alpha == 1, beta != 1
       alpha = 1.0;
-      beta  = par[fRunInfo->GetBetaParamNo()-1];
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        beta = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       break;
     case 4: // alpha != 1, beta != 1
-      alpha = par[fRunInfo->GetAlphaParamNo()-1];
-      beta  = par[fRunInfo->GetBetaParamNo()-1];
+      if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
+        alpha = par[fRunInfo->GetAlphaParamNo()-1];
+      } else { // alpha is function
+        // get function number
+        UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
+      if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
+        beta = par[fRunInfo->GetBetaParamNo()-1];
+      } else { // beta is a function
+        // get function number
+        UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
+        // evaluate function
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+      }
       break;
     default:
       break;