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

improve the doxygen docu with the help of Claude AI.
All checks were successful
Build and Deploy Documentation / build-and-deploy (push) Successful in 17s
Details

Browse Source
This commit is contained in:
suter_a
2025-11-12 20:32:20 +01:00
parent 4ab01f4faf
commit 25024aae6b
4 changed files with 623 additions and 63 deletions
Download Patch File Download Diff File Expand all files Collapse all files

123
src/include/PFindRun.h
View File

@@ -32,25 +32,134 @@
#include "PMusr.h"
//--------------------------------------------------------------------------
/**
* @brief PFindRun - Locates muSR data files using template-based path resolution.
*
* This class searches for muSR run data files across multiple paths using
* configurable templates that encode instrument naming conventions. It supports
* various file formats (ROOT, NeXus, PSI-BIN, PSI-MDU, MUD, WKM) and handles
* year/run number substitution in file paths.
*
* The template system uses placeholders:
* - %yyyy% : 4-digit year (e.g., 2023)
* - %yy% : 2-digit year (e.g., 23)
* - %rr...r% : Run number with varying digits (%rr%, %rrr%, up to %rrrrrrrrr%)
*
* @par Example Usage:
* @code
* PStringVector paths = {"/data/gps", "/data/lem"};
* PRunNameTemplateList templates;
* PRunNameTemplate gpsTemplate;
* gpsTemplate.instrument = "GPS";
* gpsTemplate.runNameTemplate = "%yyyy%/%rrrrr%.root";
* templates.push_back(gpsTemplate);
*
* PFindRun finder(paths, templates, "GPS", 2023, 2425, "MusrRoot");
* if (finder.FoundPathName()) {
* TString fullPath = finder.GetPathName();
* // fullPath = "/data/gps/2023/02425.root"
* }
* @endcode
*
* @see PRunNameTemplate
* @see PRunNameTemplateList
*/
class PFindRun {
public:
//----------------------------------------------------------------------
/**
* @brief Default constructor - Creates instance without search parameters.
*
* Initializes the finder with paths and templates but no specific run to search.
* Use the full constructor to perform automatic searches.
*
* @param path Vector of directory paths to search
* @param runNameTemplateList List of template patterns for different instruments
*/
PFindRun(const PStringVector path, const PRunNameTemplateList runNameTemplateList);
//----------------------------------------------------------------------
/**
* @brief Full constructor - Creates instance and prepares for file search.
*
* Initializes the finder with all parameters needed to locate a specific run file.
* Call FoundPathName() after construction to perform the actual search.
*
* @param path Vector of directory paths to search
* @param runNameTemplateList List of template patterns for different instruments
* @param instrument Instrument name (must match a template entry, e.g., "GPS", "LEM")
* @param year Run year (e.g., 2023)
* @param run Run number (e.g., 2425)
* @param file_format Optional file format filter: "MusrRoot"/"ROOT", "NeXus",
* "PSI-BIN", "PSI-MDU", "MUD", "WKM". Empty string matches any format.
*/
PFindRun(const PStringVector path, const PRunNameTemplateList runNameTemplateList,
const TString &instrument, const UInt_t year, const UInt_t run, const TString file_format="");
//----------------------------------------------------------------------
/**
* @brief Searches for the run file using configured templates and paths.
*
* Iterates through all paths containing the instrument name, applies matching
* templates, and checks filesystem for file existence. If a file format is
* specified, only files with matching extensions are considered.
*
* @return true if file was found, false otherwise
*
* @par Search Algorithm:
* 1. Filter paths containing instrument name
* 2. For each matching path, try all templates for that instrument
* 3. Substitute year/run placeholders to create full path
* 4. Check if file exists on filesystem
* 5. If file_format specified, verify extension matches
*
* @note After successful search, use GetPathName() to retrieve the full path.
*/
Bool_t FoundPathName();
//----------------------------------------------------------------------
/**
* @brief Returns the full path to the found run file.
*
* @return Full filesystem path including filename and extension, or empty
* string if no file was found (call FoundPathName() first).
*/
TString GetPathName() { return fPathName; }
//----------------------------------------------------------------------
/**
* @brief Debug utility - Prints current search configuration to stdout.
*
* Outputs instrument name, year, run number, and all available templates
* with their patterns. Useful for troubleshooting path resolution issues.
*/
void DumpTemplateList();
private:
const PStringVector fPath;
const PRunNameTemplateList fRunNameTemplateList;
TString fInstrument{""};
Int_t fYear{-1};
Int_t fRun{-1};
TString fFileFormat{""};
TString fPathName{""};
const PStringVector fPath; ///< Search paths for data files
const PRunNameTemplateList fRunNameTemplateList; ///< Template patterns per instrument
TString fInstrument{""}; ///< Target instrument name (e.g., "GPS", "LEM")
Int_t fYear{-1}; ///< Run year (-1 if not specified)
Int_t fRun{-1}; ///< Run number (-1 if not specified)
TString fFileFormat{""}; ///< Optional file format filter (empty = any)
TString fPathName{""}; ///< Resolved full path (empty until found)
//----------------------------------------------------------------------
/**
* @brief Generates full file path by substituting template placeholders.
*
* Internal helper that replaces year and run number placeholders in a template
* with actual values. Supports variable-length run number formatting (2-9 digits).
*
* @param path Base directory path
* @param runNameTemplate Template string with placeholders (%yyyy%, %yy%, %rr...r%)
* @return Full path with placeholders substituted (e.g., "/data/gps/2023/02425.root")
*
* @par Template Examples:
* - "%yyyy%/%rrrrr%.root" with year=2023, run=42 → "2023/00042.root"
* - "run_%yy%_%rrr%.nxs" with year=2023, run=425 → "run_23_425.nxs"
*/
TString CreatePathName(const TString path, const TString runNameTemplate);
};

