diff --git a/CMakeLists.txt b/CMakeLists.txt
index 685d9fe5..26f04c9c 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.1 LANGUAGES C CXX)
+project(musrfit VERSION 1.6.2 LANGUAGES C CXX)
 
 #--- musrfit specific options -------------------------------------------------
 option(dks "build musrfit with DKS (GPU/MIC) support" ON)
diff --git a/ChangeLog b/ChangeLog
index 4ace7dca..a8b1968f 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -12,6 +12,13 @@ or
 
 https://bitbucket.org/muonspin/musrfit/commits/all
 
+Release of V1.6.2, 2020/06/17
+=============================
+
+Functions not only can use parameters and maps but now also meta data from the 
+data files like, field (B or b), energy (En or en), or temperature (T<x> pr t<x>
+where <x> is a number starting from 0, e.g. t0).
+
 Release of V1.6.1, 2020/05/25
 =============================
 
diff --git a/src/classes/PFunction.cpp b/src/classes/PFunction.cpp
index a8233589..2bdcc935 100644
--- a/src/classes/PFunction.cpp
+++ b/src/classes/PFunction.cpp
@@ -187,6 +187,16 @@ void PFunction::FillFuncEvalTree(iter_t const& i, PFuncTreeNode &node)
   } else if (i->value.id() == PFunctionGrammar::constGammaMuID) { // handle constant gamma_mu
     node.fID = PFunctionGrammar::constGammaMuID; // keep the ID
     node.fDvalue = GAMMA_BAR_MUON; // keep the value
+  } else if (i->value.id() == PFunctionGrammar::constFieldID) { // handle constant field from meta data
+    node.fID = PFunctionGrammar::constFieldID; // keep the ID
+  } else if (i->value.id() == PFunctionGrammar::constEnergyID) { // handle constant energy from meta data
+    node.fID = PFunctionGrammar::constEnergyID; // keep the ID
+  } else if (i->value.id() == PFunctionGrammar::constTempID) { // handle constant temperature from meta data
+    str = std::string(i->value.begin(), i->value.end()); // get string
+    boost::algorithm::trim(str);
+    status = sscanf(str.c_str(), "T%d", &ivalue); // convert string to temperature index
+    node.fID = PFunctionGrammar::constTempID; // keep the ID
+    node.fIvalue = ivalue; // Temp idx
   } else if (i->value.id() == PFunctionGrammar::parameterID) { // handle parameter number
     str = std::string(i->value.begin(), i->value.end()); // get string
     boost::algorithm::trim(str);
@@ -343,6 +353,12 @@ Bool_t PFunction::FindAndCheckMapAndParamRange(PFuncTreeNode &node, UInt_t mapSi
     return true;
   } else if (node.fID == PFunctionGrammar::constGammaMuID) {
     return true;
+  } else if (node.fID == PFunctionGrammar::constFieldID) {
+    return true;
+  } else if (node.fID == PFunctionGrammar::constEnergyID) {
+    return true;
+  } else if (node.fID == PFunctionGrammar::constTempID) {
+    return true;
   } else if (node.fID == PFunctionGrammar::parameterID) {
     if (node.fIvalue <= static_cast<Int_t>(paramSize))
       return true;
@@ -391,9 +407,10 @@ Bool_t PFunction::FindAndCheckMapAndParamRange(PFuncTreeNode &node, UInt_t mapSi
  *
  * \param param fit parameter vector
  */
-Double_t PFunction::Eval(std::vector<Double_t> param)
+Double_t PFunction::Eval(std::vector<Double_t> param, PMetaData metaData)
 {
   fParam = param;
+  fMetaData = metaData;
 
   return EvalNode(fFunc);
 }
@@ -414,6 +431,22 @@ Double_t PFunction::EvalNode(PFuncTreeNode &node)
     return node.fDvalue;
   } else if (node.fID == PFunctionGrammar::constGammaMuID) {
     return node.fDvalue;
+  } else if (node.fID == PFunctionGrammar::constFieldID) {
+    return fMetaData.fField;
+  } else if (node.fID == PFunctionGrammar::constEnergyID) {
+    if (fMetaData.fEnergy == PMUSR_UNDEFINED) {
+      std::cerr << std::endl << "**PANIC ERROR**: PFunction::EvalNode: energy meta data not available." << std::endl;
+      std::cerr << std::endl;
+      exit(0);
+    }
+    return fMetaData.fEnergy;
+  } else if (node.fID == PFunctionGrammar::constTempID) {
+    if (node.fIvalue >= fMetaData.fTemp.size()) {
+      std::cerr << std::endl << "**PANIC ERROR**: PFunction::EvalNode: Temp idx=" << node.fIvalue << " requested which is >= #Temp(s)=" << fMetaData.fTemp.size() << " available." << std::endl;
+      std::cerr << std::endl;
+      exit(0);
+    }
+    return fMetaData.fTemp[node.fIvalue];
   } else if (node.fID == PFunctionGrammar::parameterID) {
     Double_t dval;
     if (node.fSign)
@@ -576,6 +609,13 @@ void PFunction::EvalTreeForStringExpression(iter_t const& i)
     fFuncString += "Pi";
   } else if (i->value.id() == PFunctionGrammar::constGammaMuID) {
     fFuncString += "gamma_mu";
+  } else if (i->value.id() == PFunctionGrammar::constFieldID) {
+    fFuncString += "B";
+  } else if (i->value.id() == PFunctionGrammar::constEnergyID) {
+    fFuncString += "EN";
+  } else if (i->value.id() == PFunctionGrammar::constTempID) {
+    assert(i->children.size() == 0);
+    fFuncString += boost::algorithm::trim_copy(std::string(i->value.begin(), i->value.end())).c_str();
   } else if (i->value.id() == PFunctionGrammar::funLabelID) {
     assert(i->children.size() == 0);
     //SetFuncNo(i);
diff --git a/src/classes/PFunctionHandler.cpp b/src/classes/PFunctionHandler.cpp
index 2d4ede62..fb2a30d6 100644
--- a/src/classes/PFunctionHandler.cpp
+++ b/src/classes/PFunctionHandler.cpp
@@ -167,7 +167,7 @@ Bool_t PFunctionHandler::CheckMapAndParamRange(UInt_t mapSize, UInt_t paramSize)
  * \param map map vector
  * \param param fit parameter vector
  */
-Double_t PFunctionHandler::Eval(Int_t funNo, std::vector<Int_t> map, std::vector<double> param)
+Double_t PFunctionHandler::Eval(Int_t funNo, std::vector<Int_t> map, std::vector<double> param, PMetaData metaData)
 {
   if (GetFuncIndex(funNo) == -1) {
     std::cerr << std::endl << "**ERROR**: Couldn't find FUN" << funNo << " for evaluation";
@@ -179,7 +179,7 @@ Double_t PFunctionHandler::Eval(Int_t funNo, std::vector<Int_t> map, std::vector
   fFuncs[GetFuncIndex(funNo)].SetMap(map);
 
   // return evaluated function
-  return fFuncs[GetFuncIndex(funNo)].Eval(param);
+  return fFuncs[GetFuncIndex(funNo)].Eval(param, metaData);
 }
 
 //-------------------------------------------------------------
diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp
index c21577d1..a5626121 100644
--- a/src/classes/PRunAsymmetry.cpp
+++ b/src/classes/PRunAsymmetry.cpp
@@ -199,7 +199,7 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate chi square
@@ -219,7 +219,7 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       b = 1.0;
       break;
@@ -231,7 +231,7 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     case 4: // alpha != 1, beta != 1
@@ -241,7 +241,7 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
         b = par[fRunInfo->GetBetaParamNo()-1];
@@ -249,7 +249,7 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     default:
@@ -455,7 +455,7 @@ void PRunAsymmetry::CalcTheory()
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate asymmetry
@@ -475,7 +475,7 @@ void PRunAsymmetry::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0));
@@ -487,7 +487,7 @@ void PRunAsymmetry::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0));
@@ -499,7 +499,7 @@ void PRunAsymmetry::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
           b = par[fRunInfo->GetBetaParamNo()-1];
@@ -507,7 +507,7 @@ void PRunAsymmetry::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
@@ -560,11 +560,14 @@ Bool_t PRunAsymmetry::PrepareData()
   }
 
   // keep the field from the meta-data from the data-file
