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

  PFourier.cpp

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

#include <cmath>

#include <iostream>
#include <iomanip>

#include "TF1.h"
#include "TAxis.h"

#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnMinimize.h"

#include "PMusr.h"
#include "PFourier.h"

#define PI      3.14159265358979312
#define PI_HALF 1.57079632679489656

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// PFTPhaseCorrection
//--------------------------------------------------------------------------
// 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)
{
  Init();
}

//--------------------------------------------------------------------------
// 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)
{
  Init();

  Int_t realSize = static_cast<Int_t>(fReal.size());
  if (fMinBin == -1)
    fMinBin = 0;
  if (fMaxBin == -1)
    fMaxBin = realSize;
  if (fMaxBin > realSize)
    fMaxBin = realSize;

  fRealPh.resize(fReal.size());
  fImagPh.resize(fReal.size());
}

//--------------------------------------------------------------------------
// 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()
{
  // create Minuit2 parameters
  ROOT::Minuit2::MnUserParameters upar;

  upar.Add("c0", fPh_c0, 2.0);
  upar.Add("c1", fPh_c1, 2.0);

  // create minimizer
  ROOT::Minuit2::MnMinimize mn_min(*this, upar);

  // minimize
  ROOT::Minuit2::FunctionMinimum min = mn_min();

  if (min.IsValid()) {
    fPh_c0 = min.UserState().Value("c0");
    fPh_c1 = min.UserState().Value("c1");
    fMin = min.Fval();
  } else {
    fMin = -1.0;
    fValid = false;
    std::cout << std::endl << ">> **WARNING** minimize failed to find a minimum for the real FT phase correction ..." << std::endl;
    return;
  }
}

//--------------------------------------------------------------------------
// 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)
{
  Double_t result=0.0;

  if (idx == 0)
    result = fPh_c0;
  else if (idx == 1)
    result = fPh_c1;
  else
    std::cerr << ">> **ERROR** requested phase correction parameter with index=" << idx << " does not exist!" << std::endl;

  return result;
}

//--------------------------------------------------------------------------
// 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()
{
  if (!fValid) {
    std::cerr << ">> **ERROR** requested minimum is invalid!" << std::endl;
    return -1.0;
  }

  return fMin;
}

//--------------------------------------------------------------------------
// 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()
{
  fValid = true;
  fPh_c0 = 0.0;
  fPh_c1 = 0.0;
  fGamma = 1.0;
  fMin   = -1.0;
}

//--------------------------------------------------------------------------
// 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
{
  Double_t phi=0.0;
  Double_t w=0.0;

  for (UInt_t i=0; i<fRealPh.size(); i++) {
    w = static_cast<Double_t>(i) / static_cast<Double_t>(fReal.size());
    phi = fPh_c0 + fPh_c1 * w;
    fRealPh[i] = fReal[i]*cos(phi) - fImag[i]*sin(phi);
    fImagPh[i] = fReal[i]*sin(phi) + fImag[i]*cos(phi);
  }
}

//--------------------------------------------------------------------------
// 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
{
  fRealPhD.resize(fRealPh.size());

  fRealPhD[0] = 1.0;
  fRealPhD[fRealPh.size()-1] = 1.0;

  for (UInt_t i=1; i<fRealPh.size()-1; i++)
    fRealPhD[i] = fRealPh[i+1]-fRealPh[i];
}

//--------------------------------------------------------------------------
// 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
{
  Double_t penalty = 0.0;

  for (UInt_t i=fMinBin; i<fMaxBin; i++) {
    if (fRealPh[i] < 0.0)
      penalty += fRealPh[i]*fRealPh[i];
  }

  return fGamma*penalty;
}

//--------------------------------------------------------------------------
// 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
{
  Double_t norm = 0.0;

  for (UInt_t i=fMinBin; i<fMaxBin; i++)
    norm += fabs(fRealPhD[i]);

  Double_t entropy = 0.0;
  Double_t dval = 0.0, hh = 0.0;

  for (UInt_t i=fMinBin; i<fMaxBin; i++) {
    dval = fabs(fRealPhD[i]);
    if (dval > 1.0e-15) {
      hh = dval / norm;
      entropy -= hh * log(hh);
    }
  }

  return entropy;
}

