From 0fc7f8c4bf63f3596efba05a2a2636fb88bb5b2c Mon Sep 17 00:00:00 2001 From: nemu Date: Fri, 27 Feb 2009 09:45:52 +0000 Subject: [PATCH] added dynamic Gauss KT --- src/classes/PFitter.cpp | 19 ++-- src/classes/PTheory.cpp | 193 +++++++++++++++++++++++++++++++- src/include/PTheory.h | 4 + src/tests/dynKT_LF/dynKT_LF.cpp | 4 +- 4 files changed, 207 insertions(+), 13 deletions(-) diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index 74823be3..eb6c7f46 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -65,17 +65,18 @@ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, bo { /* PRunData *data=runListCollection->GetSingleHisto(0); -fstream fout( "__test.dat", ios_base::out ); -fout << data->fDataTimeStart << endl; -fout << data->fDataTimeStep << endl; -fout << "------" << endl; -fout << data->fValue.size() << endl; -fout << "------" << endl; -for (unsigned int i=0; ifValue.size(); i++) { - fout << data->fValue[i] << ", " << data->fError[i] << endl; -} +fstream fout( "__test.dat", ios_base::out | ios_base::binary); +fout.write(reinterpret_cast(&data->fDataTimeStart),sizeof(double)); +fout.write(reinterpret_cast(&data->fDataTimeStep),sizeof(double)); +unsigned int ss=data->fValue.size(); +fout.write(reinterpret_cast(&ss),sizeof(ss)); +for (unsigned int i=0; i(&data->fValue[i]),sizeof(double)); +for (unsigned int i=0; i(&data->fError[i]),sizeof(double)); fout.close(); */ + fUseChi2 = true; // chi^2 is the default fParams = *(runInfo->GetMsrParamList()); diff --git a/src/classes/PTheory.cpp b/src/classes/PTheory.cpp index 03befa34..10d1a234 100644 --- a/src/classes/PTheory.cpp +++ b/src/classes/PTheory.cpp @@ -93,6 +93,7 @@ PTheory::PTheory(PMsrHandler *msrInfo, unsigned int runNo, const bool hasParent) fMul = 0; fUserFcnClassName = TString(""); fUserFcn = 0; + fDynLFdt = 0.0; static unsigned int lineNo = 1; // lineNo static unsigned int depth = 0; // needed to handle '+' properly @@ -298,6 +299,7 @@ PTheory::~PTheory() fUserParam.clear(); fLFIntegral.clear(); + fDynLFFuncValue.clear(); if (fMul) { delete fMul; @@ -1090,7 +1092,68 @@ double PTheory::StaticGaussKTLF(register double t, const PDoubleVector& paramVal */ double PTheory::DynamicGaussKTLF(register double t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const { - return 0.0; + // expected parameters: frequency damping hopping [tshift] + + double val[4]; + double result = 0.0; + bool useKeren = false; + + assert(fParamNo.size() <= 4); + + // check if FUNCTIONS are used + for (unsigned int i=0; i 5.0) // nu/Delta > 5.0 + useKeren = true; + + if (!useKeren) { + // check if the parameter values have changed, and if yes recalculate the non-analytic integral + bool newParam = false; + for (unsigned int i=0; i<4; i++) { + if (val[i] != fPrevParam[i]) { + newParam = true; + break; + } + } + + if (newParam) { // new parameters found + for (unsigned int i=0; i<4; i++) + fPrevParam[i] = val[i]; + CalculateDynKTLF(val, 0); // 0 means Gauss + } + } + + double tt; + if (fParamNo.size() == 3) // no tshift + tt = t; + else // tshift present + tt = t-val[3]; + + if (tt < 0.0) // for times < 0 return a function value of 0.0 + return 0.0; + + + if (useKeren) { // see PRB50, 10039 (1994) + double wL = TWO_PI * val[0]; + double wL2 = wL*wL; + double nu2 = val[2]*val[2]; + double Gamma_t = 2.0*val[1]*val[1]/((wL2+nu2)*(wL2+nu2))* + ((wL2+nu2)*val[2]*t + + (wL2-nu2)*(1.0 - TMath::Exp(-val[2]*t)*TMath::Cos(wL*t)) + - 2.0*val[2]*wL*TMath::Exp(-val[2]*t)*TMath::Sin(wL*t)); + result = TMath::Exp(-Gamma_t); + } else { // from Voltera + result = GetDynKTLFValue(tt); + } + + return result; } //-------------------------------------------------------------------------- @@ -1150,7 +1213,7 @@ double PTheory::StaticLorentzKTLF(register double t, const PDoubleVector& paramV double at = a*tt; double w0 = 2.0*TMath::Pi()*val[0]; double a_w0 = a/w0; - double w0t = w0*t; + double w0t = w0*tt; double j1, j0; if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x) @@ -1696,6 +1759,132 @@ double PTheory::GetLFIntegralValue(const double t) const return fLFIntegral[idx]+df; } +//-------------------------------------------------------------------------- +/** + *

Number of sampling points is estimated according to w0: w0 T = 2 pi + * t_sampling = T/16 (i.e. 16 points on 1 period) + * N = 8/pi w0 Tmax = 16 nu0 Tmax + * + * \param val parameters needed to solve the voltera integral equation + * \param tag 0=Gauss, 1=Lorentz + */ +void PTheory::CalculateDynKTLF(const double *val, int tag) const +{ + // val: 0=nu0, 1=Delta (Gauss) / a (Lorentz), 2=nu + const double Tmax = 20.0; // 20 usec + unsigned int N = static_cast(16.0*Tmax*val[0]); + + if (N < 300) // if too view points, i.e. nu0 very small, take 300 points + N = 300; + + // allocate memory for dyn KT LF function vector + fDynLFFuncValue.clear(); // get rid of a possible previous vector + fDynLFFuncValue.resize(N); + + // calculate the non-analytic integral of the static KT LF function + switch (tag) { + case 0: // Gauss + CalculateGaussLFIntegral(val); + break; + case 1: // Lorentz + CalculateLorentzLFIntegral(val); + break; + default: + cout << endl << "**FATAL ERROR** in PTheory::CalculateDynKTLF: You should never have reached this point." << endl; + assert(false); + break; + } + + // calculate the P^(0)(t) exp(-nu t) vector + PDoubleVector p0exp(N); + double t = 0.0; + double dt = Tmax/N; + fDynLFdt = dt; // keep it since it is needed in GetDynKTLFValue() + for (unsigned int i=0; i + * + * \param t time in (usec) + */ +double PTheory::GetDynKTLFValue(const double t) const +{ + unsigned int idx = static_cast(t/fDynLFdt); + + if (idx < 0) + return 0.0; + + if (idx > fDynLFFuncValue.size()-1) + return fDynLFFuncValue[fDynLFFuncValue.size()-1]; + + // liniarly interpolate between the two relvant function bins + double df = (fDynLFFuncValue[idx+1]-fDynLFFuncValue[idx])*(t-idx*fDynLFdt)/fDynLFdt; + + return fDynLFFuncValue[idx]+df; +} + //-------------------------------------------------------------------------- // END //-------------------------------------------------------------------------- diff --git a/src/include/PTheory.h b/src/include/PTheory.h index c380d579..55fc0065 100644 --- a/src/include/PTheory.h +++ b/src/include/PTheory.h @@ -223,6 +223,8 @@ class PTheory virtual void CalculateGaussLFIntegral(const double *val) const; virtual void CalculateLorentzLFIntegral(const double *val) const; virtual double GetLFIntegralValue(const double t) const; + virtual void CalculateDynKTLF(const double *val, int tag) const; + virtual double GetDynKTLFValue(const double t) const; // variables bool fValid; @@ -240,6 +242,8 @@ class PTheory mutable double fPrevParam[THEORY_MAX_PARAM]; ///< needed for LF-stuff mutable PDoubleVector fLFIntegral; ///< needed for LF-stuff. Keeps the non-analytic integral values + mutable double fDynLFdt; + mutable PDoubleVector fDynLFFuncValue; ///< needed for LF-stuff. Keeps the dynamic LF KT function values }; #endif // _PTHEORY_H_ diff --git a/src/tests/dynKT_LF/dynKT_LF.cpp b/src/tests/dynKT_LF/dynKT_LF.cpp index 8cbed311..b974873b 100644 --- a/src/tests/dynKT_LF/dynKT_LF.cpp +++ b/src/tests/dynKT_LF/dynKT_LF.cpp @@ -145,8 +145,8 @@ int main(int argc, char *argv[]) } else { // w0 criteria, i.e. w0 T = 2 pi, ts = T/16, N = Tmax/ts, if N < 300, N == 300 double val = 8.0/PI*Tmax*param[0]; - if (val < 250) - N = 250; + if (val < 300) + N = 300; else N = static_cast(val);