PTheory.cpp

Go to the documentation of this file.

/***************************************************************************
 
  PTheory.cpp
 
  Author: Andreas Suter
  e-mail: andreas.suter@psi.ch
 
***************************************************************************/
 
/***************************************************************************
 *   Copyright (C) 2007-2026 by Andreas Suter                              *
 *   andreas.suter@psi.ch                                                  *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
 
#include <iostream>
#include <vector>
#include <string>
#include <cmath>
 
#include <TObject.h>
#include <TString.h>
#include <TF1.h>
#include <TClass.h>
#include <TMath.h>
 
#include <Math/SpecFuncMathMore.h>
 
#include "PMsrHandler.h"
#include "PStringUtils.h"
#include "PTheory.h"
 
#define SQRT_TWO 1.41421356237
#define SQRT_PI  1.77245385091
 
extern std::vector<void*> gGlobalUserFcn;
 
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
PTheory::PTheory(PMsrHandler *msrInfo, UInt_t runNo, const Bool_t hasParent) : fMsrInfo(msrInfo)
{
  // init stuff
  fValid = true;
  fAdd = nullptr;
  fMul = nullptr;
  fUserFcnClassName = TString("");
  fUserFcn = nullptr;
  fDynLFdt = 0.0;
  fSamplingTime = 0.001; // default = 1ns (units in us)
 
  static UInt_t lineNo = 1; // lineNo
  static UInt_t depth  = 0; // needed to handle '+' properly
 
  if (hasParent == false) { // reset static counters if root object
    lineNo = 1; // the lineNo counter and the depth counter need to be
    depth  = 0; // reset for every root object (new run).
  }
 
  for (UInt_t i=0; i<THEORY_MAX_PARAM; i++)
    fPrevParam[i] = 0.0;
 
  // keep the number of user functions found up to this point
  fUserFcnIdx = GetUserFcnIdx(lineNo);
 
  // get the input to be analyzed from the msr handler
  PMsrLines *fullTheoryBlock = msrInfo->GetMsrTheory();
  if (lineNo > fullTheoryBlock->size()-1) {
    return;
  }
  // get the line to be parsed
  PMsrLineStructure *line = &(*fullTheoryBlock)[lineNo];
 
  // copy line content to str in order to remove comments
  TString str = line->fLine.Copy();
 
  // remove theory line comment if present, i.e. something starting with '('
  Int_t index = str.Index("(");
  if (index > 0) // theory line comment present
    str.Resize(index);
 
  // remove msr-file comment if present, i.e. something starting with '#'
  index = str.Index("#");
  if (index > 0) // theory line comment present
    str.Resize(index);
 
  // tokenize line
  std::vector<std::string> tokens = PStringUtils::Split(str.Data(), " \t");
  if (tokens.empty()) {
    std::cerr << std::endl << ">> PTheory::PTheory: **SEVERE ERROR** Couldn't tokenize theory block line " << line->fLineNo << ".";
    std::cerr << std::endl << ">>  line content: " << line->fLine.Data();
    std::cerr << std::endl;
    exit(0);
  }
  str = tokens[0];
 
  // search the theory function
  UInt_t idx = SearchDataBase(str);
 
  // function found is not defined
  if (idx == static_cast<UInt_t>(THEORY_UNDEFINED)) {
    std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** Theory line '" << line->fLine.Data() << "'";
    std::cerr << std::endl << ">>  in line no " << line->fLineNo << " is undefined!";
    std::cerr << std::endl;
    fValid = false;
    return;
  }
 
  // line is a valid function, hence analyze parameters
  if ((static_cast<UInt_t>(tokens.size()-1) < fNoOfParam) &&
      ((idx != THEORY_USER_FCN) && (idx != THEORY_POLYNOM))) {
    std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** Theory line '" << line->fLine.Data() << "'";
    std::cerr << std::endl << ">>  in line no " << line->fLineNo;
    std::cerr << std::endl << ">>  expecting " << fgTheoDataBase[idx].fNoOfParam << ", but found " << tokens.size()-1;
    std::cerr << std::endl;
    fValid = false;
  }
  // keep function index
  fType = idx;
  // filter out the parameters
  Int_t status;
  UInt_t value;
  Bool_t ok = false;
  for (UInt_t i=1; i<tokens.size(); i++) {
    str = tokens[i];
 
    // if userFcn, the first entry is the function name and needs to be handled specially
    if ((fType == THEORY_USER_FCN) && ((i == 1) || (i == 2))) {
      if (i == 1) {
        fUserFcnSharedLibName = str;
      }
      if (i == 2) {
        fUserFcnClassName = str;
      }
      continue;
    }
 
    // check if str is map
    if (str.Contains("map")) {
      status = sscanf(str.Data(), "map%u", &value);
      if (status == 1) { // everthing ok
        ok = true;
        // get parameter from map
        PIntVector maps = *(*msrInfo->GetMsrRunList())[runNo].GetMap();
        if ((value <= maps.size()) && (value > 0)) { // everything fine
          fParamNo.push_back(maps[value-1]-1);
        } else { // map index out of range
          std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** map index " << value << " out of range! See line no " << line->fLineNo;
          std::cerr << std::endl;
          fValid = false;
        }
      } else { // something wrong
        std::cerr << std::endl << ">> PTheory::PTheory: **ERROR**: map '" << str.Data() << "' not allowed. See line no " << line->fLineNo;
        std::cerr << std::endl;
        fValid = false;
      }
    } else if (str.Contains("fun")) { // check if str is fun
      status = sscanf(str.Data(), "fun%u", &value);
      if (status == 1) { // everthing ok
        ok = true;
        // handle function, i.e. get, from the function number x (FUNx), the function index, 
        // add function offset and fill "parameter vector" 
        fParamNo.push_back(msrInfo->GetFuncIndex(value)+MSR_PARAM_FUN_OFFSET);
      } else { // something wrong
        fValid = false;
      }
    } else { // check if str is param no
      status = sscanf(str.Data(), "%u", &value);
      if (status == 1) { // everthing ok
        ok = true;
        fParamNo.push_back(value-1);
      }
      // check if one of the valid entries was found
      if (!ok) {
        std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** '" << str.Data() << "' not allowed. See line no " << line->fLineNo;
        std::cerr << std::endl;
        fValid = false;
      }
    }
  }
 
  // call the next line (only if valid so far and not the last line)
  // handle '*'
  if (fValid && (lineNo < fullTheoryBlock->size()-1)) {
    line = &(*fullTheoryBlock)[lineNo+1];
    if (!line->fLine.Contains("+")) { // make sure next line is not a '+'
      depth++;
      lineNo++;
      fMul = new PTheory(msrInfo, runNo, true);
      depth--;
    }
  }
  // call the next line (only if valid so far and not the last line)
  // handle '+'
  if (fValid && (lineNo < fullTheoryBlock->size()-1)) {
    line = &(*fullTheoryBlock)[lineNo+1];
    if ((depth == 0) && line->fLine.Contains("+")) {
      lineNo += 2; // go to the next theory function line
      fAdd = new PTheory(msrInfo, runNo, true);
    }
  }
 
  // make clean and tidy theory block for the msr-file
  if (fValid && !hasParent) { // parent theory object
    MakeCleanAndTidyTheoryBlock(fullTheoryBlock);
  }
 
  // check if user function, if so, check if it is reachable (root) and if yes invoke object
  if (!fUserFcnClassName.IsWhitespace()) {
    std::cout << std::endl << ">> user function class name: " << fUserFcnClassName.Data() << std::endl;
    if (!TClass::GetDict(fUserFcnClassName.Data())) {
      if (gSystem->Load(fUserFcnSharedLibName.Data()) < 0) {
        std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** user function class '" << fUserFcnClassName.Data() << "' not found.";
        std::cerr << std::endl << ">>           Tried to load " << fUserFcnSharedLibName.Data() << " but failed.";
        std::cerr << std::endl << ">>           See line no " << line->fLineNo;
        std::cerr << std::endl;
        fValid = false;
        return;
      } else if (!TClass::GetDict(fUserFcnClassName.Data())) {
        std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** user function class '" << fUserFcnClassName.Data() << "' not found.";
        std::cerr << std::endl << ">>           " << fUserFcnSharedLibName.Data() << " loaded successfully, but no dictionary present.";
        std::cerr << std::endl << ">>           See line no " << line->fLineNo;
        std::cerr << std::endl;
        fValid = false;
        return;
      }
    }
 
    // invoke user function object
    fUserFcn = nullptr;
    fUserFcn = static_cast<PUserFcnBase*>(TClass::GetClass(fUserFcnClassName.Data())->New());
    if (fUserFcn == nullptr) {
      std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** user function object could not be invoked. See line no " << line->fLineNo;
      std::cerr << std::endl;
      fValid = false;
      return;
    } else { // user function valid, hence expand the fUserParam vector to the proper size
      fUserParam.resize(fParamNo.size());
    }
 
    // check if the global part of the user function is needed
    if (fUserFcn->NeedGlobalPart()) {
      fUserFcn->SetGlobalPart(gGlobalUserFcn, fUserFcnIdx); // invoke or retrieve global user function object
      if (!fUserFcn->GlobalPartIsValid()) {
        std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** global user function object could not be invoked/retrived. See line no " << line->fLineNo;
        std::cerr << std::endl;
        fValid = false;
      }
    }
  }
}
 
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
PTheory::~PTheory()
{
  fParamNo.clear();
  fUserParam.clear();
 
  fLFIntegral.clear();
  fDynLFFuncValue.clear();
 
  // recursive clean up
  CleanUp(this);
 
  if (fUserFcn) {
    delete fUserFcn;
    fUserFcn = nullptr;
  }
 
  gGlobalUserFcn.clear();
}
 
//--------------------------------------------------------------------------
// IsValid
//--------------------------------------------------------------------------
Bool_t PTheory::IsValid()
{
 
  if (fMul) {
    if (fAdd) {
      return (fValid && fMul->IsValid() && fAdd->IsValid());
    } else {
      return (fValid && fMul->IsValid());
    }
  } else {
    if (fAdd) {
      return (fValid && fAdd->IsValid());
    } else {
      return fValid;
    }
  }
}
 
//--------------------------------------------------------------------------
Double_t PTheory::Func(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  if (fMul) {
    if (fAdd) { // fMul != 0 && fAdd != 0
      switch (fType) {
        case THEORY_CONST:
          return Constant(paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
              fAdd->Func(t, paramValues, funcValues);
        case THEORY_ASYMMETRY:
          return Asymmetry(paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_SIMPLE_EXP:
          return SimpleExp(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_GENERAL_EXP:
          return GeneralExp(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_SIMPLE_GAUSS:
          return SimpleGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_STATIC_GAUSS_KT:
          return StaticGaussKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_STATIC_GAUSS_KT_LF:
          return StaticGaussKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAUSS_KT_LF:
          return DynamicGaussKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_STATIC_LORENTZ_KT:
          return StaticLorentzKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_STATIC_LORENTZ_KT_LF:
          return StaticLorentzKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_LORENTZ_KT_LF:
          return DynamicLorentzKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_FAST_KT_ZF:
          return DynamicGauLorKTZFFast(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_FAST_KT_LF:
          return DynamicGauLorKTLFFast(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_KT_LF:
          return DynamicGauLorKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_COMBI_LGKT:
          return CombiLGKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_STR_KT:
          return StrKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_SPIN_GLASS:
          return SpinGlass(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_RANDOM_ANISOTROPIC_HYPERFINE:
          return RandomAnisotropicHyperfine(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_ABRAGAM:
          return Abragam(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) + 
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_TF_COS:
          return TFCos(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD:
          return InternalField(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD_KORNILOV:
          return InternalFieldGK(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
              fAdd->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD_LARKIN:
          return InternalFieldLL(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
              fAdd->Func(t, paramValues, funcValues);
        case THEORY_BESSEL:
          return Bessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) + 
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_BESSEL:
          return InternalBessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_SKEWED_GAUSS:
          return SkewedGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_STATIC_ZF_NK:
          return StaticNKZF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);          
        case THEORY_STATIC_TF_NK:
          return StaticNKTF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_ZF_NK:
          return DynamicNKZF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_TF_NK:
          return DynamicNKTF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_F_MU_F:
          return FmuF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_POLYNOM:
          return Polynom(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_MU_MINUS_EXP:
          return MuMinusExpTF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        case THEORY_USER_FCN:
          return UserFcn(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                 fAdd->Func(t, paramValues, funcValues);
        default:
          std::cerr << std::endl << ">> PTheory::Func: **PANIC ERROR** You never should have reached this line?!?! (" << fType << ")";
          std::cerr << std::endl;
          exit(0);
      }
    } else { // fMul !=0 && fAdd == 0
      switch (fType) {
        case THEORY_CONST:
          return Constant(paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_ASYMMETRY:
          return Asymmetry(paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_SIMPLE_EXP:
          return SimpleExp(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_GENERAL_EXP:
          return GeneralExp(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_SIMPLE_GAUSS:
          return SimpleGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_STATIC_GAUSS_KT:
          return StaticGaussKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_STATIC_GAUSS_KT_LF:
          return StaticGaussKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAUSS_KT_LF:
          return DynamicGaussKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_STATIC_LORENTZ_KT:
          return StaticLorentzKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_STATIC_LORENTZ_KT_LF:
          return StaticLorentzKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_LORENTZ_KT_LF:
          return DynamicLorentzKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_FAST_KT_ZF:
          return DynamicGauLorKTZFFast(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_FAST_KT_LF:
          return DynamicGauLorKTLFFast(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_KT_LF:
          return DynamicGauLorKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_COMBI_LGKT:
          return CombiLGKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_STR_KT:
          return StrKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_SPIN_GLASS:
          return SpinGlass(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_RANDOM_ANISOTROPIC_HYPERFINE:
          return RandomAnisotropicHyperfine(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_ABRAGAM:
          return Abragam(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_TF_COS:
          return TFCos(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD:
          return InternalField(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD_KORNILOV:
          return InternalFieldGK(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD_LARKIN:
          return InternalFieldLL(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_BESSEL:
          return Bessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_BESSEL:
          return InternalBessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_SKEWED_GAUSS:
          return SkewedGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_STATIC_ZF_NK:
          return StaticNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_STATIC_TF_NK:
          return StaticNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_ZF_NK:
          return DynamicNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_TF_NK:
          return DynamicNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_MU_MINUS_EXP:
          return MuMinusExpTF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_F_MU_F:
          return FmuF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_POLYNOM:
          return Polynom(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        case THEORY_USER_FCN:
          return UserFcn(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
        default:
          std::cerr << std::endl << ">> PTheory::Func: **PANIC ERROR** You never should have reached this line?!?! (" << fType << ")";
          std::cerr << std::endl;
          exit(0);
      }
    }
  } else { // fMul == 0 && fAdd != 0
    if (fAdd) {
      switch (fType) {
        case THEORY_CONST:
          return Constant(paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_ASYMMETRY:
          return Asymmetry(paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_SIMPLE_EXP:
          return SimpleExp(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_GENERAL_EXP:
          return GeneralExp(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_SIMPLE_GAUSS:
          return SimpleGauss(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_STATIC_GAUSS_KT:
          return StaticGaussKT(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_STATIC_GAUSS_KT_LF:
          return StaticGaussKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAUSS_KT_LF:
          return DynamicGaussKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_STATIC_LORENTZ_KT:
          return StaticLorentzKT(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_STATIC_LORENTZ_KT_LF:
          return StaticLorentzKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_LORENTZ_KT_LF:
          return DynamicLorentzKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_FAST_KT_ZF:
          return DynamicGauLorKTZFFast(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_FAST_KT_LF:
          return DynamicGauLorKTLFFast(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_KT_LF:
          return DynamicGauLorKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_COMBI_LGKT:
          return CombiLGKT(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_STR_KT:
          return StrKT(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_SPIN_GLASS:
          return SpinGlass(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_RANDOM_ANISOTROPIC_HYPERFINE:
          return RandomAnisotropicHyperfine(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_ABRAGAM:
          return Abragam(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_TF_COS:
          return TFCos(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD:
          return InternalField(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD_KORNILOV:
          return InternalFieldGK(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD_LARKIN:
          return InternalFieldLL(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_BESSEL:
          return Bessel(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_INTERNAL_BESSEL:
          return InternalBessel(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_SKEWED_GAUSS:
          return SkewedGauss(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_STATIC_ZF_NK:
          return StaticNKZF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_STATIC_TF_NK:
          return StaticNKTF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_ZF_NK:
          return DynamicNKZF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_DYNAMIC_TF_NK:
          return DynamicNKTF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_MU_MINUS_EXP:
          return MuMinusExpTF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_F_MU_F:
            return FmuF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_POLYNOM:
          return Polynom(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        case THEORY_USER_FCN:
          return UserFcn(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
        default:
          std::cerr << std::endl << ">> PTheory::Func: **PANIC ERROR** You never should have reached this line?!?! (" << fType << ")";
          std::cerr << std::endl;
          exit(0);
      }
    } else { // fMul == 0 && fAdd == 0
      switch (fType) {
        case THEORY_CONST:
          return Constant(paramValues, funcValues);
        case THEORY_ASYMMETRY:
          return Asymmetry(paramValues, funcValues);
        case THEORY_SIMPLE_EXP:
          return SimpleExp(t, paramValues, funcValues);
        case THEORY_GENERAL_EXP:
          return GeneralExp(t, paramValues, funcValues);
        case THEORY_SIMPLE_GAUSS:
          return SimpleGauss(t, paramValues, funcValues);
        case THEORY_STATIC_GAUSS_KT:
          return StaticGaussKT(t, paramValues, funcValues);
        case THEORY_STATIC_GAUSS_KT_LF:
          return StaticGaussKTLF(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAUSS_KT_LF:
          return DynamicGaussKTLF(t, paramValues, funcValues);
        case THEORY_STATIC_LORENTZ_KT:
          return StaticLorentzKT(t, paramValues, funcValues);
        case THEORY_STATIC_LORENTZ_KT_LF:
          return StaticLorentzKTLF(t, paramValues, funcValues);
        case THEORY_DYNAMIC_LORENTZ_KT_LF:
          return DynamicLorentzKTLF(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_FAST_KT_ZF:
          return DynamicGauLorKTZFFast(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_FAST_KT_LF:
          return DynamicGauLorKTLFFast(t, paramValues, funcValues);
        case THEORY_DYNAMIC_GAULOR_KT_LF:
          return DynamicGauLorKTLF(t, paramValues, funcValues);
        case THEORY_COMBI_LGKT:
          return CombiLGKT(t, paramValues, funcValues);
        case THEORY_STR_KT:
          return StrKT(t, paramValues, funcValues);
        case THEORY_SPIN_GLASS:
          return SpinGlass(t, paramValues, funcValues);
        case THEORY_RANDOM_ANISOTROPIC_HYPERFINE:
          return RandomAnisotropicHyperfine(t, paramValues, funcValues);
        case THEORY_ABRAGAM:
          return Abragam(t, paramValues, funcValues);
        case THEORY_TF_COS:
          return TFCos(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD:
          return InternalField(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD_KORNILOV:
          return InternalFieldGK(t, paramValues, funcValues);
        case THEORY_INTERNAL_FIELD_LARKIN:
          return InternalFieldLL(t, paramValues, funcValues);
        case THEORY_BESSEL:
          return Bessel(t, paramValues, funcValues);
        case THEORY_INTERNAL_BESSEL:
          return InternalBessel(t, paramValues, funcValues);
        case THEORY_SKEWED_GAUSS:
          return SkewedGauss(t, paramValues, funcValues);
        case THEORY_STATIC_ZF_NK:
          return StaticNKZF(t, paramValues, funcValues);
        case THEORY_STATIC_TF_NK:
          return StaticNKTF(t, paramValues, funcValues);
        case THEORY_DYNAMIC_ZF_NK:
          return DynamicNKZF(t, paramValues, funcValues);
        case THEORY_DYNAMIC_TF_NK:
          return DynamicNKTF(t, paramValues, funcValues);
        case THEORY_MU_MINUS_EXP:
          return MuMinusExpTF(t, paramValues, funcValues);
        case THEORY_F_MU_F:
          return FmuF(t, paramValues, funcValues);
        case THEORY_POLYNOM:
          return Polynom(t, paramValues, funcValues);
        case THEORY_USER_FCN:
          return UserFcn(t, paramValues, funcValues);
        default:
          std::cerr << std::endl << ">> PTheory::Func: **PANIC ERROR** You never should have reached this line?!?! (" << fType << ")";
          std::cerr << std::endl;
          exit(0);
      }
    }
  }
}
 
//--------------------------------------------------------------------------
void PTheory::CleanUp(PTheory *theo)
{
  if (theo->fMul) { // '*' present
     delete theo->fMul;
     theo->fMul = nullptr;
  }
 
  if (theo->fAdd) {
     delete theo->fAdd;
     theo->fAdd = nullptr;
  }
}
 
//--------------------------------------------------------------------------
Int_t PTheory::SearchDataBase(TString name)
{
  Int_t idx = THEORY_UNDEFINED;
 
  for (UInt_t i=0; i<THEORY_MAX; i++) {
    if (!name.CompareTo(fgTheoDataBase[i].fName, TString::kIgnoreCase) ||
        !name.CompareTo(fgTheoDataBase[i].fAbbrev, TString::kIgnoreCase)) {
      idx = static_cast<Int_t>(fgTheoDataBase[i].fType);
      fType = fgTheoDataBase[i].fType;
      fNoOfParam = fgTheoDataBase[i].fNoOfParam;
    }
  }
 
  return idx;
}
 
//--------------------------------------------------------------------------
// GetUserFcnIdx (private)
//--------------------------------------------------------------------------
Int_t PTheory::GetUserFcnIdx(UInt_t lineNo) const
{
  Int_t userFcnIdx = -1;
 
  // retrieve the theory block from the msr-file handler
  PMsrLines *fullTheoryBlock = fMsrInfo->GetMsrTheory();
 
  // make sure that lineNo is within proper bounds
  if (lineNo > fullTheoryBlock->size())
    return userFcnIdx;
 
  // count the number of user function present up to the lineNo
  for (UInt_t i=1; i<=lineNo; i++) {
    if (fullTheoryBlock->at(i).fLine.Contains("userFcn", TString::kIgnoreCase))
      userFcnIdx++;
  }
 
  return userFcnIdx;
}
 
//--------------------------------------------------------------------------
// MakeCleanAndTidyTheoryBlock private
//--------------------------------------------------------------------------
void PTheory::MakeCleanAndTidyTheoryBlock(PMsrLines *fullTheoryBlock)
{
  PMsrLineStructure *line;
  TString str, tidy;
  Char_t substr[256];
  Int_t idx = THEORY_UNDEFINED;
 
  for (UInt_t i=1; i<fullTheoryBlock->size(); i++) {
    // get the line to be prettyfied
    line = &(*fullTheoryBlock)[i];
    // copy line content to str in order to remove comments
    str = line->fLine.Copy();
    // remove theory line comment if present, i.e. something starting with '('
    Int_t index = str.Index("(");
    if (index > 0) // theory line comment present
      str.Resize(index);
    // tokenize line
    std::vector<std::string> tokens = PStringUtils::Split(str.Data(), " \t");
    if (tokens.empty())
      continue;
    // make a handable string out of the asymmetry token
    str = tokens[0];
    // check if the line is just a '+' if so nothing to be done
    if (str.Contains("+"))
      continue;
    // check if the function is a polynom
    if (!str.CompareTo("p") || str.Contains("polynom")) {
      MakeCleanAndTidyPolynom(i, fullTheoryBlock);
      continue;
    }
    // check if the function is a userFcn
    if (!str.CompareTo("u") || str.Contains("userFcn")) {
      MakeCleanAndTidyUserFcn(i, fullTheoryBlock);
      continue;
    }
    // search the theory function
    for (UInt_t j=0; j<THEORY_MAX; j++) {
      if (!str.CompareTo(fgTheoDataBase[j].fName, TString::kIgnoreCase) ||
          !str.CompareTo(fgTheoDataBase[j].fAbbrev, TString::kIgnoreCase)) {
        idx = static_cast<Int_t>(fgTheoDataBase[j].fType);
      }
    }
    // check if theory is indeed defined. This should not be necessay at this point but ...
    if (idx == THEORY_UNDEFINED)
      return;
    // check that there enough tokens. This should not be necessay at this point but ...
    if (static_cast<UInt_t>(tokens.size()) < fgTheoDataBase[idx].fNoOfParam + 1)
      return;
    // make tidy string
    snprintf(substr, sizeof(substr), "%-10s", fgTheoDataBase[idx].fName.Data());
    tidy = TString(substr);
    for (UInt_t j=1; j<tokens.size(); j++) {
      str = tokens[j];
      snprintf(substr, sizeof(substr), "%6s", str.Data());
      tidy += TString(substr);
    }
    if (fgTheoDataBase[idx].fComment.Length() != 0) {
      if (tidy.Length() < 35) {
        for (Int_t k=0; k<35-tidy.Length(); k++)
          tidy += TString(" ");
      } else {
        tidy += TString(" ");
      }
      if (static_cast<UInt_t>(tokens.size()) == fgTheoDataBase[idx].fNoOfParam + 1) // no tshift
        tidy += fgTheoDataBase[idx].fComment;
      else
        tidy += fgTheoDataBase[idx].fCommentTimeShift;
    }
    // write tidy string back into theory block
    (*fullTheoryBlock)[i].fLine = tidy;
  }
 
}
 
//--------------------------------------------------------------------------
// MakeCleanAndTidyPolynom private
//--------------------------------------------------------------------------
void PTheory::MakeCleanAndTidyPolynom(UInt_t i, PMsrLines *fullTheoryBlock)
{
  PMsrLineStructure *line;
  TString str, tidy;
  Char_t substr[256];
 
  // init tidy
  tidy = TString("polynom ");
  // get the line to be prettyfied
  line = &(*fullTheoryBlock)[i];
  // copy line content to str in order to remove comments
  str = line->fLine.Copy();
  // tokenize line
  std::vector<std::string> tokens = PStringUtils::Split(str.Data(), " \t");
 
  // check if comment is already present, and if yes ignore it by setting max correctly
  Int_t max = static_cast<Int_t>(tokens.size());
  for (Int_t j=1; j<max; j++) {
    str = tokens[j];
    if (str.Contains("(")) { // comment present
      max=j;
      break;
    }
  }
 
  for (Int_t j=1; j<max; j++) {
    str = tokens[j];
    snprintf(substr, sizeof(substr), "%6s", str.Data());
    tidy += TString(substr);
  }
 
  // add comment
  tidy += " (tshift p0 p1 ... pn)";
 
  // write tidy string back into theory block
  (*fullTheoryBlock)[i].fLine = tidy;
}
 
//--------------------------------------------------------------------------
// MakeCleanAndTidyUserFcn private
//--------------------------------------------------------------------------
void PTheory::MakeCleanAndTidyUserFcn(UInt_t i, PMsrLines *fullTheoryBlock)
{
  PMsrLineStructure *line;
  TString str, tidy;
 
  // init tidy
  tidy = TString("userFcn ");
  // get the line to be prettyfied
  line = &(*fullTheoryBlock)[i];
  // copy line content to str in order to remove comments
  str = line->fLine.Copy();
  // tokenize line
  std::vector<std::string> tokens = PStringUtils::Split(str.Data(), " \t");
 
  for (UInt_t j=1; j<tokens.size(); j++) {
    str = tokens[j];
    tidy += TString(" ") + str;
  }
 
  // write tidy string back into theory block
  (*fullTheoryBlock)[i].fLine = tidy;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::Constant(const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: const
 
  Double_t constant;
 
  // check if FUNCTIONS are used
  if (fParamNo[0] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
    constant = paramValues[fParamNo[0]];
  } else {
    constant = funcValues[fParamNo[0]-MSR_PARAM_FUN_OFFSET];
  }
 
  return constant;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::Asymmetry(const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: asym
 
  Double_t asym;
 
  // check if FUNCTIONS are used
  if (fParamNo[0] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
    asym = paramValues[fParamNo[0]];
  } else {
    asym = funcValues[fParamNo[0]-MSR_PARAM_FUN_OFFSET];
  }
 
  return asym;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::SimpleExp(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: lambda [tshift]
 
  Double_t val[2];
 
  assert(fParamNo.size() <= 2);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 1) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[1];
 
  return TMath::Exp(-tt*val[0]);
}
 
//--------------------------------------------------------------------------
Double_t PTheory::GeneralExp(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: lambda beta [tshift]
 
  Double_t val[3];
  Double_t result;
 
  assert(fParamNo.size() <= 3);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 2) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[2];
 
  // check if tt*val[0] < 0 and 
  if ((tt*val[0] < 0) && (trunc(val[1])-val[1] != 0.0)) {
    result = 0.0;
  } else {
    result = TMath::Exp(-TMath::Power(tt*val[0], val[1]));
  }
 
  return result;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::SimpleGauss(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: sigma [tshift]
 
  Double_t val[2];
 
  assert(fParamNo.size() <= 2);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 1) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[1];
 
  return TMath::Exp(-0.5*TMath::Power(tt*val[0], 2.0));
}
 
//--------------------------------------------------------------------------
Double_t PTheory::StaticGaussKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: sigma [tshift]
 
  Double_t val[2];
 
  assert(fParamNo.size() <= 2);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t sigma_t_2;
  if (fParamNo.size() == 1) // no tshift
    sigma_t_2 = t*t*val[0]*val[0];
  else // tshift present
    sigma_t_2 = (t-val[1])*(t-val[1])*val[0]*val[0];
 
  return 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*TMath::Exp(-0.5*sigma_t_2));
}
 
//--------------------------------------------------------------------------
Double_t PTheory::StaticGaussKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
 
  // expected parameters: frequency damping [tshift]
 
  Double_t val[3];
  Double_t result;
 
  assert(fParamNo.size() <= 3);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  // check if all parameters == 0
  if ((val[0] == 0.0) && (val[1] == 0.0))
    return 1.0;
 
  // check if the parameter values have changed, and if yes recalculate the non-analytic integral
  // check only the first two parameters since the tshift is irrelevant for the LF-integral calculation!!
  Bool_t newParam = false;
  for (UInt_t i=0; i<2; i++) {
    if (val[i] != fPrevParam[i]) {
      newParam = true;
      break;
    }
  }
 
  if (newParam) { // new parameters found
    {
      for (UInt_t i=0; i<2; i++)
        fPrevParam[i] = val[i];
      CalculateGaussLFIntegral(val);
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 2) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[2];
 
  if (tt < 0.0) // for times < 0 return a function value of 1.0
    return 1.0;
 
  Double_t sigma_t_2 = 0.0;
  if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use the ZF formula
    sigma_t_2 = tt*tt*val[1]*val[1];
    result = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*TMath::Exp(-0.5*sigma_t_2));
  } else if (val[1]/val[0] > 79.5775) { // check if Delta/w0 > 500.0, in which case the ZF formula is used
    sigma_t_2 = tt*tt*val[1]*val[1];
    result = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*TMath::Exp(-0.5*sigma_t_2));
  } else {
    Double_t delta = val[1];
    Double_t w0    = 2.0*TMath::Pi()*val[0];
 
    result = 1.0 - 2.0*TMath::Power(delta/w0,2.0)*(1.0 - 
                TMath::Exp(-0.5*TMath::Power(delta*tt, 2.0))*TMath::Cos(w0*tt)) +
                GetLFIntegralValue(tt);
  }
 
  return result;
 
}
 
//--------------------------------------------------------------------------
Double_t PTheory::DynamicGaussKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: frequency damping hopping [tshift]
 
  Double_t val[4];
  Double_t result = 0.0;
  Bool_t useKeren = false;
 
  assert(fParamNo.size() <= 4);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  // check if all parameters == 0
  if ((val[0] == 0.0) && (val[1] == 0.0) && (val[2] == 0.0))
    return 1.0;
 
  // make sure that damping and hopping are positive definite
  if (val[1] < 0.0)
    val[1] = -val[1];
  if (val[2] < 0.0)
    val[2] = -val[2];
 
  // check that Delta != 0, if not (i.e. stupid parameter) return 1, which is the correct limit
  if (fabs(val[1]) < 1.0e-6) {
    return 1.0;
  }
 
  // check if Keren approximation can be used
  if (val[2]/val[1] > 5.0) // nu/Delta > 5.0
    useKeren = true;
 
  if (!useKeren) {
    // check if the parameter values have changed, and if yes recalculate the non-analytic integral
    // check only the first three parameters since the tshift is irrelevant for the LF-integral calculation!!
    Bool_t newParam = false;
    for (UInt_t i=0; i<3; i++) {
      if (val[i] != fPrevParam[i]) {
        newParam = true;
        break;
      }
    }
 
    if (newParam) { // new parameters found
      for (UInt_t i=0; i<3; i++)
        fPrevParam[i] = val[i];
      CalculateDynKTLF(val, 0); // 0 means Gauss
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 3) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[3];
 
  if (tt < 0.0) // for times < 0 return a function value of 0.0
    return 0.0;
 
 
  if (useKeren) { // see PRB50, 10039 (1994)
    Double_t wL  = TWO_PI * val[0];
    Double_t wL2 = wL*wL;
    Double_t nu2 = val[2]*val[2];
    Double_t Gamma_t = 2.0*val[1]*val[1]/((wL2+nu2)*(wL2+nu2))*
                       ((wL2+nu2)*val[2]*t
                        + (wL2-nu2)*(1.0 - TMath::Exp(-val[2]*t)*TMath::Cos(wL*t))
                        - 2.0*val[2]*wL*TMath::Exp(-val[2]*t)*TMath::Sin(wL*t));
    result = TMath::Exp(-Gamma_t);
  } else { // from Voltera
    result = GetDynKTLFValue(tt);
  }
 
  return result;
 
}
 
//--------------------------------------------------------------------------
Double_t PTheory::StaticLorentzKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: lambda [tshift]
 
  Double_t val[2];
 
  assert(fParamNo.size() <= 2);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t a_t;
  if (fParamNo.size() == 1) // no tshift
    a_t = t*val[0];
  else // tshift present
    a_t = (t-val[1])*val[0];
 
  return 0.333333333333333 * (1.0 + 2.0*(1.0 - a_t)*TMath::Exp(-a_t));
}
 
//--------------------------------------------------------------------------
Double_t PTheory::StaticLorentzKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: frequency damping [tshift]
 
  Double_t val[3];
  Double_t result;
 
  assert(fParamNo.size() <= 3);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  // check if all parameters == 0
  if ((val[0] == 0.0) && (val[1] == 0.0))
    return 1.0;
 
  // check if the parameter values have changed, and if yes recalculate the non-analytic integral
  // check only the first two parameters since the tshift is irrelevant for the LF-integral calculation!!
  Bool_t newParam = false;
  for (UInt_t i=0; i<2; i++) {
    if (val[i] != fPrevParam[i]) {
      newParam = true;
      break;
    }
  }
 
  if (newParam) { // new parameters found
    for (UInt_t i=0; i<2; i++)
      fPrevParam[i] = val[i];
    CalculateLorentzLFIntegral(val);
  }
 
  Double_t tt;
  if (fParamNo.size() == 2) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[2];
 
  if (tt < 0.0) // for times < 0 return a function value of 1.0
    return 1.0;
 
 if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use the ZF formula
    Double_t at = tt*val[1];
    result = 0.333333333333333 * (1.0 + 2.0*(1.0 - at)*TMath::Exp(-at));
  } else if (val[1]/val[0] > 159.1549) { // check if a/w0 > 1000.0, in which case the ZF formula is used
    Double_t at = tt*val[1];
    result = 0.333333333333333 * (1.0 + 2.0*(1.0 - at)*TMath::Exp(-at));
  } else {
    Double_t a    = val[1];
    Double_t at   = a*tt;
    Double_t w0   = 2.0*TMath::Pi()*val[0];
    Double_t a_w0 = a/w0;
    Double_t w0t  = w0*tt;
 
    Double_t j1, j0;
    if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
      j0 = 1.0;
      j1 = 0.0;
    } else {
      j0 = sin(w0t)/w0t;
      j1 = (sin(w0t)-w0t*cos(w0t))/(w0t*w0t);
    }
 
    result = 1.0 - a_w0*j1*exp(-at) - a_w0*a_w0*(j0*exp(-at) - 1.0) - GetLFIntegralValue(tt);
  }
 
  return result;
 
}
 
//--------------------------------------------------------------------------
Double_t PTheory::DynamicLorentzKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: frequency damping hopping [tshift]
 
  Double_t val[4];
  Double_t result = 0.0;
 
  assert(fParamNo.size() <= 4);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  // check if all parameters == 0
  if ((val[0] == 0.0) && (val[1] == 0.0) && (val[2] == 0.0))
    return 1.0;
 
  // make sure that damping and hopping are positive definite
  if (val[1] < 0.0)
    val[1] = -val[1];
  if (val[2] < 0.0)
    val[2] = -val[2];
 
 
  Double_t tt;
  if (fParamNo.size() == 3) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[3];
 
  if (tt < 0.0) // for times < 0 return a function value of 1.0
    return 1.0;
 
  // check if hopping > 5 * damping, of Larmor angular frequency is > 30 * damping (BMW limit)
  Double_t w0 = 2.0*TMath::Pi()*val[0];
  Double_t a  = val[1];
  Double_t nu = val[2];
  if ((nu > 5.0 * a) || (w0 >= 30.0 * a)) {
    // 'c' and 'd' are parameters BMW obtained by fitting large parameter space LF-curves to the model below
    const Double_t c[7] = {1.15331, 1.64826, -0.71763, 3.0, 0.386683, -5.01876, 2.41854};
    const Double_t d[4] = {2.44056, 2.92063, 1.69581, 0.667277};
    Double_t w0N[4];
    Double_t nuN[4];
    w0N[0] = w0;
    nuN[0] = nu;
    for (UInt_t i=1; i<4; i++) {
      w0N[i] = w0 * w0N[i-1];
      nuN[i] = nu * nuN[i-1];
    }
    Double_t denom = w0N[3]+d[0]*w0N[2]*nuN[0]+d[1]*w0N[1]*nuN[1]+d[2]*w0N[0]*nuN[2]+d[3]*nuN[3];
    Double_t b1 = (c[0]*w0N[2]+c[1]*w0N[1]*nuN[0]+c[2]*w0N[0]*nuN[1])/denom;
    Double_t b2 = (c[3]*w0N[2]+c[4]*w0N[1]*nuN[0]+c[5]*w0N[0]*nuN[1]+c[6]*nuN[2])/denom;
 
    Double_t w0t = w0*tt;
    Double_t j1, j0;
    if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
      j0 = 1.0;
      j1 = 0.0;
    } else {
      j0 = sin(w0t)/w0t;
      j1 = (sin(w0t)-w0t*cos(w0t))/(w0t*w0t);
    }
 
    Double_t Gamma_t = -4.0/3.0*a*(b1*(1.0-j0*TMath::Exp(-nu*tt))+b2*j1*TMath::Exp(-nu*tt)+(1.0-b2*w0/3.0-b1*nu)*tt);
 
    return TMath::Exp(Gamma_t);
  }
 
  // check if the parameter values have changed, and if yes recalculate the non-analytic integral
  // check only the first three parameters since the tshift is irrelevant for the LF-integral calculation!!
  Bool_t newParam = false;
  for (UInt_t i=0; i<3; i++) {
    if (val[i] != fPrevParam[i]) {
      newParam = true;
      break;
    }
  }
 
  if (newParam) { // new parameters found
    for (UInt_t i=0; i<3; i++)
      fPrevParam[i] = val[i];
    CalculateDynKTLF(val, 1); // 1 means Lorentz
  }
 
  result = GetDynKTLFValue(tt);
 
  return result;
 
}
 
//--------------------------------------------------------------------------
Double_t PTheory::DynamicGauLorKTZFFast(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: damping hopping [tshift]
 
  Double_t val[3];
 
  assert(fParamNo.size() <= 3);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 2) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[2];
 
  Double_t nut = val[1]*tt;
  return exp(-sqrt(4.0*pow(val[0]/val[1], 2.0)*(exp(-nut)-1.0+nut)));
}
 
//--------------------------------------------------------------------------
Double_t PTheory::DynamicGauLorKTLFFast(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: frequency damping hopping [tshift]
 
  Double_t val[4];
 
  assert(fParamNo.size() <= 4);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 3) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[3];
 
  Double_t w0 = TMath::TwoPi()*val[0];
  Double_t w0_2 = w0*w0;
  Double_t nu_2 = val[2]*val[2];
  Double_t nu_t = val[2]*tt;
  Double_t w0_t = w0*tt;
  Double_t Gamma_t = ((w0_2+nu_2)*nu_t+(w0_2-nu_2)*(1.0-exp(-nu_t)*cos(w0_t))-2.0*val[2]*w0*exp(-nu_t)*sin(w0_t))/pow(w0_2+nu_2,2.0);
  if (Gamma_t < 0.0)
    Gamma_t = 0.0;
 
  return exp(-sqrt(4.0*val[1]*val[1]*Gamma_t));
}
 
//--------------------------------------------------------------------------
Double_t PTheory::DynamicGauLorKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: frequency damping hopping [tshift]
 
  Double_t val[4];
 
  assert(fParamNo.size() <= 4);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 3) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[3];
 
  // check if the parameter values have changed, and if yes recalculate DynamicGaussKTLF
  Bool_t newParam = false;
  for (UInt_t i=0; i<3; i++) {
    if (val[i] != fPrevParam[i]) {
      newParam = true;
      break;
    }
  }
 
  if (newParam) { // new parameters found, hence calculate DynamicGauLorKTLF
    // keep parameters
    for (UInt_t i=0; i<3; i++)
      fPrevParam[i] = val[i];
 
    // reset GL LF polarzation vector
    fDyn_GL_LFFuncValue.clear();
    fDyn_GL_LFFuncValue.resize(20000); // Tmax=20us, dt=1ns
 
    PDoubleVector rr={0.2, 0.4, 0.6, 0.8, 1.0, 1.25, 1.5, 1.75, 2.0, 2.5, 3.0, 4.0, 5.0, 7.5, 10.0,
                      12.8125, 15.625, 18.4375, 21.25, 26.875, 32.5, 43.75, 55.0, 77.5, 100.0};
    Double_t par[3] = {val[0], val[1], val[2]};
    Double_t sqrtTwoInv = 1.0/sqrt(2.0);
    Bool_t isOneVec{false};
    Bool_t useKeren{false};
    Double_t scale, up{0.0}, low{-1.0};
 
    for (UInt_t i=0; i<rr.size(); i++) {
      useKeren = false;
      isOneVec = false;
 
      // Delta_G = rr * Delta_L
      par[1] = rr[i] * val[1];
 
      // check if all parameters == 0
      if ((par[0] == 0.0) && (par[1] == 0.0) && (par[2] == 0.0)) {
        isOneVec = true;
      }
 
      // make sure that damping and hopping are positive definite
      if (par[1] < 0.0)
        par[1] = -par[1];
      if (par[2] < 0.0)
        par[2] = -par[2];
 
      // check that Delta != 0, if not (i.e. stupid parameter) return 1, which is the correct limit
      if (fabs(par[1]) < 1.0e-6) {
        isOneVec = true;
      }
 
      // check if Keren approximation can be used
      if (par[2]/par[1] > 5.0) // nu/Delta > 5.0
        useKeren = true;
 
      if (!useKeren && !isOneVec) {
        CalculateDynKTLF(par, 0); // 0 means Gauss
      }
 
      // calculate polarization vector for the given parameters
      up = -std::erf(sqrtTwoInv/rr[i]);
      scale = up - low;
      low = up;
 
      const Double_t dt=0.001;
      for (UInt_t i=0; i<20000; i++) {
        if (isOneVec) {
          fDyn_GL_LFFuncValue[i] += scale;
        } else if (useKeren && !isOneVec) {// see PRB50, 10039 (1994)
          Double_t wL  = TWO_PI * par[0];
          Double_t wL2 = wL*wL;
          Double_t nu2 = par[2]*par[2];
          Double_t Gamma_t = 2.0*par[1]*par[1]/((wL2+nu2)*(wL2+nu2))*
                             ((wL2+nu2)*val[2]*i*dt
                              + (wL2-nu2)*(1.0 - TMath::Exp(-val[2]*i*dt)*TMath::Cos(wL*i*dt))
                              - 2.0*val[2]*wL*TMath::Exp(-val[2]*i*dt)*TMath::Sin(wL*i*dt));
          fDyn_GL_LFFuncValue[i] += scale*TMath::Exp(-Gamma_t);
        } else if (!useKeren && !isOneVec) {
          fDyn_GL_LFFuncValue[i] += scale*GetDynKTLFValue(i*dt);
        }
      }
    }
  }
 
 
  // get the proper value from the look-up table
  Double_t result{1.0};
  if (tt>=0)
    result=GetDyn_GL_KTLFValue(tt);
 
  return result;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::CombiLGKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: lambdaL lambdaG [tshift]
 
  Double_t val[3];
 
  assert(fParamNo.size() <= 3);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 2) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[2];
 
  Double_t lambdaL_t   = tt*val[0];
  Double_t lambdaG_t_2 = tt*tt*val[1]*val[1];
 
  return 0.333333333333333 * 
         (1.0 + 2.0*(1.0-lambdaL_t-lambdaG_t_2)*TMath::Exp(-(lambdaL_t+0.5*lambdaG_t_2)));
}
 
//--------------------------------------------------------------------------
Double_t PTheory::StrKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: sigma beta [tshift]
 
  Double_t val[3];
 
  assert(fParamNo.size() <= 3);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  // check for beta too small (beta < 0.1) in which case numerical problems could arise and the function is anyhow
  // almost identical to a constant of 1/3.
  if (val[1] < 0.1)
    return 0.333333333333333;
 
  Double_t tt;
  if (fParamNo.size() == 2) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[2];
 
  Double_t sigma_t   = TMath::Power(tt*val[0],val[1]);
 
  return 0.333333333333333 * 
    (1.0 + 2.0*(1.0-sigma_t)*TMath::Exp(-sigma_t/val[1]));
}
 
//--------------------------------------------------------------------------
Double_t PTheory::SpinGlass(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: lambda gamma q [tshift]
 
  if (paramValues[fParamNo[0]] == 0.0)
    return 1.0;
 
  Double_t val[4];
 
  assert(fParamNo.size() <= 4);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 3) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[3];
 
  Double_t lambda_2     = val[0]*val[0];
  Double_t lambda_t_2_q = tt*tt*lambda_2*val[2];
  Double_t rate_2       = 4.0*lambda_2*(1.0-val[2])*tt/val[1];
 
  Double_t rateL = TMath::Sqrt(fabs(rate_2));
  Double_t rateT = TMath::Sqrt(fabs(rate_2)+lambda_t_2_q);
 
  return 0.333333333333333*(TMath::Exp(-rateL) + 2.0*(1.0-lambda_t_2_q/rateT)*TMath::Exp(-rateT));
}
 
//--------------------------------------------------------------------------
Double_t PTheory::RandomAnisotropicHyperfine(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: nu lambda [tshift]
 
  Double_t val[3];
 
  assert(fParamNo.size() <= 3);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 2) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[2];
 
  Double_t nu_t     = tt*val[0];
  Double_t lambda_t = tt*val[1];
 
  return 0.166666666666667*(1.0-0.5*nu_t)*TMath::Exp(-0.5*nu_t) +
         0.333333333333333*(1.0-0.25*nu_t)*TMath::Exp(-0.25*(nu_t+2.44949*lambda_t));
}
 
//--------------------------------------------------------------------------
Double_t PTheory::Abragam(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: sigma gamma [tshift]
 
  Double_t val[3];
 
  assert(fParamNo.size() <= 3);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 2) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[2];
 
  Double_t gamma_t = tt*val[1];
 
  return TMath::Exp(-TMath::Power(val[0]/val[1],2.0)*
                    (TMath::Exp(-gamma_t)-1.0+gamma_t));
}
 
//--------------------------------------------------------------------------
Double_t PTheory::TFCos(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: phase frequency [tshift]
 
  Double_t val[3];
 
  assert(fParamNo.size() <= 3);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 2) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[2];
 
  return TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt);
}
 
//--------------------------------------------------------------------------
Double_t PTheory::InternalField(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: fraction phase frequency rateT rateL [tshift]
 
  Double_t val[6];
 
  assert(fParamNo.size() <= 6);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 5) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[5];
 
  return val[0]*TMath::Cos(DEG_TO_RAD*val[1]+TWO_PI*val[2]*tt)*TMath::Exp(-val[3]*tt) +
         (1-val[0])*TMath::Exp(-val[4]*tt);
}
 
//--------------------------------------------------------------------------
Double_t PTheory::InternalFieldGK(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: [0]:fraction [1]:frequency [2]:sigma [3]:lambda [4]:beta [[5]:tshift]
 
  Double_t val[6];
 
  assert(fParamNo.size() <= 6);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 5) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[5];
 
  Double_t result = 0.0;
  Double_t w_t = TWO_PI*val[1]*tt;
  Double_t rateLF = TMath::Power(val[3]*tt, val[4]);
  Double_t rate2 = val[2]*val[2]*tt*tt; // (sigma t)^2
 
  if (val[1] < 0.01) { // internal field frequency is approaching zero
    result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(1.0-rate2)*TMath::Exp(-0.5*rate2);
  } else {
    result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(TMath::Cos(w_t)-val[2]*val[2]*tt/(TWO_PI*val[1])*TMath::Sin(w_t))*TMath::Exp(-0.5*rate2);
  }
 
  return result;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::InternalFieldLL(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: [0]:fraction [1]:frequency [2]:a [3]:lambda [4]:beta [[5]:tshift]
 
  Double_t val[6];
 
  assert(fParamNo.size() <= 6);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 5) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[5];
 
  Double_t result = 0.0;
  Double_t w_t = TWO_PI*val[1]*tt;
  Double_t rateLF = TMath::Power(val[3]*tt, val[4]);
  Double_t a_t = val[2]*tt; // a t
 
  if (val[1] < 0.01) { // internal field frequency is approaching zero
    result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(1.0-a_t)*TMath::Exp(-a_t);
  } else {
    result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(TMath::Cos(w_t)-val[3]/(TWO_PI*val[1])*TMath::Sin(w_t))*TMath::Exp(-a_t);
  }
 
  return result;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::Bessel(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: phase frequency [tshift]
 
  Double_t val[3];
 
  assert(fParamNo.size() <= 3);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 2) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[2];
 
  return TMath::BesselJ0(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt);
}
 
//--------------------------------------------------------------------------
Double_t PTheory::InternalBessel(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: fraction phase frequency rateT rateL [tshift]
 
  Double_t val[6];
 
  assert(fParamNo.size() <= 6);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 5) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[5];
 
  return val[0]* TMath::BesselJ0(DEG_TO_RAD*val[1]+TWO_PI*val[2]*tt)*
                 TMath::Exp(-val[3]*tt) +
         (1.0-val[0])*TMath::Exp(-val[4]*tt);
}
 
//--------------------------------------------------------------------------
Double_t PTheory::SkewedGauss(Double_t t, const PDoubleVector &paramValues,
                              const PDoubleVector &funcValues) const {
  // Expected parameters: phase, frequency, sigma-, sigma+, [tshift].
  // To be stored in the array "val" as:
  // val[0] = phase
  // val[1] = frequency
  // val[2] = sigma-
  // val[3] = sigma+
  // val[4] = tshift [optional]
  Double_t val[5];
 
  // Check that we have the correct number of fit parameters.
  assert(fParamNo.size() <= 5);
 
  // Check if FUNCTIONS are used.
  for (UInt_t i = 0; i < fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i] - MSR_PARAM_FUN_OFFSET];
    }
  }
 
  // Apply the tshift (if required).
  Double_t tt = t;
  if (fParamNo.size() == 5) {
    tt = t - val[4];
  }
 
  // Evaluate the skewed Gaussian!
 
  // First, calculate some "helper" terms.
  Double_t sigma_p = std::abs(val[3]);
  Double_t sigma_m = std::abs(val[2]);
  Double_t arg_p = sigma_p * tt;
  Double_t arg_m = sigma_m * tt;
  Double_t z_p = arg_p / SQRT_TWO;                 // sigma+
  Double_t z_m = arg_m / SQRT_TWO;                 // sigma-
  Double_t g_p = TMath::Exp(-0.5 * arg_p * arg_p); // gauss sigma+
  Double_t g_m = TMath::Exp(-0.5 * arg_m * arg_m); // gauss sigma-
  Double_t w_p = sigma_p / (sigma_p + sigma_m);
  Double_t w_m = 1.0 - w_p;
  Double_t phase = DEG_TO_RAD * val[0];
  Double_t freq = TWO_PI * val[1];
 
  // Evalute the EVEN frequency component of the skewed Gaussian.
  Double_t skg_cos = TMath::Cos(phase + freq * tt) * (w_m * g_m + w_p * g_p);
 
  // Evalute the ODD frequency component of the skewed Gaussian.
  constexpr Double_t z_max = 26.7776;
  // Note: the check against z_max is needed to prevent floating-point overflow
  // in the return value of ROOT::Math::conf_hyperg(1/2, 3/2, z * z)
  // (i.e., confluent hypergeometric function of the first kind, 1F1).
  // In the case that z > z_max, return zero (otherwise there is some
  // numeric discontinuity at later times).
  Double_t skg_sin =
      TMath::Sin(phase + freq * tt) *
      ((z_m > z_max) or (z_p > z_max)
           ? 0.0
           : (w_m * g_m * 2.0 * z_m / SQRT_PI) *
                     ROOT::Math::conf_hyperg(0.5, 1.5, z_m * z_m) -
                 (w_p * g_p * 2.0 * z_p / SQRT_PI) *
                     ROOT::Math::conf_hyperg(0.5, 1.5, z_p * z_p));
 
  // Return the skewed Gaussian: skg = skg_cos + skg_sin.
  // Also check that skg_sin is finite!
  return skg_cos + (std::isfinite(skg_sin) ? skg_sin : 0.0);
}
 
//--------------------------------------------------------------------------
Double_t PTheory::StaticNKZF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected paramters: damping_D0 [0], R_b tshift [1]
 
  Double_t val[3];
  Double_t result = 1.0;
 
  assert(fParamNo.size() <= 3);
 
  if (t < 0.0)
    return result;
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 2) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[2];
 
  Double_t D2_t2 = val[0]*val[0]*tt*tt;
  Double_t denom = 1.0+val[1]*val[1]*D2_t2;
 
  result = 0.333333333333333 + 0.666666666666666667 * TMath::Power(1.0/denom, 1.5) * (1.0 - (D2_t2/denom)) * exp(-0.5*D2_t2/denom);
 
  return result;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::StaticNKTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected paramters: phase [0], frequency [1], damping_D0 [2], R_b [3],  [tshift [4]]
 
  Double_t val[5];
  Double_t result = 1.0;
 
  assert(fParamNo.size() <= 5);
 
  if (t < 0.0)
    return result;
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 4) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[4];
 
  Double_t D2_t2 = val[2]*val[2]*tt*tt;
  Double_t denom = 1.0+val[3]*val[3]*D2_t2;
 
  result = sqrt(1.0/denom)*exp(-0.5*D2_t2/denom)*TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt);
 
  return result;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::DynamicNKZF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected paramters: damping_D0 [0], R_b [1], nu_c [2], [tshift [3]]
 
  Double_t val[4];
  Double_t result = 1.0;
 
  assert(fParamNo.size() <= 4);
 
  if (t < 0.0)
    return result;
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 3) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[3];
 
  Double_t theta;
  if (val[2] < 1.0e-6) { // nu_c -> 0
    theta = 0.5*tt*tt;
  } else {
    theta = (exp(-val[2]*tt) - 1.0 + val[2]*tt)/(val[2]*val[2]);
  }
  Double_t denom = 1.0 + 4.0*val[0]*val[0]*val[1]*val[1]*theta;
 
  result = sqrt(1.0/denom)*exp(-2.0*val[0]*val[0]*theta/denom);
 
  return result;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::DynamicNKTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected paramters: phase [0], frequency [1], damping_D0 [2], R_b [3], nu_c [4], [tshift [5]]
 
  Double_t val[6];
  Double_t result = 1.0;
 
  assert(fParamNo.size() <= 6);
 
  if (t < 0.0)
    return result;
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 5) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[5];
 
  Double_t theta;
  if (val[4] < 1.0e-6) { // nu_c -> 0
    theta = 0.5*tt*tt;
  } else {
    theta = (exp(-val[4]*tt) - 1.0 + val[4]*tt)/(val[4]*val[4]);
  }
  Double_t denom = 1.0 + 2.0*val[2]*val[2]*val[3]*val[3]*theta;
 
  result = sqrt(1.0/denom)*exp(-val[2]*val[2]*theta/denom)*TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt);
 
  return result;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::FmuF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
    // expected paramters: w_d [0], [tshift [1]]
 
    Double_t val[2];
 
    assert(fParamNo.size() <= 2);
 
    if (t < 0.0)
      return 1.0;
 
    // check if FUNCTIONS are used
    for (UInt_t i=0; i<fParamNo.size(); i++) {
      if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
        val[i] = paramValues[fParamNo[i]];
      } else { // function
        val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
      }
    }
 
    Double_t tt;
    if (fParamNo.size() == 1) // no tshift
      tt = t;
    else // tshift present
      tt = t-val[1];
 
    const Double_t sqrt3 = sqrt(3.0);
    const Double_t wd_t = val[0]*tt;
 
    return (3.0+cos(sqrt3*wd_t)+(1.0-1.0/sqrt3)*cos(((3.0-sqrt3)/2.0)*wd_t)+(1.0+1.0/sqrt3)*cos(((3.0 + sqrt3)/2.0)*wd_t))/6.0;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::Polynom(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: tshift p0 p1 p2 ...
 
  Double_t result = 0.0;
  Double_t tshift = 0.0;
  Double_t val;
  Double_t expo = 0.0;
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val = paramValues[fParamNo[i]];
    } else { // function
      val = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
    if (i==0) { // tshift
      tshift = val;
      continue;
    }
    result += val*pow(t-tshift, expo);
    expo++;
  }
 
  return result;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::UserFcn(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fUserParam.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      fUserParam[i] = paramValues[fParamNo[i]];
    } else { // function
      fUserParam[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  return (*fUserFcn)(t, fUserParam);
}
 
//--------------------------------------------------------------------------
void PTheory::CalculateGaussLFIntegral(const Double_t *val) const
{
  // val[0] = nu (field), val[1] = Delta
 
  if (val[0] == 0.0) { // field == 0.0, hence nothing to be done
    return;
  } else if (val[1]/val[0] > 79.5775) { // check if a/w0 > 500.0, in which case the ZF formula is used and here nothing has to be done
    return;
  }
 
 
  Double_t dt=0.001; // all times in usec
  Double_t t, ft;
  Double_t w0 = TMath::TwoPi()*val[0];
  Double_t Delta = val[1];
  Double_t preFactor = 2.0*TMath::Power(Delta, 4.0) / TMath::Power(w0, 3.0);
 
  // check if the time resolution needs to be increased
  const Int_t samplingPerPeriod = 20;
  const Int_t samplingOnExp = 3000;
  if ((Delta <= w0) && (1.0/val[0] < 20.0)) { // makes sure that the frequency sampling is fine enough
    if (1.0/val[0]/samplingPerPeriod < 0.001) {
      dt = 1.0/val[0]/samplingPerPeriod;
    }
  } else if ((Delta > w0) && (Delta <= 10.0)) {
    if (Delta/w0 > 10.0) {
      dt = 0.00005;
    }
  } else if ((Delta > w0) && (Delta > 10.0)) { // makes sure there is a fine enough sampling for large Delta's
    if (1.0/Delta/samplingOnExp < 0.001) {
      dt = 1.0/Delta/samplingOnExp;
    }
  }
 
  // keep sampling time
  fSamplingTime = dt;
 
  // clear previously allocated vector
  fLFIntegral.clear();
 
  // calculate integral
  t  = 0.0;
  fLFIntegral.push_back(0.0); // start value of the integral
 
  ft = 0.0;
  Double_t step = 0.0, lastft = 1.0, diff = 0.0;
  do {
    t  += dt;
    step = 0.5*dt*preFactor*(exp(-0.5*pow(Delta * (t-dt), 2.0))*sin(w0*(t-dt))+
                             exp(-0.5*pow(Delta * t, 2.0))*sin(w0*t));
    ft += step;
    diff = fabs(fabs(lastft)-fabs(ft));
    lastft = ft;
    fLFIntegral.push_back(ft);
  } while ((t <= 20.0) && (diff > 1.0e-10));
}
 
//--------------------------------------------------------------------------
void PTheory::CalculateLorentzLFIntegral(const Double_t *val) const
{
  // val[0] = nu, val[1] = a
 
  // a few checks if the integral actually needs to be calculated
  if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use the ZF formula and here nothing has to be done
    return;
  } else if (val[1]/val[0] > 159.1549) { // check if a/w0 > 1000.0, in which case the ZF formula is used and here nothing has to be done
    return;
  }
 
  Double_t dt=0.001; // all times in usec
  Double_t t, ft;
  Double_t w0 = TMath::TwoPi()*val[0];
  Double_t a = val[1];
  Double_t preFactor = a*(1+pow(a/w0,2.0));
 
  // check if the time resolution needs to be increased
  const Int_t samplingPerPeriod = 20;
  const Int_t samplingOnExp = 3000;
  if ((a <= w0) && (1.0/val[0] < 20.0)) { // makes sure that the frequency sampling is fine enough
    if (1.0/val[0]/samplingPerPeriod < 0.001) {
      dt = 1.0/val[0]/samplingPerPeriod;
    }
  } else if ((a > w0) && (a <= 10.0)) {
    if (a/w0 > 10.0) {
      dt = 0.00005;
    }
  } else if ((a > w0) && (a > 10.0)) { // makes sure there is a fine enough sampling for large a's
    if (1.0/a/samplingOnExp < 0.001) {
      dt = 1.0/a/samplingOnExp;
    }
  }
 
  // keep sampling time
  fSamplingTime = dt;
 
  // clear previously allocated vector
  fLFIntegral.clear();
 
  // calculate integral
  t  = 0.0;
  fLFIntegral.push_back(0.0); // start value of the integral
 
  ft = 0.0;
  // calculate first integral bin value (needed bcause of sin(x)/x x->0)
  t  += dt;
  ft += 0.5*dt*preFactor*(1.0+sin(w0*t)/(w0*t)*exp(-a*t));
  fLFIntegral.push_back(ft);
  // calculate all the other integral bin values
  Double_t step = 0.0, lastft = 1.0, diff = 0.0;
  do {
    t  += dt;
    step = 0.5*dt*preFactor*(sin(w0*(t-dt))/(w0*(t-dt))*exp(-a*(t-dt))+sin(w0*t)/(w0*t)*exp(-a*t));
    ft += step;
    diff = fabs(fabs(lastft)-fabs(ft));
    lastft = ft;
    fLFIntegral.push_back(ft);
  } while ((t <= 20.0) && (diff > 1.0e-10));
}
 
 
//--------------------------------------------------------------------------
Double_t PTheory::GetLFIntegralValue(const Double_t t) const
{
  if (t < 0.0)
    return 0.0;
 
  UInt_t idx = static_cast<UInt_t>(t/fSamplingTime);
 
  if (idx + 2 > fLFIntegral.size())
    return fLFIntegral.back();
 
  // linearly interpolate between the two relevant function bins
  Double_t df = (fLFIntegral[idx+1]-fLFIntegral[idx])*(t/fSamplingTime-static_cast<Double_t>(idx));
 
  return fLFIntegral[idx]+df;
}
 
//--------------------------------------------------------------------------
void PTheory::CalculateDynKTLF(const Double_t *val, Int_t tag) const
{
  // val: 0=nu0, 1=Delta (Gauss) / a (Lorentz), 2=nu
  const Double_t Tmax = 20.0; // 20 usec
  UInt_t N = static_cast<UInt_t>(16.0*Tmax*val[0]);
 
  // check if rate (Delta or a) is very high
  if (fabs(val[1]) > 0.1) {
    Double_t tmin = 20.0;
    switch (tag) {
      case 0: // Gauss
        tmin = fabs(sqrt(3.0)/val[1]);
        break;
      case 1: // Lorentz
        tmin = fabs(2.0/val[1]);
        break;
      default:
        break;
    }
    UInt_t Nrate = static_cast<UInt_t>(25.0 * Tmax / tmin);
    if (Nrate > N) {
      N = Nrate;
    }
  }
 
  if (N < 300) // if too few points, i.e. nu0 very small, take 300 points
    N = 300;
 
  if (N>1e6) // make sure that N is not too large to prevent memory overflow
    N = 1e6;
 
  // allocate memory for dyn KT LF function vector
  fDynLFFuncValue.clear(); // get rid of a possible previous vector
  fDynLFFuncValue.resize(N);
 
  // calculate the non-analytic integral of the static KT LF function
  switch (tag) {
    case 0: // Gauss
      CalculateGaussLFIntegral(val);
      break;
    case 1: // Lorentz
      CalculateLorentzLFIntegral(val);
      break;
    default:
      std::cerr << std::endl << ">> PTheory::CalculateDynKTLF: **FATAL ERROR** You should never have reached this point." << std::endl;
      assert(false);
      break;
  }
 
  // calculate the P^(0)(t) exp(-nu t) vector
  PDoubleVector p0exp(N);
  Double_t t = 0.0;
  Double_t dt = Tmax/N;
  fDynLFdt = dt; // keep it since it is needed in GetDynKTLFValue()
  for (UInt_t i=0; i<N; i++) {
    switch (tag) {
      case 0: // Gauss
        if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use zero field formula
          Double_t sigma_t_2 = t*t*val[1]*val[1];
          p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*TMath::Exp(-0.5*sigma_t_2));
        } else if (val[1]/val[0] > 79.5775) { // check if Delta/w0 > 500.0, in which case the ZF formula is used
          Double_t sigma_t_2 = t*t*val[1]*val[1];
          p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*TMath::Exp(-0.5*sigma_t_2));
        } else {
          Double_t delta = val[1];
          Double_t w0    = TWO_PI*val[0];
 
          p0exp[i] = 1.0 - 2.0*TMath::Power(delta/w0,2.0)*(1.0 -
                    TMath::Exp(-0.5*TMath::Power(delta*t, 2.0))*TMath::Cos(w0*t)) +
                    GetLFIntegralValue(t);
        }
        break;
      case 1: // Lorentz
        if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use zero field formula
          Double_t at = t*val[1];
          p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - at)*TMath::Exp(-at));
        } else if (val[1]/val[0] > 159.1549) { // check if a/w0 > 1000.0, in which case the ZF formula is used
          Double_t at = t*val[1];
          p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - at)*TMath::Exp(-at));
        } else {
          Double_t a    = val[1];
          Double_t at   = a*t;
          Double_t w0   = TWO_PI*val[0];
          Double_t a_w0 = a/w0;
          Double_t w0t  = w0*t;
 
          Double_t j1, j0;
          if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
            j0 = 1.0;
            j1 = 0.0;
          } else {
            j0 = sin(w0t)/w0t;
            j1 = (sin(w0t)-w0t*cos(w0t))/(w0t*w0t);
          }
 
          p0exp[i] = 1.0 - a_w0*j1*exp(-at) - a_w0*a_w0*(j0*exp(-at) - 1.0) - GetLFIntegralValue(t);
        }
        break;
      default:
        break;
    }
    p0exp[i] *= TMath::Exp(-val[2]*t);
    t += dt;
  }
 
  // solve the volterra equation (trapezoid integration)
  fDynLFFuncValue[0]=p0exp[0];
 
  Double_t sum;
  Double_t a;
  Double_t preFactor = dt*val[2];
  for (UInt_t i=1; i<N; i++) {
    sum = p0exp[i];
    sum += 0.5*preFactor*p0exp[i]*fDynLFFuncValue[0];
    for (UInt_t j=1; j<i; j++) {
      sum += preFactor*p0exp[i-j]*fDynLFFuncValue[j];
    }
    a = 1.0-0.5*preFactor*p0exp[0];
 
    fDynLFFuncValue[i]=sum/a;
  }
 
  // clean up
  p0exp.clear();
}
 
//--------------------------------------------------------------------------
Double_t PTheory::GetDynKTLFValue(const Double_t t) const
{
  if (t < 0.0)
    return 1.0;
 
  UInt_t idx = static_cast<UInt_t>(t/fDynLFdt);
 
  if (idx + 2 > fDynLFFuncValue.size())
    return fDynLFFuncValue.back();
 
  // linearly interpolate between the two relevant function bins
  Double_t df = (fDynLFFuncValue[idx+1]-fDynLFFuncValue[idx])*(t/fDynLFdt-static_cast<Double_t>(idx));
 
  return fDynLFFuncValue[idx]+df;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::GetDyn_GL_KTLFValue(const Double_t t) const
{
  if (t < 0.0)
    return 1.0;
 
  const Double_t dt=0.001; // 1ns
  UInt_t idx = static_cast<UInt_t>(t/dt);
 
  if (idx + 2 > fDyn_GL_LFFuncValue.size())
    return fDyn_GL_LFFuncValue.back();
 
  // linearly interpolate between the two relevant function bins
  Double_t df = (fDyn_GL_LFFuncValue[idx+1]-fDyn_GL_LFFuncValue[idx])*(t/dt-static_cast<Double_t>(idx));
 
  return fDyn_GL_LFFuncValue[idx]+df;
}
 
//--------------------------------------------------------------------------
Double_t PTheory::MuMinusExpTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
  // expected parameters: N0 tau A lambda phase frequency [tshift]
 
  Double_t val[7];
 
  assert(fParamNo.size() <= 7);
 
  // check if FUNCTIONS are used
  for (UInt_t i=0; i<fParamNo.size(); i++) {
    if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
      val[i] = paramValues[fParamNo[i]];
    } else { // function
      val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
    }
  }
 
  Double_t tt;
  if (fParamNo.size() == 6) // no tshift
    tt = t;
  else // tshift present
    tt = t-val[6];
 
  return val[0]*exp(-tt/val[1])*(1.0+val[2]*exp(-val[3]*tt)*cos(TWO_PI*val[5]*tt+DEG_TO_RAD*val[4]));
}
 
//--------------------------------------------------------------------------
// END
//--------------------------------------------------------------------------