285
src/include/PFitter.h
View File

@@ -291,73 +291,294 @@ class PFitter
Bool_t DoFit();
private:
Bool_t fIsValid; ///< flag. true: the fit is valid.
Bool_t fIsScanOnly; ///< flag. true: scan along some parameters (no fitting).
Bool_t fConverged; ///< flag. true: the fit has converged.
Bool_t fChisqOnly; ///< flag. true: calculate chi^2 only (no fitting).
Bool_t fYamlOut; ///< flag. true: generate yaml output file of the fit results (MINUIT2.OUTPUT -> yaml)
Bool_t fUseChi2; ///< flag. true: chi^2 fit. false: log-max-likelihood
UInt_t fPrintLevel; ///< tag, showing the level of messages whished. 0=minimum, 1=standard, 2=maximum
// State flags
Bool_t fIsValid; ///< Overall validity flag: true if fitter initialized successfully
Bool_t fIsScanOnly; ///< Scan mode flag: true if only parameter scans requested (no minimization)
Bool_t fConverged; ///< Convergence flag: true if fit converged to a valid minimum
Bool_t fChisqOnly; ///< Evaluation-only flag: true to calculate χ² without fitting
Bool_t fYamlOut; ///< Output flag: true to generate YAML output file (MINUIT2.OUTPUT → yaml)
Bool_t fUseChi2; ///< Fit mode: true = χ² minimization, false = log-max-likelihood
UInt_t fPrintLevel; ///< Verbosity level: 0=quiet, 1=normal, 2=verbose (Minuit output)
UInt_t fStrategy; ///< fitting strategy (see minuit2 manual).
UInt_t fStrategy; ///< Minuit2 strategy: 0=fast/low-accuracy, 1=default, 2=careful/high-accuracy
PMsrHandler *fRunInfo; ///< pointer to the msr-file handler
PRunListCollection *fRunListCollection; ///< pointer to the run list collection
// Core data structures
PMsrHandler *fRunInfo; ///< Pointer to MSR file handler (parameters, theory, commands)
PRunListCollection *fRunListCollection; ///< Pointer to preprocessed run data collection
PMsrParamList fParams; ///< msr-file parameters
PMsrParamList fParams; ///< Copy of parameter list from MSR file
PMsrLines fCmdLines; ///< all the Minuit commands from the msr-file
PIntPairVector fCmdList; ///< command list, first=cmd, second=cmd line index
PMsrLines fCmdLines; ///< Raw command lines from MSR COMMANDS block
PIntPairVector fCmdList; ///< Parsed commands: first=command ID, second=line number
std::unique_ptr<PFitterFcn> fFitterFcn; ///< pointer to the fitter function object
std::unique_ptr<PFitterFcn> fFitterFcn; ///< Objective function for Minuit2 minimization
ROOT::Minuit2::MnUserParameters fMnUserParams; ///< minuit2 input parameter list
std::unique_ptr<ROOT::Minuit2::FunctionMinimum> fFcnMin; ///< function minimum object
ROOT::Minuit2::MnUserParameters fMnUserParams; ///< Minuit2 parameter state (values, errors, limits)
std::unique_ptr<ROOT::Minuit2::FunctionMinimum> fFcnMin; ///< Minuit2 function minimum result
// minuit2 scan/contours command relate variables (see MnScan/MnContours in the minuit2 user manual)
Bool_t fScanAll; ///< flag. false: single parameter scan, true: not implemented yet (see MnScan/MnContours in the minuit2 user manual)
UInt_t fScanParameter[2]; ///< scan parameter. idx=0: used for scan and contour, idx=1: used for contour (see MnScan/MnContours in the minuit2 user manual)
UInt_t fScanNoPoints; ///< number of points in a scan/contour (see MnScan/MnContours in the minuit2 user manual)
Double_t fScanLow; ///< scan range low. default=0.0 which means 2 std dev. (see MnScan/MnContours in the minuit2 user manual)
Double_t fScanHigh; ///< scan range high. default=0.0 which means 2 std dev. (see MnScan/MnContours in the minuit2 user manual)
PDoublePairVector fScanData; ///< keeps the scan/contour data
// Scan and contour analysis
Bool_t fScanAll; ///< Multi-parameter scan flag: false=1D scan, true=2D scan (not fully implemented)
UInt_t fScanParameter[2]; ///< Parameter indices: [0]=primary scan/contour, [1]=secondary (contours only)
UInt_t fScanNoPoints; ///< Number of scan/contour evaluation points (default=41)
Double_t fScanLow; ///< Scan lower bound: 0.0 = auto (2σ below current value)
Double_t fScanHigh; ///< Scan upper bound: 0.0 = auto (2σ above current value)
PDoublePairVector fScanData; ///< Scan results: (parameter_value, χ²) pairs
PDoublePairVector fOriginalFitRange; ///< keeps the original fit range in case there is a range command in the COMMAND block
PDoublePairVector fOriginalFitRange; ///< Original fit ranges per run (saved for FIT_RANGE command)
PStringVector fElapsedTime;
PStringVector fElapsedTime; ///< Timing information for each fit command
Bool_t fSectorFlag; ///< sector command present flag
std::vector<PSectorChisq> fSector; ///< stores all chisq/maxLH sector information
// Sector χ² analysis
Bool_t fSectorFlag; ///< SECTOR command present flag
std::vector<PSectorChisq> fSector; ///< Sector analysis results (χ² vs. time windows)
std::vector<bool> fPhase; ///< flag array in which an entry is true if the related parameter value is a phase
std::vector<bool> fPhase; ///< Phase parameter flags: true if parameter is a phase angle
// phase related functions
//----------------------------------------------------------------------
// Phase parameter identification (private helpers)
//----------------------------------------------------------------------
/**
* @brief Identifies which parameters represent phase angles.
*
* Scans the THEORY block to detect parameters used as phases in
* standard functions (TFieldCos, bessel, etc.). Phase parameters
* are constrained to [-360°, +360°] during fitting.
*/
void GetPhaseParams();
/**
* @brief Extracts parameter numbers from a FUNCTIONS block entry.
*
* Parses "funX" references in theory lines to find all parameters
* used in the function definition.
*
* @param funStr Function identifier string (e.g., "fun1", "fun23")
* @return Vector of parameter numbers (1-indexed) used in the function
*/
PIntVector GetParFromFun(const TString funStr);
/**
* @brief Extracts parameter numbers from a map reference.
*
* Parses "mapX" references to find mapped parameters across all runs.
* Maps allow different runs to use different parameters for the same
* theoretical component.
*
* @param mapStr Map identifier string (e.g., "map1", "map5")
* @return Vector of parameter numbers (1-indexed) referenced by the map
*/
PIntVector GetParFromMap(const TString mapStr);
// commands
Bool_t CheckCommands();
Bool_t SetParameters();
//----------------------------------------------------------------------
// Command validation and execution (private methods)
//----------------------------------------------------------------------
/**
* @brief Validates COMMANDS block syntax and builds execution queue.
*
* Parses all command lines, checks for syntax errors, extracts parameters,
* and populates fCmdList for sequential execution.
*
* @return true if all commands are valid, false on syntax errors
*/
Bool_t CheckCommands();
/**
* @brief Transfers MSR parameters to Minuit2 parameter state.
*
* Initializes fMnUserParams with values, errors, and bounds from the
* MSR file's PARAMETERS block.
*
* @return true if parameters set successfully
*/
Bool_t SetParameters();
/**
* @brief Executes CONTOURS command (2D error contours).
*
* Calculates confidence regions in 2D parameter space by evaluating
* χ² on a grid around the minimum.
*
* @return true if contour calculation succeeded
*/
Bool_t ExecuteContours();
/**
* @brief Executes FIT_RANGE command (optimal time-window search).
*
* Scans fit quality vs. fit start time to find the optimal first-good-bin.
* Useful for determining when background subtraction is adequate.
*
* @param lineNo Command line number in MSR file
* @return true if range scan succeeded
*/
Bool_t ExecuteFitRange(UInt_t lineNo);
/**
* @brief Executes FIX command (freeze parameters).
*
* Prevents specified parameters from varying during subsequent minimization.
*
* @param lineNo Command line number in MSR file
* @return true if parameters fixed successfully
*/
Bool_t ExecuteFix(UInt_t lineNo);
/**
* @brief Executes HESSE command (calculate error matrix).
*
* Computes the covariance matrix by evaluating second derivatives at
* the current minimum. Provides symmetric (parabolic) parameter errors.
*
* @return true if Hessian calculation succeeded
*/
Bool_t ExecuteHesse();
/**
* @brief Executes MIGRAD command (gradient descent minimization).
*
* Runs Minuit2's MIGRAD algorithm, the recommended robust minimizer
* using first derivatives and approximate Hessian updates.
*
* @return true if MIGRAD converged to a valid minimum
*/
Bool_t ExecuteMigrad();
/**
* @brief Executes MINIMIZE command (automatic algorithm selection).
*
* Lets Minuit2 choose the best minimization strategy. Usually equivalent
* to MIGRAD for well-behaved problems.
*
* @return true if minimization converged
*/
Bool_t ExecuteMinimize();
/**
* @brief Executes MINOS command (asymmetric error analysis).
*
* Computes accurate asymmetric confidence intervals by scanning χ²
* along each parameter axis. Slower but more accurate than HESSE.
*
* @return true if MINOS analysis completed
*/
Bool_t ExecuteMinos();
/**
* @brief Executes PLOT command (visualize scan/contour results).
*
* Displays scan or contour data from previous SCAN/CONTOURS commands.
*
* @return true if plot generated successfully
*/
Bool_t ExecutePlot();
/**
* @brief Executes PRINT command (set verbosity level).
*
* Controls Minuit2 output detail: 0=minimal, 1=normal, 2=debug.
*
* @param lineNo Command line number in MSR file
* @return true if print level set successfully
*/
Bool_t ExecutePrintLevel(UInt_t lineNo);
/**
* @brief Executes RELEASE command (unfreeze parameters).
*
* Allows previously fixed parameters to vary in subsequent fits.
*
* @param lineNo Command line number in MSR file
* @return true if parameters released successfully
*/
Bool_t ExecuteRelease(UInt_t lineNo);
/**
* @brief Executes RESTORE command (reload saved parameters).
*
* Restores parameter values from the last SAVE command.
*
* @return true if parameters restored successfully
*/
Bool_t ExecuteRestore();
/**
* @brief Executes SCAN command (1D parameter space scan).
*
* Evaluates χ² along one or two parameter axes to visualize the
* objective function landscape near the minimum.
*
* @return true if scan completed
*/
Bool_t ExecuteScan();
/**
* @brief Executes SAVE command (store current parameters).
*
* Saves current parameter state for later RESTORE. Updates MSR file
* statistics on first save (after final fit).
*
* @param first True if this is the first SAVE command in the session
* @return true if parameters saved successfully
*/
Bool_t ExecuteSave(Bool_t first);
/**
* @brief Executes SIMPLEX command (non-gradient minimization).
*
* Runs the Nelder-Mead simplex algorithm. Robust for rough objective
* functions but slow to converge. Often used before MIGRAD for difficult fits.
*
* @return true if SIMPLEX found a minimum
*/
Bool_t ExecuteSimplex();
/**
* @brief Prepares sector χ² analysis data structures.
*
* Initializes sector time windows and allocates storage for sector results.
*
* @param param Current parameter values
* @param error Current parameter errors
*/
void PrepareSector(PDoubleVector &param, PDoubleVector &error);
/**
* @brief Executes SECTOR command (time-dependent χ² analysis).
*
* Calculates χ² for progressively wider time windows to identify
* optimal fit ranges and systematic time-dependent effects.
*
* @param fout Output stream for sector analysis results
* @return true if sector analysis completed
*/
Bool_t ExecuteSector(std::ofstream &fout);
//----------------------------------------------------------------------
// Utility functions (private)
//----------------------------------------------------------------------
/**
* @brief Returns current time in milliseconds.
*
* Used for timing fit commands and generating performance statistics.
*
* @return Timestamp in milliseconds since epoch
*/
Double_t MilliTime();
/**
* @brief Rounds parameters for output with appropriate precision.
*
* Determines significant figures based on errors and formats parameters
* for display in MSR file output.
*
* @param par Parameter values
* @param err Parameter errors
* @param ok Output flag: false if rounding failed
* @return Rounded parameter values
*/
PDoubleVector ParamRound(const PDoubleVector &par, const PDoubleVector &err, Bool_t &ok);
};

