diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index e1944d779..0b23549c4 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -260,11 +260,32 @@ UInt_t PSectorChisq::GetNDF(UInt_t idx) // Constructor //-------------------------------------------------------------------------- /** - *

Constructor. + * @brief Constructor for the fitting engine. * - * \param runInfo pointer of the msr-file handler - * \param runListCollection pointer of the run list collection (pre-processed historgrams) - * \param chisq_only flag: true=calculate chisq only (no fitting) + * Initializes the fitter with MSR configuration and preprocessed data. + * Sets up the fitting environment including: + * - Parameter lists and command queues + * - Fit mode (χ² vs. likelihood) + * - Strategy and print level defaults + * - Original fit ranges for RANGE commands + * - Phase parameter identification + * + * The constructor validates the COMMANDS block and creates the objective + * function (PFitterFcn) but does not start the fit. Call DoFit() to execute. + * + * @param runInfo Pointer to MSR file handler containing fit configuration + * @param runListCollection Pointer to preprocessed data collection + * @param chisq_only If true, only evaluate χ² without fitting (useful for validation) + * @param yaml_out If true, generate YAML output file with fit results + * + * @note Sets fIsValid = false if command validation fails. Check IsValid() before DoFit(). + * + * @par Initialization sequence: + * 1. Copy parameters and commands from MSR handler + * 2. Store original fit ranges (for FIT_RANGE command) + * 3. Parse and validate COMMANDS block + * 4. Identify phase parameters (for angle wrapping) + * 5. Create objective function interface */ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bool_t chisq_only, Bool_t yaml_out) : fChisqOnly(chisq_only), fYamlOut(yaml_out), fRunInfo(runInfo), fRunListCollection(runListCollection) @@ -323,7 +344,10 @@ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bo // Destructor //-------------------------------------------------------------------------- /** - *

