improve the doxygen docu.

src/classes/PFourier.cpp

@@ -51,7 +51,18 @@
 // Constructor
 //--------------------------------------------------------------------------
 /**
+ * <p>Default constructor for phase correction optimizer.
  *
+ * <p>Initializes a phase correction object with default Fourier data.
+ * The actual Fourier data must be set externally before minimization.
+ *
+ * @param minBin Minimum frequency bin index for optimization region.
+ *               Use -1 to optimize over all bins. Default: -1.
+ * @param maxBin Maximum frequency bin index for optimization region.
+ *               Use -1 to optimize over all bins. Default: -1.
+ *
+ * <p><b>Note:</b> Restricting the optimization range to signal-containing
+ * regions improves results by excluding noisy baseline regions.
  */
 PFTPhaseCorrection::PFTPhaseCorrection(const Int_t minBin, const Int_t maxBin) :
   fMinBin(minBin), fMaxBin(maxBin)
@@ -63,7 +74,28 @@ PFTPhaseCorrection::PFTPhaseCorrection(const Int_t minBin, const Int_t maxBin) :
 // Constructor
 //--------------------------------------------------------------------------
 /**
+ * <p>Constructor with explicit Fourier transform data.
  *
+ * <p>Creates a phase correction optimizer with pre-computed Fourier data.
+ * The optimizer will find phase parameters φ(ω) = c₀ + c₁·ω that maximize
+ * the real component while minimizing the imaginary component.
+ *
+ * @param reFT Real part of the Fourier transform (vector of amplitudes).
+ * @param imFT Imaginary part of the Fourier transform (vector of amplitudes).
+ * @param minBin Minimum frequency bin index for optimization region.
+ *               Use -1 to optimize over all bins. Default: -1.
+ * @param maxBin Maximum frequency bin index for optimization region.
+ *               Use -1 to optimize over all bins. Default: -1.
+ *
+ * <p><b>Usage example:</b>
+ * @code
+ * std::vector<Double_t> real, imag;
+ * // ... fill real and imag with FFT results ...
+ * PFTPhaseCorrection phCorr(real, imag, 10, 500);
+ * phCorr.Minimize();
+ * Double_t c0 = phCorr.GetPhaseCorrectionParam(0);
+ * Double_t c1 = phCorr.GetPhaseCorrectionParam(1);
+ * @endcode
  */
 PFTPhaseCorrection::PFTPhaseCorrection(std::vector<Double_t> &reFT, std::vector<Double_t> &imFT, const Int_t minBin, const Int_t maxBin) :
   fReal(reFT), fImag(imFT), fMinBin(minBin), fMaxBin(maxBin)