118
src/include/PFitterFcn.h
View File

@@ -36,26 +36,136 @@
#include "PRunListCollection.h"
//--------------------------------------------------------------------------
/**
* <p>This is the minuit2 interface function class porviding the function to be optimized (chisq or log max-likelihood).
* @brief Objective function interface for ROOT Minuit2 minimization.
*
* This class implements the FCNBase interface required by ROOT's Minuit2
* minimizer. It provides the objective function (χ² or log-likelihood)
* that Minuit2 minimizes during parameter optimization.
*
* The class serves as a bridge between musrfit's data structures
* (PRunListCollection) and Minuit2's optimization algorithms, calculating
* the goodness-of-fit measure for any given parameter set.
*
* @par Fitting modes:
* - **χ² minimization:** Standard least-squares fitting for Gaussian errors
* - **Maximum likelihood:** Poisson statistics, better for low-count data
*
* @par Usage in fitting workflow:
* 1. PFitter creates a PFitterFcn instance with data and fit mode
* 2. Minuit2 calls operator()() repeatedly with trial parameter sets
* 3. operator()() calculates χ²/maxLH by evaluating theory vs. data
* 4. Minuit2 searches parameter space to minimize the returned value
* 5. Up() defines the error criterion (Δχ²=1 or ΔmaxLH=0.5 for 1σ)
*
* @see PFitter, PRunListCollection
* @see ROOT::Minuit2::FCNBase in ROOT Minuit2 documentation
*/
class PFitterFcn : public ROOT::Minuit2::FCNBase
{
public:
//----------------------------------------------------------------------
/**
* @brief Constructor for objective function.
*
* Initializes the function evaluator with preprocessed data and
* configures the error definition based on the fitting mode.
*
* @param runList Pointer to collection of preprocessed run data
* @param useChi2 If true, use χ² minimization; if false, use maximum likelihood
*
* @note The runList pointer must remain valid for the lifetime of this object.
*/
PFitterFcn(PRunListCollection *runList, Bool_t useChi2);
/**
* @brief Destructor.
*/
~PFitterFcn();
//----------------------------------------------------------------------
/**
* @brief Returns error definition for Minuit2 (Up value).
*
* The "Up" value defines what change in the objective function
* corresponds to 1σ error bars on parameters:
* - For χ² fits: Up = 1.0 (parabolic errors, Δχ²=1)
* - For max likelihood: Up = 0.5 (asymmetric errors, ΔmaxLH=0.5)
*
* This value is used by Minuit2's error analysis algorithms (HESSE, MINOS).
*
* @return Error definition value (1.0 for χ², 0.5 for likelihood)
*
* @see ROOT::Minuit2::FCNBase::Up() in Minuit2 manual
*/
Double_t Up() const { return fUp; }
//----------------------------------------------------------------------
/**
* @brief Evaluates objective function for given parameters.
*
* This is the core function called by Minuit2 during minimization.
* It computes either χ² or negative log-likelihood by:
* 1. Passing parameters to PRunListCollection
* 2. Calculating theory predictions for all runs
* 3. Comparing theory vs. data across all fitted bins
* 4. Returning the total χ²/maxLH value
*
* @param par Parameter vector with current trial values
* @return χ² value (if fUseChi2=true) or -2×log-likelihood (if fUseChi2=false)
*
* @note This function must be const as required by FCNBase interface.
* @note For likelihood fits, returns -2×ln(L) so minimization is equivalent to maximizing L.
*
* @par Performance:
* This function is called hundreds to thousands of times during
* a fit, so it's optimized for speed (parallel evaluation if OpenMP enabled).
*/
Double_t operator()(const std::vector<Double_t> &par) const;
//----------------------------------------------------------------------
/**
* @brief Returns total number of bins used in the fit across all runs.
*
* @return Total count of fitted bins (summed over all runs)
*/
UInt_t GetTotalNoOfFittedBins() { return fRunListCollection->GetTotalNoOfBinsFitted(); }
//----------------------------------------------------------------------
/**
* @brief Returns number of fitted bins for a specific run.
*
* @param idx Run index (0-based)
* @return Number of bins fitted for the specified run
*/
UInt_t GetNoOfFittedBins(const UInt_t idx) { return fRunListCollection->GetNoOfBinsFitted(idx); }
//----------------------------------------------------------------------
/**
* @brief Calculates expected χ² (or maxLH) for quality assessment.
*
* Computes the theoretical expected value of χ² assuming the model
* is correct. This is used to assess goodness-of-fit:
* - If observed χ² ≈ expected χ²: fit is consistent with data quality
* - If observed χ² >> expected χ²: systematic deviations present
* - If observed χ² << expected χ²: possible overestimated errors
*
* For single histogram fits, expected χ² = NDF. For asymmetry fits,
* the calculation is more complex due to error propagation.
*
* @param par Parameter vector for evaluation
* @param totalExpectedChisq Returns total expected χ²/maxLH (output)
* @param expectedChisqPerRun Returns expected χ²/maxLH for each run (output)
*
* @note The expectedChisqPerRun vector is resized and filled by this method.
*/
void CalcExpectedChiSquare(const std::vector<Double_t> &par, Double_t &totalExpectedChisq, std::vector<Double_t> &expectedChisqPerRun);
private:
Double_t fUp; ///< for chisq == 1.0, i.e. errors are 1 std. deviation errors. for log max-likelihood == 0.5, i.e. errors are 1 std. deviation errors (for details see the minuit2 user manual).
Bool_t fUseChi2; ///< true = chisq fit, false = log max-likelihood fit
PRunListCollection *fRunListCollection; ///< pre-processed data to be fitted
Double_t fUp; ///< Error definition: 1.0 for χ² (1σ = Δχ²=1), 0.5 for maxLH (1σ = ΔmaxLH=0.5)
Bool_t fUseChi2; ///< Fit mode flag: true = χ² minimization, false = max log-likelihood
PRunListCollection *fRunListCollection; ///< Pointer to preprocessed muSR data collection
};
#endif // _PFITTERFCN_H_