/***************************************************************************

  PRunSingleHisto.h

  Author: Andreas Suter
  e-mail: andreas.suter@psi.ch

***************************************************************************/

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

#ifndef _PRUNSINGLEHISTO_H_
#define _PRUNSINGLEHISTO_H_

#include "PRunBase.h"

/**
 * \brief Class for fitting single detector histograms (basic time-differential μSR).
 *
 * PRunSingleHisto implements the most fundamental μSR analysis: fitting a single
 * positron detector histogram to extract relaxation parameters. This is the basis
 * for all μSR measurements and is used when asymmetry analysis is not appropriate
 * or desired.
 *
 * \section singlehisto_physics Physics and Applications
 *
 * <b>Single histogram measurements are used for:</b>
 * - <b>Time-differential μSR:</b> Measure μ⁺ decay positron time spectrum
 * - <b>Detector calibration:</b> Characterize individual detector response
 * - <b>Transverse field (TF):</b> When forward/backward separation unnecessary
 * - <b>High statistics:</b> Use all positrons (no F-B discrimination)
 * - <b>Specialized geometries:</b> Non-standard detector arrangements
 * - <b>Method development:</b> Test fitting algorithms on single detector
 *
 * \section singlehisto_data Data Structure and Analysis
 *
 * <b>Histogram content:</b>
 * - Raw counts vs. time from a single positron detector
 * - Time-zero (t0): Muon arrival time marking start of decay
 * - Background: Constant or estimated from pre-t0 bins
 * - Signal: N(t) = N₀·exp(-t/τ_μ)·P(t) + B
 *
 * where:
 * - N₀ = initial count rate (normalization parameter)
 * - τ_μ = 2.197 μs (muon lifetime)
 * - P(t) = polarization function (contains physics: relaxation, oscillation, etc.)
 * - B = background (random coincidences, accidentals)
 *
 * \section singlehisto_workflow Analysis Workflow
 *
 * 1. <b>Load Histogram:</b> Read raw detector counts from data file
 * 2. <b>Determine t0:</b> Identify muon arrival time
 * 3. <b>Background Subtraction:</b>
 *    - Fixed: Subtract specified constant
 *    - Estimated: Calculate from pre-t0 bins, subtract with error propagation
 * 4. <b>Bin Packing:</b> Rebin to improve statistics (optional)
 * 5. <b>Fit Range:</b> Select time window for parameter extraction
 * 6. <b>Theory Evaluation:</b> N_theory(t) = N₀·exp(-t/τ_μ)·P_theory(t) + B
 * 7. <b>Minimization:</b> χ² or maximum likelihood via MINUIT
 *
 * \section singlehisto_theory Theory Function
 *
 * The fitted function is typically:
 * \f[ N(t) = N_0 e^{-t/\tau_\mu} P(t) + B \f]
 *
 * Common polarization functions P(t):
 * - <b>Static Gaussian:</b> P(t) = exp(-σ²t²/2)
 * - <b>Static Lorentzian:</b> P(t) = exp(-λt)
 * - <b>Dynamic relaxation:</b> P(t) = exp(-(λt)^β)  (stretched exponential)
 * - <b>Oscillating:</b> P(t) = cos(ωt + φ) · exp(-λt)
 * - <b>Kubo-Toyabe:</b> Complex relaxation functions for spin systems
 *
 * \section singlehisto_parameters Key Parameters
 *
 * <b>Normalization (N₀):</b>
 * - Can be fit parameter or fixed value
 * - Can be derived from FUNCTIONS block
 * - Automatically scaled to 1/ns or 1/bin (fScaleN0AndBkg)
 *
 * <b>Background (B):</b>
 * - Fixed: From "background" entry in RUN block
 * - Estimated: Calculated from pre-t0 bin range
 * - Units: counts/bin
 *
 * <b>Packing:</b>
 * - Number of consecutive bins to combine
 * - REQUIRED parameter (RUN or GLOBAL block)
 * - Higher packing → better statistics, worse time resolution
 *
 * \section singlehisto_msr MSR File Example
 *
 * \code
 * RUN data/run2425 PSI MUE4 PSI MUSR-ROOT   (name beamline)
 *   fittype         0       (SingleHisto)
 *   map             1       (forward histogram number)
 *   forward         1
 *   packing         50      (combine 50 bins → one packed bin)
 *   background      50 150  (estimate from bins 50-150)
 *   data            200 2000 (use bins 200-2000 for analysis)
 *   t0              210.5   (muon arrival bin)
 *   fit             0.1 10.0 (fit from 0.1 to 10.0 μs after t0)
 *
 * THEORY
 *   asymmetry      1
 *   simpleGss      2       (Gaussian relaxation σ)
 *   +              3       (constant background offset)
 * \endcode
 *
 * \section singlehisto_vs_asymmetry Single Histo vs. Asymmetry
 *
 * <table>
 * <tr><th>Feature</th><th>Single Histogram</th><th>Asymmetry</th></tr>
 * <tr><td>Detectors</td><td>One (forward)</td><td>Two (forward + backward)</td></tr>
 * <tr><td>Quantity fitted</td><td>N(t) counts</td><td>A(t) = (F-αB)/(F+αB)</td></tr>
 * <tr><td>Statistics</td><td>All positrons</td><td>F and B separately</td></tr>
 * <tr><td>Background</td><td>Additive B</td><td>Cancels in asymmetry</td></tr>
 * <tr><td>α parameter</td><td>N/A</td><td>Required (F/B asymmetry)</td></tr>
 * <tr><td>Use cases</td><td>TF, calibration</td><td>LF, ZF, weak TF</td></tr>
 * </table>
 *
 * \see PRunAsymmetry for forward-backward asymmetry analysis
 * \see PRunSingleHistoRRF for rotating reference frame (high-TF)
 * \see PRunBase for base class interface and common functionality
 */