@@ -86,7 +118,28 @@ PFTPhaseCorrection::PFTPhaseCorrection(std::vector<Double_t> &reFT, std::vector<
 // Minimize (public)
 //--------------------------------------------------------------------------
 /**
+ * <p>Performs phase correction minimization using Minuit2.
  *
+ * <p>This method finds optimal phase parameters (c₀, c₁) by minimizing
+ * a combined entropy-penalty functional:
+ * - <b>Entropy term:</b> Promotes smooth, concentrated real spectra
+ * - <b>Penalty term:</b> Penalizes negative values in real spectrum
+ *
+ * <p>The phase correction is: φ(ω) = c₀ + c₁·ω where:
+ * - c₀ corrects constant phase offset (time-zero uncertainty)
+ * - c₁ corrects linear phase dispersion (detector timing effects)
+ *
+ * <p><b>Algorithm:</b>
+ * 1. Initialize Minuit2 parameters with starting values (c₀=0, c₁=0)
+ * 2. Minimize entropy + γ×penalty using MnMinimize
+ * 3. Store optimal parameters in fPh_c0, fPh_c1
+ *
+ * <p>After calling this method, use GetPhaseCorrectionParam() to retrieve
+ * the optimal phase parameters and IsValid() to check success.
+ *
+ * @see GetPhaseCorrectionParam()
+ * @see SetGamma()
+ * @see IsValid()
  */
 void PFTPhaseCorrection::Minimize()
 {
@@ -118,7 +171,26 @@ void PFTPhaseCorrection::Minimize()
 // GetPhaseCorrectionParam (public)
 //--------------------------------------------------------------------------
 /**
+ * <p>Retrieves optimized phase correction parameter by index.
  *
+ * <p>Returns the phase parameter determined by Minimize(). The phase
+ * correction formula is: φ(ω) = c₀ + c₁·(ω/ω_max)
+ *
+ * @param idx Parameter index:
+ *            - 0: Returns c₀ (constant phase offset in radians)
+ *            - 1: Returns c₁ (linear phase dispersion coefficient)
+ *
+ * @return Phase correction parameter value, or 0.0 if idx is invalid.
+ *
+ * <p><b>Example:</b>
+ * @code
+ * phCorr.Minimize();
+ * if (phCorr.IsValid()) {
+ *   Double_t c0 = phCorr.GetPhaseCorrectionParam(0);
+ *   Double_t c1 = phCorr.GetPhaseCorrectionParam(1);
+ *   std::cout << "Phase: " << c0 << " + " << c1 << "*w" << std::endl;
+ * }
+ * @endcode
  */
 Double_t PFTPhaseCorrection::GetPhaseCorrectionParam(UInt_t idx)
 {
@@ -138,7 +210,15 @@ Double_t PFTPhaseCorrection::GetPhaseCorrectionParam(UInt_t idx)
 // GetMinimum (public)
 //--------------------------------------------------------------------------
 /**
+ * <p>Returns the minimum value of the optimization functional.
  *
+ * <p>This value represents the combined entropy-penalty functional at
+ * the optimal phase parameters. Lower values indicate better phase
+ * correction. A value of -1.0 indicates failed minimization.
+ *
+ * @return Minimum functional value, or -1.0 if minimization failed.
+ *
+ * <p><b>Note:</b> Always check IsValid() before using this value.
  */
 Double_t PFTPhaseCorrection::GetMinimum()
 {
@@ -154,7 +234,16 @@ Double_t PFTPhaseCorrection::GetMinimum()
 // Init (private)
 //--------------------------------------------------------------------------
 /**
+ * <p>Initializes phase correction object with default values.
  *
+ * <p>Sets initial state:
+ * - fValid = true (object ready for use)
+ * - fPh_c0 = 0.0 (no constant phase offset)
+ * - fPh_c1 = 0.0 (no linear phase dispersion)
+ * - fGamma = 1.0 (equal weighting of entropy and penalty)
+ * - fMin = -1.0 (no minimization performed yet)
+ *
+ * <p>Called by constructors to establish consistent initial state.
  */
 void PFTPhaseCorrection::Init()
 {
@@ -169,7 +258,20 @@ void PFTPhaseCorrection::Init()
 // CalcPhasedFT (private)
 //--------------------------------------------------------------------------
 /**
+ * <p>Calculates phase-corrected Fourier transform.
  *
+ * <p>Applies the phase correction φ(ω) = c₀ + c₁·(i/N) to the complex
+ * Fourier spectrum using rotation in the complex plane:
+ * - F'_re(ω) = F_re(ω)·cos(φ) - F_im(ω)·sin(φ)
+ * - F'_im(ω) = F_re(ω)·sin(φ) + F_im(ω)·cos(φ)
+ *
+ * <p>The corrected spectra are stored in fRealPh and fImagPh.
+ *
+ * <p><b>Physics:</b> This rotation compensates for:
+ * - Time-zero offset (constant phase c₀)
+ * - Timing dispersion (linear phase c₁)
+ *
+ * @see CalcRealPhFTDerivative()
  */
 void PFTPhaseCorrection::CalcPhasedFT() const
 {
@@ -188,7 +290,19 @@ void PFTPhaseCorrection::CalcPhasedFT() const
 // CalcRealPhFTDerivative (private)
 //--------------------------------------------------------------------------
 /**
+ * <p>Calculates first derivative of phase-corrected real Fourier spectrum.
  *
+ * <p>Computes the finite difference derivative: dF/dω ≈ F(i+1) - F(i)
+ * for all interior points. Boundary points are set to 1.0.
+ *
+ * <p>The derivative is used in the entropy calculation, where smooth
+ * spectra (small derivatives) have lower entropy and are preferred by
+ * the optimization algorithm.
+ *
+ * <p><b>Purpose:</b> Entropy minimization favors phase corrections that
+ * produce smooth, well-resolved real spectra with minimal oscillations.
+ *
+ * @see Entropy()
  */
 void PFTPhaseCorrection::CalcRealPhFTDerivative() const
 {
@@ -205,7 +319,22 @@ void PFTPhaseCorrection::CalcRealPhFTDerivative() const
 // Penalty (private)
 //--------------------------------------------------------------------------
 /**
+ * <p>Calculates penalty term for negative real spectrum values.
  *
+ * <p>Computes: P = γ·Σ[F_re(ω)² for all F_re(ω) < 0]
+ *
+ * <p>The penalty term penalizes negative values in the real Fourier
+ * spectrum, since physically meaningful absorption-mode spectra should
+ * be predominantly positive. The γ parameter controls the relative
+ * importance of this constraint.
+ *
+ * @return Weighted penalty value (γ×sum of squared negative amplitudes)
+ *
+ * <p><b>Note:</b> Only bins in range [fMinBin, fMaxBin) are considered,
+ * allowing focus on signal-containing regions.
+ *
+ * @see SetGamma()
+ * @see Entropy()
  */
 Double_t PFTPhaseCorrection::Penalty() const
 {
@@ -223,7 +352,26 @@ Double_t PFTPhaseCorrection::Penalty() const
 // Entropy (private)
 //--------------------------------------------------------------------------
 /**
+ * <p>Calculates Shannon entropy of the real spectrum derivative.
  *
+ * <p>Computes: S = -Σ[p_i·ln(p_i)] where p_i = |dF/dω|_i / Σ|dF/dω|
+ *
+ * <p>This entropy measure quantifies the "smoothness" of the real
+ * Fourier spectrum:
+ * - <b>Low entropy:</b> Concentrated, smooth spectrum (desired)
+ * - <b>High entropy:</b> Dispersed, noisy spectrum (undesired)
+ *
+ * <p><b>Physical interpretation:</b> Correct phase alignment produces
+ * sharp, well-defined peaks with smooth baselines, resulting in
+ * concentrated derivatives and low entropy.
+ *
+ * @return Shannon entropy of normalized derivative distribution
+ *
+ * <p><b>Note:</b> Small values (< 10⁻¹⁵) are skipped to avoid
+ * numerical issues with log(0).
+ *
+ * @see CalcRealPhFTDerivative()
+ * @see Penalty()
  */
 Double_t PFTPhaseCorrection::Entropy() const
 {
@@ -250,7 +398,29 @@ Double_t PFTPhaseCorrection::Entropy() const
 // operator() (private)
 //--------------------------------------------------------------------------
 /**
+ * <p>Objective function for Minuit2 minimization (FCNBase interface).
  *
+ * <p>Evaluates the combined entropy-penalty functional:
+ * f(c₀, c₁) = S + γ·P
+ * where:
+ * - S = entropy of phase-corrected real spectrum derivative
+ * - P = penalty for negative values in real spectrum
+ * - γ = balancing parameter (fGamma)
+ *
+ * @param par Parameter vector: [0]=c₀ (constant phase), [1]=c₁ (linear phase)
+ * @return Combined functional value to be minimized
+ *
+ * <p><b>Algorithm flow:</b>
+ * 1. Apply phase correction with parameters par
+ * 2. Calculate phase-corrected Fourier transform
+ * 3. Compute spectrum derivative
+ * 4. Return entropy + penalty
+ *
+ * <p>Called repeatedly by Minuit2 during optimization.
+ *
+ * @see Entropy()
+ * @see Penalty()
+ * @see Minimize()
  */
 double PFTPhaseCorrection::operator()(const std::vector<double> &par) const
 {
@@ -271,15 +441,43 @@ double PFTPhaseCorrection::operator()(const std::vector<double> &par) const
 // Constructor
 //--------------------------------------------------------------------------
 /**
- * <p>Constructor.
+ * <p>Constructor for Fourier transform engine.
  *
- * \param data data histogram
+ * <p>Creates a PFourier object configured to transform μSR time-domain
- * \param unitTag tag telling in which units the Fourier transform shall be represented. Possible tags are:
+ * data to frequency/field domain. The constructor validates inputs,
- *                FOURIER_UNIT_GAUSS, FOURIER_UNIT_TESLA, FOURIER_UNIT_FREQ, FOURIER_UNIT_CYCLES
+ * allocates FFTW resources, and prepares for FFT computation.
- * \param startTime start time of the data time window
+ *
- * \param endTime end time of the data time window
+ * @param data TH1F histogram containing time-domain μSR signal.
- * \param dcCorrected if true, removed DC offset from signal before Fourier transformation, otherwise not
+ *             X-axis: time in microseconds, Y-axis: asymmetry or counts.
- * \param zeroPaddingPower if set to values > 0, there will be zero padding up to 2^zeroPaddingPower
+ * @param unitTag Output units for frequency axis:
+ *        - FOURIER_UNIT_GAUSS (1): Magnetic field in Gauss
+ *        - FOURIER_UNIT_TESLA (2): Magnetic field in Tesla
+ *        - FOURIER_UNIT_FREQ (3): Frequency in MHz
+ *        - FOURIER_UNIT_CYCLES (4): Angular frequency in Mc/s (Mrad/s)
+ * @param startTime Start of transform window in μs (0.0 = from t₀). Default: 0.0
+ * @param endTime End of transform window in μs (0.0 = to end). Default: 0.0
+ * @param dcCorrected If true, remove DC offset before FFT to eliminate
+ *                    zero-frequency artifact. Default: false
+ * @param zeroPaddingPower Zero-pad to 2^N points where N=zeroPaddingPower.
+ *                         Value 0 disables padding. Padding improves frequency
+ *                         interpolation but doesn't add information. Default: 0
+ *
+ * <p><b>Unit conversions (using γ_μ/2π = 0.01355 MHz/G):</b>
+ * - Gauss: B(G) = ω(MHz) / 0.01355
+ * - Tesla: B(T) = ω(MHz) / 135.54
+ *
+ * <p><b>Example usage:</b>
+ * @code
+ * TH1F *timeHisto = ...; // μSR time spectrum
+ * PFourier ft(timeHisto, FOURIER_UNIT_GAUSS, 0.1, 10.0, true, 12);
+ * ft.Transform(F_APODIZATION_WEAK);
+ * TH1F *powerSpec = ft.GetPowerFourier();
+ * @endcode
+ *
+ * <p>After construction, check IsValid() before calling Transform().
+ *
+ * @see Transform()
+ * @see IsValid()
  */
 PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTime, Bool_t dcCorrected, UInt_t zeroPaddingPower, Bool_t useFFTW) :
                    fData(data), fUnitTag(unitTag), fStartTime(startTime), fEndTime(endTime),
@@ -429,7 +627,16 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
 // Destructor
 //--------------------------------------------------------------------------
 /**
- * <p>Destructor
+ * <p>Destructor - cleans up FFTW resources.
+ *
+ * <p>Releases all FFTW-allocated memory and destroys the FFT plan.
+ * Proper cleanup is essential to avoid memory leaks, as FFTW uses
+ * special aligned memory allocation.
+ *
+ * <p><b>Resources freed:</b>
+ * - fFFTwPlan: FFTW execution plan
+ * - fIn: Input array (time-domain data)
+ * - fOut: Output array (frequency-domain data)
  */
 PFourier::~PFourier()
 {
@@ -458,10 +665,37 @@ PFourier::~PFourier()
 // Transform (public)
 //--------------------------------------------------------------------------
 /**
- * <p>Carries out the Fourier transform. It is assumed that fStartTime is the time zero
+ * <p>Executes the Fourier transform with optional apodization.
- * for the Fourier frame. Hence if fStartTime != 0.0 the phase shift will be corrected.
  *
- * \param apodizationTag 0=no apod., 1=weak apod., 2=medium apod., 3=strong apod.
+ * <p>This is the main workhorse method that performs the complete FFT
+ * pipeline:
+ * 1. Prepare input data (DC correction, zero padding)
+ * 2. Apply apodization window if requested
+ * 3. Execute FFTW plan (compute FFT)
+ * 4. Correct phase for non-zero start time
+ *
+ * <p><b>Phase correction:</b> If fStartTime ≠ 0, applies phase shift
+ * e^(i·ω·t₀) to account for the time-zero offset. This ensures proper
+ * phase relationships in the frequency spectrum.
+ *
+ * <p><b>Apodization:</b> Window functions reduce spectral leakage
+ * (Gibbs phenomenon) at the cost of slight peak broadening:
+ * - NONE: Best amplitude accuracy, worst leakage
+ * - WEAK: Minimal smoothing
+ * - MEDIUM: Balanced (recommended for most cases)
+ * - STRONG: Maximum leakage suppression
+ *
+ * @param apodizationTag Window function strength:
+ *        - 0 or F_APODIZATION_NONE: No windowing
+ *        - 1 or F_APODIZATION_WEAK: Gentle taper
+ *        - 2 or F_APODIZATION_MEDIUM: Moderate taper
+ *        - 3 or F_APODIZATION_STRONG: Heavy taper
+ *
+ * <p>After calling Transform(), retrieve results using:
+ * GetRealFourier(), GetImaginaryFourier(), GetPowerFourier(), GetPhaseFourier()
+ *
+ * @see PrepareFFTwInputData()
+ * @see GetPowerFourier()
  */
 void PFourier::Transform(UInt_t apodizationTag)
 {
@@ -542,7 +776,25 @@ void PFourier::SetUseFFTW(const Bool_t flag)
 // GetMaxFreq (public)
 //--------------------------------------------------------------------------
 /**
- * <p>returns the maximal frequency in units choosen, i.e. Gauss, Tesla, MHz, Mc/s
+ * <p>Returns the maximum frequency (Nyquist frequency) of the spectrum.
+ *
+ * <p>The Nyquist frequency is the highest frequency that can be
+ * unambiguously represented given the time resolution:
+ * f_Nyquist = 1 / (2·Δt)
+ *
+ * <p>Higher frequencies are aliased back into the [0, f_Nyquist] range,
+ * so the useful spectrum extends from 0 to f_Nyquist.
+ *
+ * @return Maximum frequency in units specified by fUnitTag:
+ *         - Gauss (if FOURIER_UNIT_GAUSS)
+ *         - Tesla (if FOURIER_UNIT_TESLA)
+ *         - MHz (if FOURIER_UNIT_FREQ)
+ *         - Mc/s (if FOURIER_UNIT_CYCLES)
+ *
+ * <p><b>Example:</b> For time resolution Δt = 0.01 μs:
+ * f_Nyquist = 1/(2×0.01) = 50 MHz ≈ 3690 Gauss
+ *
+ * @see GetResolution()
  */
 Double_t PFourier::GetMaxFreq()
 {
@@ -555,9 +807,36 @@ Double_t PFourier::GetMaxFreq()
 // GetRealFourier (public)
 //--------------------------------------------------------------------------
 /**
- * <p>returns the real part Fourier as a histogram.
+ * <p>Returns real part of Fourier transform as ROOT histogram.
  *
- * \param scale normalisation factor
+ * <p>The real component represents the in-phase (cosine) projection
+ * of the signal. Without phase correction, both positive and negative
+ * values typically appear due to arbitrary phase offsets.
+ *
+ * <p><b>Interpretation:</b>
+ * - Phase-corrected: Absorption-mode spectrum (predominantly positive)
+ * - Uncorrected: Mixed-phase spectrum (positive + negative regions)
+ *
+ * <p>Histogram covers [0, f_Nyquist] with N/2 bins where N is the
+ * FFT size (including zero padding if applied).
+ *
+ * @param scale Multiplication factor applied to all amplitudes.
+ *              Use for normalization or unit conversion. Default: 1.0
+ *
+ * @return Pointer to newly allocated TH1F histogram. <b>Caller owns
+ *         the histogram and must delete it.</b> Returns nullptr on error.
+ *
+ * <p><b>Example:</b>
+ * @code
+ * ft.Transform(F_APODIZATION_WEAK);
+ * TH1F *realSpec = ft.GetRealFourier(1.0);
+ * realSpec->Draw();
+ * // ... use histogram ...
+ * delete realSpec;
+ * @endcode
+ *
+ * @see GetImaginaryFourier()
+ * @see GetPhaseOptRealFourier()
  */
 TH1F* PFourier::GetRealFourier(const Double_t scale)
 {
@@ -602,16 +881,54 @@ TH1F* PFourier::GetRealFourier(const Double_t scale)
 // GetPhaseOptRealFourier (public, static)
 //--------------------------------------------------------------------------
 /**
- * <p>returns the phase corrected real Fourier transform.
+ * <p>Returns phase-optimized real Fourier spectrum (static method).
  *
- * \return the TH1F histo of the phase 'optimzed' real Fourier transform.
+ * <p>This static method performs automatic phase correction on a complex
+ * Fourier spectrum to maximize the real component (absorption mode) while
+ * minimizing the imaginary component. The algorithm finds optimal phase
+ * parameters φ(ω) = c₀ + c₁·ω that produce the cleanest real spectrum.
  *
- * \param re real part Fourier histogram
+ * <p><b>Why phase correction?</b> Raw FFT produces mixed-phase spectra
- * \param im imaginary part Fourier histogram
+ * with arbitrary phase offsets due to:
- * \param phase return value of the optimal phase dispersion phase[0]+phase[1]*i/N
+ * - Uncertain time-zero (t₀) determination
- * \param scale normalisation factor
+ * - Detector timing offsets
- * \param min minimal freq / field from which to optimise. Given in the choosen unit.
+ * - Signal propagation delays
- * \param max maximal freq / field up to which to optimise. Given in the choosen unit.
+ *
+ * <p><b>Algorithm:</b>
+ * 1. Extract Re/Im data in specified frequency range [min, max]
+ * 2. Use PFTPhaseCorrection to minimize entropy+penalty functional
+ * 3. Apply optimal phase rotation to entire spectrum
+ * 4. Return phase-corrected real part
+ *
+ * <p><b>Optimization range [min, max]:</b> Restrict to signal-containing
+ * regions for best results. Including noisy baseline degrades optimization.
+ *
+ * @param re Real part of Fourier transform (TH1F histogram)
+ * @param im Imaginary part of Fourier transform (TH1F histogram)
+ * @param phase Output parameter: phase[0]=c₀, phase[1]=c₁ (vector resized to 2)
+ * @param scale Amplitude scaling factor. Default: 1.0
+ * @param min Minimum frequency/field for optimization window.
+ *            Use -1.0 to include all frequencies. Default: -1.0
+ * @param max Maximum frequency/field for optimization window.
+ *            Use -1.0 to include all frequencies. Default: -1.0
+ *
+ * @return Pointer to newly allocated phase-corrected TH1F histogram.
+ *         <b>Caller owns and must delete.</b> Returns nullptr on error.
+ *
+ * <p><b>Example:</b>
+ * @code
+ * TH1F *real = ft.GetRealFourier();
+ * TH1F *imag = ft.GetImaginaryFourier();
+ * std::vector<Double_t> phaseParams;
+ * TH1F *phCorrected = PFourier::GetPhaseOptRealFourier(real, imag, phaseParams, 1.0, 5.0, 50.0);
+ * std::cout << "Phase: " << phaseParams[0] << " + " << phaseParams[1] << "*w" << std::endl;
+ * phCorrected->Draw();
+ * delete real; delete imag; delete phCorrected;
+ * @endcode
+ *
+ * @see PFTPhaseCorrection
+ * @see GetRealFourier()
+ * @see GetImaginaryFourier()
  */
 TH1F* PFourier::GetPhaseOptRealFourier(const TH1F *re, const TH1F *im, std::vector<Double_t> &phase,
                                        const Double_t scale, const Double_t min, const Double_t max)
@@ -703,9 +1020,29 @@ TH1F* PFourier::GetPhaseOptRealFourier(const TH1F *re, const TH1F *im, std::vect
 // GetImaginaryFourier (public)
 //--------------------------------------------------------------------------
 /**
- * <p>returns the imaginary part Fourier as a histogram.
+ * <p>Returns imaginary part of Fourier transform as ROOT histogram.
  *
- * \param scale normalisation factor
+ * <p>The imaginary component represents the out-of-phase (sine)
+ * projection of the signal. Ideally zero for perfect phase correction,
+ * but contains signal information when phase is uncorrected.
+ *
+ * <p><b>Uses:</b>
+ * - Phase correction optimization (minimize imaginary component)
+ * - Dispersion-mode spectroscopy
+ * - Computing power spectrum: |F|² = Re² + Im²
+ * - Phase spectrum calculation: φ = atan(Im/Re)
+ *
+ * <p>Histogram covers [0, f_Nyquist] with N/2 bins where N is the
+ * FFT size (including zero padding if applied).
+ *
+ * @param scale Multiplication factor applied to all amplitudes.
+ *              Use for normalization or unit conversion. Default: 1.0
+ *
+ * @return Pointer to newly allocated TH1F histogram. <b>Caller owns
+ *         the histogram and must delete it.</b> Returns nullptr on error.
+ *
+ * @see GetRealFourier()
+ * @see GetPhaseOptRealFourier()
  */
 TH1F* PFourier::GetImaginaryFourier(const Double_t scale)
 {
@@ -750,9 +1087,37 @@ TH1F* PFourier::GetImaginaryFourier(const Double_t scale)
 // GetPowerFourier (public)
 //--------------------------------------------------------------------------
 /**
- * <p>returns the Fourier power spectrum as a histogram.
+ * <p>Returns power spectrum |F(ω)| as ROOT histogram.
  *
- * \param scale normalisation factor
+ * <p>The power spectrum shows signal magnitude at each frequency,
+ * computed as: |F(ω)| = √(Re² + Im²)
+ *
+ * <p><b>Advantages:</b>
+ * - Always positive (easier interpretation)
+ * - Phase-independent (no phase correction needed)
+ * - Directly shows signal strength vs. frequency
+ * - Ideal for identifying all frequency components
+ *
+ * <p><b>Applications:</b>
+ * - Quick survey of frequency content
+ * - Field distribution analysis in superconductors
+ * - Multiple muon site identification
+ * - TF-μSR field measurements
+ *
+ * <p>Histogram covers [0, f_Nyquist] with N/2 bins where N is the
+ * FFT size (including zero padding if applied).
+ *
+ * @param scale Multiplication factor applied to all amplitudes.
+ *              Use for normalization or unit conversion. Default: 1.0
+ *
+ * @return Pointer to newly allocated TH1F histogram. <b>Caller owns
+ *         the histogram and must delete it.</b> Returns nullptr on error.
+ *
+ * <p><b>Note:</b> This returns the amplitude |F|, not the power |F|².
+ * For power values, square the histogram contents.
+ *
+ * @see GetRealFourier()
+ * @see GetPhaseFourier()
  */
 TH1F* PFourier::GetPowerFourier(const Double_t scale)
 {
@@ -797,9 +1162,39 @@ TH1F* PFourier::GetPowerFourier(const Double_t scale)
 // GetPhaseFourier (public)
 //--------------------------------------------------------------------------
 /**
- * <p>returns the Fourier phase spectrum as a histogram.
+ * <p>Returns phase spectrum φ(ω) = arg(F) as ROOT histogram.
  *
- * \param scale normalisation factor
+ * <p>The phase spectrum shows the phase angle of the complex Fourier
+ * transform at each frequency: φ(ω) = atan2(Im, Re)
+ *
+ * <p><b>Range:</b> Phase values are in [-π, +π] radians, with proper
+ * quadrant handling using the two-argument arctangent.
+ *
+ * <p><b>Interpretation:</b>
+ * - Constant phase: All frequencies in phase (simple exponential decay)
+ * - Linear phase: Time-zero offset (shift theorem)
+ * - Quadratic phase: Dispersion or chirp
+ * - Random phase: Noise
+ *
+ * <p><b>Applications:</b>
+ * - Diagnosing time-zero problems
+ * - Verifying phase correction quality
+ * - Identifying signal vs. noise regions
+ *
+ * <p>Histogram covers [0, f_Nyquist] with N/2 bins where N is the
+ * FFT size (including zero padding if applied).
+ *
+ * @param scale Multiplication factor applied to phase values.
+ *              Typically 1.0 (radians) or 180/π (degrees). Default: 1.0
+ *
+ * @return Pointer to newly allocated TH1F histogram. <b>Caller owns
+ *         the histogram and must delete it.</b> Returns nullptr on error.
+ *
+ * <p><b>Note:</b> Phase is undefined where amplitude is zero. In practice,
+ * phase values in noise regions are meaningless.
+ *
+ * @see GetRealFourier()
+ * @see GetImaginaryFourier()
  */
 TH1F* PFourier::GetPhaseFourier(const Double_t scale)
 {
@@ -865,11 +1260,32 @@ TH1F* PFourier::GetPhaseFourier(const Double_t scale)
 // PrepareFFTwInputData (private)
 //--------------------------------------------------------------------------
 /**
- * <p>Feeds the Fourier data and apply the apodization.
+ * <p>Prepares time-domain data for FFT with optional DC correction and apodization.
  *
- * \param apodizationTag apodization tag. Possible are currently: F_APODIZATION_NONE = no apodization,
+ * <p>This method performs essential preprocessing steps:
- *                       F_APODIZATION_WEAK = weak apodization, F_APODIZATION_MEDIUM = intermediate apodization,
+ * 1. Locates time-zero bin in the histogram
- *                       F_APODIZATION_STRONG = strong apodization
+ * 2. Extracts data from [fStartTime, fEndTime] window
+ * 3. Removes DC offset (if fDCCorrected=true) for baseline correction
+ * 4. Zero-pads to fNoOfBins (power of 2 for optimal FFT performance)
+ * 5. Applies apodization window to reduce spectral leakage
+ *
+ * <p><b>DC Correction:</b> Subtracts mean value from time signal to
+ * remove constant offset, which would otherwise create a large artifact
+ * at zero frequency.
+ *
+ * <p><b>Zero Padding:</b> Increases frequency resolution by interpolation
+ * without adding information content. Often used for smoother spectra.
+ *
+ * @param apodizationTag Window function strength:
+ *        - F_APODIZATION_NONE (1): Rectangular window (no smoothing)
+ *        - F_APODIZATION_WEAK (2): Gentle roll-off at edges
+ *        - F_APODIZATION_MEDIUM (3): Moderate smoothing
+ *        - F_APODIZATION_STRONG (4): Heavy smoothing for maximum leakage suppression
+ *
+ * <p>Result is stored in fIn[] array ready for FFTW execution.
+ *
+ * @see ApodizeData()
+ * @see Transform()
  */
 void PFourier::PrepareFFTwInputData(UInt_t apodizationTag)
 {
@@ -928,11 +1344,33 @@ void PFourier::PrepareFFTwInputData(UInt_t apodizationTag)
 // ApodizeData (private)
 //--------------------------------------------------------------------------
 /**
- * <p>Carries out the appodization of the data.
+ * <p>Applies apodization (windowing) to time-domain data before FFT.
  *
- * \param apodizationTag apodization tag. Possible are currently: F_APODIZATION_NONE = no apodization,
+ * <p><b>Purpose:</b> Apodization multiplies the time signal by a smooth
- *                       F_APODIZATION_WEAK = weak apodization, F_APODIZATION_MEDIUM = intermediate apodization,
+ * window function that tapers to zero at the edges, reducing the Gibbs
- *                       F_APODIZATION_STRONG = strong apodization
+ * phenomenon (spectral leakage) caused by finite observation windows.
+ *
+ * <p><b>Trade-off:</b>
+ * - <b>Advantage:</b> Sharper, cleaner frequency peaks with less sidelobes
+ * - <b>Disadvantage:</b> Slightly broader main peaks, reduced amplitude accuracy
+ *
+ * <p><b>Window functions:</b> Polynomial windows of the form:
+ * w(t) = Σ c_j·(t/T)^(2j) where t ∈ [0,T]
+ *
+ * <p>Three predefined window strengths with optimized coefficients:
+ * - <b>Weak:</b> Minimal smoothing, preserves amplitude
+ * - <b>Medium:</b> Balanced smoothing for general use
+ * - <b>Strong:</b> Maximum leakage suppression for crowded spectra
+ *
+ * @param apodizationTag Window function type:
+ *        - F_APODIZATION_NONE: No windowing (rectangular)
+ *        - F_APODIZATION_WEAK: Gentle taper
+ *        - F_APODIZATION_MEDIUM: Moderate taper
+ *        - F_APODIZATION_STRONG: Heavy taper
+ *
+ * <p>Window is applied in-place to fIn[] array.
+ *
+ * @see PrepareFFTwInputData()
  */
 void PFourier::ApodizeData(Int_t apodizationTag) {