From c713a367d6733ee8db799c8021b5ae369ddf299e Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Thu, 31 May 2012 09:02:05 +0000 Subject: [PATCH] added Noakes-Kalvius function --- ChangeLog | 2 + doc/musrfit.dox | 15 +- doc/musrfit_dox.cfg | 4 +- src/classes/PRunAsymmetry.cpp | 2 +- src/classes/PTheory.cpp | 282 +++++++++++++++++++++++++++++++++- src/include/PTheory.h | 30 +++- 6 files changed, 317 insertions(+), 18 deletions(-) diff --git a/ChangeLog b/ChangeLog index cd8dc402..4b19022c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -6,6 +6,8 @@ changes since 0.11.0 =================================== +NEW 2012-05-31 added Noakes-Kalvius function (see A. Yaouanc and P. Dalmas de Reotiers, + "Muon Spin Rotation, Relaxation, and Resonance" Oxford, Section 6.4.1.3). NEW 2012-05-25 musredit/musrgui: added a dump muSR data header file information. NEW 2012-05-22 added spin rotation angle to the LEM MusrRoot file. NEW 2012-05-12 added dump_header. This is a little program which dumps the diff --git a/doc/musrfit.dox b/doc/musrfit.dox index 34325042..8b0af484 100644 --- a/doc/musrfit.dox +++ b/doc/musrfit.dox @@ -50,9 +50,8 @@ How to setup musrfit on different platforms: \texttt{http://lmu.web.psi.ch/facil on >= Qt4.6. - \ref MuSRFit A graphical user interface based on PerlQt (written by Z. Salman) for an easy to use interface to the musrfit framework. Compared to the more general approach of writting msr-files, it has some limitations, though it might be easier for a first user of the musrfit framework. - \ref any2many Should be a "universal" muSR data-file-format converter. -- \ref nexus_dump Is a small program to dump NeXus file information (mainly run header info) to the standard output. +- \ref dump_header Is a small program to dump the header information of a muSR data file to the standard output. - \ref musrRootValidation This is a program to validate MusrRoot files. -- \ref read_musrRoot_runHeader Is a small program to dump MusrRoot file information (mainly run header info) to the standard output. - \ref write_musrRoot_runHeader Is a little example program showing how to write MusrRoot files. \section roadmap Road map and missing features @@ -62,7 +61,6 @@ How to setup musrfit on different platforms: \texttt{http://lmu.web.psi.ch/facil

The following features should eventually be implemented, but are still missing: - there are still issues with MUD files on 64bit systems which should eventually be fixed. - non-muSR: The plan is to add an option to fit/plot \f$f(x_1,\ldots,x_n)\f$ versus \f$g(x_1,\ldots,x_n)\f$, where \f$x_i\f$ is a given data set element. -- mu-Minus: The handling for \f$\mu_{-}\f$ fits is still missing and should be implemented. - as soon as ROOT will properly support MS Windows platforms, some better support for MS Windows will be added. Currently only the cygwin version will be supported. - check if it is possible to add FIR filtering for muSR data - add an interface to maxent @@ -141,9 +139,9 @@ PSI firewall required).