//--------------------------------------------------------------------------
// 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
{
  // par : [0]: c0, [1]: c1

  fPh_c0 = par[0];
  fPh_c1 = par[1];

  CalcPhasedFT();
  CalcRealPhFTDerivative();

  return Entropy()+Penalty();
}

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// PFourier
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
 * <p>Constructor for Fourier transform engine.
 *
 * <p>Creates a PFourier object configured to transform μSR time-domain
 * data to frequency/field domain. The constructor validates inputs,
 * allocates FFTW resources, and prepares for FFT computation.
 *
 * @param data TH1F histogram containing time-domain μSR signal.
 *             X-axis: time in microseconds, Y-axis: asymmetry or counts.
 * @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) :
                   fData(data), fUnitTag(unitTag), fStartTime(startTime), fEndTime(endTime),
                   fDCCorrected(dcCorrected), fZeroPaddingPower(zeroPaddingPower)
{
  // some necessary checks and initialization
  if (fData == nullptr) {
    std::cerr << std::endl << "**ERROR** PFourier::PFourier: no valid data" << std::endl << std::endl;
    fValid = false;
    return;
  }

  fValid = true;
  fIn  = nullptr;
  fOut = nullptr;
//as  fPhCorrectedReFT = 0;

  fApodization = F_APODIZATION_NONE;

  // calculate time resolution in (us)
  fTimeResolution = fData->GetBinWidth(1);

  // if endTime == 0 set it to the last time slot
  if (fEndTime == 0.0) {
    Int_t last = fData->GetNbinsX()-1;
    fEndTime = fData->GetBinCenter(last);
  }

  // swap start and end time if necessary
  if (fStartTime > fEndTime) {
    Double_t keep = fStartTime;
    fStartTime = fEndTime;
    fEndTime = keep;
  }

  // calculate start and end bin
  fNoOfData = static_cast<UInt_t>((fEndTime-fStartTime)/fTimeResolution);

  // check if zero padding is whished
  if (fZeroPaddingPower > 0) {
    UInt_t noOfBins = static_cast<UInt_t>(pow(2.0, static_cast<Double_t>(fZeroPaddingPower)));
    if (noOfBins > fNoOfData)
      fNoOfBins = noOfBins;
    else
      fNoOfBins = fNoOfData;
  } else {
    fNoOfBins = fNoOfData;
  }

  // calculate fourier resolution, depending on the units
  Double_t resolution = 1.0/(fTimeResolution*fNoOfBins); // in MHz
  switch (fUnitTag) {
    case FOURIER_UNIT_GAUSS:
      fResolution = resolution/GAMMA_BAR_MUON;
      break;
    case FOURIER_UNIT_TESLA:
      fResolution = 1e-4*resolution/GAMMA_BAR_MUON;
      break;
    case FOURIER_UNIT_FREQ:
      fResolution = resolution;
      break;
    case FOURIER_UNIT_CYCLES:
      fResolution = 2.0*PI*resolution;
      break;
    default:
      fValid = false;
      return;
  }

  // allocate necessary memory
  fIn  = static_cast<fftw_complex *>(fftw_malloc(sizeof(fftw_complex)*fNoOfBins));
  fOut = static_cast<fftw_complex *>(fftw_malloc(sizeof(fftw_complex)*fNoOfBins));

  // check if memory allocation has been successful
  if ((fIn == nullptr) || (fOut == nullptr)) {
    fValid = false;
    return;
  }

  // get the FFTW3 plan (see FFTW3 manual)
  fFFTwPlan = fftw_plan_dft_1d(static_cast<Int_t>(fNoOfBins), fIn, fOut, FFTW_FORWARD, FFTW_ESTIMATE);

  // check if a valid plan has been generated
  if (!fFFTwPlan) {
    fValid = false;
  }
}

