implemented static Gauss and Lorentz KT LF functions. Removed some unnecessary output from PMusrCanvas

src/classes/PMusrCanvas.cpp

@@ -266,8 +266,6 @@ void PMusrCanvas::UpdateParamTheoryPad()
   else
     yoffset = 0.05;
 
-cout << endl << ">> yoffset parameter = " << yoffset;
-
   // add parameters to the pad
   for (unsigned int i=0; i<param.size(); i++) {
     str = "";
@@ -317,7 +315,6 @@ cout << endl << ">> yoffset parameter = " << yoffset;
       str += cnum;
     }
     ypos = 0.98-i*yoffset;
-cout << endl << ">> ypos = " << ypos;
     fParameterPad->AddText(0.03, ypos, str.Data());
   }
 
@@ -327,7 +324,6 @@ cout << endl << ">> ypos = " << ypos;
     yoffset = 1.0/(theory.size()+1);
   else
     yoffset = 0.05;
-cout << endl << ">> yoffset theory = " << yoffset << endl;
   for (unsigned int i=1; i<theory.size(); i++) {
     // remove comment if present
     str = theory[i].fLine;
@@ -2156,7 +2152,7 @@ void PMusrCanvas::PlotData()
       }
       // plot all the theory
       for (unsigned int i=0; i<fData.size(); i++) {
-        fData[i].theory->Draw("csame");
+        fData[i].theory->Draw("lsame");
       }
     }
   } else { // fPlotType == MSR_PLOT_NO_MUSR

src/classes/PTheory.cpp

@@ -91,7 +91,6 @@ PTheory::PTheory(PMsrHandler *msrInfo, unsigned int runNo, const bool hasParent)
   fValid = true;
   fAdd = 0;
   fMul = 0;
-  fStaticKTLFFunc = 0;
   fUserFcnClassName = TString("");
   fUserFcn = 0;
 
@@ -103,6 +102,9 @@ PTheory::PTheory(PMsrHandler *msrInfo, unsigned int runNo, const bool hasParent)
     depth  = 0; // reset for every root object (new run).
   }
 
+  for (unsigned int i=0; i<THEORY_MAX_PARAM; i++)
+    fPrevParam[i] = 0.0;
+
   // get the input to be analyzed from the msr handler
   PMsrLines *fullTheoryBlock = msrInfo->GetMsrTheory();
   if (lineNo > fullTheoryBlock->size()-1) {
@@ -149,22 +151,6 @@ PTheory::PTheory(PMsrHandler *msrInfo, unsigned int runNo, const bool hasParent)
     return;
   }
 
-  // for static KT LF some TF1 object needs to be invoked 
-  if (idx == (unsigned int) THEORY_STATIC_GAUSS_KT_LF) {
-    // check if it is necessary to create fStaticKTLFFunc
-    if (fStaticKTLFFunc == 0) {
-      fStaticKTLFFunc = new TF1("fStaticKTLFFunc",
-                                "exp(-0.5*pow([0]*x,2.0))*sin(2.0*pi*[1]*x)",
-                                0.0, 13.0);
-      if (!fStaticKTLFFunc) {
-        cout << endl << "**WARNING** in StaticKTLF: Couldn't invoke necessary function";
-        fValid = false;
-      } else {
-        fStaticKTLFFunc->SetNpx(1000);
-      }
-    }
-  }
-
   // line is a valid function, hence analyze parameters
   if (((unsigned int)(tokens->GetEntries()-1) < fNoOfParam) &&
       ((idx != THEORY_USER_FCN) && (idx != THEORY_POLYNOM))) {
@@ -311,6 +297,8 @@ PTheory::~PTheory()
   fParamNo.clear();
   fUserParam.clear();
 
+  fLFIntegral.clear();
+
   if (fMul) {
     delete fMul;
     fMul = 0;
@@ -321,11 +309,6 @@ PTheory::~PTheory()
     fAdd = 0;
   }
 
-  if (fStaticKTLFFunc) {
-    delete fStaticKTLFFunc;
-    fStaticKTLFFunc = 0;
-  }
-
   if (fUserFcn) {
     delete fUserFcn;
     fUserFcn = 0;
@@ -1059,25 +1042,40 @@ double PTheory::StaticGaussKTLF(register double t, const PDoubleVector& paramVal
     }
   }
 
+  // check if the parameter values have changed, and if yes recalculate the non-analytic integral
+  bool newParam = false;
+  for (unsigned int i=0; i<3; i++) {
+    if (val[i] != fPrevParam[i]) {
+      newParam = true;
+      break;
+    }
+  }
+
+  if (newParam) { // new parameters found
+    for (unsigned int i=0; i<3; i++)
+      fPrevParam[i] = val[i];
+    CalculateGaussLFIntegral(val);
+  }
+
   double tt;
   if (fParamNo.size() == 2) // no tshift
     tt = t;
   else // tshift present
     tt = t-val[2];
 
-  if (val[0] == 0.0) {
+  if (tt < 0.0) // for times < 0 return a function value of 0.0
+    return 0.0;
+
+  if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use zero field formula
     double 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 {
-    fStaticKTLFFunc->SetParameters(val[1], val[0]); // damping, frequency
-
     double delta = val[1];
     double 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)) +
-                2.0*TMath::Power(delta, 4.0)/TMath::Power(w0, 3.0)*
+                GetLFIntegralValue(tt);
-                fStaticKTLFFunc->Integral(0.0, tt);
   }
 
   return result;
