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

use Claude ai to generate doxygen documentation.
All checks were successful
Build and Deploy Documentation / build-and-deploy (push) Successful in 17s
Details

Browse Source
This commit is contained in:
suter_a
2025-11-10 15:14:08 +01:00
parent 262b5a36aa
commit 3d07894b2d
8 changed files with 1749 additions and 129 deletions
Download Patch File Download Diff File Expand all files Collapse all files

179
src/include/PFitter.h
View File

@@ -40,56 +40,170 @@
#include "PRunListCollection.h"
#include "PFitterFcn.h"
//-------------------------------------------------------------
/**
* <p>Minuit2 command identifiers for COMMANDS block.
*
* <p>These constants map MSR file COMMANDS block keywords to internal
* command identifiers. Commands control the fitting process, including
* minimization algorithms, parameter fixing/releasing, error analysis,
* and output options.
*
* <p><b>Categories:</b>
* - <b>Minimizers:</b> MIGRAD (gradient descent), SIMPLEX (non-gradient), MINIMIZE (automatic)
* - <b>Error analysis:</b> HESSE (covariance matrix), MINOS (asymmetric errors)
* - <b>Parameter control:</b> FIX, RELEASE, RESTORE, SAVE
* - <b>Exploration:</b> SCAN (1D/2D parameter space), CONTOURS (error ellipses)
* - <b>Settings:</b> STRATEGY, MACHINE_PRECISION, PRINT
* - <b>Analysis:</b> FIT_RANGE (time-window scan), SECTOR (chisq vs. time)
*/
/// Interactive mode (not implemented in musrfit)
#define PMN_INTERACTIVE 0
/// Contour plot command (2D error ellipses)
#define PMN_CONTOURS 1
/// Eigenvalue analysis (not implemented)
#define PMN_EIGEN 2
/// Fit range scan to find optimal fitting window
#define PMN_FIT_RANGE 3
/// Fix parameters at current values
#define PMN_FIX 4
/// Calculate Hessian matrix for error estimation
#define PMN_HESSE 5
/// Set machine precision for numerical derivatives
#define PMN_MACHINE_PRECISION 6
/// MIGRAD minimizer (gradient-based, recommended)
#define PMN_MIGRAD 7
/// MINIMIZE (automatic algorithm selection)
#define PMN_MINIMIZE 8
/// MINOS error analysis (asymmetric confidence intervals)
#define PMN_MINOS 9
/// Plot command (for scan/contour results)
#define PMN_PLOT 10
/// Release previously fixed parameters
#define PMN_RELEASE 11
/// Restore parameters from previous SAVE
#define PMN_RESTORE 12
/// Save current parameter state
#define PMN_SAVE 13
/// Parameter scan (1D or 2D)
#define PMN_SCAN 14
/// SIMPLEX minimizer (robust but slow)
#define PMN_SIMPLEX 15
/// Set minimization strategy (0=fast, 1=default, 2=careful)
#define PMN_STRATEGY 16
/// Set user covariance matrix (not implemented)
#define PMN_USER_COVARIANCE 17
/// Set user parameter state (not implemented)
#define PMN_USER_PARAM_STATE 18
/// Set print level for fit output (0=quiet, 1=normal, 2=verbose)
#define PMN_PRINT 19
/// Calculate χ² vs. time in sectors (time-window analysis)
#define PMN_SECTOR 20
//-----------------------------------------------------------------------------
/**
* <p>The PSectorChisq class is needed to store the chisq/maxLH of a sector.
* A sector is a time window from fgb to fLast (fLast < lgb). It also stores
* these information for each run of the msr-file.
* <p>Container for sector χ² (or max likelihood) analysis results.
*
* <p>The SECTOR command analyzes how χ² changes as the fit range is extended
* from the first good bin (fgb) to progressively later times. This helps
* identify:
* - Optimal fitting windows
* - Time regions with systematic deviations
* - Data quality issues
* - Relaxation component contributions
*
* <p>A "sector" is a time window [fgb, fLast] where fLast < lgb. The class
* stores χ²/maxLH and NDF for each sector, both for the total fit and for
* individual runs.
*
* <p><b>Use case:</b> If early-time data shows poor χ², the sector analysis
* reveals whether to adjust t0, fgb, or background estimation.
*/
class PSectorChisq
{
public:
/**
* <p>Constructor.
*
* @param noOfRuns Number of runs in the fit (allocates storage for per-run statistics)
*/
PSectorChisq(UInt_t noOfRuns);
/// Sets the first good bin time for a specific run
/// @param first First good bin time in microseconds
/// @param idx Run index
void SetRunFirstTime(Double_t first, UInt_t idx);
/// Sets the sector end time (last bin included in this sector)
/// @param last Sector end time in microseconds
void SetSectorTime(Double_t last) { fLast = last; }
/// Sets the total χ² (or max likelihood) for this sector
/// @param chisq Chi-squared or maximum likelihood value
void SetChisq(Double_t chisq) { fChisq = chisq; }
/// Sets the χ² for a specific run in this sector
/// @param chisq Chi-squared value for run
/// @param idx Run index
void SetChisq(Double_t chisq, UInt_t idx);
/// Sets the expected total χ² for this sector
/// @param expChisq Expected chi-squared value
void SetExpectedChisq(Double_t expChisq) { fExpectedChisq = expChisq; }
/// Sets the expected χ² for a specific run
/// @param chisq Expected chi-squared value
/// @param idx Run index
void SetExpectedChisq(Double_t chisq, UInt_t idx);
/// Sets the total number of degrees of freedom for this sector
/// @param ndf Degrees of freedom
void SetNDF(UInt_t ndf) { fNDF = ndf; }
/// Sets the NDF for a specific run
/// @param ndf Degrees of freedom for run
/// @param idx Run index
void SetNDF(UInt_t ndf, UInt_t idx);
/// Gets the first good bin time for a specific run
/// @param idx Run index
/// @return First good bin time in microseconds
Double_t GetTimeRangeFirst(UInt_t idx);
/// Returns the sector end time
/// @return Last bin time in microseconds
Double_t GetTimeRangeLast() { return fLast; }
/// Returns the total χ² for this sector
/// @return Chi-squared or max likelihood value
Double_t GetChisq() { return fChisq; }
/// Returns the χ² for a specific run
/// @param idx Run index
/// @return Chi-squared value
Double_t GetChisq(UInt_t idx);
/// Returns the expected total χ²
/// @return Expected chi-squared value
Double_t GetExpectedChisq() { return fExpectedChisq; }
/// Returns the expected χ² for a specific run
/// @param idx Run index
/// @return Expected chi-squared value
Double_t GetExpectedChisq(UInt_t idx);
/// Returns the total NDF for this sector
/// @return Degrees of freedom
UInt_t GetNDF() { return fNDF; }
/// Returns the NDF for a specific run
/// @param idx Run index
/// @return Degrees of freedom
UInt_t GetNDF(UInt_t idx);
/// Returns the number of runs
/// @return Number of runs in the analysis
UInt_t GetNoRuns() { return fNoOfRuns; }
private:
@@ -106,17 +220,74 @@ class PSectorChisq
//-----------------------------------------------------------------------------
/**
* <p>Interface class to minuit2.
* <p>Main fitting engine interfacing with ROOT Minuit2.
*
* <p>PFitter orchestrates the entire fitting process for musrfit:
* - Initializes Minuit2 with parameters from MSR file
* - Executes COMMANDS block directives (MIGRAD, HESSE, MINOS, etc.)
* - Manages parameter fixing, releasing, and boundaries
* - Performs error analysis (Hessian, MINOS)
* - Conducts parameter scans and contour plots
* - Calculates sector χ² for time-window analysis
* - Generates fit output and statistics
*
* <p><b>Fitting workflow:</b>
* 1. Initialize parameters and set boundaries
* 2. Execute COMMANDS block sequentially:
* - MIGRAD: Find minimum using gradient descent
* - HESSE: Calculate symmetric errors from covariance matrix
* - MINOS: Calculate asymmetric errors (optional, slower)
* - SAVE: Save parameter state
* 3. Update MSR file with fitted parameters and statistics
*
* <p><b>Minimization modes:</b>
* - <b>χ² minimization:</b> Standard least-squares fitting
* - <b>Maximum likelihood:</b> Poisson statistics (better for low counts)
*
* <p><b>Example COMMANDS block:</b>
* @code
* COMMANDS
* SET PRINT 1
* MIGRAD
* MINOS
* SAVE
* @endcode
*/
class PFitter
{
public:
/**
* <p>Constructor for fitting engine.
*
* @param runInfo Pointer to MSR handler containing fit configuration
* @param runListCollection Pointer to collection of data runs to fit
* @param chisq_only If true, only calculate χ² without fitting
* @param yaml_out If true, generate YAML output file with fit results
*/
PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bool_t chisq_only = false, Bool_t yaml_out = false);
virtual ~PFitter();
/// Returns true if fitter initialized successfully
/// @return Validity status
Bool_t IsValid() { return fIsValid; }
/// Returns true if only parameter scan requested (no minimization)
/// @return Scan-only flag
Bool_t IsScanOnly() { return fIsScanOnly; }
/// Returns true if fit converged successfully
/// @return Convergence status
Bool_t HasConverged() { return fConverged; }
/**
* <p>Executes the complete fitting procedure.
*
* <p>Processes all commands from the COMMANDS block sequentially,
* performs the fit, calculates errors, and prepares output statistics.
*
* @return true if fit completed successfully, false on error
*/
Bool_t DoFit();
private:

189
src/include/PFourier.h
View File

@@ -40,29 +40,105 @@
#include "PMusr.h"
//-------------------------------------------------------------
/**
* <p>Apodization (windowing) strength constants.
*
* <p>Apodization applies a window function to time-domain data before
* Fourier transformation to reduce spectral leakage (Gibbs phenomenon).
* Stronger apodization improves frequency resolution but reduces amplitude
* accuracy.
*/
/// No apodization (rectangular window)
#define F_APODIZATION_NONE 1
/// Weak apodization (gentle roll-off at edges)
#define F_APODIZATION_WEAK 2
/// Medium apodization (moderate roll-off)
#define F_APODIZATION_MEDIUM 3
/// Strong apodization (heavy roll-off for best frequency resolution)
#define F_APODIZATION_STRONG 4
//-------------------------------------------------------------
/**
* Re Fourier phase correction
* <p>Phase correction optimizer for Fourier transforms.
*
* <p>This class performs automatic phase correction on complex Fourier spectra
* to maximize the real component and minimize the imaginary component. Phase
* errors arise from:
* - Uncertain time-zero determination
* - Detector time offsets
* - Signal dispersion
*
* <p><b>Algorithm:</b> Minimizes a combined entropy-penalty functional using
* Minuit2, finding optimal phase parameters (constant + linear dispersion):
* φ(ω) = c₀ + c₁·ω
*
* <p><b>Applications:</b>
* - Improving signal clarity in real Fourier spectra
* - Identifying field distributions in vortex lattices
* - Resolving closely-spaced frequency components
*
* <p><b>Usage:</b> Specify frequency range for optimization to focus on
* signal peaks while avoiding noise regions.
*/
class PFTPhaseCorrection : public ROOT::Minuit2::FCNBase
{
public:
/**
* <p>Constructor for phase correction with default Fourier data.
*
* @param minBin Minimum frequency bin for optimization (-1 = use all)
* @param maxBin Maximum frequency bin for optimization (-1 = use all)
*/
PFTPhaseCorrection(const Int_t minBin=-1, const Int_t maxBin=-1);
/**
* <p>Constructor with explicit Fourier data.
*
* @param reFT Real part of Fourier transform
* @param imFT Imaginary part of Fourier transform
* @param minBin Minimum frequency bin for optimization
* @param maxBin Maximum frequency bin for optimization
*/
PFTPhaseCorrection(std::vector<Double_t> &reFT, std::vector<Double_t> &imFT, const Int_t minBin=-1, const Int_t maxBin=-1);
virtual ~PFTPhaseCorrection() {}
/// Returns true if phase correction initialized successfully
/// @return Validity status
virtual Bool_t IsValid() { return fValid; }
/**
* <p>Performs phase correction minimization.
*
* <p>Uses Minuit2 to find optimal phase parameters that maximize
* the real spectrum while minimizing imaginary components.
*/
virtual void Minimize();
/// Sets the gamma balancing parameter between entropy and penalty
/// @param gamma Balancing factor (typical range: 0.1 to 10)
virtual void SetGamma(const Double_t gamma) { fGamma = gamma; }
/// Sets phase correction parameters manually
/// @param c0 Constant phase offset in degrees
/// @param c1 Linear phase dispersion coefficient
virtual void SetPh(const Double_t c0, const Double_t c1) { fPh_c0 = c0; fPh_c1 = c1; CalcPhasedFT(); CalcRealPhFTDerivative(); }
/// Returns the gamma parameter
/// @return Balancing factor between entropy and penalty
virtual Double_t GetGamma() { return fGamma; }
/**
* <p>Gets phase correction parameter.
*
* @param idx Parameter index (0=c₀, 1=c₁)
* @return Phase parameter value
*/
virtual Double_t GetPhaseCorrectionParam(UInt_t idx);
/// Returns the minimum value of the optimization functional
/// @return Minimum value achieved
virtual Double_t GetMinimum();
private:
@@ -91,32 +167,139 @@ class PFTPhaseCorrection : public ROOT::Minuit2::FCNBase
virtual Double_t operator()(const std::vector<Double_t>&) const;
};
//-------------------------------------------------------------
/**
* muSR Fourier class.
* <p>Fourier transform engine for μSR time-domain data.
*
* <p>PFourier converts time-domain μSR signals to frequency domain,
* revealing:
* - Muon precession frequencies (field measurements)
* - Internal field distributions (superconductors, magnets)
* - Multiple muon stopping sites
* - Dynamic frequency fluctuations
*
* <p><b>Key features:</b>
* - Uses FFTW3 library for efficient FFT computation
* - DC offset removal (for baseline correction)
* - Zero-padding (improves frequency interpolation)
* - Apodization/windowing (reduces spectral leakage)
* - Multiple output formats (real, imaginary, power, phase)
* - Unit conversion (field ↔ frequency)
*
* <p><b>Workflow:</b>
* 1. Create PFourier with time histogram and settings
* 2. Call Transform() with desired apodization
* 3. Retrieve results: GetRealFourier(), GetPowerFourier(), etc.
*
* <p><b>Unit conversions:</b>
* - Gauss: ω(MHz) = γ_μ/(2π) × B(G) = 0.01355 × B(G)
* - Tesla: ω(MHz) = γ_μ/(2π) × B(T) = 135.54 × B(T)
*
* <p><b>Example:</b> TF-μSR measurement at 100 G produces a peak at
* ~1.36 MHz in the Fourier spectrum.
*/
class PFourier
{
public:
/**
* <p>Constructor for Fourier transformation.
*
* @param data Time histogram to transform
* @param unitTag Output units (1=Gauss, 2=Tesla, 3=MHz, 4=Mc/s)
* @param startTime Start time for transform in microseconds (0=from t0)
* @param endTime End time for transform in microseconds (0=to end)
* @param dcCorrected If true, remove DC offset before FFT
* @param zeroPaddingPower Zero-pad to 2^N points (0=no padding)
*/
PFourier(TH1F *data, Int_t unitTag,
Double_t startTime = 0.0, Double_t endTime = 0.0,
Bool_t dcCorrected = false, UInt_t zeroPaddingPower = 0);
virtual ~PFourier();
/**
* <p>Performs the Fourier transformation.
*
* <p>Applies optional apodization, computes FFT using FFTW3,
* and prepares output histograms in requested units.
*
* @param apodizationTag Apodization strength (0/1=none, 2=weak, 3=medium, 4=strong)
*/
virtual void Transform(UInt_t apodizationTag = 0);
/// Returns the original data histogram title
/// @return Title string
virtual const char* GetDataTitle() { return fData->GetTitle(); }
/// Returns the output unit tag (1=G, 2=T, 3=MHz, 4=Mc/s)
/// @return Unit identifier
virtual const Int_t GetUnitTag() { return fUnitTag; }
/// Returns the frequency resolution (bin width in output units)
/// @return Frequency resolution
virtual Double_t GetResolution() { return fResolution; }
/**
* <p>Returns the maximum frequency (Nyquist frequency).
*
* @return Maximum frequency in output units
*/
virtual Double_t GetMaxFreq();
/**
* <p>Gets real part of Fourier transform as histogram.
*
* @param scale Scaling factor for amplitudes (default=1.0)
* @return Pointer to TH1F histogram (caller must delete)
*/
virtual TH1F* GetRealFourier(const Double_t scale = 1.0);
//as virtual TH1F* GetPhaseOptRealFourier(std::vector<Double_t> &phase, const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0);
/**
* <p>Gets imaginary part of Fourier transform as histogram.
*
* @param scale Scaling factor for amplitudes (default=1.0)
* @return Pointer to TH1F histogram (caller must delete)
*/
virtual TH1F* GetImaginaryFourier(const Double_t scale = 1.0);
/**
* <p>Gets power spectrum |F(ω)|² as histogram.
*
* <p>Power spectrum is always positive and shows signal strength
* at each frequency, useful for identifying dominant frequencies.
*
* @param scale Scaling factor for power (default=1.0)
* @return Pointer to TH1F histogram (caller must delete)
*/
virtual TH1F* GetPowerFourier(const Double_t scale = 1.0);
/**
* <p>Gets phase spectrum arg(F(ω)) as histogram.
*
* @param scale Scaling factor (default=1.0)
* @return Pointer to TH1F histogram (caller must delete)
*/
virtual TH1F* GetPhaseFourier(const Double_t scale = 1.0);
/**
* <p>Static method for phase-optimized real Fourier spectrum.
*
* <p>Applies phase correction to maximize real component using
* provided phase parameters.
*
* @param re Real part of Fourier transform
* @param im Imaginary part of Fourier transform
* @param phase Phase correction parameters [c₀, c₁]
* @param scale Scaling factor (default=1.0)
* @param min Minimum frequency for correction (-1=all)
* @param max Maximum frequency for correction (-1=all)
* @return Pointer to phase-corrected TH1F histogram
*/
static TH1F* GetPhaseOptRealFourier(const TH1F *re, const TH1F *im, std::vector<Double_t> &phase,
const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0);
/// Returns true if Fourier transform is ready
/// @return Validity status
virtual Bool_t IsValid() { return fValid; }
private:

229
src/include/PMsrHandler.h
View File

@@ -44,67 +44,292 @@
//-------------------------------------------------------------
/**
* <p>This class provides the routines needed to handle msr-files, i.e. reading, writing, parsing, etc.
* <p>MSR file parser and manager for musrfit framework.
*
* <p>This class handles all MSR (Muon Spin Rotation) file operations including:
* - Reading and parsing MSR files
* - Writing MSR files and log files
* - Validating MSR file content (parameter ranges, theory consistency, etc.)
* - Managing fit parameters, theory definitions, run configurations
* - Evaluating user-defined functions
* - Handling parameter mapping between runs
*
* <p>The PMsrHandler performs comprehensive validation during parsing, checking for:
* - Parameter name uniqueness
* - Valid theory-to-parameter mappings
* - Histogram grouping consistency
* - RRF (Rotating Reference Frame) settings
* - Function syntax and parameter usage
*
* <p><b>Usage pattern:</b>
* @code
* PMsrHandler handler("myfile.msr");
* if (handler.ReadMsrFile() == PMUSR_SUCCESS) {
* PMsrParamList *params = handler.GetMsrParamList();
* // Process parameters, run fits, etc.
* handler.WriteMsrFile("output.msr");
* }
* @endcode
*/
class PMsrHandler
{
public:
/**
* <p>Constructor for PMsrHandler.
*
* @param fileName Path to MSR file to read/write
* @param startupOptions Optional startup configuration (from musrfit_startup.xml)
* @param fourierOnly If true, only parse Fourier-related blocks (for musrFT)
*/
PMsrHandler(const Char_t *fileName, PStartupOptions *startupOptions=0, const Bool_t fourierOnly=false);
virtual ~PMsrHandler();
/**
* <p>Reads and parses the MSR file.
*
* <p>Performs comprehensive parsing of all MSR file blocks including
* TITLE, FITPARAMETER, THEORY, FUNCTIONS, GLOBAL, RUN, COMMANDS,
* FOURIER, PLOT, and STATISTIC blocks. Validates consistency and
* reports detailed error messages on failure.
*
* @return PMUSR_SUCCESS on success, negative error code on failure
*
* @see PMUSR_MSR_FILE_NOT_FOUND
* @see PMUSR_MSR_SYNTAX_ERROR
*/
virtual Int_t ReadMsrFile();
/**
* <p>Writes a log file with MSR file content and parsing information.
*
* <p>Creates a .mlog file containing the parsed MSR structure,
* useful for debugging and verifying parameter interpretation.
*
* @param messages If true, includes informational messages in log
* @return PMUSR_SUCCESS on success, negative error code on failure
*/
virtual Int_t WriteMsrLogFile(const Bool_t messages = true);
/**
* <p>Writes MSR file with updated parameters and results.
*
* <p>Writes a complete MSR file, optionally preserving user comments
* from specific blocks. Typically called after fitting to save
* fitted parameters and statistics.
*
* @param filename Output MSR file path
* @param commentsPAR Optional comments for FITPARAMETER block (line number → comment)
* @param commentsTHE Optional comments for THEORY block
* @param commentsFUN Optional comments for FUNCTIONS block
* @param commentsRUN Optional comments for RUN block
* @return PMUSR_SUCCESS on success, negative error code on failure
*/
virtual Int_t WriteMsrFile(const Char_t *filename, std::map<UInt_t, TString> *commentsPAR = 0, std::map<UInt_t, TString> *commentsTHE = 0, \
std::map<UInt_t, TString> *commentsFUN = 0, std::map<UInt_t, TString> *commentsRUN = 0);
/// Returns pointer to MSR file title string
virtual TString* GetMsrTitle() { return &fTitle; }
/// Returns pointer to fit parameter list
virtual PMsrParamList* GetMsrParamList() { return &fParam; }
/// Returns pointer to THEORY block lines
virtual PMsrLines* GetMsrTheory() { return &fTheory; }
/// Returns pointer to FUNCTIONS block lines
virtual PMsrLines* GetMsrFunctions() { return &fFunctions; }
/// Returns pointer to GLOBAL block settings
virtual PMsrGlobalBlock* GetMsrGlobal() { return &fGlobal; }
/// Returns pointer to list of RUN blocks
virtual PMsrRunList* GetMsrRunList() { return &fRuns; }
/// Returns pointer to COMMANDS block lines
virtual PMsrLines* GetMsrCommands() { return &fCommands; }
/// Returns pointer to FOURIER block settings
virtual PMsrFourierStructure* GetMsrFourierList() { return &fFourier; }
/// Returns pointer to list of PLOT blocks
virtual PMsrPlotList* GetMsrPlotList() { return &fPlots; }
/// Returns pointer to STATISTIC block
virtual PMsrStatisticStructure* GetMsrStatistic() { return &fStatistic; }
/// Returns pointer to MSR file directory path
virtual TString* GetMsrFileDirectoryPath() { return &fMsrFileDirectoryPath; }
/// Returns the number of RUN blocks in MSR file
virtual UInt_t GetNoOfRuns() { return fRuns.size(); }
/// Returns the number of fit parameters in FITPARAMETER block
virtual UInt_t GetNoOfParams() { return fParam.size(); }
/// Returns the MSR file name
virtual const TString& GetFileName() const { return fFileName; }
/// Sets the MSR file title
/// @param title New title string
virtual void SetMsrTitle(const TString &title) { fTitle = title; }
/**
* <p>Sets the value of a fit parameter.
*
* @param i Parameter index (0-based)
* @param value New parameter value
* @return true on success, false if index out of range
*/
virtual Bool_t SetMsrParamValue(UInt_t i, Double_t value);
/**
* <p>Sets the step size (or error) of a fit parameter.
*
* @param i Parameter index (0-based)
* @param value New step/error value
* @return true on success, false if index out of range
*/
virtual Bool_t SetMsrParamStep(UInt_t i, Double_t value);
/**
* <p>Sets whether positive error is present for a parameter.
*
* @param i Parameter index (0-based)
* @param value True if positive error is defined
* @return true on success, false if index out of range
*/
virtual Bool_t SetMsrParamPosErrorPresent(UInt_t i, Bool_t value);
/**
* <p>Sets the positive error value for a parameter (asymmetric errors).
*
* @param i Parameter index (0-based)
* @param value Positive error value
* @return true on success, false if index out of range
*/
virtual Bool_t SetMsrParamPosError(UInt_t i, Double_t value);
/**
* <p>Sets a time-zero bin entry for a specific run.
*
* @param runNo Run block number (0-based)
* @param idx Histogram index within t0 list
* @param bin Time-zero bin value
*/
virtual void SetMsrT0Entry(UInt_t runNo, UInt_t idx, Double_t bin);
/**
* <p>Sets a time-zero bin for an addrun histogram.
*
* @param runNo Run block number (0-based)
* @param addRunIdx Index of addrun entry
* @param histoIdx Histogram index within addrun
* @param bin Time-zero bin value
*/
virtual void SetMsrAddT0Entry(UInt_t runNo, UInt_t addRunIdx, UInt_t histoIdx, Double_t bin);
/**
* <p>Sets a data range bin entry for a specific run.
*
* @param runNo Run block number (0-based)
* @param idx Data range index (0=start, 1=end, etc.)
* @param bin Data range bin value
*/
virtual void SetMsrDataRangeEntry(UInt_t runNo, UInt_t idx, Int_t bin);
/**
* <p>Sets a background range bin entry for a specific run.
*
* @param runNo Run block number (0-based)
* @param idx Background range index (0=start, 1=end, etc.)
* @param bin Background range bin value
*/
virtual void SetMsrBkgRangeEntry(UInt_t runNo, UInt_t idx, Int_t bin);
/// Flags that STATISTIC block should be copied as-is (for musrt0)
virtual void CopyMsrStatisticBlock() { fCopyStatisticsBlock = true; }
/// Sets whether fit converged in STATISTIC block
/// @param converged True if fit converged successfully
virtual void SetMsrStatisticConverged(Bool_t converged) { fStatistic.fValid = converged; }
/// Sets the minimum χ² (or max likelihood) in STATISTIC block
/// @param min Minimum value
virtual void SetMsrStatisticMin(Double_t min) { fStatistic.fMin = min; }
/// Sets the number of degrees of freedom in STATISTIC block
/// @param ndf Degrees of freedom
virtual void SetMsrStatisticNdf(UInt_t ndf) { fStatistic.fNdf = ndf; }
/// Returns the number of user-defined functions in FUNCTIONS block
virtual Int_t GetNoOfFuncs() { return fFuncHandler->GetNoOfFuncs(); }
/**
* <p>Gets function number by index.
*
* @param idx Function index (0-based)
* @return Function number as defined in FUNCTIONS block
*/
virtual UInt_t GetFuncNo(Int_t idx) { return fFuncHandler->GetFuncNo(idx); }
/**
* <p>Gets function index from function number.
*
* @param funNo Function number
* @return Function index (0-based)
*/
virtual UInt_t GetFuncIndex(Int_t funNo) { return fFuncHandler->GetFuncIndex(funNo); }
/**
* <p>Checks if map and parameter ranges are valid for functions.
*
* @param mapSize Size of map vector
* @param paramSize Number of available parameters
* @return true if ranges are valid
*/
virtual Bool_t CheckMapAndParamRange(UInt_t mapSize, UInt_t paramSize)
{ return fFuncHandler->CheckMapAndParamRange(mapSize, paramSize); }
/**
* <p>Evaluates a user-defined function.
*
* @param i Function index
* @param map Parameter mapping vector
* @param param Parameter value vector
* @param metaData Experimental metadata (field, temperature, etc.)
* @return Evaluated function value
*/
virtual Double_t EvalFunc(UInt_t i, std::vector<Int_t> map, std::vector<Double_t> param, PMetaData metaData)
{ return fFuncHandler->Eval(i, map, param, metaData); }
/**
* <p>Gets the number of fit parameters used in a specific theory line.
*
* @param idx Theory line index
* @return Number of parameters used
*/
virtual UInt_t GetNoOfFitParameters(UInt_t idx);
/**
* <p>Checks if a parameter is used in theory or functions.
*
* @param paramNo Parameter number (1-based as in MSR file)
* @return 1 if used, 0 if unused, -1 on error
*/
virtual Int_t ParameterInUse(UInt_t paramNo);
/**
* <p>Generates a grouping string for histogram display.
*
* @param runNo Run block number
* @param detector Detector identifier ("forward" or "backward")
* @param groupingStr Output grouping string
*/
virtual void GetGroupingString(Int_t runNo, TString detector, TString &groupingStr);
/**
* <p>Estimates N0 parameter for single histogram fits.
*
* <p>Uses data amplitude at t=0 to provide initial N0 estimate,
* improving fit convergence for single histogram fits.
*
* @return true on success
*/
virtual Bool_t EstimateN0();
/// Returns the last error message as a string
/// @return Error message string
virtual std::string GetLastErrorMsg() { return fLastErrorMsg.str(); }
private:

684
src/include/PMusr.h
View File

File diff suppressed because it is too large Load Diff

161
src/include/PRunBase.h
View File

@@ -42,53 +42,166 @@
//------------------------------------------------------------------------------------------
/**
* <p>The run base class is enforcing a common interface to all supported fit-types.
* <p>Abstract base class defining the interface for all μSR fit types.
*
* <p>PRunBase establishes a common API for processing and fitting different
* types of μSR data (single histogram, asymmetry, RRF, etc.). Derived classes
* implement specific data processing and χ² calculation for each fit type:
* - <b>PRunSingleHisto:</b> Single detector histogram fits
* - <b>PRunAsymmetry:</b> Asymmetry fits (forward - backward)
* - <b>PRunSingleHistoRRF:</b> Single histogram in rotating reference frame
* - <b>PRunAsymmetryRRF:</b> Asymmetry in rotating reference frame
* - <b>PRunMuMinus:</b> Negative muon fits
* - <b>PRunNonMusr:</b> General x-y data fits
*
* <p><b>Key responsibilities:</b>
* - Loading and preprocessing raw histogram data
* - Background subtraction and normalization
* - Time bin packing/rebinning
* - Theory evaluation at data points
* - χ² or maximum likelihood calculation
* - RRF transformations (if applicable)
*
* <p><b>Processing workflow:</b>
* 1. <tt>PrepareData()</tt> - Load and preprocess raw data
* 2. <tt>CalcTheory()</tt> - Evaluate theory function
* 3. <tt>CalcChiSquare()</tt> or <tt>CalcMaxLikelihood()</tt> - Compute fit metric
*
* <p><b>Design pattern:</b> Template Method - base class defines workflow,
* derived classes implement fit-type-specific algorithms.
*/
class PRunBase
{
public:
/// Default constructor
PRunBase();
/**
* <p>Constructor initializing run from MSR file and raw data.
*
* @param msrInfo Pointer to MSR file handler
* @param rawData Pointer to raw data handler
* @param runNo Run number (0-based index in MSR file)
* @param tag Operation mode (kFit, kView)
*/
PRunBase(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag);
virtual ~PRunBase();
virtual Double_t CalcChiSquare(const std::vector<Double_t>& par) = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!!
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par) = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!!
/**
* <p>Calculates χ² between data and theory (pure virtual).
*
* <p>χ² = Σ[(data_i - theory_i)² / σ_i²] summed over all data points.
* This is the standard least-squares metric for fitting.
*
* @param par Vector of fit parameter values
* @return Chi-squared value
*/
virtual Double_t CalcChiSquare(const std::vector<Double_t>& par) = 0;
/**
* <p>Calculates maximum likelihood (pure virtual).
*
* <p>For Poisson statistics: -2ln(L) = 2Σ[theory_i - data_i·ln(theory_i)]
* Better than χ² for low-count data where Gaussian approximation fails.
*
* @param par Vector of fit parameter values
* @return Negative 2 times log-likelihood
*/
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par) = 0;
/**
* <p>Sets the fit time range for this run.
*
* <p>Updates the fitting window, useful for FIT_RANGE command which
* scans different time windows to find the optimal range.
*
* @param fitRange Vector of (start, end) time pairs in microseconds
*/
virtual void SetFitRange(PDoublePairVector fitRange);
virtual void CalcTheory() = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!!
/**
* <p>Evaluates theory function at all data points (pure virtual).
*
* <p>Uses current parameter values to calculate expected signal
* at each time point. Called during each fit iteration.
*/
virtual void CalcTheory() = 0;
virtual UInt_t GetRunNo() { return fRunNo; } ///< returns the number of runs of the msr-file
virtual PRunData* GetData() { return &fData; } ///< returns the data to be fitted
/// Returns the run number (0-based index in MSR file)
/// @return Run number
virtual UInt_t GetRunNo() { return fRunNo; }
/// Returns pointer to processed data (background-corrected, binned)
/// @return Pointer to PRunData object
virtual PRunData* GetData() { return &fData; }
/// Cleans up internal data structures
virtual void CleanUp();
virtual Bool_t IsValid() { return fValid; } ///< returns if the state is valid
/// Returns validity status of this run object
/// @return true if run initialized successfully
virtual Bool_t IsValid() { return fValid; }
protected:
Bool_t fValid; ///< flag showing if the state of the class is valid
Bool_t fValid; ///< Flag showing if the run object initialized successfully
EPMusrHandleTag fHandleTag; ///< tag telling whether this is used for fit, view, ...
EPMusrHandleTag fHandleTag; ///< Operation mode: kFit (fitting), kView (display only)
Int_t fRunNo; ///< number of the run within the msr-file
PMsrHandler *fMsrInfo; ///< msr-file handler
PMsrRunBlock *fRunInfo; ///< run info used to filter out needed infos of a run
PRunDataHandler *fRawData; ///< holds the raw run data
Int_t fRunNo; ///< Run number (0-based index in MSR file RUN blocks)
PMsrHandler *fMsrInfo; ///< Pointer to MSR file handler
PMsrRunBlock *fRunInfo; ///< Pointer to this run's RUN block settings
PRunDataHandler *fRawData; ///< Pointer to raw data handler for this run
PRunData fData; ///< data to be fitted, viewed, i.e. binned data
Double_t fTimeResolution; ///< time resolution in (us)
PMetaData fMetaData; ///< keeps the meta data from the data file like field, temperature, energy, ...
PDoubleVector fT0s; ///< all t0 bins of a run! The derived classes will handle it.
std::vector<PDoubleVector> fAddT0s; ///< all t0 bins of all addrun's of a run! The derived classes will handle it.
PRunData fData; ///< Processed data ready for fitting (background-corrected, packed, etc.)
Double_t fTimeResolution; ///< Time resolution of raw data in microseconds (μs)
PMetaData fMetaData; ///< Experimental metadata (field, temperature, energy) from data file
PDoubleVector fT0s; ///< Time-zero bins for all histograms in this run
std::vector<PDoubleVector> fAddT0s; ///< Time-zero bins for all addrun histograms
Double_t fFitStartTime; ///< fit start time
Double_t fFitEndTime; ///< fit end time
Double_t fFitStartTime; ///< Fit range start time in microseconds (μs)
Double_t fFitEndTime; ///< Fit range end time in microseconds (μs)
PDoubleVector fFuncValues; ///< is keeping the values of the functions from the FUNCTIONS block
std::unique_ptr<PTheory> fTheory; ///< theory needed to calculate chi-square
PDoubleVector fFuncValues; ///< Values of user-defined functions from FUNCTIONS block
std::unique_ptr<PTheory> fTheory; ///< Theory function evaluator for χ² calculation
PDoubleVector fKaiserFilter; ///< stores the Kaiser filter vector (needed for the RRF).
PDoubleVector fKaiserFilter; ///< Kaiser window coefficients for RRF filtering
virtual Bool_t PrepareData() = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!!
/**
* <p>Prepares raw data for fitting (pure virtual).
*
* <p>Performs fit-type-specific preprocessing:
* - Background subtraction
* - Asymmetry calculation (if applicable)
* - Time bin packing/rebinning
* - RRF transformation (if applicable)
* - Error propagation
*
* <p>Called during initialization before fitting begins.
*
* @return true on success, false on error
*/
virtual Bool_t PrepareData() = 0;
/**
* <p>Calculates Kaiser window filter coefficients for RRF.
*
* <p>The Kaiser window reduces spectral leakage when transforming
* to the rotating reference frame. Parameters control the trade-off
* between main lobe width and side lobe suppression.
*
* @param wc Cutoff frequency (normalized, 0 to π)
* @param A Attenuation in dB (controls side lobe level)
* @param dw Transition width (normalized)
*/
virtual void CalculateKaiserFilterCoeff(Double_t wc, Double_t A, Double_t dw);
/**
* <p>Applies Kaiser filter to theory values for RRF.
*
* <p>Filters the theory function in the same way data is filtered,
* ensuring consistent comparison between data and theory in RRF fits.
*/
virtual void FilterTheo();
};

148
src/include/PRunDataHandler.h
View File

@@ -35,42 +35,188 @@
#include "PMusr.h"
#include "PMsrHandler.h"
//-------------------------------------------------------------
/**
* <p>Data file format identifiers for any2many converter.
*
* <p>These constants identify different μSR data file formats supported
* by musrfit for reading and conversion.
*/
/// Undefined or unknown format
#define A2M_UNDEFINED 0
/// Generic ROOT file
#define A2M_ROOT 1
/// MusrRoot format (PSI-specific ROOT structure)
#define A2M_MUSR_ROOT 2
/// MusrRoot with directory structure
#define A2M_MUSR_ROOT_DIR 3
/// PSI binary format
#define A2M_PSIBIN 4
/// PSI MDU ASCII format
#define A2M_PSIMDU 5
/// TRIUMF MUD (Muon Data) format
#define A2M_MUD 6
/// NeXus HDF5 format (ISIS, JPARC)
#define A2M_NEXUS 7
/// WKM format (older PSI format)
#define A2M_WKM 8
/// Generic ASCII format
#define A2M_ASCII 9
//-------------------------------------------------------------
/**
* <p>Handler class needed to read/handle raw data files.
* <p>Raw data file reader and format converter.
*
* <p>PRunDataHandler is the I/O layer for musrfit, responsible for:
* - Loading raw histogram data from multiple file formats
* - Converting between different μSR data formats
* - Searching multiple data paths for run files
* - Reading run metadata (field, temperature, time, etc.)
* - Managing collections of runs for simultaneous fits
*
* <p><b>Supported file formats:</b>
* - <b>MusrRoot:</b> PSI ROOT-based format with comprehensive metadata
* - <b>NeXus:</b> HDF5-based format (ISIS, JPARC)
* - <b>MUD:</b> TRIUMF's binary format
* - <b>PSI-BIN:</b> Legacy PSI binary format
* - <b>WKM:</b> Older PSI format
* - <b>ASCII:</b> Text-based formats (for non-μSR data)
* - <b>DB/DAT:</b> Database formats for general x-y data
*
* <p><b>Key features:</b>
* - Automatic format detection
* - Run name template expansion (e.g., "run%r.root" → "run2425.root")
* - Multi-path search (search common data directories)
* - Batch file conversion (any2many tool)
* - Metadata extraction and validation
*
* <p><b>Example usage:</b>
* @code
* PRunDataHandler handler(msrInfo, dataPaths);
* handler.ReadData(); // Reads all runs specified in MSR file
* PRawRunData *data = handler.GetRunData(0);
* @endcode
*/
class PRunDataHandler
{
public:
/// Default constructor
PRunDataHandler();
/**
* <p>Constructor for single file with explicit format.
*
* @param fileName Path to data file
* @param fileFormat Format string ("root", "nexus", "mud", etc.)
*/
PRunDataHandler(TString fileName, const TString fileFormat);
/**
* <p>Constructor with search paths.
*
* @param fileName File name or template
* @param fileFormat Format string
* @param dataPath Vector of directories to search for data files
*/
PRunDataHandler(TString fileName, const TString fileFormat, const PStringVector dataPath);
/**
* <p>Constructor for reading single file into provided structure.
*
* @param fileName File name
* @param fileFormat Format string
* @param dataPath Data search path
* @param runData Reference to PRawRunData to fill
*/
PRunDataHandler(TString fileName, const TString fileFormat, const TString dataPath, PRawRunData &runData);
/**
* <p>Constructor for format conversion (any2many).
*
* @param any2ManyInfo Conversion configuration structure
*/
PRunDataHandler(PAny2ManyInfo *any2ManyInfo);
/**
* <p>Constructor for format conversion with search paths.
*
* @param any2ManyInfo Conversion configuration
* @param dataPath Vector of search directories
*/
PRunDataHandler(PAny2ManyInfo *any2ManyInfo, const PStringVector dataPath);
/**
* <p>Constructor for MSR-based data loading.
*
* @param msrInfo MSR file handler (provides run list and paths)
*/
PRunDataHandler(PMsrHandler *msrInfo);
/**
* <p>Constructor for MSR-based loading with custom search paths.
*
* @param msrInfo MSR file handler
* @param dataPath Vector of directories to search
*/
PRunDataHandler(PMsrHandler *msrInfo, const PStringVector dataPath);
virtual ~PRunDataHandler();
/**
* <p>Reads all data files specified in MSR file or configuration.
*
* <p>Searches data paths, detects formats, loads histograms and
* metadata. Call this before attempting to access run data.
*/
virtual void ReadData();
/**
* <p>Performs format conversion (for any2many utility).
*
* <p>Converts loaded data to the target format specified in
* any2many configuration.
*/
virtual void ConvertData();
/**
* <p>Writes data to file in specified format.
*
* @param fileName Output file name (empty = use default)
* @return true on success, false on error
*/
virtual Bool_t WriteData(TString fileName="");
/// Returns true if all required data files were successfully loaded
/// @return Data availability status
virtual Bool_t IsAllDataAvailable() const { return fAllDataAvailable; }
/**
* <p>Gets run data by run name.
*
* @param runName Name of run to retrieve
* @return Pointer to PRawRunData, or nullptr if not found
*/
virtual PRawRunData* GetRunData(const TString &runName);
/**
* <p>Gets run data by index.
*
* @param idx Run index (0-based)
* @return Pointer to PRawRunData, or nullptr if out of range
*/
virtual PRawRunData* GetRunData(const UInt_t idx=0);
/// Returns the number of loaded run data sets
/// @return Number of runs
virtual Int_t GetNoOfRunData() {return fData.size(); }
/**
* <p>Sets or replaces run data at specified index.
*
* @param data Pointer to PRawRunData to store
* @param idx Index where to store data (default=0)
* @return true on success
*/
virtual Bool_t SetRunData(PRawRunData *data, UInt_t idx=0);
private:

178
src/include/PTheory.h
View File

@@ -37,49 +37,109 @@
#include "PMsrHandler.h"
#include "PUserFcnBase.h"
// --------------------------------------------------------
// function handling tags
// --------------------------------------------------------
//-------------------------------------------------------------
/**
* <p>Theory function type tags.
*
* <p>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.
*
* <p>Theory functions can be combined using addition (+) and multiplication (*):
* - <b>Addition:</b> Independent relaxation channels (e.g., 1/3 fast + 2/3 slow)
* - <b>Multiplication:</b> Multiple relaxation mechanisms (e.g., precession * damping)
*
* <p><b>Categories:</b>
* - <b>Basic:</b> Constants, asymmetries, simple exponentials
* - <b>Static relaxation:</b> Gaussian, Lorentzian (Kubo-Toyabe)
* - <b>Dynamic relaxation:</b> Motional narrowing, spin fluctuations
* - <b>Precession:</b> Cosine, Bessel functions for oscillations
* - <b>Vortex lattice:</b> Internal field distributions in superconductors
* - <b>Special:</b> Spin glass, Abragam, mu-minus, polynomials, user functions
*/
// function tags
/// 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
// function parameter tags, i.e. how many parameters has a specific function
// if there is a comment with a (tshift), the number of parameters is increased by one
//-------------------------------------------------------------
/**
* <p>Number of parameters for each theory function.
*
* <p>These constants define how many parameters each theory function requires,
* <b>excluding</b> optional time shift. If a function includes time shift,
* add 1 to the parameter count.
*
* <p>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)
@@ -113,31 +173,64 @@
#define THEORY_PARAM_F_MU_F 1 // frequency (tshift)
#define THEORY_PARAM_MU_MINUS_EXP 6 // N0, tau, A, damping, phase, frequency (tshift)
// number of available user functions
//-------------------------------------------------------------
/**
* <p>Maximum number of theory functions in database.
*
* <p>This is the total number of built-in theory functions available,
* including all relaxation, precession, and special functions.
*/
#define THEORY_MAX 34
// maximal number of parameters. Needed in the contents of LF
//-------------------------------------------------------------
/**
* <p>Maximum number of parameters for any theory function.
*
* <p>Used to allocate arrays for longitudinal field calculations where
* parameter values from previous iterations must be cached.
*/
#define THEORY_MAX_PARAM 10
// deg -> rad factor
//-------------------------------------------------------------
/**
* <p>Conversion factor from degrees to radians.
*
* <p>Value: π/180 = 0.017453292519943295
* <p>Used for phase parameters which are specified in degrees in MSR files
* but converted to radians for calculations.
*/
#define DEG_TO_RAD 0.0174532925199432955
// 2 pi
/**
* <p>Mathematical constant 2π.
*
* <p>Used extensively in frequency-to-angular-frequency conversions:
* ω = 2πν
*/
#define TWO_PI 6.28318530717958623
class PTheory;
//--------------------------------------------------------------------------------------
/**
* <p>Structure holding the infos of a the available internally defined functions.
* <p>Database entry for a theory function definition.
*
* <p>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; ///< function tag
UInt_t fNoOfParam; ///< number of parameters for this function
Bool_t fTable; ///< table flag, indicating if the function is generate from a table
TString fName; ///< name of the function as written into the msr-file
TString fAbbrev; ///< abbreviation of the function name
TString fComment; ///< comment added in the msr-file theory block to help the used
TString fCommentTimeShift; ///< comment added in the msr-file theory block if there is a time shift
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;
//--------------------------------------------------------------------------------------
@@ -251,15 +344,62 @@ static PTheoDataBase fgTheoDataBase[THEORY_MAX] = {
//--------------------------------------------------------------------------------------
/**
* <p>Class handling a theory of a msr-file.
* <p>Theory function evaluator and manager.
*
* <p>This class parses, validates, and evaluates theory functions specified
* in the THEORY block of MSR files. It handles:
* - Parsing theory syntax (function names, parameters, operators)
* - Building theory expression trees (addition/multiplication)
* - Resolving parameter mappings and user functions
* - Evaluating theory at specific time points
* - Caching longitudinal field calculations
* - Loading and managing user-defined shared library functions
*
* <p><b>Theory syntax:</b>
* - Single function: <tt>asymmetry simpleExp 1 3</tt>
* - Addition: <tt>func1 par1... + func2 par1...</tt>
* - Multiplication: <tt>func1 par1... * func2 par1...</tt>
* - Nested: <tt>(func1 + func2) * func3</tt>
*
* <p><b>Parameter resolution:</b>
* Parameters can be direct numbers, parameter references (1, 2, 3, ...),
* map references (map1, map2, ...), or function references (fun1, fun2, ...).
*
* <p><b>Example:</b>
* @code
* THEORY
* asymmetry 1
* TFieldCos 2 3 (phase, freq)
* +
* simpleExp 4 (rate)
* @endcode
* This evaluates to: par1 * cos(par2 + 2π*par3*t) + exp(-par4*t)
*/
class PTheory
{
public:
/**
* <p>Constructor for theory handler.
*
* @param msrInfo Pointer to MSR file handler
* @param runNo Run number for this theory (0-based)
* @param hasParent True if this is a sub-theory (for +/* operators)
*/
PTheory(PMsrHandler *msrInfo, UInt_t runNo, const Bool_t hasParent = false);
virtual ~PTheory();
/// Returns true if theory is valid and ready for evaluation
virtual Bool_t IsValid();
/**
* <p>Evaluates the theory function at time t.
*
* @param t Time in microseconds (μs)
* @param paramValues Vector of fit parameter values
* @param funcValues Vector of evaluated user function values
* @return Theory value at time t
*/
virtual Double_t Func(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const;
private:

110
src/include/PUserFcnBase.h
View File

@@ -37,26 +37,124 @@
//--------------------------------------------------------------------------------------------
/**
* <p>Interface class for the user function.
* <p>Abstract base class for user-defined theory functions.
*
* <p>PUserFcnBase enables extending musrfit with custom theory functions
* beyond the 33 built-in functions. Users create derived classes implementing
* specific physics models, compile them into shared libraries, and load them
* dynamically at runtime.
*
* <p><b>Use cases:</b>
* - Novel relaxation mechanisms not in standard library
* - Material-specific models (e.g., Skyrmion lattices)
* - Complex multi-component functions
* - Proprietary or experimental theory functions
* - Functions requiring external libraries (GSL, CUDA, etc.)
*
* <p><b>Implementation steps:</b>
* 1. Create a class deriving from PUserFcnBase
* 2. Implement operator()(t, param) with your theory
* 3. Optionally implement global part for heavy initialization
* 4. Compile to shared library (.so/.dylib/.dll)
* 5. Reference in MSR file THEORY block: "userFcn libMyFunc TMyFuncClass"
*
* <p><b>Example minimal implementation:</b>
* @code
* class TMyRelaxation : public PUserFcnBase {
* public:
* Double_t operator()(Double_t t, const std::vector<Double_t> &par) const {
* // par[0] = rate, par[1] = exponent, par[2] = time shift
* Double_t tt = t - par[2];
* if (tt < 0) return 0.0;
* return exp(-pow(par[0]*tt, par[1]));
* }
* ClassDef(TMyRelaxation, 1)
* };
* @endcode
*
* <p><b>Global part:</b> For expensive one-time computations (lookup tables,
* matrix inversions), override NeedGlobalPart(), SetGlobalPart(), and
* GlobalPartIsValid(). The global part is initialized once and shared across
* all fit iterations.
*
* <p><b>MSR file usage:</b>
* @code
* THEORY
* asymmetry 1
* userFcn libMyRelax.so TMyRelaxation map1 2 0.5 (rate, expo, tshift)
* @endcode
*/
class PUserFcnBase : public TObject
{
public:
/// Default constructor
PUserFcnBase() {}
/// Virtual destructor
virtual ~PUserFcnBase() {}
virtual Bool_t NeedGlobalPart() const { return false; } ///< if a user function needs a global part this function should return true, otherwise false (default: false)
virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) {} ///< if a user function is using a global part, this function is used to invoke and retrieve the proper global object
virtual Bool_t GlobalPartIsValid() const { return false; } ///< if a user function is using a global part, this function returns if the global object part is valid (default: false)
/**
* <p>Indicates if this function requires global initialization.
*
* <p>Override to return true if your function needs expensive one-time
* setup (e.g., calculating lookup tables, loading external data).
* The global part is computed once and reused across fit iterations.
*
* @return true if global part needed, false otherwise (default: false)
*/
virtual Bool_t NeedGlobalPart() const { return false; }
/**
* <p>Sets up the global part of the user function.
*
* <p>Called once during initialization if NeedGlobalPart() returns true.
* Use this to allocate and initialize shared data structures.
*
* @param globalPart Vector of global objects (one per run)
* @param idx Index of this run's global object in the vector
*/
virtual void SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) {}
/**
* <p>Checks if the global part initialized correctly.
*
* <p>Override to validate your global data structure after SetGlobalPart().
*
* @return true if global part is valid and ready, false otherwise (default: false)
*/
virtual Bool_t GlobalPartIsValid() const { return false; }
/**
* <p>Evaluates the user function at time t (pure virtual).
*
* <p>This is the core evaluation method called for each data point
* during fitting. Must be implemented by derived classes.
*
* <p><b>Parameter convention:</b> The last parameter is typically the
* time shift. Earlier parameters are model-specific (rates, amplitudes,
* exponents, etc.).
*
* @param t Time in microseconds (μs)
* @param param Vector of function parameters (from MSR file + maps)
* @return Function value at time t
*/
virtual Double_t operator()(Double_t t, const std::vector<Double_t> &param) const = 0;
ClassDef(PUserFcnBase, 1)
};
//--------------------------------------------------------------------------
// This function is a replacement for the ParseFile method of TSAXParser.
//--------------------------------------------------------------------------
/**
* <p>XML file parser for user function configurations.
*
* <p>This function provides a replacement for TSAXParser::ParseFile with
* better error handling for XML configuration files used by some user
* functions to load parameters or settings.
*
* @param parser Pointer to TSAXParser object
* @param fileName Path to XML file to parse
* @return 0 on success, error code on failure
*/
Int_t parseXmlFile(TSAXParser*, const Char_t*);
#endif // _PUSERFCNBASE_H_