LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 4 Activity

added Noakes-Kalvius function

Browse Source
This commit is contained in:
suter_a
2012-05-31 09:02:05 +00:00
parent e6e4bc19da
commit c713a367d6
6 changed files with 317 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files

282
src/classes/PTheory.cpp
View File

@ -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;
}
//--------------------------------------------------------------------------
/**
* <p> 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]
*
* <b>meaning of paramValues:</b> \f$\Delta_0\f$, \f$R_{\rm b} = \Delta_{\rm GbG}/\Delta_0\f$ [,\f$t_{\rm shift}\f$]
*
* <b>return:</b> 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<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 t2 = tt*tt;
Double_t Rb2 = val[1]*val[1];
Double_t Rb2p = 1.0+Rb2;
Double_t Deff2_t2 = val[0]*val[0]*(1.0+Rb2)*t2;
Double_t denom = (Rb2p+Rb2*Deff2_t2);
result = 0.333333333333333 + 0.666666666666666667 * TMath::Power(Rb2p/denom, 1.5) * (1.0 - (Deff2_t2/denom)) * exp(-0.5*Deff2_t2/denom);
return result;
}
//--------------------------------------------------------------------------
/**
* <p> 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]
*
* <b>meaning of paramValues:</b> \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$]
*
* <b>return:</b> 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<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 D0t_2 = val[2]*val[2]*tt*tt;
Double_t DGt_2p = 1.0 + val[2]*val[2]*val[3]*val[3]*tt*tt;
result = 1.0/sqrt(DGt_2p)*exp(-0.5*D0t_2/DGt_2p)*TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt);
return result;
}
//--------------------------------------------------------------------------
/**
* <p> 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}
*
* <b>meaning of paramValues:</b> \f$\Delta_0\f$, \f$R_{\rm b} = \Delta_{\rm GbG}/\Delta_0\f$, \f$\nu_c\f$ [,\f$t_{\rm shift}\f$]
*
* <b>return:</b> 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<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/(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;
}
//--------------------------------------------------------------------------
/**
* <p> 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}
*
* <b>meaning of paramValues:</b> \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$]
*
* <b>return:</b> 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<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/(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;
}
//--------------------------------------------------------------------------
/**
* <p> theory function: polynom