-  fField = runData->GetField();
+  fMetaData.fField = runData->GetField();
+
+  // keep the energy from the meta-data from the data-file
+  fMetaData.fEnergy = runData->GetEnergy();
 
   // keep the temperature(s) from the meta-data from the data-file
   for (unsigned int i=0; i<runData->GetNoOfTemperatures(); i++)
-    fTemp.push_back(runData->GetTemperature(i));
+    fMetaData.fTemp.push_back(runData->GetTemperature(i));
 
   // collect histogram numbers
   PUIntVector forwardHistoNo;
@@ -1219,7 +1222,7 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2])
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       beta  = 1.0;
       break;
@@ -1231,7 +1234,7 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2])
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     case 4: // alpha != 1, beta != 1
@@ -1241,7 +1244,7 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2])
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
         beta = par[fRunInfo->GetBetaParamNo()-1];
@@ -1249,7 +1252,7 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2])
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     default:
@@ -1287,7 +1290,7 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2])
   // fill theory vector for kView
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate theory
@@ -1458,7 +1461,7 @@ Bool_t PRunAsymmetry::PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       beta  = 1.0;
       break;
@@ -1470,7 +1473,7 @@ Bool_t PRunAsymmetry::PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     case 4: // alpha != 1, beta != 1
@@ -1480,7 +1483,7 @@ Bool_t PRunAsymmetry::PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
         beta = par[fRunInfo->GetBetaParamNo()-1];
