From e891f0489adf92c49f5f61c339c8f90f785f5679 Mon Sep 17 00:00:00 2001
From: Andreas Suter <andreas.suter@psi.ch>
Date: Fri, 6 Jun 2025 10:21:09 +0200
Subject: [PATCH] 1st full Gauss/Lorentz LF. Still room for optimization, and
 further testing.

---
 src/classes/PTheory.cpp | 172 ++++++++++++++++++++++++++++++++++++----
 src/include/PTheory.h   |   2 +
 2 files changed, 157 insertions(+), 17 deletions(-)

diff --git a/src/classes/PTheory.cpp b/src/classes/PTheory.cpp
index d1be50fd..3fd45c08 100644
--- a/src/classes/PTheory.cpp
+++ b/src/classes/PTheory.cpp
@@ -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;
 }
 
 //--------------------------------------------------------------------------
@@ -2619,22 +2732,22 @@ Double_t PTheory::FmuF(Double_t t, const PDoubleVector& paramValues, const PDoub
     assert(fParamNo.size() <= 2);
 
     if (t < 0.0)
-        return 1.0;
+      return 1.0;
 
     // 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];
-        }
+      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() == 1) // no tshift
-        tt = t;
+      tt = t;
     else // tshift present
-        tt = t-val[1];
+      tt = t-val[1];
 
     const Double_t sqrt3 = sqrt(3.0);
     const Double_t wd_t = val[0]*tt;
@@ -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
diff --git a/src/include/PTheory.h b/src/include/PTheory.h
index ddbccf9c..d873c62d 100644
--- a/src/include/PTheory.h
+++ b/src/include/PTheory.h
@@ -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_