1st full Gauss/Lorentz LF. Still room for optimization, and further testing.
This commit is contained in:
parent
66af890157
commit
ecec59c8a3
@ -29,6 +29,7 @@
|
||||
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
|
||||
#include <TObject.h>
|
||||
#include <TString.h>
|
||||
@ -1759,16 +1760,128 @@ Double_t PTheory::DynamicGauLorKTLFFast(Double_t t, const PDoubleVector& paramVa
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* @brief PTheory::DynamicGauLorKTLF
|
||||
* @param t
|
||||
* @param paramValues
|
||||
* @param funcValues
|
||||
* @return
|
||||
* <p>Local Gaussian, global Lorentzian in LF.
|
||||
* For details see "Muon Spin Rotation, Relaxation, and Resonance",
|
||||
* A. Yaouanc and P. Dalmas Sec. 6.4, Eq.(6.86).
|
||||
*
|
||||
* @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)
|
||||
*
|
||||
* @return Polarization value of this function
|
||||
*/
|
||||
Double_t PTheory::DynamicGauLorKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
|
||||
{
|
||||
// NOT YET IMPLEMENTED. Will be Eq.(6.86)
|
||||
return 0.0;
|
||||
// expected parameters: frequency damping hopping [tshift]
|
||||
|
||||
Double_t val[4];
|
||||
|
||||
assert(fParamNo.size() <= 4);
|
||||
|
||||
// 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];
|
||||
|
||||
// check if the parameter values have changed, and if yes recalculate DynamicGaussKTLF
|
||||
Bool_t newParam = false;
|
||||
for (UInt_t i=0; i<3; i++) {
|
||||
if (val[i] != fPrevParam[i]) {
|
||||
newParam = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (newParam) { // new parameters found, hence calculate DynamicGauLorKTLF
|
||||
// keep parameters
|
||||
for (UInt_t i=0; i<3; i++)
|
||||
fPrevParam[i] = val[i];
|
||||
|
||||
// reset GL LF polarzation vector
|
||||
fDyn_GL_LFFuncValue.clear();
|
||||
fDyn_GL_LFFuncValue.resize(20000); // Tmax=20us, dt=1ns
|
||||
|
||||
PDoubleVector rr={0.2, 0.4, 0.6, 0.8, 1.0, 1.25, 1.5, 1.75, 2.0, 2.5, 3.0, 4.0, 5.0, 7.5, 10.0,
|
||||
12.8125, 15.625, 18.4375, 21.25, 26.875, 32.5, 43.75, 55.0, 77.5, 100.0};
|
||||
Double_t par[3] = {val[0], val[1], val[2]};
|
||||
Double_t sqrtTwoInv = 1.0/sqrt(2.0);
|
||||
Bool_t isOneVec{false};
|
||||
Bool_t useKeren{false};
|
||||
Double_t scale, up{0.0}, low{-1.0};
|
||||
|
||||
for (UInt_t i=0; i<rr.size(); i++) {
|
||||
useKeren = false;
|
||||
isOneVec = false;
|
||||
|
||||
// Delta_G = rr * Delta_L
|
||||
par[1] = rr[i] * val[1];
|
||||
|
||||
// check if all parameters == 0
|
||||
if ((par[0] == 0.0) && (par[1] == 0.0) && (par[2] == 0.0)) {
|
||||
isOneVec = true;
|
||||
}
|
||||
|
||||
// make sure that damping and hopping are positive definite
|
||||
if (par[1] < 0.0)
|
||||
par[1] = -par[1];
|
||||
if (par[2] < 0.0)
|
||||
par[2] = -par[2];
|
||||
|
||||
// check that Delta != 0, if not (i.e. stupid parameter) return 1, which is the correct limit
|
||||
if (fabs(par[1]) < 1.0e-6) {
|
||||
isOneVec = true;
|
||||
}
|
||||
|
||||
// check if Keren approximation can be used
|
||||
if (par[2]/par[1] > 5.0) // nu/Delta > 5.0
|
||||
useKeren = true;
|
||||
|
||||
if (!useKeren && !isOneVec) {
|
||||
CalculateDynKTLF(par, 0); // 0 means Gauss
|
||||
}
|
||||
|
||||
// calculate polarization vector for the given parameters
|
||||
up = -std::erf(sqrtTwoInv/rr[i]);
|
||||
scale = up - low;
|
||||
low = up;
|
||||
|
||||
const Double_t dt=0.001;
|
||||
for (UInt_t i=0; i<20000; i++) {
|
||||
if (isOneVec) {
|
||||
fDyn_GL_LFFuncValue[i] += scale;
|
||||
} else if (useKeren && !isOneVec) {// see PRB50, 10039 (1994)
|
||||
Double_t wL = TWO_PI * par[0];
|
||||
Double_t wL2 = wL*wL;
|
||||
Double_t nu2 = par[2]*par[2];
|
||||
Double_t Gamma_t = 2.0*par[1]*par[1]/((wL2+nu2)*(wL2+nu2))*
|
||||
((wL2+nu2)*val[2]*i*dt
|
||||
+ (wL2-nu2)*(1.0 - TMath::Exp(-val[2]*i*dt)*TMath::Cos(wL*i*dt))
|
||||
- 2.0*val[2]*wL*TMath::Exp(-val[2]*i*dt)*TMath::Sin(wL*i*dt));
|
||||
fDyn_GL_LFFuncValue[i] += scale*TMath::Exp(-Gamma_t);
|
||||
} else if (!useKeren && !isOneVec) {
|
||||
fDyn_GL_LFFuncValue[i] += scale*GetDynKTLFValue(i*dt);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// get the proper value from the look-up table
|
||||
Double_t result{1.0};
|
||||
if (tt>=0)
|
||||
result=GetDyn_GL_KTLFValue(tt);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
@ -3024,19 +3137,44 @@ void PTheory::CalculateDynKTLF(const Double_t *val, Int_t tag) const
|
||||
Double_t PTheory::GetDynKTLFValue(const Double_t t) const
|
||||
{
|
||||
if (t < 0.0)
|
||||
return 0.0;
|
||||
return 1.0;
|
||||
|
||||
UInt_t idx = static_cast<UInt_t>(t/fDynLFdt);
|
||||
|
||||
if (idx + 2 > fDynLFFuncValue.size())
|
||||
return fDynLFFuncValue.back();
|
||||
|
||||
// linearly interpolate between the two relvant function bins
|
||||
// linearly interpolate between the two relevant function bins
|
||||
Double_t df = (fDynLFFuncValue[idx+1]-fDynLFFuncValue[idx])*(t/fDynLFdt-static_cast<Double_t>(idx));
|
||||
|
||||
return fDynLFFuncValue[idx]+df;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Gets value of the dynamic local Gauss / global Lorentzian Kubo-Toyabe LF function.
|
||||
*
|
||||
* <b>return:</b> function value
|
||||
*
|
||||
* \param t time in (usec)
|
||||
*/
|
||||
Double_t PTheory::GetDyn_GL_KTLFValue(const Double_t t) const
|
||||
{
|
||||
if (t < 0.0)
|
||||
return 1.0;
|
||||
|
||||
const Double_t dt=0.001; // 1ns
|
||||
UInt_t idx = static_cast<UInt_t>(t/dt);
|
||||
|
||||
if (idx + 2 > fDyn_GL_LFFuncValue.size())
|
||||
return fDyn_GL_LFFuncValue.back();
|
||||
|
||||
// linearly interpolate between the two relevant function bins
|
||||
Double_t df = (fDyn_GL_LFFuncValue[idx+1]-fDyn_GL_LFFuncValue[idx])*(t/dt-static_cast<Double_t>(idx));
|
||||
|
||||
return fDyn_GL_LFFuncValue[idx]+df;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p> theory function: MuMinusExpTF
|
||||
|
||||
@ -310,6 +310,7 @@ class PTheory
|
||||
virtual Double_t GetLFIntegralValue(const Double_t t) const;
|
||||
virtual void CalculateDynKTLF(const Double_t *val, Int_t tag) const;
|
||||
virtual Double_t GetDynKTLFValue(const Double_t t) const;
|
||||
virtual Double_t GetDyn_GL_KTLFValue(const Double_t t) const;
|
||||
|
||||
// variables
|
||||
Bool_t fValid; ///< flag to tell if the theory is valid
|
||||
@ -331,6 +332,7 @@ class PTheory
|
||||
mutable PDoubleVector fLFIntegral; ///< needed for LF. Keeps the non-analytic integral values
|
||||
mutable Double_t fDynLFdt; ///< needed for LF. Keeps the time step for the dynamic LF calculation, used in the integral equation approach
|
||||
mutable PDoubleVector fDynLFFuncValue; ///< needed for LF. Keeps the dynamic LF KT function values
|
||||
mutable PDoubleVector fDyn_GL_LFFuncValue; ///< needed for LF. Keeps the dynamic LF local Gauss/global Lorentzian KT function values
|
||||
};
|
||||
|
||||
#endif // _PTHEORY_H_
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user