@@ -1488,7 +1491,7 @@ Bool_t PRunAsymmetry::PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     default:
@@ -1591,7 +1594,7 @@ Bool_t PRunAsymmetry::PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   Double_t theoryValue;
diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp
index bab19fe5..7f6b6f75 100644
--- a/src/classes/PRunAsymmetryBNMR.cpp
+++ b/src/classes/PRunAsymmetryBNMR.cpp
@@ -199,7 +199,7 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector<Double_t>& par)
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate chi square
@@ -220,7 +220,7 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       b = 1.0;
       break;
@@ -232,7 +232,7 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     case 4: // alpha != 1, beta != 1
@@ -242,7 +242,7 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
         b = par[fRunInfo->GetBetaParamNo()-1];
@@ -250,7 +250,7 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     case 5: // alpha ?? , beta == 1
@@ -265,7 +265,7 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     default:
@@ -471,7 +471,7 @@ void PRunAsymmetryBNMR::CalcTheory()
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate asymmetry
@@ -495,7 +495,7 @@ void PRunAsymmetryBNMR::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         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));
@@ -507,7 +507,7 @@ void PRunAsymmetryBNMR::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         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));
@@ -519,7 +519,7 @@ void PRunAsymmetryBNMR::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
           b = par[fRunInfo->GetBetaParamNo()-1];
@@ -527,7 +527,7 @@ void PRunAsymmetryBNMR::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         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));
@@ -545,7 +545,7 @@ void PRunAsymmetryBNMR::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         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));
@@ -598,11 +598,14 @@ Bool_t PRunAsymmetryBNMR::PrepareData()
   }
 
   // keep the field from the meta-data from the data-file
-  fField = runData->GetField();
+  fMetaData.fField = runData->GetField();
+
+  // keep the energy from the meta-data from the data-file
+  fMetaData.fEnergy = runData->GetEnergy();
 
   // keep the temperature(s) from the meta-data from the data-file
   for (unsigned int i=0; i<runData->GetNoOfTemperatures(); i++)
-    fTemp.push_back(runData->GetTemperature(i));
+    fMetaData.fTemp.push_back(runData->GetTemperature(i));
 
   // collect histogram numbers
   PUIntVector forwardHistoNo;
@@ -1369,7 +1372,7 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       beta  = 1.0;
       break;
@@ -1381,7 +1384,7 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     case 4: // alpha != 1, beta != 1
@@ -1391,7 +1394,7 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
         beta = par[fRunInfo->GetBetaParamNo()-1];
