/*************************************************************************** PTheory.h Author: Andreas Suter e-mail: andreas.suter@psi.ch ***************************************************************************/ /*************************************************************************** * Copyright (C) 2007-2025 by Andreas Suter * * andreas.suter@psi.ch * * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the * * Free Software Foundation, Inc., * * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * ***************************************************************************/ #ifndef _PTHEORY_H_ #define _PTHEORY_H_ #include #include #include "PMusr.h" #include "PMsrHandler.h" #include "PUserFcnBase.h" //------------------------------------------------------------- /** *

Theory function type tags. * *

These constants identify the built-in theory functions available in * the THEORY block of MSR files. Each function represents a specific physical * model for muon spin relaxation, precession, or depolarization. * *

Theory functions can be combined using addition (+) and multiplication (*): * - Addition: Independent relaxation channels (e.g., 1/3 fast + 2/3 slow) * - Multiplication: Multiple relaxation mechanisms (e.g., precession * damping) * *

Categories: * - Basic: Constants, asymmetries, simple exponentials * - Static relaxation: Gaussian, Lorentzian (Kubo-Toyabe) * - Dynamic relaxation: Motional narrowing, spin fluctuations * - Precession: Cosine, Bessel functions for oscillations * - Vortex lattice: Internal field distributions in superconductors * - Special: Spin glass, Abragam, mu-minus, polynomials, user functions */ /// Undefined or invalid theory function #define THEORY_UNDEFINED -1 /// Constant value (baseline, background) #define THEORY_CONST 0 /// Initial asymmetry (multiplicative factor) #define THEORY_ASYMMETRY 1 /// Simple exponential relaxation: exp(-λt) #define THEORY_SIMPLE_EXP 2 /// General exponential relaxation: exp(-(λt)^β) #define THEORY_GENERAL_EXP 3 /// Simple Gaussian relaxation: exp(-σ²t²/2) #define THEORY_SIMPLE_GAUSS 4 /// Static Gaussian Kubo-Toyabe (zero-field) #define THEORY_STATIC_GAUSS_KT 5 /// Static Gaussian Kubo-Toyabe in longitudinal field #define THEORY_STATIC_GAUSS_KT_LF 6 /// Dynamic Gaussian Kubo-Toyabe in longitudinal field #define THEORY_DYNAMIC_GAUSS_KT_LF 7 /// Static Lorentzian Kubo-Toyabe (zero-field) #define THEORY_STATIC_LORENTZ_KT 8 /// Static Lorentzian Kubo-Toyabe in longitudinal field #define THEORY_STATIC_LORENTZ_KT_LF 9 /// Dynamic Lorentzian Kubo-Toyabe in longitudinal field #define THEORY_DYNAMIC_LORENTZ_KT_LF 10 /// Fast dynamic Gauss-Lorentz Kubo-Toyabe (zero-field) #define THEORY_DYNAMIC_GAULOR_FAST_KT_ZF 11 /// Fast dynamic Gauss-Lorentz Kubo-Toyabe in longitudinal field #define THEORY_DYNAMIC_GAULOR_FAST_KT_LF 12 /// Dynamic Gauss-Lorentz Kubo-Toyabe in longitudinal field #define THEORY_DYNAMIC_GAULOR_KT_LF 13 /// Combined Lorentzian-Gaussian Kubo-Toyabe #define THEORY_COMBI_LGKT 14 /// Stretched Kubo-Toyabe relaxation #define THEORY_STR_KT 15 /// Spin glass order parameter function #define THEORY_SPIN_GLASS 16 /// Random anisotropic hyperfine coupling #define THEORY_RANDOM_ANISOTROPIC_HYPERFINE 17 /// Abragam relaxation function (diffusion) #define THEORY_ABRAGAM 18 /// Transverse field cosine precession #define THEORY_TF_COS 19 /// Internal magnetic field distribution (superconductors) #define THEORY_INTERNAL_FIELD 20 /// Internal field (Kornilov vortex lattice model) #define THEORY_INTERNAL_FIELD_KORNILOV 21 /// Internal field (Larkin-Ovchinnikov model) #define THEORY_INTERNAL_FIELD_LARKIN 22 /// Bessel function (modulated precession) #define THEORY_BESSEL 23 /// Internal Bessel (field distribution with Bessel) #define THEORY_INTERNAL_BESSEL 24 /// Skewed Gaussian relaxation (asymmetric rates) #define THEORY_SKEWED_GAUSS 25 /// Static Nakajima zero-field function #define THEORY_STATIC_ZF_NK 26 /// Static Nakajima transverse field function #define THEORY_STATIC_TF_NK 27 /// Dynamic Nakajima zero-field function #define THEORY_DYNAMIC_ZF_NK 28 /// Dynamic Nakajima transverse field function #define THEORY_DYNAMIC_TF_NK 29 /// F-μ-F (μ-fluorine) oscillation #define THEORY_F_MU_F 30 /// Negative muon (μ-) exponential TF decay #define THEORY_MU_MINUS_EXP 31 /// Polynomial function (arbitrary order) #define THEORY_POLYNOM 32 /// User-defined external function (shared library) #define THEORY_USER_FCN 33 //------------------------------------------------------------- /** *