class PRunSingleHisto : public PRunBase
{
  public:
    /**
     * \brief Default constructor creating an empty, invalid single histogram run object.
     *
     * Initializes all member variables to default/invalid states:
     * - fScaleN0AndBkg = true (scale to 1/ns by default)
     * - fNoOfFitBins = 0 (no bins to fit)
     * - fBackground = 0 (no background)
     * - fPacking = -1 (unspecified - will cause error if not set)
     * - fTheoAsData = false (high-resolution theory grid)
     * - Good bins markers = -1 (unset)
     * - Fit range bins = -1 (unset)
     *
     * This constructor is needed for creating vectors of PRunSingleHisto objects.
     * The resulting object cannot be used until properly initialized via the main constructor.
     */
    PRunSingleHisto();

    /**
     * \brief Main constructor initializing a single histogram run from MSR file and data.
     *
     * Performs comprehensive initialization:
     *
     * 1. <b>Base Class Initialization:</b>
     *    - Calls PRunBase constructor
     *    - Initializes theory engine, parameter mappings, FUNCTIONS block
     *
     * 2. <b>N₀/Background Scaling:</b>
     *    - Calls IsScaleN0AndBkg() to determine scaling mode
     *    - Sets fScaleN0AndBkg flag (true = scale to 1/ns, false = leave as 1/bin)
     *
     * 3. <b>Packing Validation (CRITICAL):</b>
     *    - Attempts to read packing from RUN block
     *    - Falls back to GLOBAL block if not in RUN
     *    - SEVERE ERROR if packing == -1 (mandatory parameter)
     *    - Marks run invalid and returns if packing not found
     *
     * 4. <b>Member Initialization:</b>
     *    - Good bin markers, fit range bins set to -1 (determined later)
     *    - Background initialized to 0 (set during PrepareData if specified)
     *
     * 5. <b>Data Preparation:</b>
     *    - Calls PrepareData() to load and preprocess histogram
     *    - If PrepareData() fails → marks run invalid
     *
     * The object is marked as invalid (fValid=false) if:
     * - Packing parameter is missing from both RUN and GLOBAL blocks
     * - PrepareData() fails (file not found, invalid t0, etc.)
     *
     * \param msrInfo Pointer to MSR file handler (must remain valid)
     * \param rawData Pointer to raw data handler for histogram loading
     * \param runNo Run number (0-based index in MSR file RUN blocks)
     * \param tag Operation mode: kFit (fitting), kView (display/plotting)
     * \param theoAsData Theory mode: true = at data points, false = high-resolution
     *
     * \warning Always check IsValid() after construction. Packing is MANDATORY.
     *
     * \see PrepareData() for data preprocessing details
     * \see IsScaleN0AndBkg() for scaling determination
     */
    PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData);

    /**
     * \brief Virtual destructor cleaning up allocated resources.
     *
     * Releases memory used by the forward histogram vector (fForward).
     * Base class destructor handles cleanup of theory objects and other shared resources.
     */
    virtual ~PRunSingleHisto();

    /**
     * \brief Calculates χ² between histogram data and theory.
     *
     * Computes chi-squared for single histogram fitting:
     * \f[ \chi^2 = \sum_{i} \frac{(N_i^{\rm data} - N_i^{\rm theory})^2}{\sigma_i^2} \f]
     *
     * where N_theory(t) = N₀·exp(-t/τ_μ)·P(t) + B
     *
     * Uses OpenMP parallelization when available. N₀ can be a fit parameter
     * or derived from FUNCTIONS block.
     *
     * \param par Parameter vector from MINUIT
     * \return χ² value (minimize during fitting)
     */
    virtual Double_t CalcChiSquare(const std::vector<Double_t>& par);

    /**
     * \brief Calculates expected χ² based on theory predictions.
     * \param par Parameter vector from MINUIT
     * \return Expected χ² (for statistical diagnostics)
     */
    virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par);

    /**
     * \brief Calculates maximum likelihood for Poisson-distributed histogram counts.
     *
     * Computes -2ln(L) for low-count data (< 10-20 counts/bin).
     * Superior to χ² when Gaussian approximation invalid.
     *
     * \param par Parameter vector from MINUIT
     * \return -2×ln(L) value
     */
    virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);

    /**
     * \brief Calculates expected maximum likelihood.
     * \param par Parameter vector from MINUIT
     * \return Expected -2×ln(L)
     */
    virtual Double_t CalcMaxLikelihoodExpected(const std::vector<Double_t>& par);

    /**
     * \brief Evaluates theory function at all data points or high-resolution grid.
     *
     * Calculates N_theory(t) = N₀·exp(-t/τ_μ)·P_theory(t) + B
     * using THEORY block functions. Stores results in fData for χ² calculation.
     */
    virtual void CalcTheory();

    /**
     * \brief Returns the number of bins included in the fit range.
     *
     * Used for degrees of freedom: ν = N_bins - N_params
     *
     * \return Number of bins in fit range
     */
    virtual UInt_t GetNoOfFitBins();

    /**
     * \brief Sets fit range using bin-offset specification (COMMANDS block syntax).
     *
     * Format: "fit_range fgb+n0 lgb-n1"
     *
     * \param fitRange String with bin offsets from good bin markers
     */
    virtual void SetFitRangeBin(const TString fitRange);

    /**
     * \brief Returns the estimated background level.
     * \return Background in counts/bin
     */
    virtual Double_t GetBackground() { return fBackground; }

    /**
     * \brief Returns the first bin index in the fit range.
     * \return Start bin index (0-based, after packing)
     */
    virtual Int_t GetStartTimeBin() { return fStartTimeBin; }

    /**
     * \brief Returns the last bin index in the fit range (exclusive).
     * \return End bin index (loop condition: i < fEndTimeBin)
     */
    virtual Int_t GetEndTimeBin() { return fEndTimeBin; }

    /**
     * \brief Returns the bin packing factor.
     * \return Number of raw bins combined into one packed bin
     */
    virtual Int_t GetPacking() { return fPacking; }

    /**
     * \brief Returns the N₀/background scaling mode.
     * \return true = scaled to 1/ns, false = left as 1/bin
     */
    virtual Bool_t GetScaleN0AndBkg() { return fScaleN0AndBkg; }

    /**
     * \brief Calculates start/end bin indices from fit time range.
     *
     * Converts fit range (μs) to bin indices, accounting for t0, time resolution,
     * and packing. Updates fStartTimeBin, fEndTimeBin, fNoOfFitBins.
     */
    virtual void CalcNoOfFitBins();

  protected:
    /**
     * \brief Main data preparation orchestrator.
     *
     * Coordinates histogram loading and preprocessing: determines operation mode,
     * calls PrepareFitData() or PrepareViewData(), validates success.
     *
     * \return True if data preparation succeeds, false on error
     */
    virtual Bool_t PrepareData();

    /**
     * \brief Prepares histogram data for fitting.
     *
     * Loads forward histogram, extracts metadata, determines t0, subtracts background,
     * packs bins, propagates errors, sets up time grid and fit ranges.
     *
     * \param runData Raw run data handler
     * \param histoNo Histogram index in data file
     * \return True on success, false if preprocessing fails
     */
    virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo);

    /**
     * \brief Prepares raw histogram data for viewing (minimal processing).
     *
     * Lighter-weight preprocessing for raw histogram visualization without
     * background subtraction or full fitting infrastructure.
     *
     * \param runData Raw run data handler
     * \param histoNo Histogram index
     * \return True on success
     */
    virtual Bool_t PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo);

    /**
     * \brief Prepares processed histogram data for viewing/plotting.
     *
     * Similar to PrepareFitData() but optimized for visualization with potentially
     * wider time range for context.
     *
     * \param runData Raw run data handler
     * \param histoNo Histogram index
     * \return True on success
     */
    virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo);

  private:
    Bool_t fScaleN0AndBkg;  ///< Scaling mode: true = scale N₀ and B to 1/ns, false = leave as 1/bin (determined by IsScaleN0AndBkg())
    UInt_t fNoOfFitBins;    ///< Number of bins within fit range (fStartTimeBin to fEndTimeBin)
    Double_t fBackground;   ///< Background level in counts/bin (estimated from pre-t0 bins or fixed value from RUN block)
    Int_t fPacking;         ///< Bin packing factor (REQUIRED: from RUN or GLOBAL block)
    Bool_t fTheoAsData;     ///< Theory mode: true = at data points, false = high-resolution grid for smooth Fourier transforms

    Int_t fGoodBins[2];     ///< Good bin markers for COMMANDS block: [0]=fgb (first good bin/t0), [1]=lgb (last good bin)

    PDoubleVector fForward; ///< Forward detector histogram (background-corrected, packed)

    Int_t fStartTimeBin;    ///< First bin index in fit range (inclusive, 0-based after packing)
    Int_t fEndTimeBin;      ///< Last bin index in fit range (exclusive: loop as i < fEndTimeBin)

    /**
     * \brief Determines and validates t0 values for histogram.
     *
     * Extracts time-zero from RUN block, data file header, GLOBAL block, or
     * automatic determination. Validates t0 is within histogram bounds.
     *
     * \param runData Raw run data
     * \param globalBlock GLOBAL block settings
     * \param histoNo Vector of histogram indices
     * \return True if valid t0 found
     */
    virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo);

    /**
     * \brief Determines data range (region of valid histogram data).
     *
     * Establishes start/end bins for analysis from RUN block "data" entry.
     * Data range is typically wider than fit range.
     *
     * \return True if valid data range determined
     */
    virtual Bool_t GetProperDataRange();

    /**
     * \brief Determines fit range from MSR file settings.
     *
     * Extracts fit time window from RUN or GLOBAL block "fit" entry.
     * Format: time-based (μs) or bin-based (fgb+n0 lgb-n1).
     *
     * \param globalBlock GLOBAL block with default fit settings
     */
    virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock);

    /**
     * \brief Estimates initial normalization N₀ from histogram data.
     *
     * Calculates N₀ estimate from histogram amplitude, used as starting
     * value if N₀ is a fit parameter.
     */
    virtual void EstimateN0();

    /**
     * \brief Estimates background from pre-t0 bins.
     *
     * Calculates background average and error from specified bin range
     * (typically before t0). Sets fBackground member.
     *
     * \param histoNo Histogram index
     * \return True on success, false if background range invalid
     */
    virtual Bool_t EstimateBkg(UInt_t histoNo);

    /**
     * \brief Determines if N₀ and background should be scaled to 1/ns.
     *
     * Checks time resolution and fitting preferences to decide scaling mode.
     * Returns true for standard time bins (scale to 1/ns), false otherwise.
     *
     * \return True if scaling should be applied
     */
    virtual Bool_t IsScaleN0AndBkg();
};

#endif // _PRUNSINGLEHISTO_H_