@@ -1104,7 +1102,69 @@ double PTheory::DynamicGaussKTLF(register double t, const PDoubleVector& paramVa
  */
 double PTheory::StaticLorentzKTLF(register double t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
 {
-  return 0.0;
+  // expected parameters: frequency damping [tshift]
+
+  double val[3];
+  double result;
+
+  assert(fParamNo.size() <= 3);
+
+  // check if FUNCTIONS are used
+  for (unsigned int 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 the parameter values have changed, and if yes recalculate the non-analytic integral
+  bool newParam = false;
+  for (unsigned int i=0; i<3; i++) {
+    if (val[i] != fPrevParam[i]) {
+      newParam = true;
+      break;
+    }
+  }
+
+  if (newParam) { // new parameters found
+    for (unsigned int i=0; i<3; i++)
+      fPrevParam[i] = val[i];
+    CalculateLorentzLFIntegral(val);
+  }
+
+  double 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 0.0
+   return 0.0;
+
+ if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use zero field formula
+    double at = tt*val[1];
+    result = 0.333333333333333 * (1.0 + 2.0*(1.0 - at)*TMath::Exp(-at));
+  } else {
+    double a    = val[1];
+    double at   = a*tt;
+    double w0   = 2.0*TMath::Pi()*val[0];
+    double a_w0 = a/w0;
+    double w0t  = w0*t;
+
+    double 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;
 }
 
 //--------------------------------------------------------------------------
@@ -1539,3 +1599,103 @@ if (first) {
 */
   return (*fUserFcn)(t, fUserParam);
 }
+
+//--------------------------------------------------------------------------
+/**
+ * <p>The fLFIntegral is given in steps of 1 ns, e.g. index i=100 coresponds to t=100ns
+ *
+ * \param val parameters needed to calculate the non-analytic integral of the static Gauss LF function.
+ */
+void PTheory::CalculateGaussLFIntegral(const double *val) const
+{
+  // val[0] = nu, val[1] = Delta
+
+  if (val[0] == 0.0) // field == 0.0, hence nothing to be done
+    return;
+
+  const double dt=0.001; // all times in usec
+  double t, ft;
+  double w0 = TMath::TwoPi()*val[0];
+  double Delta = val[1];
+  double preFactor = 2.0*TMath::Power(Delta, 4.0) / TMath::Power(w0, 3.0);
+
+  // clear previously allocated vector
+  fLFIntegral.clear();
+
+  // calculate integral
+  t  = 0.0;
+  fLFIntegral.push_back(0.0); // start value of the integral
+
+  ft = 0.0;
+  for (unsigned int i=1; i<2e4; i++) {
+    t  += dt;
+    ft += 0.5*dt*preFactor*(TMath::Exp(-0.5*TMath::Power(Delta * (t-dt), 2.0))*TMath::Sin(w0*(t-dt))+
+                           TMath::Exp(-0.5*TMath::Power(Delta * t, 2.0))*TMath::Sin(w0*t));
+    fLFIntegral.push_back(ft);
+  }
+}
+
+//--------------------------------------------------------------------------
+/**
+ * <p>
+ *
+ * \param val parameters needed to calculate the non-analytic integral of the static Lorentz LF function.
+ */
+void PTheory::CalculateLorentzLFIntegral(const double *val) const
+{
+  // val[0] = nu, val[1] = a
+
+  if (val[0] == 0.0) // field == 0.0, hence nothing to be done
+    return;
+
+  const double dt=0.001; // all times in usec
+  double t, ft;
+  double w0 = TMath::TwoPi()*val[0];
+  double a = val[1];
+  double preFactor = a*(1+pow(a/w0,2.0));
+
+  // 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
+  for (unsigned int i=2; i<2e4; i++) {
+    t  += dt;
+    ft += 0.5*dt*preFactor*(sin(w0*(t-dt))/(w0*(t-dt))*exp(-a*(t-dt))+sin(w0*t)/(w0*t)*exp(-a*t));
+    fLFIntegral.push_back(ft);
+  }
+}
+
+//--------------------------------------------------------------------------
+/**
+ * <p>
+ *
+ * \param t time in (usec)
+ */
+double PTheory::GetLFIntegralValue(const double t) const
+{
+  unsigned int idx = static_cast<unsigned int>(t/0.001); // since LF-integral is calculated in nsec
+
+  if (idx < 0)
+    return 0.0;
+
+  if (idx > fLFIntegral.size()-1)
+    return fLFIntegral[fLFIntegral.size()-1];
+
+  // liniarly interpolate between the two relvant function bins
+  double df = (fLFIntegral[idx+1]-fLFIntegral[idx])/0.001*(t-idx*0.001);
+
+  return fLFIntegral[idx]+df;
+}
+
+//--------------------------------------------------------------------------
+// END
+//--------------------------------------------------------------------------

src/include/PTheory.h

@@ -34,18 +34,11 @@
 
 #include <TSystem.h>
 #include <TString.h>
-#include <TF1.h>
 
 #include "PMusr.h"
 #include "PMsrHandler.h"
 #include "PUserFcnBase.h"
 
-// #include <gsl_sf_hyperg.h>
-// 
-// extern "C" {
-//   double gsl_sf_hyperg_1F1(double a, double b, double x);
-// }
-
 // --------------------------------------------------------
 // function handling tags
 // --------------------------------------------------------
@@ -97,6 +90,9 @@
 // number of available user functions
 #define THEORY_MAX 20
 
+// maximal number of parameters. Needed in the contents of LF
+#define THEORY_MAX_PARAM 10
+
 // deg -> rad factor
 #define DEG_TO_RAD 0.0174532925199432955
 // 2 pi
@@ -224,13 +220,16 @@ class PTheory
     virtual double Polynom(register double t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
     virtual double UserFcn(register double t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
 
+    virtual void CalculateGaussLFIntegral(const double *val) const;
+    virtual void CalculateLorentzLFIntegral(const double *val) const;
+    virtual double GetLFIntegralValue(const double t) const;
+
     // variables
     bool fValid;
     unsigned int fType;
     vector<unsigned int> fParamNo; ///< holds the parameter numbers for the theory (including maps and functions, see constructor desciption)
     unsigned int fNoOfParam;
     PTheory *fAdd, *fMul;
-    TF1 *fStaticKTLFFunc;
 
     TString fUserFcnClassName; ///< name of the user function class for within root
     TString fUserFcnSharedLibName; ///< name of the shared lib to which the user function belongs
@@ -238,6 +237,9 @@ class PTheory
     mutable PDoubleVector fUserParam;  ///< vector holding the resolved user function parameters, i.e. map and function resolved.
 
     PMsrHandler *fMsrInfo;
+
+    mutable double fPrevParam[THEORY_MAX_PARAM]; ///< needed for LF-stuff
+    mutable PDoubleVector fLFIntegral;           ///< needed for LF-stuff. Keeps the non-analytic integral values
 };
 
 #endif //  _PTHEORY_H_