//--------------------------------------------------------------------------
// 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()
{
  if (fFFTwPlan)
    fftw_destroy_plan(fFFTwPlan);
  if (fIn)
    fftw_free(fIn);
  if (fOut)
    fftw_free(fOut);
//as  if (fPhCorrectedReFT)
//as    delete fPhCorrectedReFT;
}

//--------------------------------------------------------------------------
// Transform (public)
//--------------------------------------------------------------------------
/**
 * <p>Executes the Fourier transform with optional apodization.
 *
 * <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)
{
  if (!fValid)
    return;

  PrepareFFTwInputData(apodizationTag);

  fftw_execute(fFFTwPlan);

  // correct the phase for tstart != 0.0
  // find the first bin >= fStartTime
  Double_t shiftTime = 0.0;
  for (Int_t i=1; i<fData->GetXaxis()->GetNbins(); i++) {
    if (fData->GetXaxis()->GetBinCenter(i) >= fStartTime) {
      shiftTime = fData->GetXaxis()->GetBinCenter(i);
      break;
    }
  }

  Double_t phase, re, im;
  for (UInt_t i=0; i<fNoOfBins; i++) {
    phase = 2.0*PI/(fTimeResolution*fNoOfBins) * i * shiftTime;
    re =  fOut[i][0] * cos(phase) + fOut[i][1] * sin(phase);
    im = -fOut[i][0] * sin(phase) + fOut[i][1] * cos(phase);
    fOut[i][0] = re;
    fOut[i][1] = im;
  }
}

//--------------------------------------------------------------------------
// GetMaxFreq (public)
//--------------------------------------------------------------------------
/**
 * <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()
{
  UInt_t noOfFourierBins = 0;
  if (fNoOfBins % 2 == 0)
    noOfFourierBins = fNoOfBins/2;
  else
    noOfFourierBins = (fNoOfBins+1)/2;

  return fResolution*noOfFourierBins;
}

//--------------------------------------------------------------------------
// GetRealFourier (public)
//--------------------------------------------------------------------------
/**
 * <p>Returns real part of Fourier transform as ROOT histogram.
 *
 * <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)
{
  // check if valid flag is set
  if (!fValid)
    return nullptr;

  // invoke realFourier
  Char_t name[256];
  Char_t title[256];
  snprintf(name, sizeof(name), "%s_Fourier_Re", fData->GetName());
  snprintf(title, sizeof(title), "%s_Fourier_Re", fData->GetTitle());

  UInt_t noOfFourierBins = 0;
  if (fNoOfBins % 2 == 0)
    noOfFourierBins = fNoOfBins/2;
  else
    noOfFourierBins = (fNoOfBins+1)/2;

  TH1F *realFourier = new TH1F(name, title, static_cast<Int_t>(noOfFourierBins), -fResolution/2.0, static_cast<Double_t>(noOfFourierBins-1)*fResolution+fResolution/2.0);
  if (realFourier == nullptr) {
    fValid = false;
    std::cerr << std::endl << "**SEVERE ERROR** couldn't allocate memory for the real part of the Fourier transform." << std::endl;
    return nullptr;
  }

  // fill realFourier vector
  for (Int_t i=0; i<static_cast<Int_t>(noOfFourierBins); i++) {
    realFourier->SetBinContent(i+1, scale*fOut[i][0]);
    realFourier->SetBinError(i+1, 0.0);
  }
  return realFourier;
}

//--------------------------------------------------------------------------
// GetPhaseOptRealFourier (public, static)
//--------------------------------------------------------------------------
/**
 * <p>Returns phase-optimized real Fourier spectrum (static method).
 *
 * <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.
 *
 * <p><b>Why phase correction?</b> Raw FFT produces mixed-phase spectra
 * with arbitrary phase offsets due to:
 * - Uncertain time-zero (t₀) determination
 * - Detector timing offsets
 * - Signal propagation delays
 *
 * <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)
{
  if ((re == nullptr) || (im == nullptr))
    return nullptr;

  phase.resize(2); // c0, c1

  const TAxis *axis = re->GetXaxis();

  Int_t minBin = 1;
  Int_t maxBin = axis->GetNbins();
  Int_t noOfBins = axis->GetNbins();
  Double_t res = axis->GetBinWidth(1);

  // check if minimum frequency is given. If yes, get the proper minBin
  if (min != -1.0) {
    minBin = axis->FindFixBin(min);
    if ((minBin == 0) || (minBin > maxBin)) {
      minBin = 1;
      std::cerr << "**WARNING** minimum frequency/field out of range. Will adopted it." << std::endl;
    }
  }

  // check if maximum frequency is given. If yes, get the proper maxBin
  if (max != -1.0) {
    maxBin = axis->FindFixBin(max);
    if ((maxBin == 0) || (maxBin > axis->GetNbins())) {
      maxBin = axis->GetNbins();
      std::cerr << "**WARNING** maximum frequency/field out of range. Will adopted it." << std::endl;
    }
  }

  // copy the real/imag Fourier from min to max
  std::vector<Double_t> realF, imagF;
  for (Int_t i=minBin; i<=maxBin; i++) {
    realF.push_back(re->GetBinContent(i));
    imagF.push_back(im->GetBinContent(i));
  }

  // optimize real FT phase
  PFTPhaseCorrection *phCorrectedReFT = new PFTPhaseCorrection(realF, imagF);
  if (phCorrectedReFT == nullptr) {
    std::cerr << std::endl << "**SEVERE ERROR** couldn't invoke PFTPhaseCorrection object." << std::endl;
    return nullptr;
  }

  phCorrectedReFT->Minimize();
  if (!phCorrectedReFT->IsValid()) {
    std::cerr << std::endl << "**ERROR** could not find a valid phase correction minimum." << std::endl;
    return nullptr;
  }
  phase[0] = phCorrectedReFT->GetPhaseCorrectionParam(0);
  phase[1] = phCorrectedReFT->GetPhaseCorrectionParam(1);

  // clean up
  if (phCorrectedReFT) {
    delete phCorrectedReFT;
    phCorrectedReFT = nullptr;
  }
  realF.clear();
  imagF.clear();

  // invoke the real phase optimised histo to be filled. Caller is the owner!
  Char_t name[256];
  Char_t title[256];
  snprintf(name, sizeof(name), "%s_Fourier_PhOptRe", re->GetName());
  snprintf(title, sizeof(title), "%s_Fourier_PhOptRe", re->GetTitle());

  TH1F *realPhaseOptFourier = new TH1F(name, title, noOfBins, -res/2.0, static_cast<Double_t>(noOfBins-1)*res+res/2.0);
  if (realPhaseOptFourier == nullptr) {
    std::cerr << std::endl << "**SEVERE ERROR** couldn't allocate memory for the real part of the Fourier transform." << std::endl;
    return nullptr;
  }

  // fill realFourier vector
  Double_t ph;
  for (Int_t i=0; i<noOfBins; i++) {
    ph = phase[0] + phase[1] * static_cast<Double_t>(i-static_cast<Int_t>(minBin)) / static_cast<Double_t>(maxBin-minBin);
    realPhaseOptFourier->SetBinContent(i+1, scale*(re->GetBinContent(i+1)*cos(ph) - im->GetBinContent(i+1)*sin(ph)));
    realPhaseOptFourier->SetBinError(i+1, 0.0);
  }

  return realPhaseOptFourier;
}

//--------------------------------------------------------------------------
// GetImaginaryFourier (public)
//--------------------------------------------------------------------------
/**
 * <p>Returns imaginary part of Fourier transform as ROOT histogram.
 *
 * <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)
{
  // check if valid flag is set
  if (!fValid)
    return nullptr;

  // invoke imaginaryFourier
  Char_t name[256];
  Char_t title[256];
  snprintf(name, sizeof(name), "%s_Fourier_Im", fData->GetName());
  snprintf(title, sizeof(title), "%s_Fourier_Im", fData->GetTitle());

  UInt_t noOfFourierBins = 0;
  if (fNoOfBins % 2 == 0)
    noOfFourierBins = fNoOfBins/2;
  else
    noOfFourierBins = (fNoOfBins+1)/2;

  TH1F* imaginaryFourier = new TH1F(name, title, static_cast<Int_t>(noOfFourierBins), -fResolution/2.0, static_cast<Double_t>(noOfFourierBins-1)*fResolution+fResolution/2.0);
  if (imaginaryFourier == nullptr) {
    fValid = false;
    std::cerr << std::endl << "**SEVERE ERROR** couldn't allocate memory for the imaginary part of the Fourier transform." << std::endl;
    return nullptr;
  }

  // fill imaginaryFourier vector
  for (Int_t i=0; i<static_cast<Int_t>(noOfFourierBins); i++) {
    imaginaryFourier->SetBinContent(i+1, scale*fOut[i][1]);
    imaginaryFourier->SetBinError(i+1, 0.0);
  }

  return imaginaryFourier;
}

//--------------------------------------------------------------------------
// GetPowerFourier (public)
//--------------------------------------------------------------------------
/**
 * <p>Returns power spectrum |F(ω)| as ROOT histogram.
 *
 * <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)
{
  // check if valid flag is set
  if (!fValid)
    return nullptr;

  // invoke powerFourier
  Char_t name[256];
  Char_t title[256];
  snprintf(name, sizeof(name), "%s_Fourier_Pwr", fData->GetName());
  snprintf(title, sizeof(title), "%s_Fourier_Pwr", fData->GetTitle());

  UInt_t noOfFourierBins = 0;
  if (fNoOfBins % 2 == 0)
    noOfFourierBins = fNoOfBins/2;
  else
    noOfFourierBins = (fNoOfBins+1)/2;

  TH1F* pwrFourier = new TH1F(name, title, static_cast<Int_t>(noOfFourierBins), -fResolution/2.0, static_cast<Double_t>(noOfFourierBins-1)*fResolution+fResolution/2.0);
  if (pwrFourier == nullptr) {
    fValid = false;
    std::cerr << std::endl << "**SEVERE ERROR** couldn't allocate memory for the power part of the Fourier transform." << std::endl;
    return nullptr;
  }

  // fill powerFourier vector
  for (Int_t i=0; i<static_cast<Int_t>(noOfFourierBins); i++) {
    pwrFourier->SetBinContent(i+1, scale*sqrt(fOut[i][0]*fOut[i][0]+fOut[i][1]*fOut[i][1]));
    pwrFourier->SetBinError(i+1, 0.0);
  }

  return pwrFourier;
}

//--------------------------------------------------------------------------
// GetPhaseFourier (public)
//--------------------------------------------------------------------------
/**
 * <p>Returns phase spectrum φ(ω) = arg(F) as ROOT histogram.
 *
 * <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)
{
  // check if valid flag is set
  if (!fValid)
    return nullptr;

  // invoke phaseFourier
  Char_t name[256];
  Char_t title[256];
  snprintf(name, sizeof(name), "%s_Fourier_Phase", fData->GetName());
  snprintf(title, sizeof(title), "%s_Fourier_Phase", fData->GetTitle());

  UInt_t noOfFourierBins = 0;
  if (fNoOfBins % 2 == 0)
    noOfFourierBins = fNoOfBins/2;
  else
    noOfFourierBins = (fNoOfBins+1)/2;

  TH1F* phaseFourier = new TH1F(name, title, static_cast<Int_t>(noOfFourierBins), -fResolution/2.0, static_cast<Double_t>(noOfFourierBins-1)*fResolution+fResolution/2.0);
  if (phaseFourier == nullptr) {
    fValid = false;
    std::cerr << std::endl << "**SEVERE ERROR** couldn't allocate memory for the phase part of the Fourier transform." << std::endl;
    return nullptr;
  }

  // fill phaseFourier vector
  Double_t value = 0.0;
  for (Int_t i=0; i<static_cast<Int_t>(noOfFourierBins); i++) {
    // calculate the phase
    if (fOut[i][0] == 0.0) {
      if (fOut[i][1] >= 0.0)
        value = PI_HALF;
      else
        value = -PI_HALF;
    } else {
      value = atan(fOut[i][1]/fOut[i][0]);
      // check sector
      if (fOut[i][0] < 0.0) {
        if (fOut[i][1] > 0.0)
          value = PI + value;
        else
          value = PI - value;
      }
    }
    phaseFourier->SetBinContent(i+1, scale*value);
    phaseFourier->SetBinError(i+1, 0.0);
  }

  return phaseFourier;
}

//--------------------------------------------------------------------------
// PrepareFFTwInputData (private)
//--------------------------------------------------------------------------
/**
 * <p>Prepares time-domain data for FFT with optional DC correction and apodization.
 *
 * <p>This method performs essential preprocessing steps:
 * 1. Locates time-zero bin in the histogram
 * 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)
{
  // 1st find t==0. fData start at times t<0!!
  Int_t t0bin = -1;
  for (Int_t i=1; i<fData->GetNbinsX(); i++) {
    if (fData->GetBinCenter(i) >= 0.0) {
      t0bin = i;
      break;
    }
  }

  Int_t ival = static_cast<Int_t>(fStartTime/fTimeResolution) + t0bin;
  UInt_t start = 0;
  if (ival >= 0) {
    start = static_cast<UInt_t>(ival);
  }

  Double_t mean = 0.0;
  if (fDCCorrected) {
    for (UInt_t i=0; i<fNoOfData; i++) {
      mean += fData->GetBinContent(static_cast<Int_t>(i+start));
    }
    mean /= static_cast<Double_t>(fNoOfData);
  }

  // 2nd fill fIn
  for (UInt_t i=0; i<fNoOfData; i++) {
    fIn[i][0] = fData->GetBinContent(static_cast<Int_t>(i+start)) - mean;
    fIn[i][1] = 0.0;
  }
  for (UInt_t i=fNoOfData; i<fNoOfBins; i++) {
    fIn[i][0] = 0.0;
    fIn[i][1] = 0.0;
  }

  // 3rd apodize data (if wished)
  ApodizeData(static_cast<Int_t>(apodizationTag));
}

//--------------------------------------------------------------------------
// ApodizeData (private)
//--------------------------------------------------------------------------
/**
 * <p>Applies apodization (windowing) to time-domain data before FFT.
 *
 * <p><b>Purpose:</b> Apodization multiplies the time signal by a smooth
 * window function that tapers to zero at the edges, reducing the Gibbs
 * 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) {

  const Double_t cweak[3]   = { 0.384093, -0.087577, 0.703484 };
  const Double_t cmedium[3] = { 0.152442, -0.136176, 0.983734 };
  const Double_t cstrong[3] = { 0.045335,  0.554883, 0.399782 };

  Double_t c[5];
  for (UInt_t i=0; i<5; i++) {
    c[i] = 0.0;
  }

  switch (apodizationTag) {
    case F_APODIZATION_NONE:
      return;
    case F_APODIZATION_WEAK:
      c[0] = cweak[0]+cweak[1]+cweak[2];
      c[1] = -(cweak[1]+2.0*cweak[2]);
      c[2] = cweak[2];
      break;
    case F_APODIZATION_MEDIUM:
      c[0] = cmedium[0]+cmedium[1]+cmedium[2];
      c[1] = -(cmedium[1]+2.0*cmedium[2]);
      c[2] = cmedium[2];
      break;
    case F_APODIZATION_STRONG:
      c[0] = cstrong[0]+cstrong[1]+cstrong[2];
      c[1] = -2.0*(cstrong[1]+2.0*cstrong[2]);
      c[2] = cstrong[1]+6.0*cstrong[2];
      c[3] = -4.0*cstrong[2];
      c[4] = cstrong[2];
      break;
    default:
      std::cerr << std::endl << ">> **ERROR** User Apodization tag " << apodizationTag << " unknown, sorry ..." << std::endl;
      break;
  }

  Double_t q;
  for (UInt_t i=0; i<fNoOfData; i++) {
    q = c[0];
    for (UInt_t j=1; j<5; j++) {
      q += c[j] * pow(static_cast<Double_t>(i)/static_cast<Double_t>(fNoOfData), 2.0*static_cast<Double_t>(j));
    }
    fIn[i][0] *= q;
  }
}