@@ -1399,7 +1402,7 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     case 5: // alpha ?? , beta == 1
@@ -1414,7 +1417,7 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     default:
@@ -1465,7 +1468,7 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2
   // fill theory vector for kView
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate theory
diff --git a/src/classes/PRunAsymmetryRRF.cpp b/src/classes/PRunAsymmetryRRF.cpp
index ace43f52..4f366ed7 100644
--- a/src/classes/PRunAsymmetryRRF.cpp
+++ b/src/classes/PRunAsymmetryRRF.cpp
@@ -181,7 +181,7 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate chi square
@@ -201,7 +201,7 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       b = 1.0;
       break;
@@ -213,7 +213,7 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     case 4: // alpha != 1, beta != 1
@@ -223,7 +223,7 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
         b = par[fRunInfo->GetBetaParamNo()-1];
@@ -231,7 +231,7 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     default:
@@ -437,7 +437,7 @@ void PRunAsymmetryRRF::CalcTheory()
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate asymmetry
@@ -457,7 +457,7 @@ void PRunAsymmetryRRF::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0));
@@ -469,7 +469,7 @@ void PRunAsymmetryRRF::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0));
@@ -481,7 +481,7 @@ void PRunAsymmetryRRF::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
           b = par[fRunInfo->GetBetaParamNo()-1];
@@ -489,7 +489,7 @@ void PRunAsymmetryRRF::CalcTheory()
           // get function number
           UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
           // evaluate function
-          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+          b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
         }
         f = fTheory->Func(time, par, fFuncValues);
         asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
@@ -542,11 +542,14 @@ Bool_t PRunAsymmetryRRF::PrepareData()
   }
 
   // keep the field from the meta-data from the data-file
-  fField = runData->GetField();
+  fMetaData.fField = runData->GetField();
+
+  // keep the energy from the meta-data from the data-file
+  fMetaData.fEnergy = runData->GetEnergy();
 
   // keep the temperature(s) from the meta-data from the data-file
   for (unsigned int i=0; i<runData->GetNoOfTemperatures(); i++)
-    fTemp.push_back(runData->GetTemperature(i));
+    fMetaData.fTemp.push_back(runData->GetTemperature(i));
 
   // collect histogram numbers
   PUIntVector forwardHistoNo;
@@ -1105,7 +1108,7 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       beta  = 1.0;
       break;
@@ -1117,7 +1120,7 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     case 4: // alpha != 1, beta != 1
@@ -1127,7 +1130,7 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]
         // get function number
         UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
         beta = par[fRunInfo->GetBetaParamNo()-1];
@@ -1135,7 +1138,7 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]
         // get function number
         UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
         // evaluate function
-        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+        beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
       }
       break;
     default:
@@ -1220,7 +1223,7 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]
   // fill theory vector for kView
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate theory
diff --git a/src/classes/PRunBase.cpp b/src/classes/PRunBase.cpp
index a55707d2..4951d563 100644
--- a/src/classes/PRunBase.cpp
+++ b/src/classes/PRunBase.cpp
@@ -49,7 +49,8 @@ PRunBase::PRunBase()
   fRunNo = -1;
   fRunInfo = nullptr;
   fRawData = nullptr;
-  fField = -1.0;
+  fMetaData.fField  = PMUSR_UNDEFINED;
+  fMetaData.fEnergy = PMUSR_UNDEFINED;
   fTimeResolution = -1.0;
 
   fFitStartTime = PMUSR_UNDEFINED;
@@ -91,7 +92,8 @@ PRunBase::PRunBase(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo,
   }
 
   // init private variables
-  fField = -1.0;
+  fMetaData.fField  = PMUSR_UNDEFINED;
+  fMetaData.fEnergy = PMUSR_UNDEFINED;
   fTimeResolution = -1.0;
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++)
     fFuncValues.push_back(0.0);
diff --git a/src/classes/PRunMuMinus.cpp b/src/classes/PRunMuMinus.cpp
index ec923faf..aff00181 100644
--- a/src/classes/PRunMuMinus.cpp
+++ b/src/classes/PRunMuMinus.cpp
@@ -138,7 +138,7 @@ Double_t PRunMuMinus::CalcChiSquare(const std::vector<Double_t>& par)
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
     Int_t funcNo = fMsrInfo->GetFuncNo(i);
