added skewed Gaussian to PTheory, and other minor stuff

src/.kdbgrc.musrfit

@@ -1,18 +0,0 @@
-[Breakpoint 0]
-Enabled=true
-File=musrfit.cpp
-Line=464
-Temporary=false
-
-[General]
-DebuggerCmdStr=
-DriverName=GDB
-FileVersion=1
-OptionsSelected=
-ProgramArgs=
-TTYLevel=7
-WorkingDirectory=
-
-[Memory]
-ColumnWidths=80,0
-NumExprs=0

src/Makefile.musrfit

@@ -37,7 +37,8 @@ CXX           = g++
 CXXFLAGS      = -g -Wall -fPIC
 PMUSRPATH     = ./include
 MNPATH        = $(ROOTSYS)/include
-INCLUDES      = -I $(PMUSRPATH) -I $(MNPATH)
+GSLPATH       = /usr/include/gsl
+INCLUDES      = -I $(PMUSRPATH) -I $(MNPATH) -I $(GSLPATH)
 LD            = g++
 LDFLAGS       = -g
 endif
@@ -63,6 +64,8 @@ GLIBS         = $(ROOTGLIBS) -lXMLParser
 PSILIBS       = -lTLemRunHeader -lPMusr
 # Minuit2 lib
 MNLIB         = -L$(ROOTSYS)/lib -lMinuit2
+# GSL lib
+GSLLIB        = -lgslcblas -lgsl
 
 
 EXEC = musrfit
@@ -78,7 +81,7 @@ all:  $(EXEC)
 $(EXEC): $(OBJS)
 	@echo "---> Building $(EXEC) ..."
 	/bin/rm -f $(SHLIB)
-	$(LD) $(OBJS) -o $(EXEC) $(GLIBS) $(PSILIBS) $(MNLIB)
+	$(LD) $(OBJS) -o $(EXEC) $(GLIBS) $(PSILIBS) $(MNLIB) $(GSLLIB)
 	@echo "done"
 
 # clean up: remove all object file (and core files)

src/classes/Makefile.PMusr

@@ -43,7 +43,8 @@ CXX           = g++
 CXXFLAGS      = -g -Wall -fPIC
 PMUSRPATH     = ../include
 MNPATH        = $(ROOTSYS)/include
-INCLUDES      = -I $(PMUSRPATH) -I $(MNPATH)
+GSLPATH       = /usr/include/gsl
+INCLUDES      = -I $(PMUSRPATH) -I $(MNPATH) -I $(GSLPATH)
 LD            = g++
 LDFLAGS       = -g
 SOFLAGS       = -O -shared
@@ -71,6 +72,8 @@ GLIBS         = $(ROOTGLIBS) -lXMLParser
 PSILIBS       = -lTLemRunHeader
 # Minuit2 lib
 MNLIB         = -L$(ROOTSYS)/lib -lMinuit2
+# GSL lib
+GSLLIB        = -lgslcblas -lgsl
 
 
 # some definitions: headers, sources, objects,...
@@ -116,7 +119,7 @@ all:  $(SHLIB)
 $(SHLIB): $(OBJS)
 	@echo "---> Building shared library $(SHLIB) ..."
 	/bin/rm -f $(SHLIB)
-	$(LD) $(OBJS) $(SOFLAGS) -o $(SHLIB) $(GLIBS) $(PSILIBS) $(MNLIB)
+	$(LD) $(OBJS) $(SOFLAGS) -o $(SHLIB) $(GLIBS) $(PSILIBS) $(MNLIB) $(GSLLIB)
 	@echo "done"
 
 # clean up: remove all object file (and core files)

src/classes/PTheory.cpp

@@ -117,7 +117,7 @@ PTheory::PTheory(PMsrHandler *msrInfo, unsigned int runNo, const bool hasParent)
   TObjArray *tokens;
   TObjString *ostr;
 
