diff --git a/src/.kdbgrc.musrfit b/src/.kdbgrc.musrfit deleted file mode 100755 index 5b22ad34..00000000 --- a/src/.kdbgrc.musrfit +++ /dev/null @@ -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 diff --git a/src/Makefile.musrfit b/src/Makefile.musrfit index 89cd546c..33492a91 100644 --- a/src/Makefile.musrfit +++ b/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) diff --git a/src/classes/Makefile.PMusr b/src/classes/Makefile.PMusr index 3e23593b..30a54bbc 100644 --- a/src/classes/Makefile.PMusr +++ b/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) diff --git a/src/classes/PTheory.cpp b/src/classes/PTheory.cpp index 35b20167..9fefa68d 100644 --- a/src/classes/PTheory.cpp +++ b/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& paramValues, const vector& funcValues) const { + if (fMul) { if (fAdd) { // fMul != 0 && fAdd != 0 switch (fType) { @@ -351,6 +352,10 @@ double PTheory::Func(register double t, const vector& 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& 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& 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& 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(tokens->At(0)); str = ostr->GetString(); @@ -968,3 +982,43 @@ double PTheory::InternalBessel(register double t, const vector& paramVal 0.333333333333333*TMath::Exp(-val[3]*t); } +//-------------------------------------------------------------------------- +/** + *

+ * + * \param t time in \f$(\mu\mathrm{s})\f$ + * \param paramValues + */ +double PTheory::SkewedGauss(register double t, const vector& paramValues, const vector& 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; +} diff --git a/src/include/PTheory.h b/src/include/PTheory.h index ae42a6dc..be70d9d3 100644 --- a/src/include/PTheory.h +++ b/src/include/PTheory.h @@ -38,6 +38,12 @@ #include "PMsrHandler.h" +#include + +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& paramValues, const vector& funcValues) const; virtual double Bessel(register double t, const vector& paramValues, const vector& funcValues) const; virtual double InternalBessel(register double t, const vector& paramValues, const vector& funcValues) const; + virtual double SkewedGauss(register double t, const vector& paramValues, const vector& funcValues) const; // variables bool fValid;