-    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate chi square
@@ -186,7 +186,7 @@ Double_t PRunMuMinus::CalcChiSquareExpected(const std::vector<Double_t>& par)
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
     Int_t funcNo = fMsrInfo->GetFuncNo(i);
-    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate chi square
@@ -233,7 +233,7 @@ Double_t PRunMuMinus::CalcMaxLikelihood(const std::vector<Double_t>& par)
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
     Int_t funcNo = fMsrInfo->GetFuncNo(i);
-    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate maximum log likelihood
@@ -419,7 +419,7 @@ void PRunMuMinus::CalcTheory()
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate theory
@@ -467,6 +467,16 @@ Bool_t PRunMuMinus::PrepareData()
     return false;
   }
 
+  // keep the field from the meta-data from the data-file
+  fMetaData.fField = runData->GetField();
+
+  // keep the energy from the meta-data from the data-file
+  fMetaData.fEnergy = runData->GetEnergy();
+
+  // keep the temperature(s) from the meta-data from the data-file
+  for (unsigned int i=0; i<runData->GetNoOfTemperatures(); i++)
+    fMetaData.fTemp.push_back(runData->GetTemperature(i));
+
   // collect histogram numbers
   PUIntVector histoNo; // histoNo = msr-file forward + redGreen_offset - 1
   for (UInt_t i=0; i<fRunInfo->GetForwardHistoNoSize(); i++) {
@@ -727,7 +737,7 @@ Bool_t PRunMuMinus::PrepareRawViewData(PRawRunData* runData, const UInt_t histoN
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate theory
diff --git a/src/classes/PRunNonMusr.cpp b/src/classes/PRunNonMusr.cpp
index f5d35851..0976d241 100644
--- a/src/classes/PRunNonMusr.cpp
+++ b/src/classes/PRunNonMusr.cpp
@@ -103,7 +103,7 @@ Double_t PRunNonMusr::CalcChiSquare(const std::vector<Double_t>& par)
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate chi square
@@ -345,7 +345,7 @@ Bool_t PRunNonMusr::PrepareViewData()
     par.push_back((*paramList)[i].fValue);
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // get plot range
diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp
index 17b852b3..6d1914a2 100644
--- a/src/classes/PRunSingleHisto.cpp
+++ b/src/classes/PRunSingleHisto.cpp
@@ -149,7 +149,7 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector<Double_t>& par)
     // get function number
     UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
     // evaluate function
-    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // get tau
@@ -174,7 +174,7 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector<Double_t>& par)
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
     UInt_t funcNo = fMsrInfo->GetFuncNo(i);
-    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate chi square
@@ -234,7 +234,7 @@ Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector<Double_t>& par
     // get function number
     UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
     // evaluate function
-    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // get tau
@@ -259,7 +259,7 @@ Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector<Double_t>& par
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
     Int_t funcNo = fMsrInfo->GetFuncNo(i);
-    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate chi square
@@ -317,7 +317,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector<Double_t>& par)
     // get function number
     Int_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
     // evaluate function
-    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // get tau
@@ -342,7 +342,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector<Double_t>& par)
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
     UInt_t funcNo = fMsrInfo->GetFuncNo(i);
-    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate maximum log likelihood
@@ -415,7 +415,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector<Double_t>&
     // get function number
     Int_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
     // evaluate function
-    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // get tau
@@ -440,7 +440,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector<Double_t>&
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
     UInt_t funcNo = fMsrInfo->GetFuncNo(i);
-    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate maximum log likelihood
@@ -510,7 +510,7 @@ void PRunSingleHisto::CalcTheory()
     // get function number
     Int_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
     // evaluate function
-    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // get tau
@@ -534,7 +534,7 @@ void PRunSingleHisto::CalcTheory()
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate theory
@@ -710,11 +710,14 @@ Bool_t PRunSingleHisto::PrepareData()
   }
 
   // keep the field from the meta-data from the data-file
-  fField = runData->GetField();
+  fMetaData.fField = runData->GetField();
+
+  // keep the energy from the meta-data from the data-file
+  fMetaData.fEnergy = runData->GetEnergy();
 
   // keep the temperature(s) from the meta-data from the data-file
   for (unsigned int i=0; i<runData->GetNoOfTemperatures(); i++)