Here will eventually follow a more technical description of any2many. If you looking for a user-manual like description, please check \htmlonlymusrfit user manual\endhtmlonly \latexonly musrfit user manual: \texttt{http://lmu.web.psi.ch/facilities/software/musrfit/user/MUSR/MusrFit.html} \endlatexonly //**************************************************************************************************** -\page nexusDumpPage -\section nexus_dump nexus_dump -

This is a little help program which reads a NeXus file and dumps most of the relevant information to the standard output. +\page dumpHeaderPage +\section dump_header dump_header +

This is a little helper program which reads the header information of a muSR data file and dumps most of the relevant information to the standard output. //**************************************************************************************************** \page musrRootValidationPage @@ -151,11 +149,6 @@ PSI firewall required).

This program allows to validate a given (or supposedly given) MusrRoot file if it is consistent with the minimal required entries via XML Schema schemes. Please check \htmlonlyMusrRoot web-page\endhtmlonly \latexonly MusrRoot web-page: \texttt{http://lmu.web.psi.ch/facilities/software/musrfit/user/MUSR/MusrRoot.html}\endlatexonly -//**************************************************************************************************** -\page read_musrRoot_runHeader_Page -\section read_musrRoot_runHeader read_musrRoot_runHeader -

This is a little help program which reads a MusrRoot file and dumps the RunHeader information to the standard output. - //**************************************************************************************************** \page write_musrRoot_runHeader_Page \section write_musrRoot_runHeader write_musrRoot_runHeader diff --git a/doc/musrfit_dox.cfg b/doc/musrfit_dox.cfg index 30fbbfe5..292239bf 100644 --- a/doc/musrfit_dox.cfg +++ b/doc/musrfit_dox.cfg @@ -509,14 +509,14 @@ INPUT = musrfit.dox \ ../src/classes/PUserFcn.cpp \ ../src/external/MusrRoot/TMusrRunHeader.h \ ../src/external/MusrRoot/TMusrRunHeader.cpp \ + ../src/any2many.cpp \ + ../src/dump_header.cpp \ ../src/msr2data.cpp \ ../src/msr2msr.cpp \ ../src/musrfit.cpp \ ../src/musrt0.cpp \ ../src/musrview.cpp \ - ../src/nexus_dump.cpp \ ../src/musrRootValidation.cpp \ - ../src/read_musrRoot_runHeader.cpp \ ../src/write_musrRoot_runHeader.cpp # If the value of the INPUT tag contains directories, you can use the diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp index 6b301e90..9b9babc0 100644 --- a/src/classes/PRunAsymmetry.cpp +++ b/src/classes/PRunAsymmetry.cpp @@ -379,7 +379,7 @@ Bool_t PRunAsymmetry::PrepareData() // keep the time resolution in (us) fTimeResolution = runData->GetTimeResolution()/1.0e3; cout.precision(10); - cout << endl << ">> PRunSingleHisto::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ns)" << endl; + cout << endl << ">> PRunAsymmetry::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ns)" << endl; // collect histogram numbers PUIntVector forwardHistoNo; diff --git a/src/classes/PTheory.cpp b/src/classes/PTheory.cpp index a18296d9..97f2f4c8 100644 --- a/src/classes/PTheory.cpp +++ b/src/classes/PTheory.cpp @@ -393,7 +393,6 @@ Bool_t PTheory::IsValid() */ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const { - if (fMul) { if (fAdd) { // fMul != 0 && fAdd != 0 switch (fType) { @@ -477,6 +476,22 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co return SkewedGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues); break; + case THEORY_STATIC_ZF_NK: + return StaticNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) + + fAdd->Func(t, paramValues, funcValues); + break; + case THEORY_STATIC_TF_NK: + return StaticNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) + + fAdd->Func(t, paramValues, funcValues); + break; + case THEORY_DYNAMIC_ZF_NK: + return DynamicNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) + + fAdd->Func(t, paramValues, funcValues); + break; + case THEORY_DYNAMIC_TF_NK: + return DynamicNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) + + fAdd->Func(t, paramValues, funcValues); + break; case THEORY_POLYNOM: return Polynom(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues); @@ -552,6 +567,18 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co case THEORY_SKEWED_GAUSS: return SkewedGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues); break; + case THEORY_STATIC_ZF_NK: + return StaticNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues); + break; + case THEORY_STATIC_TF_NK: + return StaticNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues); + break; + case THEORY_DYNAMIC_ZF_NK: + return DynamicNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues); + break; + case THEORY_DYNAMIC_TF_NK: + return DynamicNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues); + break; case THEORY_POLYNOM: return Polynom(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues); break; @@ -627,6 +654,18 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co case THEORY_SKEWED_GAUSS: return SkewedGauss(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues); break; + case THEORY_STATIC_ZF_NK: + return StaticNKZF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues); + break; + case THEORY_STATIC_TF_NK: + return StaticNKTF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues); + break; + case THEORY_DYNAMIC_ZF_NK: + return DynamicNKZF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues); + break; + case THEORY_DYNAMIC_TF_NK: + return DynamicNKTF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues); + break; case THEORY_POLYNOM: return Polynom(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues); break; @@ -700,6 +739,18 @@ Double_t PTheory::Func(register Double_t t, const PDoubleVector& paramValues, co case THEORY_SKEWED_GAUSS: return SkewedGauss(t, paramValues, funcValues); break; + case THEORY_STATIC_ZF_NK: + return StaticNKZF(t, paramValues, funcValues); + break; + case THEORY_STATIC_TF_NK: + return StaticNKTF(t, paramValues, funcValues); + break; + case THEORY_DYNAMIC_ZF_NK: + return DynamicNKZF(t, paramValues, funcValues); + break; + case THEORY_DYNAMIC_TF_NK: + return DynamicNKTF(t, paramValues, funcValues); + break; case THEORY_POLYNOM: return Polynom(t, paramValues, funcValues); break; @@ -2063,6 +2114,235 @@ Double_t PTheory::SkewedGauss(register Double_t t, const PDoubleVector& paramVal return skg; } +//-------------------------------------------------------------------------- +/** + *

theory function: staticNKZF (see D.R. Noakes and G.M. Kalvius Phys. Rev. B 56, 2352 (1997) and + * A. Yaouanc and P. Dalmas de Reotiers, "Muon Spin Rotation, Relaxation, and Resonance" Oxford, Section 6.4.1.3) + * + * \f[ = \frac{1}{3} + \frac{2}{3}\,\frac{1}{\left(1+(\gamma\Delta_{\rm GbG}t)^2\right)^{3/2}}\, + * \left(1 - \frac{(\gamma\Delta_0 t)^2}{\left(1+(\gamma\Delta_{\rm GbG}t)^2\right)}\right)\, + * \exp\left[\frac{(\gamma\Delta_0 t)^2}{2\left(1+(\gamma\Delta_{\rm GbG}t)^2\right)}\right] \f] + * + * meaning of paramValues: \f$\Delta_0\f$, \f$R_{\rm b} = \Delta_{\rm GbG}/\Delta_0\f$ [,\f$t_{\rm shift}\f$] + * + * return: function value + * + * \param t time in \f$(\mu\mathrm{s})\f$, or x-axis value for non-muSR fit + * \param paramValues parameter values + * \param funcValues vector with the functions (i.e. functions of the parameters) + */ +Double_t PTheory::StaticNKZF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const +{ + // expected paramters: damping_D0 R_b tshift + + 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 theory function: staticNKTF (see D.R. Noakes and G.M. Kalvius Phys. Rev. B 56, 2352 (1997) and + * A. Yaouanc and P. Dalmas de Reotiers, "Muon Spin Rotation, Relaxation, and Resonance" Oxford, Section 6.4.1.3) + * + * \f[ = \frac{1}{\sqrt{1+(\gamma\Delta_{\rm GbG} t)^2}}\, + * \exp\left[-\frac{(\gamma\Delta_0 t)^2}{2(1+(\gamma\Delta_{\rm GbG}t)^2)}\right]\, + * \cos(\gamma B_{\rm ext} t + \varphi) \f] + * + * meaning of paramValues: \f$\varphi\f$, \f$\nu = \gamma B_{\rm ext}\f$, \f$\Delta_0\f$, \f$R_{\rm b} = \Delta_{\rm GbG}/\Delta_0\f$ [,\f$t_{\rm shift}\f$] + * + * return: function value + * + * \param t time in \f$(\mu\mathrm{s})\f$, or x-axis value for non-muSR fit + * \param paramValues parameter values + * \param funcValues vector with the functions (i.e. functions of the parameters) + */ +Double_t PTheory::StaticNKTF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const +{ + // expected paramters: phase frequency damping_D0 R_b tshift + + 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 theory function: dynamicNKZF (see D.R. Noakes and G.M. Kalvius Phys. Rev. B 56, 2352 (1997) and + * A. Yaouanc and P. Dalmas de Reotiers, "Muon Spin Rotation, Relaxation, and Resonance" Oxford, Section 6.4.1.3) + * + * \f{eqnarray*} + * \Theta(t) &=& \frac{\exp(-\nu_c t) - 1 - \nu_c t}{\nu_c^2} \\ + * \Delta_{\rm eff} &=& \sqrt{\Delta_0^2 + \Delta_{\rm GbG}^2} \\ + * P_{Z}^{\rm dyn}(t) &=& \sqrt{\frac{1+R_{\rm b}^2}{1+R_{\rm b}^2+4 (R_{\rm b}\gamma\Delta_{\rm eff})^2 \Theta(t)}}\, + * \exp\left[-\frac{2 (\gamma\Delta_{\rm eff})^2\Theta(t)}{1+R_{\rm b}^2+4 (R_{\rm b}\gamma\Delta_{\rm eff})^2 \Theta(t)}\right] \\ + * &=& \sqrt{\frac{1}{1+4 \Delta_{\rm GbG}^2 \Theta(t)}}\,\exp\left[-\frac{2 \Delta_0^2 \Theta(t)}{1+4 \Delta_{\rm GbG}^2 \Theta(t)}\right] + * \f} + * + * meaning of paramValues: \f$\Delta_0\f$, \f$R_{\rm b} = \Delta_{\rm GbG}/\Delta_0\f$, \f$\nu_c\f$ [,\f$t_{\rm shift}\f$] + * + * return: function value + * + * \param t time in \f$(\mu\mathrm{s})\f$, or x-axis value for non-muSR fit + * \param paramValues parameter values + * \param funcValues vector with the functions (i.e. functions of the parameters) + */ +Double_t PTheory::DynamicNKZF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const +{ + // expected paramters: damping_D0 R_b nu_c tshift + + 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 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/(1.0 + 4.0*val[0]*val[0]*val[1]*val[1]*theta); + + result = sqrt(denom)*exp(-2.0*val[0]*val[0]*theta*denom); + + return result; +} + +//-------------------------------------------------------------------------- +/** + *

theory function: dynamicNKTF (see D.R. Noakes and G.M. Kalvius Phys. Rev. B 56, 2352 (1997) and + * A. Yaouanc and P. Dalmas de Reotiers, "Muon Spin Rotation, Relaxation, and Resonance" Oxford, Section 6.4.1.3) + * + * \f{eqnarray*} + * \Theta(t) &=& \frac{\exp(-\nu_c t) - 1 - \nu_c t}{\nu_c^2} \\ + * \Delta_{\rm eff} &=& \sqrt{\Delta_0^2 + \Delta_{\rm GbG}^2} \\ + * P_{X}^{\rm dyn}(t) &=& \sqrt{\frac{1+R_{\rm b}^2}{1+R_{\rm b}^2+2 (R_{\rm b}\gamma\Delta_{\rm eff})^2 \Theta(t)}}\, + * \exp\left[-\frac{(\gamma\Delta_{\rm eff})^2\Theta(t)}{1+R_{\rm b}^2+2 (R_{\rm b}\gamma\Delta_{\rm eff})^2 \Theta(t)}\right]\,\cos(\gamma B_{\rm ext} t + \varphi) \\ + * &=& \sqrt{\frac{1}{1+2 \Delta_{\rm GbG}^2 \Theta(t)}}\,\exp\left[-\frac{\Delta_0^2 \Theta(t)}{1+2 \Delta_{\rm GbG}^2 \Theta(t)}\right]\,\cos(\gamma B_{\rm ext} t + \varphi) + * \f} + * + * meaning of paramValues: \f$\varphi\f$, \f$\nu = \gamma B_{\rm ext}\f$, \f$\Delta_0\f$, \f$R_{\rm b} = \Delta_{\rm GbG}/\Delta_0\f$, \f$\nu_c\f$ [,\f$t_{\rm shift}\f$] + * + * return: function value + * + * \param t time in \f$(\mu\mathrm{s})\f$, or x-axis value for non-muSR fit + * \param paramValues parameter values + * \param funcValues vector with the functions (i.e. functions of the parameters) + */ +Double_t PTheory::DynamicNKTF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const +{ + // expected paramters: phase frequency damping_D0 R_b nu_c tshift + + 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 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/(1.0 + 2.0*val[2]*val[2]*val[3]*val[3]*theta); + + result = sqrt(denom)*exp(-val[2]*val[2]*theta*denom)*TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt); + + return result; +} + //-------------------------------------------------------------------------- /** *

theory function: polynom diff --git a/src/include/PTheory.h b/src/include/PTheory.h index 5b2b8e03..94f7aae3 100644 --- a/src/include/PTheory.h +++ b/src/include/PTheory.h @@ -65,8 +65,12 @@ #define THEORY_BESSEL 17 #define THEORY_INTERNAL_BESSEL 18 #define THEORY_SKEWED_GAUSS 19 -#define THEORY_POLYNOM 20 -#define THEORY_USER_FCN 21 +#define THEORY_STATIC_ZF_NK 20 +#define THEORY_STATIC_TF_NK 21 +#define THEORY_DYNAMIC_ZF_NK 22 +#define THEORY_DYNAMIC_TF_NK 23 +#define THEORY_POLYNOM 24 +#define THEORY_USER_FCN 25 // function parameter tags, i.e. how many parameters has a specific function // if there is a comment with a (tshift), the number of parameters is increased by one @@ -90,9 +94,13 @@ #define THEORY_PARAM_BESSEL 2 // phase, frequency (tshift) #define THEORY_PARAM_INTERNAL_BESSEL 5 // fraction, phase, frequency, TF damping, LF damping (tshift) #define THEORY_PARAM_SKEWED_GAUSS 4 // phase, frequency, rate minus, rate plus (tshift) +#define THEORY_PARAM_STATIC_ZF_NK 2 // damping D0, R_b=DGbG/D0 (tshift) +#define THEORY_PARAM_STATIC_TF_NK 4 // phase, frequency, damping D0, R_b=DGbG/D0 (tshift) +#define THEORY_PARAM_DYNAMIC_ZF_NK 3 // damping D0, R_b=DGbG/D0, nu_c (tshift) +#define THEORY_PARAM_DYNAMIC_TF_NK 5 // phase, frequency, damping D0, R_b=DGbG/D0, nu_c (tshift) // number of available user functions -#define THEORY_MAX 22 +#define THEORY_MAX 26 // maximal number of parameters. Needed in the contents of LF #define THEORY_MAX_PARAM 10 @@ -184,6 +192,18 @@ static PTheoDataBase fgTheoDataBase[THEORY_MAX] = { {THEORY_SKEWED_GAUSS, THEORY_PARAM_SKEWED_GAUSS, false, "skewedGss", "skg", "(phase frequency rate_m rate_p)", "(phase frequency rate_m rate_p tshift)"}, + {THEORY_STATIC_ZF_NK, THEORY_PARAM_STATIC_ZF_NK, false, + "staticNKZF", "snkzf", "(damping_D0 R_b)", "(damping_D0 R_b tshift)"}, + + {THEORY_STATIC_TF_NK, THEORY_PARAM_STATIC_TF_NK, false, + "staticNKTF", "snktf", "(phase frequency damping_D0 R_b)", "(phase frequency damping_D0 R_b tshift)"}, + + {THEORY_DYNAMIC_ZF_NK, THEORY_PARAM_DYNAMIC_ZF_NK, false, + "dynamicNKZF", "dnkzf", "(damping_D0 R_b nu_c)", "(damping_D0 R_b nu_c tshift)"}, + + {THEORY_DYNAMIC_TF_NK, THEORY_PARAM_DYNAMIC_TF_NK, false, + "dynamicNKTF", "dnktf", "(phase frequency damping_D0 R_b nu_c)", "(phase frequency damping_D0 R_b nu_c tshift)"}, + {THEORY_POLYNOM, 0, false, "polynom", "p", "(tshift p0 p1 ... pn)", "(tshift p0 p1 ... pn)"}, @@ -232,6 +252,10 @@ class PTheory virtual Double_t Bessel(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; virtual Double_t InternalBessel(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; virtual Double_t SkewedGauss(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; + virtual Double_t StaticNKZF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; + virtual Double_t StaticNKTF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; + virtual Double_t DynamicNKZF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; + virtual Double_t DynamicNKTF(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; virtual Double_t Polynom(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; virtual Double_t UserFcn(register Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;