Destructor. + * @brief Destructor - Cleans up dynamically allocated resources. + * + * Frees memory used by command lists, scan data, and timing information. + * Smart pointers (fFitterFcn, fFcnMin) are automatically cleaned up. */ PFitter::~PFitter() { @@ -338,8 +362,30 @@ PFitter::~PFitter() // GetPhaseParams (private) //-------------------------------------------------------------------------- /** - *

Checks which parameters are phases. This information is needed to - * restrict the phases to the intervall -360 to +360 degrees. + * @brief Identifies which parameters represent phase angles. + * + * Scans the THEORY block to detect parameters used as phases in standard + * muSR functions. Phase parameters are flagged to enable automatic wrapping + * to the interval [-360°, +360°] during fitting, preventing meaningless + * phase values outside this range. + * + * Recognized phase-containing functions: + * - TFieldCos/tf: phase is 1st parameter + * - bessel/b: phase is 1st parameter + * - skewedGss/skg: phase is 1st parameter + * - staticNKTF/snktf: phase is 1st parameter + * - dynamicNKTF/dnktf: phase is 1st parameter + * - internFld/if: phase is 2nd parameter + * - internBsl/ib: phase is 2nd parameter + * - muMinusExpTF/mmsetf: phase is 5th parameter + * + * Phase references can be: + * - Direct parameter numbers (e.g., "7" → par7) + * - Function references (e.g., "fun3" → all params in function 3) + * - Map references (e.g., "map2" → parameters mapped across runs) + * + * @note Populates fPhase vector where fPhase[i] = true means parameter i+1 is a phase. + * @note User-defined functions cannot be automatically analyzed for phases. */ void PFitter::GetPhaseParams() { @@ -423,11 +469,22 @@ void PFitter::GetPhaseParams() // GetParFromFun (private) //-------------------------------------------------------------------------- /** - *

Extract from string funX the function number. Base on the function number - * the paramter numbers will be collected. + * @brief Extracts parameter numbers from a FUNCTIONS block entry. * - * @param funStr string of the form funX, where X is the function number - * @return a vector of all the parameter numbers related to funX + * Parses a function definition to find all parameters it references. + * Recursively handles nested map references within the function. + * + * @param funStr Function identifier (e.g., "fun1", "fun23") + * @return Vector of parameter numbers (1-indexed) used in the function + * + * @par Example: + * If FUNCTIONS block contains: + * @code + * fun1 = par3 * par7 + map2 + * @endcode + * Then GetParFromFun("fun1") returns [3, 7, ...params from map2...] + * + * @note Returns empty vector if function not found or parsing fails. */ PIntVector PFitter::GetParFromFun(const TString funStr) { @@ -490,11 +547,25 @@ PIntVector PFitter::GetParFromFun(const TString funStr) // GetParFromMap (private) //-------------------------------------------------------------------------- /** - *

Extract from string mapX the map number. Based on the map number the - * parameter numbers will be collected. + * @brief Extracts parameter numbers from a map reference. * - * @param mapStr string of the form mapX, where X is the map number - * @return a vector of all the parameter numbers related to mapX + * Parses "mapX" to find which parameters are mapped to the X-th position + * across all RUN blocks. Maps allow different runs to use different parameters + * for the same theoretical component, enabling multi-run fits with run-dependent + * parameters. + * + * @param mapStr Map identifier (e.g., "map1", "map5") + * @return Vector of parameter numbers (1-indexed) mapped at this position + * + * @par Example: + * If RUN blocks have: + * @code + * RUN 1: map 7 8 9 + * RUN 2: map 10 11 12 + * @endcode + * Then GetParFromMap("map1") returns [7, 10] (first map entry per run) + * + * @note Returns empty vector if map index is out of range or invalid. */ PIntVector PFitter::GetParFromMap(const TString mapStr) { @@ -541,9 +612,37 @@ PIntVector PFitter::GetParFromMap(const TString mapStr) // DoFit //-------------------------------------------------------------------------- /** - *

Main calling routine to invoke minuit2, i.e. fitting etc. + * @brief Main entry point for executing the fit. * - * return: true if all commands could be executed successfully, otherwise returns false. + * This is the primary method that orchestrates the entire fitting process: + * 1. Transfers parameters from MSR to Minuit2 + * 2. If chisq_only mode: calculates χ²/maxLH and returns + * 3. Otherwise: executes COMMANDS block sequentially + * 4. Updates MSR file with fit results and statistics + * + * @return true if fit completed successfully (or chisq calculated), false on errors + * + * @par Execution modes: + * - **chisq_only = true:** Evaluates objective function once with current parameters, + * reports χ²/maxLH and NDF, useful for validating parameter sets + * - **chisq_only = false:** Runs full fit according to COMMANDS block, which typically + * includes MIGRAD (minimization), HESSE (error matrix), optionally MINOS (asymmetric errors) + * + * @par COMMANDS execution: + * Commands are executed in the order they appear in the MSR file. Common sequence: + * @code + * SET PRINT 1 # Set verbosity + * MIGRAD # Find minimum + * HESSE # Calculate symmetric errors + * MINOS # Calculate asymmetric errors (optional) + * SAVE # Save results to MSR file + * @endcode + * + * @note Check IsValid() before calling. Check HasConverged() after completion. + * @note Updates fRunInfo with final parameters, errors, and fit statistics. + * @note Prints detailed χ² breakdown per run if multiple runs are fitted. + * + * @see ExecuteMigrad(), ExecuteHesse(), ExecuteMinos(), ExecuteSave() */ Bool_t PFitter::DoFit() { @@ -825,10 +924,31 @@ Bool_t PFitter::DoFit() // CheckCommands //-------------------------------------------------------------------------- /** - *

Check the msr-file COMMAND's, fill the command queue and make sure that given parameters (if present) - * do make any sense. + * @brief Validates COMMANDS block and builds execution queue. * - * return: true if the commands are valid, otherwise returns false. + * Parses all command lines from the MSR file's COMMANDS block, validates + * syntax and parameters, and constructs an ordered execution list. This + * ensures commands are executable before starting the fit. + * + * @return true if all commands are valid and parseable, false on any error + * + * @par Validation checks: + * - Command keywords are recognized (MIGRAD, HESSE, MINOS, etc.) + * - Numeric arguments are valid (parameter indices, ranges, point counts) + * - Parameter references are within bounds + * - SCAN/CONTOURS have required arguments + * - FIX/RELEASE specify valid parameter lists + * + * @par Side effects: + * - Populates fCmdList with (command_ID, line_number) pairs + * - Sets fIsValid = false if any command is invalid + * - Configures scan parameters (fScanParameter, fScanNoPoints, etc.) + * - Detects SECTOR command presence (sets fSectorFlag) + * + * @note Called automatically by constructor. Errors are reported to stderr. + * @note Some legacy commands (SET BATCH, END RETURN) are silently ignored. + * + * @see PMN_* command constants in PFitter.h */ Bool_t PFitter::CheckCommands() { diff --git a/src/include/PFindRun.h b/src/include/PFindRun.h index e70a4e0ba..7b236832f 100644 --- a/src/include/PFindRun.h +++ b/src/include/PFindRun.h @@ -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); }; diff --git a/src/include/PFitter.h b/src/include/PFitter.h index a7bd95d84..7dc018259 100644 --- a/src/include/PFitter.h +++ b/src/include/PFitter.h @@ -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 fFitterFcn; ///< pointer to the fitter function object + std::unique_ptr fFitterFcn; ///< Objective function for Minuit2 minimization - ROOT::Minuit2::MnUserParameters fMnUserParams; ///< minuit2 input parameter list - std::unique_ptr fFcnMin; ///< function minimum object + ROOT::Minuit2::MnUserParameters fMnUserParams; ///< Minuit2 parameter state (values, errors, limits) + std::unique_ptr 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 fSector; ///< stores all chisq/maxLH sector information + // Sector χ² analysis + Bool_t fSectorFlag; ///< SECTOR command present flag + std::vector fSector; ///< Sector analysis results (χ² vs. time windows) - std::vector fPhase; ///< flag array in which an entry is true if the related parameter value is a phase + std::vector 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 ¶m, 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); }; diff --git a/src/include/PFitterFcn.h b/src/include/PFitterFcn.h index dac342000..edea7a647 100644 --- a/src/include/PFitterFcn.h +++ b/src/include/PFitterFcn.h @@ -36,26 +36,136 @@ #include "PRunListCollection.h" +//-------------------------------------------------------------------------- /** - *

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 &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 &par, Double_t &totalExpectedChisq, std::vector &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_