-  tokens = str.Tokenize(" ");
+  tokens = str.Tokenize(" \t");
   if (!tokens) {
     cout << endl << "**SEVERE ERROR**: PTheory(): Couldn't tokenize theory block line " << line->fLineNo << ".";
     cout << endl << "  line content: " << line->fLine.Data();
@@ -288,6 +288,7 @@ bool PTheory::IsValid()
  */
 double PTheory::Func(register double t, const vector<double>& paramValues, const vector<double>& funcValues) const
 {
+
   if (fMul) {
     if (fAdd) { // fMul != 0 && fAdd != 0
       switch (fType) {
@@ -351,6 +352,10 @@ double PTheory::Func(register double t, const vector<double>& paramValues, const
           return InternalBessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
                  fAdd->Func(t, paramValues, funcValues);
           break;
+        case THEORY_SKEWED_GAUSS:
+          return SkewedGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
+                 fAdd->Func(t, paramValues, funcValues);
+          break;
         default:
           cout << endl << "**PANIC ERROR**: PTheory::Func: You never should have reached this line?!?! (" << fType << ")";
           cout << endl;
@@ -403,6 +408,9 @@ double PTheory::Func(register double t, const vector<double>& paramValues, const
         case THEORY_INTERNAL_BESSEL:
           return InternalBessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
           break;
+        case THEORY_SKEWED_GAUSS:
+          return SkewedGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
+          break;
         default:
           cout << endl << "**PANIC ERROR**: PTheory::Func: You never should have reached this line?!?! (" << fType << ")";
           cout << endl;
@@ -457,6 +465,9 @@ double PTheory::Func(register double t, const vector<double>& paramValues, const
         case THEORY_INTERNAL_BESSEL:
           return InternalBessel(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
           break;
+        case THEORY_SKEWED_GAUSS:
+          return SkewedGauss(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
+          break;
         default:
           cout << endl << "**PANIC ERROR**: PTheory::Func: You never should have reached this line?!?! (" << fType << ")";
           cout << endl;
@@ -509,6 +520,9 @@ double PTheory::Func(register double t, const vector<double>& paramValues, const
         case THEORY_INTERNAL_BESSEL:
           return InternalBessel(t, paramValues, funcValues);
           break;
+        case THEORY_SKEWED_GAUSS:
+          return SkewedGauss(t, paramValues, funcValues);
+          break;
         default:
           cout << endl << "**PANIC ERROR**: PTheory::Func: You never should have reached this line?!?! (" << fType << ")";
           cout << endl;
@@ -565,7 +579,7 @@ void PTheory::MakeCleanAndTidyTheoryBlock(PMsrLines *fullTheoryBlock)
     // copy line content to str in order to remove comments
     str = line->fLine.Copy();
     // tokenize line
-    tokens = str.Tokenize(" ");
+    tokens = str.Tokenize(" \t");
     // make a handable string out of the asymmetry token
     ostr = dynamic_cast<TObjString*>(tokens->At(0));
     str = ostr->GetString();
@@ -968,3 +982,43 @@ double PTheory::InternalBessel(register double t, const vector<double>& paramVal
          0.333333333333333*TMath::Exp(-val[3]*t);
 }
 
+//--------------------------------------------------------------------------
+/**
+ * <p>
+ *
+ * \param t time in \f$(\mu\mathrm{s})\f$
+ * \param paramValues
+ */
+double PTheory::SkewedGauss(register double t, const vector<double>& paramValues, const vector<double>& funcValues) const
+{
+  double val[4];
+  double skg;
+
+  // check if FUNCTIONS are used
+  for (unsigned int i=0; i<4; 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];
+    }
+  }
+
+
+  // val[2] = sigma-, val[3] = sigma+
+  double z1 = val[3]*t/TMath::Sqrt(2.0); // sigma+
+  double z2 = val[2]*t/TMath::Sqrt(2.0); // sigma-
+  double g1 = TMath::Exp(-0.5*TMath::Power(t*val[3], 2.0)); // gauss sigma+
+  double g2 = TMath::Exp(-0.5*TMath::Power(t*val[2], 2.0)); // gauss sigma-
+
+  if ((z1 >= 25.0) || (z2 >= 25.0)) // needed to prevent crash of 1F1
+    skg = 2.0e300;
+  else
+    skg = 0.5*TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*t) * ( g1 + g2 ) -
+          0.5*TMath::Sin(DEG_TO_RAD*val[0]+TWO_PI*val[1]*t) *
+          (
+          g1*2.0*z1/TMath::Sqrt(TMath::Pi())*gsl_sf_hyperg_1F1(0.5,1.5,z1*z1) -
+          g2*2.0*z2/TMath::Sqrt(TMath::Pi())*gsl_sf_hyperg_1F1(0.5,1.5,z2*z2)
+          );
+
+  return skg;
+}

src/include/PTheory.h

@@ -38,6 +38,12 @@
 
 #include "PMsrHandler.h"
 
+#include <gsl_sf_hyperg.h>
+
+extern "C" {
+  double gsl_sf_hyperg_1F1(double a, double b, double x);
+}
+
 // --------------------------------------------------------
 // function handling tags
 // --------------------------------------------------------
@@ -59,7 +65,8 @@
 #define THEORY_TF_COS                      12
 #define THEORY_BESSEL                      13
 #define THEORY_INTERNAL_BESSEL             14
-#define THEORY_USER                        15
+#define THEORY_SKEWED_GAUSS                15
+#define THEORY_USER                        16
 
 // function parameter tags, i.e. how many parameters has a specific function
 #define THEORY_PARAM_ASYMMETRY                    1 // asymetry
@@ -77,9 +84,10 @@
 #define THEORY_PARAM_TF_COS                       2 // phase, frequency
 #define THEORY_PARAM_BESSEL                       2 // phase, frequency
 #define THEORY_PARAM_INTERNAL_BESSEL              5 // fraction, phase, frequency, TF damping, damping
+#define THEORY_PARAM_SKEWED_GAUSS                 4 // phase, frequency, rate minus, rate plus
 
 // number of available user functions
-#define THEORY_MAX 15
+#define THEORY_MAX 16
 
 // deg -> rad factor
 #define DEG_TO_RAD 0.0174532925199432955
@@ -150,7 +158,10 @@ static PTheoDataBase fgTheoDataBase[THEORY_MAX] = {
          "bessel", "b", "(phase frequency)"},
 
         {THEORY_INTERNAL_BESSEL, THEORY_PARAM_INTERNAL_BESSEL, false,
-         "internBsl", "ib", "(phase frequency Trate Lrate)"}};
+         "internBsl", "ib", "(phase frequency Trate Lrate)"},
+
+        {THEORY_SKEWED_GAUSS, THEORY_PARAM_SKEWED_GAUSS, false,
+         "skewedGss", "skg", "(phase frequency rate_minus rate_plus)"}};
 
 //--------------------------------------------------------------------------------------
 /**
@@ -184,6 +195,7 @@ class PTheory
     virtual double TFCos(register double t, const vector<double>& paramValues, const vector<double>& funcValues) const;
     virtual double Bessel(register double t, const vector<double>& paramValues, const vector<double>& funcValues) const;
     virtual double InternalBessel(register double t, const vector<double>& paramValues, const vector<double>& funcValues) const;
+    virtual double SkewedGauss(register double t, const vector<double>& paramValues, const vector<double>& funcValues) const;
 
     // variables
     bool fValid;