Number of parameters for each theory function. * *

These constants define how many parameters each theory function requires, * excluding optional time shift. If a function includes time shift, * add 1 to the parameter count. * *

Parameters typically include: rates (λ, σ), frequencies (ω, ν), * phases (φ), fractions, exponents (β), hopping rates (ν_hop), etc. */ #define THEORY_PARAM_CONST 1 // const #define THEORY_PARAM_ASYMMETRY 1 // asymmetry #define THEORY_PARAM_SIMPLE_EXP 1 // damping (tshift) #define THEORY_PARAM_GENERAL_EXP 2 // damping, exponents (tshift) #define THEORY_PARAM_SIMPLE_GAUSS 1 // damping (tshift) #define THEORY_PARAM_STATIC_GAUSS_KT 1 // damping (tshift) #define THEORY_PARAM_STATIC_GAUSS_KT_LF 2 // frequency, damping (tshift) #define THEORY_PARAM_DYNAMIC_GAUSS_KT_LF 3 // frequency, damping, hop-rate (tshift) #define THEORY_PARAM_STATIC_LORENTZ_KT 1 // damping (tshift) #define THEORY_PARAM_STATIC_LORENTZ_KT_LF 2 // frequency, damping (tshift) #define THEORY_PARAM_DYNAMIC_LORENTZ_KT_LF 3 // frequency, damping, hop-rate (tshift) #define THEORY_PARAM_DYNAMIC_GAULOR_FAST_KT_ZF 2 // damping, hop-rate (tshift) #define THEORY_PARAM_DYNAMIC_GAULOR_FAST_KT_LF 3 // frequency, damping, hop-rate (tshift) #define THEORY_PARAM_DYNAMIC_GAULOR_KT_LF 3 // frequency, damping, hop-rate (tshift) #define THEORY_PARAM_COMBI_LGKT 2 // Lorentz rate, Gauss rate (tshift) #define THEORY_PARAM_STR_KT 2 // rate, exponent (tshift) #define THEORY_PARAM_SPIN_GLASS 3 // rate, hop-rate, order parameter (tshift) #define THEORY_PARAM_RANDOM_ANISOTROPIC_HYPERFINE 2 // frequency, rate (tshift) #define THEORY_PARAM_ABRAGAM 2 // rate, hop-rate (tshift) #define THEORY_PARAM_TF_COS 2 // phase, frequency (tshift) #define THEORY_PARAM_INTERNAL_FIELD 5 // fraction, phase, frequency, TF damping, damping (tshift) #define THEORY_PARAM_INTERNAL_FIELD_KORNILOV 5 // fraction, frequency, TF damping, damping, beta (tshift) #define THEORY_PARAM_INTERNAL_FIELD_LARKIN 4 // fraction, frequency, TF damping, damping (tshift) #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) #define THEORY_PARAM_F_MU_F 1 // frequency (tshift) #define THEORY_PARAM_MU_MINUS_EXP 6 // N0, tau, A, damping, phase, frequency (tshift) //------------------------------------------------------------- /** *

Maximum number of theory functions in database. * *

This is the total number of built-in theory functions available, * including all relaxation, precession, and special functions. */ #define THEORY_MAX 34 //------------------------------------------------------------- /** *

Maximum number of parameters for any theory function. * *

Used to allocate arrays for longitudinal field calculations where * parameter values from previous iterations must be cached. */ #define THEORY_MAX_PARAM 10 //------------------------------------------------------------- /** *