-    fTemp.push_back(runData->GetTemperature(i));
+    fMetaData.fTemp.push_back(runData->GetTemperature(i));
 
   // collect histogram numbers
   PUIntVector histoNo; // histoNo = msr-file forward + redGreen_offset - 1
@@ -1023,7 +1026,7 @@ Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t hi
     // get function number
     UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
     // evaluate function
-    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
   }
   N0 *= theoryNorm;
 
@@ -1062,7 +1065,7 @@ Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t hi
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate theory
@@ -1197,7 +1200,7 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo
     // get function number
     UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
     // evaluate function
-    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
+    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
   }
   N0 *= theoryNorm;
 
@@ -1302,7 +1305,7 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate theory
diff --git a/src/classes/PRunSingleHistoRRF.cpp b/src/classes/PRunSingleHistoRRF.cpp
index ebf3cdc0..76e9374a 100644
--- a/src/classes/PRunSingleHistoRRF.cpp
+++ b/src/classes/PRunSingleHistoRRF.cpp
@@ -154,7 +154,7 @@ Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector<Double_t>& par)
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
     UInt_t funcNo = fMsrInfo->GetFuncNo(i);
-    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate chi square
@@ -202,7 +202,7 @@ Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector<Double_t>&
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
     UInt_t funcNo = fMsrInfo->GetFuncNo(i);
-    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate chi square
@@ -265,7 +265,7 @@ void PRunSingleHistoRRF::CalcTheory()
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // calculate theory
@@ -441,11 +441,14 @@ Bool_t PRunSingleHistoRRF::PrepareData()
   }
 
   // keep the field from the meta-data from the data-file
-  fField = runData->GetField();
+  fMetaData.fField = runData->GetField();
+
+  // keep the energy from the meta-data from the data-file
+  fMetaData.fEnergy = runData->GetEnergy();
 
   // keep the temperature(s) from the meta-data from the data-file
   for (unsigned int i=0; i<runData->GetNoOfTemperatures(); i++)
-    fTemp.push_back(runData->GetTemperature(i));
+    fMetaData.fTemp.push_back(runData->GetTemperature(i));
 
   // collect histogram numbers
   PUIntVector histoNo; // histoNo = msr-file forward + redGreen_offset - 1
@@ -750,7 +753,7 @@ Bool_t PRunSingleHistoRRF::PrepareViewData(PRawRunData* runData, const UInt_t hi
 
   // calculate functions
   for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
-    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par);
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
   }
 
   // check if a finer binning for the theory is needed
diff --git a/src/include/PFunction.h b/src/include/PFunction.h
index daaa2c2b..09013f98 100644
--- a/src/include/PFunction.h
+++ b/src/include/PFunction.h
@@ -44,6 +44,7 @@
 
 #include <TString.h>
 
+#include "PMusr.h"
 #include "PFunctionGrammar.h"
 
 //----------------------------------------------------------------------------
@@ -79,7 +80,7 @@ typedef struct func_tree_node {
   Int_t    fOperatorTag; ///< tag for '+', '-', '*', '/'
   Int_t    fFunctionTag; ///< tag got "cos", "sin", ...
   Int_t    fIvalue; ///< for parameter numbers and maps
-  Bool_t fSign; ///< for sign, true means '-', false '+'
+  Bool_t   fSign; ///< for sign, true means '-', false '+'
   Double_t fDvalue; ///< for numbers
   std::vector<func_tree_node> children; ///< holding sub-tree
 } PFuncTreeNode;
