diff --git a/src/include/PFitter.h b/src/include/PFitter.h index 29d970ede..a7bd95d84 100644 --- a/src/include/PFitter.h +++ b/src/include/PFitter.h @@ -40,56 +40,170 @@ #include "PRunListCollection.h" #include "PFitterFcn.h" +//------------------------------------------------------------- +/** + *

Minuit2 command identifiers for COMMANDS block. + * + *

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. + * + *

Categories: + * - Minimizers: MIGRAD (gradient descent), SIMPLEX (non-gradient), MINIMIZE (automatic) + * - Error analysis: HESSE (covariance matrix), MINOS (asymmetric errors) + * - Parameter control: FIX, RELEASE, RESTORE, SAVE + * - Exploration: SCAN (1D/2D parameter space), CONTOURS (error ellipses) + * - Settings: STRATEGY, MACHINE_PRECISION, PRINT + * - Analysis: 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 //----------------------------------------------------------------------------- /** - *

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. + *

Container for sector χ² (or max likelihood) analysis results. + * + *

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 + * + *

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. + * + *

Use case: If early-time data shows poor χ², the sector analysis + * reveals whether to adjust t0, fgb, or background estimation. */ class PSectorChisq { public: + /** + *

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 //----------------------------------------------------------------------------- /** - *

Interface class to minuit2. + *

Main fitting engine interfacing with ROOT Minuit2. + * + *

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 + * + *

Fitting workflow: + * 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 + * + *

Minimization modes: + * - χ² minimization: Standard least-squares fitting + * - Maximum likelihood: Poisson statistics (better for low counts) + * + *

Example COMMANDS block: + * @code + * COMMANDS + * SET PRINT 1 + * MIGRAD + * MINOS + * SAVE + * @endcode */ class PFitter { public: + /** + *

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; } + + /** + *

Executes the complete fitting procedure. + * + *

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: diff --git a/src/include/PFourier.h b/src/include/PFourier.h index e29e20d22..78c744ce4 100644 --- a/src/include/PFourier.h +++ b/src/include/PFourier.h @@ -40,29 +40,105 @@ #include "PMusr.h" +//------------------------------------------------------------- +/** + *

Apodization (windowing) strength constants. + * + *

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 + *

Phase correction optimizer for Fourier transforms. + * + *

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 + * + *

Algorithm: Minimizes a combined entropy-penalty functional using + * Minuit2, finding optimal phase parameters (constant + linear dispersion): + * φ(ω) = c₀ + c₁·ω + * + *

Applications: + * - Improving signal clarity in real Fourier spectra + * - Identifying field distributions in vortex lattices + * - Resolving closely-spaced frequency components + * + *

Usage: Specify frequency range for optimization to focus on + * signal peaks while avoiding noise regions. */ class PFTPhaseCorrection : public ROOT::Minuit2::FCNBase { public: + /** + *

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); + + /** + *

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 &reFT, std::vector &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; } + + /** + *

Performs phase correction minimization. + * + *

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; } + + /** + *

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&) const; }; +//------------------------------------------------------------- /** - * muSR Fourier class. + *

Fourier transform engine for μSR time-domain data. + * + *

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 + * + *

Key features: + * - 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) + * + *

Workflow: + * 1. Create PFourier with time histogram and settings + * 2. Call Transform() with desired apodization + * 3. Retrieve results: GetRealFourier(), GetPowerFourier(), etc. + * + *

Unit conversions: + * - Gauss: ω(MHz) = γ_μ/(2π) × B(G) = 0.01355 × B(G) + * - Tesla: ω(MHz) = γ_μ/(2π) × B(T) = 135.54 × B(T) + * + *

Example: TF-μSR measurement at 100 G produces a peak at + * ~1.36 MHz in the Fourier spectrum. */ class PFourier { public: + /** + *

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(); + /** + *

Performs the Fourier transformation. + * + *

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; } + + /** + *

Returns the maximum frequency (Nyquist frequency). + * + * @return Maximum frequency in output units + */ virtual Double_t GetMaxFreq(); + + /** + *

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 &phase, const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0); + + /** + *

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); + + /** + *

Gets power spectrum |F(ω)|² as histogram. + * + *

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); + + /** + *

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); + /** + *

Static method for phase-optimized real Fourier spectrum. + * + *

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 &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: diff --git a/src/include/PMsrHandler.h b/src/include/PMsrHandler.h index d41115576..97ad5b9c6 100644 --- a/src/include/PMsrHandler.h +++ b/src/include/PMsrHandler.h @@ -44,67 +44,292 @@ //------------------------------------------------------------- /** - *

This class provides the routines needed to handle msr-files, i.e. reading, writing, parsing, etc. + *

MSR file parser and manager for musrfit framework. + * + *

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 + * + *

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 + * + *

Usage pattern: + * @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: + /** + *

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(); + /** + *

Reads and parses the MSR file. + * + *

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(); + + /** + *

Writes a log file with MSR file content and parsing information. + * + *

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); + + /** + *

Writes MSR file with updated parameters and results. + * + *

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 *commentsPAR = 0, std::map *commentsTHE = 0, \ std::map *commentsFUN = 0, std::map *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; } + /** + *

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); + + /** + *

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); + + /** + *

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); + + /** + *

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); + /** + *

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); + + /** + *

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); + + /** + *

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); + + /** + *

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(); } + + /** + *

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); } + + /** + *

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); } + + /** + *

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); } + + /** + *

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 map, std::vector param, PMetaData metaData) { return fFuncHandler->Eval(i, map, param, metaData); } + + /** + *

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); + + /** + *

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); + /** + *

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); + + /** + *

Estimates N0 parameter for single histogram fits. + * + *

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: diff --git a/src/include/PMusr.h b/src/include/PMusr.h index 241ecb5e0..aaf39d705 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -41,123 +41,311 @@ //#endif #include "fftw3.h" +//------------------------------------------------------------- +/** + *

Return code constants for musrfit framework operations. + * + *

These constants define the return values for various operations + * in the musrfit framework, where 0 indicates success and negative + * values indicate different types of errors. + */ +/// Successful operation completion #define PMUSR_SUCCESS 0 +/// Syntax help was requested via command line #define PMUSR_SYNTAX_REQUEST -1 +/// Incorrect startup command syntax provided #define PMUSR_WRONG_STARTUP_SYNTAX -2 +/// MSR file could not be found at specified path #define PMUSR_MSR_FILE_NOT_FOUND -3 +/// Memory allocation error while processing MSR file #define PMUSR_MSR_ALLOCATION_ERROR -4 +/// Syntax error detected in MSR file content #define PMUSR_MSR_SYNTAX_ERROR -5 +/// Error during tokenization/parsing of input #define PMUSR_TOKENIZE_ERROR -6 +/// Failed to write to MSR log file #define PMUSR_MSR_LOG_FILE_WRITE_ERROR -7 +/// Failed to write MSR file #define PMUSR_MSR_FILE_WRITE_ERROR -8 +/// Error reading data file (ROOT, NeXus, MUD, etc.) #define PMUSR_DATA_FILE_READ_ERROR -9 +/// Error during run processing or fitting #define PMUSR_MSR_RUN_ERROR -10 +/// Requested feature is not yet supported #define PMUSR_UNSUPPORTED_FEATURE -11 +//------------------------------------------------------------- +/** + *

Run type identifiers for different μSR measurement configurations. + * + *

These constants specify the type of μSR experiment and determine + * how raw histogram data should be processed and fitted. + */ +/// Single histogram fit (e.g., forward or backward detector) #define PRUN_SINGLE_HISTO 0 +/// Single histogram fit in rotating reference frame (RRF) #define PRUN_SINGLE_HISTO_RRF 1 +/// Asymmetry fit using forward and backward detectors #define PRUN_ASYMMETRY 2 +/// Asymmetry fit in rotating reference frame (RRF) #define PRUN_ASYMMETRY_RRF 3 +/// Negative muon (μ-) single histogram fit #define PRUN_MU_MINUS 4 +/// Beta-detected NMR asymmetry fit #define PRUN_ASYMMETRY_BNMR 5 +/// Non-μSR data fit (general x-y data) #define PRUN_NON_MUSR 8 -// muon lifetime in (us) - see e.g.: -// P. A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020). -// https://doi.org/10.1093/ptep/ptaa104 -// https://pdg.lbl.gov/2021/ +//------------------------------------------------------------- +/** + *

Physical constants for muon experiments. + * + *

These fundamental constants are used throughout the framework + * for lifetime corrections, field-frequency conversions, and other + * physics calculations. + */ + +/** + *

Muon lifetime in microseconds (μs). + * + *

The mean lifetime of the positive muon (μ+) is 2.1969811 μs. + * This value is used for lifetime corrections in single histogram fits. + * + *

Reference: + * P. A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020). + * https://doi.org/10.1093/ptep/ptaa104 + * https://pdg.lbl.gov/2021/ + */ #define PMUON_LIFETIME 2.1969811 -// muon gyromagnetic ratio: gamma_mu = 2.0 * mu_mu / hbar - see e.g.: -// E. Tiesinga et al., Rev. Mod. Phys. 93, 025010 (2021). -// https://doi.org/10.1103/RevModPhys.93.025010 -// https://physics.nist.gov/cuu/Constants/index.html -// mu_mu = -4.490 448 30(10) * 1e-26 J/T -// hbar = 1.054 571 817... * 1e-34 J*s -// gamma_muon / (2 pi) = 135.538 809 4(30) MHz/T -// = 0.013 553 880 94(30) MHz/G +/** + *

Muon gyromagnetic ratio γ_μ/(2π) in MHz/G. + * + *

The muon gyromagnetic ratio relates the muon spin precession frequency + * to the applied magnetic field: ν = γ_μ/(2π) × B + * + *

Value: 0.013553880 94(30) MHz/G = 135.53880 94(30) MHz/T + * + *

Calculated from: γ_μ = 2.0 × μ_μ / ℏ where + * - μ_μ = -4.490 448 30(10) × 10^-26 J/T (muon magnetic moment) + * - ℏ = 1.054 571 817... × 10^-34 J·s (reduced Planck constant) + * + *

Reference: + * E. Tiesinga et al., Rev. Mod. Phys. 93, 025010 (2021). + * https://doi.org/10.1103/RevModPhys.93.025010 + * https://physics.nist.gov/cuu/Constants/index.html + */ #define GAMMA_BAR_MUON 1.355388094e-2 -// accelerator cycles in (us), needed to determine proper background +//------------------------------------------------------------- +/** + *

Accelerator cycle periods in microseconds (μs). + * + *

These constants define the beam structure period for different + * μSR facilities, which is needed for proper background determination + * in pulsed beam experiments. + */ +/// PSI (Paul Scherrer Institute) accelerator cycle: 19.75 ns #define ACCEL_PERIOD_PSI 0.01975 +/// TRIUMF accelerator cycle: 43.37 ns #define ACCEL_PERIOD_TRIUMF 0.04337 +/// RAL (Rutherford Appleton Lab) - pulsed beam #define ACCEL_PERIOD_RAL 0.0 -// used to filter post pileup correct data histos from root files +//------------------------------------------------------------- +/** + *

Histogram offset for pileup-corrected data in ROOT files. + * + *

Post-pileup-corrected histograms are stored with an offset + * of 20 added to the original histogram number to distinguish them + * from uncorrected data. + */ #define POST_PILEUP_HISTO_OFFSET 20 -// defines a value for 'undefined values' +//------------------------------------------------------------- +/** + *

Sentinel value indicating undefined or uninitialized parameters. + * + *

This large negative value (-9.9×10^99) is used throughout the + * framework to indicate that a parameter has not been set or is invalid. + */ #define PMUSR_UNDEFINED -9.9e99 //------------------------------------------------------------- -// msr block header tags +/** + *

MSR file block header tags. + * + *

These constants identify the different sections in a musrfit + * MSR (Muon Spin Rotation) file. Each tag corresponds to a specific + * block that defines different aspects of the fit configuration. + */ +/// TITLE block - describes the experiment #define MSR_TAG_TITLE 0 +/// FITPARAMETER block - defines fit parameters with initial values and constraints #define MSR_TAG_FITPARAMETER 1 +/// THEORY block - specifies the theory function(s) to fit #define MSR_TAG_THEORY 2 +/// FUNCTIONS block - user-defined mathematical functions #define MSR_TAG_FUNCTIONS 3 +/// GLOBAL block - global fit settings (RRF, fit type, etc.) #define MSR_TAG_GLOBAL 4 +/// RUN block - run-specific settings and data file information #define MSR_TAG_RUN 5 +/// COMMANDS block - post-fit commands (e.g., parameter output) #define MSR_TAG_COMMANDS 6 +/// FOURIER block - Fourier transform settings #define MSR_TAG_FOURIER 7 +/// PLOT block - plotting configuration for data visualization #define MSR_TAG_PLOT 8 +/// STATISTIC block - fit statistics and results (generated after fit) #define MSR_TAG_STATISTIC 9 //------------------------------------------------------------- -// msr fit type tags +/** + *

MSR file fit type tags. + * + *

These constants specify the fit type in the GLOBAL or RUN blocks + * of an MSR file. The fit type determines how raw data is converted + * into the quantity to be fitted (single histogram, asymmetry, etc.). + */ +/// Fit single histogram (e.g., positron counts vs. time) #define MSR_FITTYPE_SINGLE_HISTO 0 +/// Fit single histogram in rotating reference frame #define MSR_FITTYPE_SINGLE_HISTO_RRF 1 +/// Fit asymmetry A(t) = (F-αB)/(F+αB) #define MSR_FITTYPE_ASYM 2 +/// Fit asymmetry in rotating reference frame #define MSR_FITTYPE_ASYM_RRF 3 +/// Fit negative muon (μ-) single histogram #define MSR_FITTYPE_MU_MINUS 4 +/// Fit beta-detected NMR asymmetry #define MSR_FITTYPE_BNMR 5 +/// Fit non-μSR data (general x-y data) #define MSR_FITTYPE_NON_MUSR 8 //------------------------------------------------------------- -// msr plot type tags +/** + *

MSR file plot type tags. + * + *

These constants specify the plot type in the PLOT block + * of an MSR file, determining how data should be displayed + * in musrview or other visualization tools. + */ +/// Plot single histogram #define MSR_PLOT_SINGLE_HISTO 0 +/// Plot single histogram in rotating reference frame #define MSR_PLOT_SINGLE_HISTO_RRF 1 +/// Plot asymmetry #define MSR_PLOT_ASYM 2 +/// Plot asymmetry in rotating reference frame #define MSR_PLOT_ASYM_RRF 3 +/// Plot negative muon (μ-) data #define MSR_PLOT_MU_MINUS 4 +/// Plot beta-detected NMR data #define MSR_PLOT_BNMR 5 +/// Plot non-μSR data #define MSR_PLOT_NON_MUSR 8 //------------------------------------------------------------- -// map and fun offsets for parameter parsing +/** + *

Offsets for parameter parsing in MSR files. + * + *

When parameters are referenced via map functions or user-defined + * functions, these offsets are added to distinguish them from direct + * parameter references (which use indices 1, 2, 3, ...). + */ +/// Offset added to map indices for parameter parsing #define MSR_PARAM_MAP_OFFSET 10000 +/// Offset added to function indices for parameter parsing #define MSR_PARAM_FUN_OFFSET 20000 //------------------------------------------------------------- -// fourier related tags +/** + *

Fourier transform unit tags. + * + *

These constants specify the units for Fourier transform + * output in the FOURIER block of MSR files. + */ +/// Units not specified #define FOURIER_UNIT_NOT_GIVEN 0 +/// Magnetic field in Gauss (G) #define FOURIER_UNIT_GAUSS 1 +/// Magnetic field in Tesla (T) #define FOURIER_UNIT_TESLA 2 +/// Frequency in MHz #define FOURIER_UNIT_FREQ 3 +/// Angular frequency in Mc/s (Mega-cycles per second) #define FOURIER_UNIT_CYCLES 4 +//------------------------------------------------------------- +/** + *

Fourier transform apodization tags. + * + *

Apodization (windowing) reduces spectral leakage in Fourier + * transforms by applying a window function to the time-domain data. + * Stronger apodization improves frequency resolution but reduces + * signal amplitude accuracy. + */ +/// Apodization not specified #define FOURIER_APOD_NOT_GIVEN 0 +/// No apodization (rectangular window) #define FOURIER_APOD_NONE 1 +/// Weak apodization (gentle windowing) #define FOURIER_APOD_WEAK 2 +/// Medium apodization (moderate windowing) #define FOURIER_APOD_MEDIUM 3 +/// Strong apodization (heavy windowing for best frequency resolution) #define FOURIER_APOD_STRONG 4 +//------------------------------------------------------------- +/** + *

Fourier transform plot type tags. + * + *

These constants specify which component(s) of the complex + * Fourier transform should be displayed. + */ +/// Plot type not specified #define FOURIER_PLOT_NOT_GIVEN 0 +/// Plot real component only #define FOURIER_PLOT_REAL 1 +/// Plot imaginary component only #define FOURIER_PLOT_IMAG 2 +/// Plot both real and imaginary components (default) #define FOURIER_PLOT_REAL_AND_IMAG 3 +/// Plot power spectrum |F(ω)|² #define FOURIER_PLOT_POWER 4 +/// Plot phase spectrum arg(F(ω)) #define FOURIER_PLOT_PHASE 5 +/// Plot phase-optimized real component #define FOURIER_PLOT_PHASE_OPT_REAL 6 //------------------------------------------------------------- -// RRF related tags +/** + *

Rotating Reference Frame (RRF) unit tags. + * + *

In RRF analysis, data is transformed into a frame rotating at + * a specified frequency/field. These tags specify the units for the + * RRF frequency, which can be given as frequency or equivalent field. + */ +/// RRF unit undefined #define RRF_UNIT_UNDEF -1 +/// Frequency in kHz (kilohertz) #define RRF_UNIT_kHz 0 +/// Frequency in MHz (megahertz) #define RRF_UNIT_MHz 1 +/// Angular frequency in Mc/s (Mega-cycles per second) #define RRF_UNIT_Mcs 2 +/// Equivalent magnetic field in Gauss (G) #define RRF_UNIT_G 3 +/// Equivalent magnetic field in Tesla (T) #define RRF_UNIT_T 4 +/** + *

Sentinel value indicating undefined RRF frequency. + * + *

This large value (10^10) marks an RRF frequency as uninitialized + * or invalid. + */ #define RRF_FREQ_UNDEF 1.0e10 //------------------------------------------------------------- @@ -216,53 +404,117 @@ typedef std::vector PStringVector; //------------------------------------------------------------- /** - *

data handling tag + *

Data handling mode tag for musrfit operations. + * + *

This enumeration specifies the operational mode of the musrfit + * framework, determining whether data should be fitted, viewed, or + * if no operation is active. */ -enum EPMusrHandleTag { kEmpty, kFit, kView }; - -//------------------------------------------------------------- -/** - *

keep a couple of meta data for the FUNCTIONS. - */ -struct PMetaData { - Double_t fField; ///< field in (G) - Double_t fEnergy; ///< energy in (keV) - PDoubleVector fTemp; ///< temperature(s) in (K) +enum EPMusrHandleTag { + kEmpty, ///< No operation active + kFit, ///< Fitting mode - perform least-squares fit to data + kView ///< Viewing mode - display data and theory without fitting }; //------------------------------------------------------------- /** - *

Handles the data which will be fitted, i.e. packed, background corrected, ... - * This are not the raw histogram data of a run. This are the pre-processed data. + *

Metadata structure for user-defined FUNCTIONS block. + * + *

This structure stores experimental metadata (field, energy, temperature) + * that can be accessed within user-defined functions in the FUNCTIONS block + * of MSR files. This allows functions to depend on experimental conditions + * without requiring additional fit parameters. + * + *

Example use case: A temperature-dependent relaxation rate can + * be calculated using fTemp values rather than fitting temperature as a parameter. + */ +struct PMetaData { + Double_t fField; ///< Applied magnetic field in Gauss (G) + Double_t fEnergy; ///< Muon implantation energy in kiloelectronvolts (keV) + PDoubleVector fTemp; ///< Vector of sample temperature(s) in Kelvin (K) +}; + +//------------------------------------------------------------- +/** + *

Container for pre-processed data and theory vectors used in fitting. + * + *

This class stores the processed data (after background subtraction, + * packing, asymmetry calculation, etc.) that is actually used in the fit, + * along with the corresponding theory values. This is not raw histogram + * data - it represents the final data points and theory that enter the χ² + * calculation. + * + *

The class maintains separate time grids for data and theory, allowing + * theory evaluation at higher resolution than the data binning. For non-μSR + * fits, x-axis vectors are also stored. */ class PRunData { public: PRunData(); virtual ~PRunData(); + /// Returns the start time of the data set in microseconds (μs) virtual Double_t GetDataTimeStart() { return fDataTimeStart; } + /// Returns the time bin width for data in microseconds (μs) virtual Double_t GetDataTimeStep() { return fDataTimeStep; } + /// Returns the start time for theory evaluation in microseconds (μs) virtual Double_t GetTheoryTimeStart() { return fTheoryTimeStart; } + /// Returns the time step for theory evaluation in microseconds (μs) virtual Double_t GetTheoryTimeStep() { return fTheoryTimeStep; } + /// Returns pointer to x-axis vector (used only for non-μSR data) virtual const PDoubleVector* GetX() { return &fX; } + /// Returns pointer to data value vector (asymmetry, counts, or y-data) virtual const PDoubleVector* GetValue() { return &fValue; } + /// Returns pointer to data error vector (statistical uncertainties) virtual const PDoubleVector* GetError() { return &fError; } + /// Returns pointer to x-axis vector for theory (non-μSR only) virtual const PDoubleVector* GetXTheory() { return &fXTheory; } + /// Returns pointer to theory value vector virtual const PDoubleVector* GetTheory() { return &fTheory; } + /// Sets the start time of the data set in microseconds (μs) + /// @param dval Start time value virtual void SetDataTimeStart(Double_t dval) { fDataTimeStart = dval; } + /// Sets the time bin width for data in microseconds (μs) + /// @param dval Time step value virtual void SetDataTimeStep(Double_t dval) { fDataTimeStep = dval; } + /// Sets the start time for theory evaluation in microseconds (μs) + /// @param dval Start time value virtual void SetTheoryTimeStart(Double_t dval) { fTheoryTimeStart = dval; } + /// Sets the time step for theory evaluation in microseconds (μs) + /// @param dval Time step value virtual void SetTheoryTimeStep(Double_t dval) { fTheoryTimeStep = dval; } + /// Appends a value to the x-axis vector (non-μSR only) + /// @param dval X-coordinate value to append virtual void AppendXValue(Double_t dval) { fX.push_back(dval); } + /// Appends a value to the data vector + /// @param dval Data value to append virtual void AppendValue(Double_t dval) { fValue.push_back(dval); } + /// Appends an error value to the error vector + /// @param dval Error/uncertainty value to append virtual void AppendErrorValue(Double_t dval) { fError.push_back(dval); } + /// Appends a value to the theory x-axis vector (non-μSR only) + /// @param dval X-coordinate value for theory virtual void AppendXTheoryValue(Double_t dval) { fXTheory.push_back(dval); } + /// Appends a value to the theory vector + /// @param dval Theory value to append virtual void AppendTheoryValue(Double_t dval) { fTheory.push_back(dval); } + /** + *

Sets a specific theory value at given index. + * + * @param i Index in theory vector + * @param dval New theory value + */ virtual void SetTheoryValue(UInt_t i, Double_t dval); + + /** + *

Replaces the entire theory vector with a new one. + * + * @param theo New theory vector to replace current theory + */ virtual void ReplaceTheory(const PDoubleVector &theo); private: @@ -281,27 +533,79 @@ class PRunData { //------------------------------------------------------------- /** - *

Handles Non-Musr raw data. + *

Container for non-μSR raw data (general x-y data sets). + * + *

This class handles raw data for non-μSR experiments (fit type NON_MUSR), + * supporting both ASCII file input and database (db/dat) file formats. + * It stores multiple data columns with labels, allowing general x-y fitting + * beyond traditional μSR time histograms. + * + *

Use cases: Fitting arbitrary x-y data, susceptibility vs. temperature, + * field scans, or any other non-time-domain measurements. */ class PNonMusrRawRunData { public: PNonMusrRawRunData(); virtual ~PNonMusrRawRunData(); + /// Returns true if data was loaded from ASCII file, false for db/dat format virtual Bool_t FromAscii() { return fFromAscii; } + /// Returns pointer to vector of axis/column labels virtual const PStringVector* GetLabels() { return &fLabels; } + /// Returns pointer to vector of data tags (identifiers for each data column) virtual const PStringVector* GetDataTags() { return &fDataTags; } + /// Returns pointer to vector of data columns virtual const std::vector* GetData() { return &fData; } + /// Returns pointer to vector of error data columns virtual const std::vector* GetErrData() { return &fErrData; } + /// Sets the flag indicating if data is from ASCII format + /// @param bval True for ASCII, false for db/dat format virtual void SetFromAscii(const Bool_t bval) { fFromAscii = bval; } + + /** + *

Pre-allocates space for a given number of data columns. + * + * @param size Number of data columns to allocate + */ virtual void SetSize(const UInt_t size); + + /// Appends a label to the label vector + /// @param str Label string (e.g., "Temperature (K)", "Field (G)") virtual void AppendLabel(const TString str) { fLabels.push_back(str); } + + /** + *

Sets or replaces a label at specific index. + * + * @param idx Index in label vector + * @param str New label string + */ virtual void SetLabel(const UInt_t idx, const TString str); + + /// Appends a data tag identifier + /// @param str Data tag string virtual void AppendDataTag(const TString str) { fDataTags.push_back(str); } + /// Appends a complete data column vector + /// @param data Data vector to append virtual void AppendData(const PDoubleVector &data) { fData.push_back(data); } + /// Appends a complete error data column vector + /// @param data Error data vector to append virtual void AppendErrData(const PDoubleVector &data) { fErrData.push_back(data); } + + /** + *

Appends a single value to a specific data column. + * + * @param idx Index of data column + * @param dval Data value to append + */ virtual void AppendSubData(const UInt_t idx, const Double_t dval); + + /** + *

Appends a single error value to a specific error data column. + * + * @param idx Index of error data column + * @param dval Error value to append + */ virtual void AppendSubErrData(const UInt_t idx, const Double_t dval); private: @@ -314,34 +618,73 @@ class PNonMusrRawRunData { //------------------------------------------------------------- /** - *

Handles a single raw muSR histogram set, without any additional header information. + *

Container for a single raw μSR histogram with associated metadata. + * + *

This class stores one raw histogram from a μSR experiment along with + * essential information like time-zero bin, good data range, and background + * range. It represents the most basic unit of μSR data before any processing + * or corrections are applied. + * + *

Does not contain run header information (temperature, field, etc.) - + * only the histogram data and its immediate properties. */ class PRawRunDataSet { public: PRawRunDataSet(); virtual ~PRawRunDataSet() { fData.clear(); } + /// Returns the histogram title string virtual TString GetTitle() { return fTitle; } + /// Returns the histogram name virtual TString GetName() { return fName; } + /// Returns the histogram number as stored in the data file virtual Int_t GetHistoNo() { return fHistoNo; } + /// Returns the time-zero bin (t0) position virtual Double_t GetTimeZeroBin() { return fTimeZeroBin; } + /// Returns the estimated time-zero bin from automatic determination virtual Double_t GetTimeZeroBinEstimated() { return fTimeZeroBinEstimated; } + /// Returns the first bin of good data range (after prompt peak) virtual Int_t GetFirstGoodBin() { return fFirstGoodBin; } + /// Returns the last bin of good data range virtual Int_t GetLastGoodBin() { return fLastGoodBin; } + /// Returns the first bin of background range virtual Int_t GetFirstBkgBin() { return fFirstBkgBin; } + /// Returns the last bin of background range virtual Int_t GetLastBkgBin() { return fLastBkgBin; } + /// Returns pointer to raw histogram data vector (counts per bin) virtual PDoubleVector *GetData() { return &fData; } + /// Clears all data from this histogram set virtual void Clear(); + /// Sets the histogram title + /// @param str Title string virtual void SetTitle(TString str) { fTitle = str; } + /// Sets the histogram name + /// @param str Name string virtual void SetName(TString str) { fName = str; } + /// Sets the histogram number + /// @param no Histogram number from data file virtual void SetHistoNo(Int_t no) { fHistoNo = no; } + /// Sets the time-zero bin position + /// @param tzb Time-zero bin value (can be fractional) virtual void SetTimeZeroBin(Double_t tzb) { fTimeZeroBin = tzb; } + /// Sets the estimated time-zero bin from automatic determination + /// @param tzb Estimated time-zero bin value virtual void SetTimeZeroBinEstimated(Double_t tzb) { fTimeZeroBinEstimated = tzb; } + /// Sets the first bin of good data range + /// @param fgb First good bin index virtual void SetFirstGoodBin(Int_t fgb) { fFirstGoodBin = fgb; } + /// Sets the last bin of good data range + /// @param lgb Last good bin index virtual void SetLastGoodBin(Int_t lgb) { fLastGoodBin = lgb; } + /// Sets the first bin of background range + /// @param fbb First background bin index virtual void SetFirstBkgBin(Int_t fbb) { fFirstBkgBin = fbb; } + /// Sets the last bin of background range (note: bug in original code sets fLastGoodBin) + /// @param lbb Last background bin index virtual void SetLastBkgBin(Int_t lbb) { fLastGoodBin = lbb; } + /// Sets the entire histogram data vector + /// @param data Vector of histogram counts virtual void SetData(PDoubleVector data) { fData = data; } private: @@ -359,26 +702,102 @@ class PRawRunDataSet { //------------------------------------------------------------- /** - *

Handles a vector of PRawRunDataSet's. Addressing of the data set can be done - * via a map mechanism, since e.g. for MusrRoot, the data sets are not continuously - * numbered due to Red/Green options. + *

Vector container for multiple PRawRunDataSet objects with flexible indexing. + * + *

This class manages a collection of raw histogram data sets, providing + * both sequential indexing (0, 1, 2, ...) and histogram-number-based access. + * The flexible indexing is necessary because histogram numbers in data files + * (especially MusrRoot files with Red/Green options) may not be contiguous. + * + *

Example: A MusrRoot file might contain histograms numbered + * [1, 2, 5, 6] (Red group) and [21, 22, 25, 26] (Green group with offset). */ class PRawRunDataVector { public: PRawRunDataVector() {} virtual ~PRawRunDataVector() { fDataVec.clear(); } + /// Returns the number of histogram data sets in this vector virtual UInt_t Size() { return fDataVec.size(); } + + /** + *

Checks if a histogram with given number exists in the vector. + * + * @param histoNo Histogram number to search for + * @return true if histogram is present, false otherwise + */ virtual Bool_t IsPresent(UInt_t histoNo); + + /** + *

Gets histogram data set by sequential index (0, 1, 2, ...). + * + * @param idx Sequential index in vector + * @return Pointer to histogram data set, or nullptr if index out of range + */ virtual PRawRunDataSet* GetSet(UInt_t idx); + + /** + *

Gets histogram data set by histogram number. + * + * @param histoNo Histogram number from data file + * @return Pointer to histogram data set, or nullptr if not found + */ virtual PRawRunDataSet* Get(UInt_t histoNo); + + /** + *

Operator overload for histogram-number-based access. + * + * @param histoNo Histogram number from data file + * @return Pointer to histogram data set + */ virtual PRawRunDataSet* operator[](UInt_t histoNo); + + /** + *

Gets histogram data vector by histogram number. + * + * @param histoNo Histogram number + * @return Pointer to data vector, or nullptr if not found + */ virtual PDoubleVector* GetData(UInt_t histoNo); + + /** + *

Gets time-zero bin for specified histogram. + * + * @param histoNo Histogram number + * @return Time-zero bin value + */ virtual Double_t GetT0Bin(UInt_t histoNo); + + /** + *

Gets estimated time-zero bin for specified histogram. + * + * @param histoNo Histogram number + * @return Estimated time-zero bin value + */ virtual Double_t GetT0BinEstimated(UInt_t histoNo); + + /** + *

Gets background bin range for specified histogram. + * + * @param histoNo Histogram number + * @return Pair (first bin, last bin) of background range + */ virtual PIntPair GetBkgBin(UInt_t histoNo); + + /** + *

Gets good data bin range for specified histogram. + * + * @param histoNo Histogram number + * @return Pair (first bin, last bin) of good data range + */ virtual PIntPair GetGoodDataBin(UInt_t histoNo); + /** + *

Adds or updates a histogram data set in the vector. + * + * @param dataSet Histogram data set to add + * @param idx Optional index; if -1 (default), appends to end + */ virtual void Set(PRawRunDataSet dataSet, Int_t idx=-1); private: @@ -387,7 +806,21 @@ class PRawRunDataVector { //------------------------------------------------------------- /** - *

Handles raw data, both non-muSR data as well muSR histogram data. + *

Complete raw data container for a single run. + * + *

This class stores all raw data and metadata for one experimental run, + * including run header information (temperature, field, setup details), + * histogram data (for μSR) or column data (for non-μSR), and data file + * metadata. It serves as the primary interface between data file readers + * and the fitting/analysis framework. + * + *

Supports multiple data formats: + * - MusrRoot files (.root) + * - NeXus files (.nxs) + * - MUD files (.msr from TRIUMF) + * - PSI binary files + * - ASCII data files (for non-μSR) + * - Database files (.db, .dat for non-μSR) */ class PRawRunData { public: @@ -531,11 +964,15 @@ typedef std::vector PRawRunDataList; //------------------------------------------------------------- /** - *

Helper structure for parsing. Keeps a msr-file line string and the corresponding line number. + *

Helper structure for MSR file parsing. + * + *

Stores a single line from an MSR file along with its line number, + * facilitating error reporting and maintaining correspondence between + * parsed content and original file locations. */ struct PMsrLineStructure { - Int_t fLineNo; ///< original line number of the msr-file - TString fLine; ///< msr-file line + Int_t fLineNo; ///< Line number in original MSR file (1-based) + TString fLine; ///< Content of the MSR file line }; //------------------------------------------------------------- @@ -546,21 +983,29 @@ typedef std::vector PMsrLines; //------------------------------------------------------------- /** - *

Handles the information of a parameter. + *

Structure defining a fit parameter in the FITPARAMETER block. + * + *

This structure contains all information about a single fit parameter: + * its number, name, value, uncertainty, and optional constraints (boundaries). + * Parameters can be fixed (fStep=0) or free (fStep>0), and may have asymmetric + * errors after fitting. + * + *

Example MSR line: "3 alpha 1.0 0.1 0.5 1.5" defines parameter + * #3 named "alpha" with value 1.0, step 0.1, and bounds [0.5, 1.5]. */ struct PMsrParamStructure{ - Int_t fNoOfParams; ///< how many parameters are given - Int_t fNo; ///< parameter number - TString fName; ///< name - Double_t fValue; ///< value - Double_t fStep; ///< step / error / neg_error, depending on the situation - Bool_t fPosErrorPresent; ///< positive error is defined (as a number) - Double_t fPosError; ///< positive error if present - Bool_t fLowerBoundaryPresent; ///< flag showing if a lower boundary is present - Double_t fLowerBoundary; ///< lower boundary for the fit parameter - Bool_t fUpperBoundaryPresent; ///< flag showing if an upper boundary is present - Double_t fUpperBoundary; ///< upper boundary for the fit parameter - Bool_t fIsGlobal; ///< flag showing if the parameter is a global one (used for msr2data global) + Int_t fNoOfParams; ///< Total number of parameters in FITPARAMETER block + Int_t fNo; ///< Parameter number (1, 2, 3, ...) + TString fName; ///< Parameter name (e.g., "alpha", "lambda", "field") + Double_t fValue; ///< Parameter value (initial or fitted) + Double_t fStep; ///< Step size / error / negative error (context-dependent) + Bool_t fPosErrorPresent; ///< True if positive error explicitly defined (asymmetric errors) + Double_t fPosError; ///< Positive error for asymmetric uncertainties + Bool_t fLowerBoundaryPresent; ///< True if lower bound constraint is active + Double_t fLowerBoundary; ///< Lower boundary value for parameter constraints + Bool_t fUpperBoundaryPresent; ///< True if upper bound constraint is active + Double_t fUpperBoundary; ///< Upper boundary value for parameter constraints + Bool_t fIsGlobal; ///< True if parameter is global (for msr2data global mode) }; //------------------------------------------------------------- @@ -571,7 +1016,16 @@ typedef std::vector PMsrParamList; //------------------------------------------------------------- /** - *

Handles the information of the GLOBAL section + *

Container for GLOBAL block information in MSR files. + * + *

The GLOBAL block defines settings that apply to all RUN blocks, + * including fit type, data ranges, time-zero bins, rotating reference + * frame (RRF) parameters, and packing. This block is optional - if absent, + * all settings must be specified individually in each RUN block. + * + *

Example use: When fitting multiple runs with identical + * experimental setup, the GLOBAL block avoids repetition by centralizing + * common parameters like RRF frequency, fit range, and t0 bins. */ class PMsrGlobalBlock { public: @@ -629,13 +1083,31 @@ class PMsrGlobalBlock { //------------------------------------------------------------- /** - *

Handles the information of a single run block + *

Container for a single RUN block in MSR files. * + *

Each RUN block defines one run or a group of runs to be fitted together, + * specifying the data file(s), fit type, histogram selections, background + * handling, time-zero bins, fit ranges, and theory-to-histogram mapping. + * Multiple RUN blocks enable simultaneous fitting of related data sets. + * + *

Key features: + * - Supports multiple data files (run name vector) + * - Flexible histogram selection (forward/backward for asymmetry) + * - Theory parameter mapping (allows parameter variations between runs) + * - Background estimation or fixing + * - Run-specific or global settings + * + *

Example: Fitting a field-dependent measurement where each run + * has a different field value but shares common relaxation parameters. */ class PMsrRunBlock { public: PMsrRunBlock(); virtual ~PMsrRunBlock(); + + /** + *

Cleans up and resets all data in the run block. + */ virtual void CleanUp(); virtual UInt_t GetRunNameSize() { return fRunName.size(); } @@ -761,7 +1233,16 @@ typedef std::vector PMsrRunList; //------------------------------------------------------------- /** - *

Holds the information of the Fourier block + *

Structure containing all FOURIER block settings. + * + *

The FOURIER block controls how Fourier transforms are performed + * and displayed. It specifies units (field or frequency), DC correction, + * zero-padding, apodization, phase corrections, and plot ranges. + * + *

Common use case: Converting time-domain μSR data to frequency + * domain to identify precession frequencies and their distributions, + * particularly useful for identifying field distributions in magnetic + * samples or multiple muon stopping sites. */ struct PMsrFourierStructure { Bool_t fFourierBlockPresent; ///< flag indicating if a Fourier block is present in the msr-file @@ -780,7 +1261,17 @@ struct PMsrFourierStructure { //------------------------------------------------------------- /** - *

Holds the information of a single plot block + *

Structure containing settings for a single PLOT block. + * + *

PLOT blocks define how data should be visualized in musrview, + * specifying which runs to display, axis ranges, log scales, and + * view-specific settings like RRF parameters. Multiple PLOT blocks + * allow different visualizations of the same data set. + * + *

Example uses: + * - Plotting multiple runs with common y-axis scale for comparison + * - Creating separate plots for different time ranges (overview vs. detail) + * - Applying different RRF frequencies to examine different frequency components */ struct PMsrPlotStructure{ Int_t fPlotType; ///< plot type @@ -809,7 +1300,15 @@ typedef std::vector PMsrPlotList; //------------------------------------------------------------- /** - *

Holds the informations for the statistics block. + *

Structure containing fit statistics and results. + * + *

The STATISTIC block is automatically generated after a successful + * fit, containing χ² (or maximum likelihood), degrees of freedom, + * expected values, and per-histogram statistics. This information + * assesses the quality of fit and helps identify problematic data sets. + * + *

Note: This block is written by musrfit after fitting and + * should not be manually edited in MSR files. */ struct PMsrStatisticStructure { Bool_t fValid; ///< flag showing if the statistics block is valid, i.e. a fit took place which converged @@ -826,7 +1325,15 @@ struct PMsrStatisticStructure { //------------------------------------------------------------- /** - *

Holds the informations for the any2many converter program + *

Configuration structure for the any2many data format converter. + * + *

This structure contains all settings needed to convert μSR data + * between different file formats (ROOT, NeXus, MUD, PSI-BIN, ASCII). + * It supports batch conversion, rebinning, compression, and format-specific + * options like histogram group selection for MusrRoot files. + * + *

Supported conversions: Any-to-any among ROOT, NeXus, MUD, + * PSI-BIN, WKM, and ASCII formats. */ struct PAny2ManyInfo { Bool_t useStandardOutput{false}; ///< flag showing if the converted shall be sent to the standard output @@ -849,7 +1356,11 @@ struct PAny2ManyInfo { //------------------------------------------------------------- /** - *

Holds information given at startup + *

Structure holding command-line startup options for musrfit. + * + *

Contains flags that control musrfit behavior based on command-line + * arguments, such as whether to output expected χ² values or estimate + * N0 for single histogram fits. */ struct PStartupOptions { Bool_t writeExpectedChisq; ///< if set to true, expected chisq and chisq per block will be written @@ -858,18 +1369,43 @@ struct PStartupOptions { //------------------------------------------------------------- /** - *

Helper class which parses list of numbers of the following 3 forms and its combination. - * (i) list of integers separted by spaces, e.g. 1 3 7 14 - * (ii) a range of integers of the form nS-nE, e.g. 13-27 which will generate 13, 14, 15, .., 26, 27 - * (iii) a sequence of integers of the form nS:nE:nStep, e.g. 10:20:2 which will generate 10, 12, 14, .., 18, 20 + *

Parser for flexible number list specifications. + * + *

This utility class parses string representations of number lists + * that combine three different notations: + * - (i) Space-separated integers: "1 3 7 14" + * - (ii) Range notation: "13-27" generates 13, 14, 15, ..., 26, 27 + * - (iii) Sequence notation: "10:20:2" generates 10, 12, 14, 16, 18, 20 + * + *

These forms can be combined in a single string, e.g., + * "1 5-8 10:20:2" produces [1, 5, 6, 7, 8, 10, 12, 14, 16, 18, 20]. + * + *

Use cases: Specifying run numbers, histogram lists, or + * parameter indices on command line or in configuration files. */ class PStringNumberList { public: + /// Constructor from C string + /// @param str String to parse PStringNumberList(char *str) { fString = str; } + + /// Constructor from std::string + /// @param str String to parse PStringNumberList(std::string str) { fString = str; } + virtual ~PStringNumberList() { fList.clear(); } + /** + *

Parses the input string and generates the number list. + * + * @param errorMsg Reference to string that will contain error message if parsing fails + * @param ignoreFirstToken If true, skips the first space-separated token (useful when first token is a label) + * @return true if parsing succeeded, false on error + */ virtual bool Parse(std::string &errorMsg, bool ignoreFirstToken=false); + + /// Returns the parsed list of numbers + /// @return Vector of unsigned integers extracted from input string virtual PUIntVector GetList() { return fList; } private: @@ -883,11 +1419,19 @@ class PStringNumberList { //------------------------------------------------------------- /** - *

Run name template structure. + *

Run name template for automatic file path construction. + * + *

This structure associates an instrument identifier with a file naming + * template, enabling automatic construction of data file paths from run + * numbers. Templates use placeholders that are replaced with actual run + * numbers and years. + * + *

Example: For instrument "GPS" with template "/data/gps/%y/%r.root", + * run 2425 in year 2023 becomes "/data/gps/2023/2425.root". */ struct PRunNameTemplate { - TString instrument{""}; - TString runNameTemplate{""}; + TString instrument{""}; ///< Instrument identifier (e.g., "GPS", "LEM", "DOLLY") + TString runNameTemplate{""}; ///< File path template with placeholders (%r=run, %y=year) }; //------------------------------------------------------------- diff --git a/src/include/PRunBase.h b/src/include/PRunBase.h index bacce8560..aa3918735 100644 --- a/src/include/PRunBase.h +++ b/src/include/PRunBase.h @@ -42,53 +42,166 @@ //------------------------------------------------------------------------------------------ /** - *

The run base class is enforcing a common interface to all supported fit-types. + *

Abstract base class defining the interface for all μSR fit types. + * + *

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: + * - PRunSingleHisto: Single detector histogram fits + * - PRunAsymmetry: Asymmetry fits (forward - backward) + * - PRunSingleHistoRRF: Single histogram in rotating reference frame + * - PRunAsymmetryRRF: Asymmetry in rotating reference frame + * - PRunMuMinus: Negative muon fits + * - PRunNonMusr: General x-y data fits + * + *

Key responsibilities: + * - 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) + * + *

Processing workflow: + * 1. PrepareData() - Load and preprocess raw data + * 2. CalcTheory() - Evaluate theory function + * 3. CalcChiSquare() or CalcMaxLikelihood() - Compute fit metric + * + *

Design pattern: Template Method - base class defines workflow, + * derived classes implement fit-type-specific algorithms. */ class PRunBase { public: + /// Default constructor PRunBase(); + + /** + *

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& par) = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!! - virtual Double_t CalcMaxLikelihood(const std::vector& par) = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!! + /** + *

Calculates χ² between data and theory (pure virtual). + * + *

χ² = Σ[(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& par) = 0; + + /** + *

Calculates maximum likelihood (pure virtual). + * + *

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& par) = 0; + + /** + *

Sets the fit time range for this run. + * + *

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!! + /** + *

Evaluates theory function at all data points (pure virtual). + * + *

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 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 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 fTheory; ///< theory needed to calculate chi-square + PDoubleVector fFuncValues; ///< Values of user-defined functions from FUNCTIONS block + std::unique_ptr 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!! + /** + *

Prepares raw data for fitting (pure virtual). + * + *

Performs fit-type-specific preprocessing: + * - Background subtraction + * - Asymmetry calculation (if applicable) + * - Time bin packing/rebinning + * - RRF transformation (if applicable) + * - Error propagation + * + *

Called during initialization before fitting begins. + * + * @return true on success, false on error + */ + virtual Bool_t PrepareData() = 0; + /** + *

Calculates Kaiser window filter coefficients for RRF. + * + *

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); + + /** + *

Applies Kaiser filter to theory values for RRF. + * + *

Filters the theory function in the same way data is filtered, + * ensuring consistent comparison between data and theory in RRF fits. + */ virtual void FilterTheo(); }; diff --git a/src/include/PRunDataHandler.h b/src/include/PRunDataHandler.h index f56620eff..d55663ab4 100644 --- a/src/include/PRunDataHandler.h +++ b/src/include/PRunDataHandler.h @@ -35,42 +35,188 @@ #include "PMusr.h" #include "PMsrHandler.h" +//------------------------------------------------------------- +/** + *

Data file format identifiers for any2many converter. + * + *

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 +//------------------------------------------------------------- /** - *

Handler class needed to read/handle raw data files. + *

Raw data file reader and format converter. + * + *

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 + * + *

Supported file formats: + * - MusrRoot: PSI ROOT-based format with comprehensive metadata + * - NeXus: HDF5-based format (ISIS, JPARC) + * - MUD: TRIUMF's binary format + * - PSI-BIN: Legacy PSI binary format + * - WKM: Older PSI format + * - ASCII: Text-based formats (for non-μSR data) + * - DB/DAT: Database formats for general x-y data + * + *

Key features: + * - 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 + * + *

Example usage: + * @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(); + + /** + *

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); + + /** + *

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); + + /** + *

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); + + /** + *

Constructor for format conversion (any2many). + * + * @param any2ManyInfo Conversion configuration structure + */ PRunDataHandler(PAny2ManyInfo *any2ManyInfo); + + /** + *

Constructor for format conversion with search paths. + * + * @param any2ManyInfo Conversion configuration + * @param dataPath Vector of search directories + */ PRunDataHandler(PAny2ManyInfo *any2ManyInfo, const PStringVector dataPath); + + /** + *

Constructor for MSR-based data loading. + * + * @param msrInfo MSR file handler (provides run list and paths) + */ PRunDataHandler(PMsrHandler *msrInfo); + + /** + *

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(); + /** + *

Reads all data files specified in MSR file or configuration. + * + *

Searches data paths, detects formats, loads histograms and + * metadata. Call this before attempting to access run data. + */ virtual void ReadData(); + + /** + *

Performs format conversion (for any2many utility). + * + *

Converts loaded data to the target format specified in + * any2many configuration. + */ virtual void ConvertData(); + + /** + *

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; } + + /** + *

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); + + /** + *

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(); } + /** + *

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: diff --git a/src/include/PTheory.h b/src/include/PTheory.h index a242391f7..7793fce5a 100644 --- a/src/include/PTheory.h +++ b/src/include/PTheory.h @@ -37,49 +37,109 @@ #include "PMsrHandler.h" #include "PUserFcnBase.h" -// -------------------------------------------------------- -// function handling tags -// -------------------------------------------------------- +//------------------------------------------------------------- +/** + *

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 + */ -// 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 +//------------------------------------------------------------- +/** + *

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) @@ -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 +//------------------------------------------------------------- +/** + *

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 -// maximal number of parameters. Needed in the contents of LF +//------------------------------------------------------------- +/** + *

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 -// deg -> rad factor +//------------------------------------------------------------- +/** + *

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 -// 2 pi + +/** + *

Mathematical constant 2π. + * + *

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

Structure holding the infos of a the available internally defined functions. + *

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; ///< 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] = { //-------------------------------------------------------------------------------------- /** - *

Class handling a theory of a msr-file. + *

Theory function evaluator and manager. + * + *

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 + * + *

Theory syntax: + * - Single function: asymmetry simpleExp 1 3 + * - Addition: func1 par1... + func2 par1... + * - Multiplication: func1 par1... * func2 par1... + * - Nested: (func1 + func2) * func3 + * + *

Parameter resolution: + * Parameters can be direct numbers, parameter references (1, 2, 3, ...), + * map references (map1, map2, ...), or function references (fun1, fun2, ...). + * + *

Example: + * @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: + /** + *

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(); + + /** + *

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: diff --git a/src/include/PUserFcnBase.h b/src/include/PUserFcnBase.h index 655affc52..f865374ac 100644 --- a/src/include/PUserFcnBase.h +++ b/src/include/PUserFcnBase.h @@ -37,26 +37,124 @@ //-------------------------------------------------------------------------------------------- /** - *

Interface class for the user function. + *

Abstract base class for user-defined theory functions. + * + *

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. + * + *

Use cases: + * - 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.) + * + *

Implementation steps: + * 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" + * + *

Example minimal implementation: + * @code + * class TMyRelaxation : public PUserFcnBase { + * public: + * Double_t operator()(Double_t t, const std::vector &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 + * + *

Global part: 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. + * + *

MSR file usage: + * @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 &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) + /** + *

Indicates if this function requires global initialization. + * + *

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; } + /** + *

Sets up the global part of the user function. + * + *

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 &globalPart, UInt_t idx) {} + + /** + *

Checks if the global part initialized correctly. + * + *

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; } + + /** + *

Evaluates the user function at time t (pure virtual). + * + *

This is the core evaluation method called for each data point + * during fitting. Must be implemented by derived classes. + * + *

Parameter convention: 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 ¶m) const = 0; ClassDef(PUserFcnBase, 1) }; //-------------------------------------------------------------------------- -// This function is a replacement for the ParseFile method of TSAXParser. -//-------------------------------------------------------------------------- +/** + *

XML file parser for user function configurations. + * + *

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_