Conversion factor from degrees to radians. * *

Value: π/180 = 0.017453292519943295 *

Used for phase parameters which are specified in degrees in MSR files * but converted to radians for calculations. */ #define DEG_TO_RAD 0.0174532925199432955 /** *

Mathematical constant 2π. * *

Used extensively in frequency-to-angular-frequency conversions: * ω = 2πν */ #define TWO_PI 6.28318530717958623 class PTheory; //-------------------------------------------------------------------------------------- /** *

Database entry for a theory function definition. * *

This structure stores metadata about each built-in theory function, * including its identifier, parameter count, name, abbreviation, and * help text. The database is used for: * - Parsing THEORY block entries in MSR files * - Validating parameter counts * - Generating help text and documentation * - Writing theory blocks with correct syntax */ typedef struct theo_data_base { UInt_t fType; ///< Theory function type tag (THEORY_CONST, THEORY_SIMPLE_EXP, etc.) UInt_t fNoOfParam; ///< Number of parameters (excluding optional time shift) Bool_t fTable; ///< True if function requires pre-calculated lookup table (e.g., LF Kubo-Toyabe) TString fName; ///< Full function name for MSR files (e.g., "simplExpo", "statGssKt") TString fAbbrev; ///< Short abbreviation (e.g., "se", "stg") TString fComment; ///< Parameter list shown as help text in MSR file TString fCommentTimeShift; ///< Parameter list with time shift included } PTheoDataBase; //-------------------------------------------------------------------------------------- /** *

Holds the functions available for the user. */ static PTheoDataBase fgTheoDataBase[THEORY_MAX] = { {THEORY_CONST, THEORY_PARAM_CONST, false, "const", "c", "", ""}, {THEORY_ASYMMETRY, THEORY_PARAM_ASYMMETRY, false, "asymmetry", "a", "", ""}, {THEORY_SIMPLE_EXP, THEORY_PARAM_SIMPLE_EXP, false, "simplExpo", "se", "(rate)", "(rate tshift)"}, {THEORY_GENERAL_EXP, THEORY_PARAM_GENERAL_EXP, false, "generExpo", "ge", "(rate exponent)", "(rate exponent tshift)"}, {THEORY_SIMPLE_GAUSS, THEORY_PARAM_SIMPLE_GAUSS, false, "simpleGss", "sg", "(rate)", "(rate tshift)"}, {THEORY_STATIC_GAUSS_KT, THEORY_PARAM_STATIC_GAUSS_KT, false, "statGssKt", "stg", "(rate)", "(rate tshift)"}, {THEORY_STATIC_GAUSS_KT_LF, THEORY_PARAM_STATIC_GAUSS_KT_LF, true, "statGssKTLF", "sgktlf", "(frequency damping)", "(frequency damping tshift)"}, {THEORY_DYNAMIC_GAUSS_KT_LF, THEORY_PARAM_DYNAMIC_GAUSS_KT_LF, true, "dynGssKTLF", "dgktlf", "(frequency damping hopping-rate)", "(frequency damping hopping-rate tshift)"}, {THEORY_STATIC_LORENTZ_KT, THEORY_PARAM_STATIC_LORENTZ_KT, true, "statExpKT", "sekt", "(rate)", "(rate tshift)"}, {THEORY_STATIC_LORENTZ_KT_LF, THEORY_PARAM_STATIC_LORENTZ_KT_LF, true, "statExpKTLF", "sektlf", "(frequency damping)", "(frequency damping tshift)"}, {THEORY_DYNAMIC_LORENTZ_KT_LF, THEORY_PARAM_DYNAMIC_LORENTZ_KT_LF, true, "dynExpKTLF", "dektlf", "(frequency damping hopping-rate)", "(frequency damping hopping-rate tshift)"}, {THEORY_DYNAMIC_GAULOR_FAST_KT_ZF, THEORY_PARAM_DYNAMIC_GAULOR_FAST_KT_ZF, true, "dynGLKT_F_ZF", "dglktfzf", "(damping hopping-rate)", "(damping hopping-rate tshift)"}, {THEORY_DYNAMIC_GAULOR_FAST_KT_LF, THEORY_PARAM_DYNAMIC_GAULOR_FAST_KT_LF, true, "dynGLKT_F_LF", "dglktflf", "(frequency damping hopping-rate)", "(frequency damping hopping-rate tshift)"}, {THEORY_DYNAMIC_GAULOR_KT_LF, THEORY_PARAM_DYNAMIC_GAULOR_KT_LF, true, "dynGLKT_LF", "dglktlf", "(frequency damping hopping-rate)", "(frequency damping hopping-rate tshift)"}, {THEORY_COMBI_LGKT, THEORY_PARAM_COMBI_LGKT, false, "combiLGKT", "lgkt", "(lorentzRate gaussRate)", "(lorentzRate gaussRate tshift)"}, {THEORY_STR_KT, THEORY_PARAM_STR_KT, false, "strKT", "skt", "(rate beta)", "(rate beta tshift)"}, {THEORY_SPIN_GLASS, THEORY_PARAM_SPIN_GLASS, false, "spinGlass", "spg", "(rate hopprate order)", "(rate hopprate order tshift)"}, {THEORY_RANDOM_ANISOTROPIC_HYPERFINE, THEORY_PARAM_RANDOM_ANISOTROPIC_HYPERFINE, false, "rdAnisoHf", "rahf", "(frequency rate)", "(frequency rate tshift)"}, {THEORY_ABRAGAM, THEORY_PARAM_ABRAGAM, false, "abragam", "ab", "(rate hopprate)", "(rate hopprate tshift)"}, {THEORY_TF_COS, THEORY_PARAM_TF_COS, false, "TFieldCos", "tf", "(phase frequency)", "(phase frequency tshift)"}, {THEORY_INTERNAL_FIELD, THEORY_PARAM_INTERNAL_FIELD, false, "internFld", "ifld", "(fraction phase frequency Trate Lrate)", "(fraction phase frequency Trate Lrate tshift)"}, {THEORY_INTERNAL_FIELD_KORNILOV, THEORY_PARAM_INTERNAL_FIELD_KORNILOV, false, "internFldGK", "ifgk", "(fraction frequency sigma lambda beta)", "(fraction frequency sigma lambda beta tshift)"}, {THEORY_INTERNAL_FIELD_LARKIN, THEORY_PARAM_INTERNAL_FIELD_LARKIN, false, "internFldLL", "ifll", "(fraction frequency sigma lambda beta)", "(fraction frequency sigma lambda beta tshift)"}, {THEORY_BESSEL, THEORY_PARAM_BESSEL, false, "bessel", "b", "(phase frequency)", "(phase frequency tshift)"}, {THEORY_INTERNAL_BESSEL, THEORY_PARAM_INTERNAL_BESSEL, false, "internBsl", "ib", "(fraction phase frequency Trate Lrate)", "(fraction phase frequency Trate Lrate tshift)"}, {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_F_MU_F, THEORY_PARAM_F_MU_F, false, "F_mu_F", "fmuf", "(frequency)", "(frequency tshift)"}, {THEORY_MU_MINUS_EXP, THEORY_PARAM_MU_MINUS_EXP, false, "muMinusExpTF", "mmsetf", "(N0 tau A lambda phase nu)", "(N0 tau A lambda phase nu tshift)"}, {THEORY_POLYNOM, 0, false, "polynom", "p", "(tshift p0 p1 ... pn)", "(tshift p0 p1 ... pn)"}, {THEORY_USER_FCN, 0, false, "userFcn", "u", "", ""} }; //-------------------------------------------------------------------------------------- /** * \brief Theory function evaluator and expression tree manager. * * PTheory is the core class responsible for parsing, validating, and evaluating * theory functions specified in the THEORY block of MSR files. It implements * a binary expression tree to handle complex combinations of theory functions. * * \section theory_overview Overview * * The theory describes the expected muon polarization as a function of time. * Different physical phenomena (relaxation, precession, field distributions) * are represented by different theory functions that can be combined. * * \section theory_tree Expression Tree Structure * * Theory expressions are parsed into a binary tree where: * - Each node represents a theory function * - Left child (fAdd) represents addition (+) * - Right child (fMul) represents multiplication (*) * * Example MSR theory block: * \verbatim * THEORY * a 1 # asymmetry * tf 2 3 # TFieldCos * se 4 # simpleExp * + * a 5 * tf 6 7 * \endverbatim * * Becomes: (par1 * cos(φ₂ + 2πν₃t) * exp(-λ₄t)) + (par5 * cos(φ₆ + 2πν₇t)) * * \section theory_syntax Syntax Details * * Function specification: * \code * function_name param1 param2 ... paramN (optional comment) * \endcode * * Parameter types: * - Direct number: \c 1, \c 2, \c 3 → parameter index (1-based) * - Map reference: \c map1, \c map2 → indirection via RUN block map * - Function reference: \c fun1, \c fun2 → evaluated FUNCTIONS block entry * * Operators: * - \c + on separate line: Addition of following terms * - Consecutive lines without +: Implicit multiplication * * \section theory_functions Available Functions * * Basic: * - \c const (c): Constant value * - \c asymmetry (a): Initial asymmetry * - \c simplExpo (se): exp(-λt) * - \c generExpo (ge): exp(-(λt)^β) * - \c simpleGss (sg): exp(-σ²t²/2) * * Kubo-Toyabe (Gaussian): * - \c statGssKt (stg): Static ZF Gaussian KT * - \c statGssKTLF (sgktlf): Static LF Gaussian KT * - \c dynGssKTLF (dgktlf): Dynamic LF Gaussian KT * * Kubo-Toyabe (Lorentzian): * - \c statExpKT (sekt): Static ZF Lorentzian KT * - \c statExpKTLF (sektlf): Static LF Lorentzian KT * - \c dynExpKTLF (dektlf): Dynamic LF Lorentzian KT * * Precession: * - \c TFieldCos (tf): cos(φ + 2πνt) * - \c bessel (b): J₀(2πνt + φ) * - \c internFld (ifld): Internal field distribution * * Special: * - \c spinGlass (spg): Spin glass relaxation * - \c abragam (ab): Motional narrowing * - \c userFcn (u): User-defined external function * * \section theory_lf Longitudinal Field Calculations * * LF Kubo-Toyabe functions require numerical integration which is cached * for efficiency. The cache is invalidated when parameters change. * Caching variables: * - fPrevParam: Previous parameter values for change detection * - fLFIntegral: Cached integral values * - fDynLFFuncValue: Dynamic KT function cache * * \section theory_user User Functions * * External functions can be loaded from shared libraries: * \code * userFcn libMyFunctions.so MyFunctionClass param1 param2 ... * \endcode * * User functions must inherit from PUserFcnBase and implement: * - operator()(Double_t t, const std::vector& par) * * \see PUserFcnBase, PMsrHandler, fgTheoDataBase */ class PTheory { public: /** * \brief Constructor that parses the THEORY block and builds the expression tree. * * Parses the theory block from the MSR file, validates function names and * parameter counts, resolves parameter references, and recursively builds * the expression tree for operators (+ and *). * * Parsing algorithm: * -# Get current line from theory block * -# Remove comments (text after '(' or '#') * -# Tokenize to extract function name and parameters * -# Look up function in fgTheoDataBase * -# Validate parameter count * -# Resolve parameter references (direct, map, fun) * -# If next line is not '+', recursively create fMul child * -# If next line is '+', skip it and recursively create fAdd child * -# For user functions: load shared library and instantiate class * * \param msrInfo Pointer to MSR file handler containing theory block * \param runNo Run number (0-based) for parameter map resolution * \param hasParent False for root theory, true for child nodes. * Controls static counter reset for recursive parsing. * * \warning If parsing fails, fValid is set to false and error messages * are output to stderr. Always check IsValid() after construction. */ PTheory(PMsrHandler *msrInfo, UInt_t runNo, const Bool_t hasParent = false); /** * \brief Destructor that recursively cleans up the expression tree. * * Releases all allocated resources: * - Clears parameter and function value vectors * - Clears LF integral caches * - Recursively deletes child theory nodes (fAdd, fMul) * - Deletes user function object if present * - Clears global user function vector */ virtual ~PTheory(); /** * \brief Checks if the entire theory expression tree is valid. * * Recursively validates all nodes in the tree. A theory is valid only * if this node and all its children (fAdd and fMul) are valid. * * \return true if all nodes are valid, false if any node has errors * * \note Call this after construction to verify the theory was parsed * successfully before attempting evaluation. */ virtual Bool_t IsValid(); /** * \brief Evaluates the theory function at a given time point. * * Recursively evaluates the expression tree at time t using the provided * parameter values. The evaluation follows the tree structure: * - If both fMul and fAdd exist: this_func * fMul + fAdd * - If only fMul exists: this_func * fMul * - If only fAdd exists: this_func + fAdd * - If neither exists: this_func * * \param t Time in microseconds (μs) for μSR fits, or x-value for non-μSR * \param paramValues Vector of current fit parameter values (FITPARAMETER block) * \param funcValues Vector of evaluated function values (FUNCTIONS block) * * \return Evaluated theory value (polarization or asymmetry) at time t * * \note For thread-safety, LF calculations cache results which may need * recalculation if parameters change between calls. */ virtual Double_t Func(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; private: /** \brief Recursively deletes child theory nodes (fAdd and fMul). */ virtual void CleanUp(PTheory *theo); /** \brief Searches fgTheoDataBase for a function by name or abbreviation. * \param name Function name (e.g., "simplExpo") or abbreviation (e.g., "se") * \return Index in fgTheoDataBase, or THEORY_UNDEFINED if not found */ virtual Int_t SearchDataBase(TString name); /** \brief Returns the index of user functions up to the given line. * \param lineNo Current line number in theory block * \return Count of userFcn entries before this line */ virtual Int_t GetUserFcnIdx(UInt_t lineNo) const; /** \brief Reformats the theory block for clean MSR file output. */ virtual void MakeCleanAndTidyTheoryBlock(PMsrLines* fullTheoryBlock); /** \brief Formats a polynomial theory line with proper spacing. */ virtual void MakeCleanAndTidyPolynom(UInt_t i, PMsrLines* fullTheoryBlock); /** \brief Formats a user function theory line with proper spacing. */ virtual void MakeCleanAndTidyUserFcn(UInt_t i, PMsrLines* fullTheoryBlock); // -------------------- Theory Function Implementations -------------------- // Each function evaluates its specific physical model at time t. // Parameters are resolved from fParamNo using paramValues and funcValues. /** \brief Returns constant value. Formula: c */ virtual Double_t Constant(const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Returns asymmetry value. Formula: A */ virtual Double_t Asymmetry(const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Simple exponential relaxation. Formula: exp(-λt) */ virtual Double_t SimpleExp(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief General (stretched) exponential. Formula: exp(-(λt)^β) */ virtual Double_t GeneralExp(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Simple Gaussian relaxation. Formula: exp(-σ²t²/2) */ virtual Double_t SimpleGauss(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Static Gaussian Kubo-Toyabe (ZF). Formula: 1/3 + 2/3(1-σ²t²)exp(-σ²t²/2) */ virtual Double_t StaticGaussKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Static Gaussian Kubo-Toyabe (LF). Requires numerical integration. */ virtual Double_t StaticGaussKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Dynamic Gaussian Kubo-Toyabe (LF). Strong collision model. */ virtual Double_t DynamicGaussKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Static Lorentzian Kubo-Toyabe (ZF). Formula: 1/3 + 2/3(1-at)exp(-at) */ virtual Double_t StaticLorentzKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Static Lorentzian Kubo-Toyabe (LF). Requires numerical integration. */ virtual Double_t StaticLorentzKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Dynamic Lorentzian Kubo-Toyabe (LF). Strong collision model. */ virtual Double_t DynamicLorentzKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Fast dynamic Gaussian-Lorentzian KT (ZF). Approximate fast calculation. */ virtual Double_t DynamicGauLorKTZFFast(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Fast dynamic Gaussian-Lorentzian KT (LF). Approximate fast calculation. */ virtual Double_t DynamicGauLorKTLFFast(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Dynamic Gaussian-Lorentzian KT (LF). Full numerical calculation. */ virtual Double_t DynamicGauLorKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Combined Lorentzian-Gaussian KT. Product of both relaxation types. */ virtual Double_t CombiLGKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Stretched Kubo-Toyabe. Formula: exp(-(σt)^β) with KT-like recovery. */ virtual Double_t StrKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Spin glass relaxation function. Edwards-Anderson order parameter. */ virtual Double_t SpinGlass(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Random anisotropic hyperfine coupling. Powder average of anisotropic coupling. */ virtual Double_t RandomAnisotropicHyperfine(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Abragam relaxation. Motional narrowing formula. */ virtual Double_t Abragam(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Transverse field cosine. Formula: cos(φ + 2πνt) */ virtual Double_t TFCos(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Internal field distribution. Gaussian field distribution model. */ virtual Double_t InternalField(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Internal field (Kornilov model). Vortex lattice field distribution. */ virtual Double_t InternalFieldGK(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Internal field (Larkin-Ovchinnikov model). Vortex lattice field distribution. */ virtual Double_t InternalFieldLL(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Bessel function precession. Formula: J₀(2πνt + φ) */ virtual Double_t Bessel(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Internal Bessel field distribution. Combines Bessel with relaxation. */ virtual Double_t InternalBessel(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Skewed Gaussian. Asymmetric relaxation rates before/after zero crossing. */ virtual Double_t SkewedGauss(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Static Nakajima-Keren (ZF). Combined nuclear and electronic relaxation. */ virtual Double_t StaticNKZF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Static Nakajima-Keren (TF). Combined nuclear and electronic relaxation with precession. */ virtual Double_t StaticNKTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Dynamic Nakajima-Keren (ZF). With spin fluctuations. */ virtual Double_t DynamicNKZF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Dynamic Nakajima-Keren (TF). With spin fluctuations and precession. */ virtual Double_t DynamicNKTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief F-μ-F oscillation. Muon bound between two fluorine atoms. */ virtual Double_t FmuF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief μ⁻ exponential TF. Negative muon in transverse field. */ virtual Double_t MuMinusExpTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief Polynomial function. Formula: Σᵢ pᵢtⁱ */ virtual Double_t Polynom(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; /** \brief User-defined function. Calls external shared library function. */ virtual Double_t UserFcn(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const; // -------------------- LF Calculation Helpers -------------------- /** \brief Calculates and caches Gaussian LF integral for static KT. */ virtual void CalculateGaussLFIntegral(const Double_t *val) const; /** \brief Calculates and caches Lorentzian LF integral for static KT. */ virtual void CalculateLorentzLFIntegral(const Double_t *val) const; /** \brief Retrieves cached LF integral value at time t using interpolation. */ virtual Double_t GetLFIntegralValue(const Double_t t) const; /** \brief Calculates dynamic KT in LF using integral equation approach. */ virtual void CalculateDynKTLF(const Double_t *val, Int_t tag) const; /** \brief Retrieves cached dynamic KT LF value at time t. */ virtual Double_t GetDynKTLFValue(const Double_t t) const; /** \brief Retrieves cached dynamic Gauss-Lorentz KT LF value at time t. */ virtual Double_t GetDyn_GL_KTLFValue(const Double_t t) const; // -------------------- Member Variables -------------------- Bool_t fValid; ///< True if this theory node and its parse state are valid UInt_t fType; ///< Theory function type (THEORY_CONST, THEORY_SIMPLE_EXP, etc.) std::vector fParamNo; ///< Resolved parameter indices (0-based). Values >= MSR_PARAM_FUN_OFFSET are function references. UInt_t fNoOfParam; ///< Expected number of parameters for this function type PTheory *fAdd; ///< Pointer to addition child node (left branch of tree) PTheory *fMul; ///< Pointer to multiplication child node (right branch of tree) // User function members Int_t fUserFcnIdx; ///< Index of this user function among all userFcn entries (for global state) TString fUserFcnClassName; ///< ROOT class name for user function (e.g., "TMyFunction") TString fUserFcnSharedLibName; ///< Shared library path (e.g., "libMyFunctions.so") PUserFcnBase *fUserFcn; ///< Pointer to instantiated user function object mutable PDoubleVector fUserParam; ///< Resolved parameter values for user function calls PMsrHandler *fMsrInfo; ///< Pointer to MSR file handler (not owned) // LF calculation caching (mutable for const Func() calls) mutable Double_t fSamplingTime; ///< Time step for LF integral calculation (default 1 ns = 0.001 μs) mutable Double_t fPrevParam[THEORY_MAX_PARAM]; ///< Previous parameter values for cache invalidation check mutable PDoubleVector fLFIntegral; ///< Cached static LF KT integral values mutable Double_t fDynLFdt; ///< Time step for dynamic LF integral equation mutable PDoubleVector fDynLFFuncValue; ///< Cached dynamic Gaussian/Lorentzian LF KT values mutable PDoubleVector fDyn_GL_LFFuncValue; ///< Cached dynamic Gauss-Lorentz LF KT values }; #endif // _PTHEORY_H_