@@ -96,7 +97,7 @@ class PFunction {
     virtual Bool_t IsValid() { return fValid; }
     virtual Int_t GetFuncNo() { return fFuncNo; }
     virtual Bool_t CheckMapAndParamRange(UInt_t mapSize, UInt_t paramSize);
-    virtual Double_t Eval(std::vector<Double_t> param);
+    virtual Double_t Eval(std::vector<Double_t> param, PMetaData metaData);
     virtual void SetMap(std::vector<Int_t> map) { fMap = map; }
 
     virtual TString* GetFuncString() { return &fFuncString; }
@@ -124,6 +125,8 @@ class PFunction {
     virtual void EvalTreeForString(tree_parse_info<> info);
     virtual void EvalTreeForStringExpression(iter_t const& i);
     TString fFuncString; ///< clear text representation of the function
+
+    PMetaData fMetaData; ///< keeps meta data from data files (field, energy, temperature, ...)
 };
 
 #endif // _PFUNCTION_H_
diff --git a/src/include/PFunctionGrammar.h b/src/include/PFunctionGrammar.h
index 95d235a4..e1ed1e1c 100644
--- a/src/include/PFunctionGrammar.h
+++ b/src/include/PFunctionGrammar.h
@@ -57,15 +57,18 @@ struct PFunctionGrammar : public grammar<PFunctionGrammar>
     static const int realID         = 1;
     static const int constPiID      = 2;
     static const int constGammaMuID = 3;
-    static const int funLabelID     = 4;
-    static const int parameterID    = 5;
-    static const int mapID          = 6;
-    static const int functionID     = 7;
-    static const int powerID        = 8;
-    static const int factorID       = 9;
-    static const int termID         = 10;
-    static const int expressionID   = 11;
-    static const int assignmentID   = 12;
+    static const int constFieldID   = 4;
+    static const int constEnergyID  = 5;
+    static const int constTempID    = 6;
+    static const int funLabelID     = 7;
+    static const int parameterID    = 8;
+    static const int mapID          = 9;
+    static const int functionID     = 10;
+    static const int powerID        = 11;
+    static const int factorID       = 12;
+    static const int termID         = 13;
+    static const int expressionID   = 14;
+    static const int assignmentID   = 15;
 
     template <typename ScannerT>
     struct definition
@@ -79,6 +82,12 @@ struct PFunctionGrammar : public grammar<PFunctionGrammar>
 
             const_gamma_mu =   leaf_node_d[ str_p("GAMMA_MU") ];
 
+            const_field    =   leaf_node_d[ str_p("B") ];
+
+            const_energy   =   leaf_node_d[ str_p("EN") ];
+
+            const_temp    =    leaf_node_d[ ( lexeme_d[ "T" >> +digit_p ] ) ];
+
             fun_label      =   leaf_node_d[ ( lexeme_d[ "FUN" >> +digit_p ] ) ];
 
             parameter      =   leaf_node_d[ ( lexeme_d[ "PAR" >> +digit_p ] ) |
@@ -110,6 +119,9 @@ struct PFunctionGrammar : public grammar<PFunctionGrammar>
             factor         =   real
                            |   const_pi
                            |   const_gamma_mu
+                           |   const_field
+                           |   const_energy
+                           |   const_temp
                            |   parameter
                            |   map
                            |   function
@@ -134,6 +146,9 @@ struct PFunctionGrammar : public grammar<PFunctionGrammar>
             BOOST_SPIRIT_DEBUG_RULE(real);
             BOOST_SPIRIT_DEBUG_RULE(const_pi);
             BOOST_SPIRIT_DEBUG_RULE(const_gamma_mu);
+            BOOST_SPIRIT_DEBUG_RULE(const_field);
+            BOOST_SPIRIT_DEBUG_RULE(const_energy);
+            BOOST_SPIRIT_DEBUG_RULE(const_temp);
             BOOST_SPIRIT_DEBUG_RULE(fun_label);
             BOOST_SPIRIT_DEBUG_RULE(parameter);
             BOOST_SPIRIT_DEBUG_RULE(map);
@@ -154,6 +169,9 @@ struct PFunctionGrammar : public grammar<PFunctionGrammar>
         rule<ScannerT, parser_context<>, parser_tag<mapID> >          map;
         rule<ScannerT, parser_context<>, parser_tag<parameterID> >    parameter;
         rule<ScannerT, parser_context<>, parser_tag<funLabelID> >     fun_label;
+        rule<ScannerT, parser_context<>, parser_tag<constTempID> >    const_temp;
+        rule<ScannerT, parser_context<>, parser_tag<constEnergyID> >  const_energy;
+        rule<ScannerT, parser_context<>, parser_tag<constFieldID> >   const_field;
         rule<ScannerT, parser_context<>, parser_tag<constGammaMuID> > const_gamma_mu;
         rule<ScannerT, parser_context<>, parser_tag<constPiID> >      const_pi;
         rule<ScannerT, parser_context<>, parser_tag<realID> >         real;
diff --git a/src/include/PFunctionHandler.h b/src/include/PFunctionHandler.h
index 0bd319a7..2caabdfc 100644
--- a/src/include/PFunctionHandler.h
+++ b/src/include/PFunctionHandler.h
@@ -51,7 +51,7 @@ class PFunctionHandler
     virtual Bool_t IsValid() { return fValid; }
     virtual Bool_t DoParse();
     virtual Bool_t CheckMapAndParamRange(UInt_t mapSize, UInt_t paramSize);
-    virtual double Eval(Int_t funNo, std::vector<Int_t> map, std::vector<double> param);
+    virtual double Eval(Int_t funNo, std::vector<Int_t> map, std::vector<double> param, PMetaData metaData);
     virtual Int_t GetFuncNo(UInt_t idx);
     virtual Int_t GetFuncIndex(Int_t funcNo);
     virtual UInt_t GetNoOfFuncs() { return fFuncs.size(); }
diff --git a/src/include/PMsrHandler.h b/src/include/PMsrHandler.h
index 87b82ec4..d8c20472 100644
--- a/src/include/PMsrHandler.h
+++ b/src/include/PMsrHandler.h
@@ -95,9 +95,8 @@ class PMsrHandler
     virtual UInt_t GetFuncIndex(Int_t funNo) { return fFuncHandler->GetFuncIndex(funNo); }
     virtual Bool_t CheckMapAndParamRange(UInt_t mapSize, UInt_t paramSize)
                        { return fFuncHandler->CheckMapAndParamRange(mapSize, paramSize); }
-    virtual Double_t EvalFunc(UInt_t i, std::vector<Int_t> map, std::vector<Double_t> param)
-                       { return fFuncHandler->Eval(i,map,param); }
-
+    virtual Double_t EvalFunc(UInt_t i, std::vector<Int_t> map, std::vector<Double_t> param, PMetaData metaData)
+                       { return fFuncHandler->Eval(i, map, param, metaData); }
     virtual UInt_t GetNoOfFitParameters(UInt_t idx);
     virtual Int_t ParameterInUse(UInt_t paramNo);
     virtual Bool_t CheckRunBlockIntegrity();
diff --git a/src/include/PMusr.h b/src/include/PMusr.h
index b36b538e..250c46f0 100644
--- a/src/include/PMusr.h
+++ b/src/include/PMusr.h
@@ -219,6 +219,16 @@ typedef std::vector<TString> PStringVector;
  */
 enum EPMusrHandleTag { kEmpty, kFit, kView };
 
+//-------------------------------------------------------------
+/**
+ * <p>keep a couple of meta data for the FUNCTIONS.
+ */
+struct PMetaData {
+  Double_t fField;  ///< field in (G)
+  Double_t fEnergy; ///< energy in (keV)
+  PDoubleVector fTemp; ///< temperature(s) in (K)
+};
+
 //-------------------------------------------------------------
 /**
  * <p>Handles the data which will be fitted, i.e. packed, background corrected, ...
diff --git a/src/include/PRunBase.h b/src/include/PRunBase.h
index f3ac90a9..febfe6a7 100644
--- a/src/include/PRunBase.h
+++ b/src/include/PRunBase.h
@@ -74,8 +74,7 @@ class PRunBase
 
     PRunData fData;             ///< data to be fitted, viewed, i.e. binned data
     Double_t fTimeResolution;   ///< time resolution in (us)
-    Double_t fField;            ///< field from the meta-data of the data file
-    PDoubleVector fTemp;        ///< temperature(s) from the meta-data of the data file : first = temp, second = std. dev. of the temperature
+    PMetaData fMetaData;        ///< keeps the meta data from the data file like field, temperature, energy, ...
     PDoubleVector fT0s;         ///< all t0 bins of a run! The derived classes will handle it.
     std::vector<PDoubleVector> fAddT0s; ///< all t0 bins of all addrun's of a run! The derived classes will handle it.