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

  PRunNonMusr.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 <iostream>

#include "PRunNonMusr.h"

//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
 * \brief Default constructor creating an empty, invalid non-μSR run object.
 *
 * Initializes all member variables to default/safe values:
 * - fNoOfFitBins = 0 (no bins to fit)
 * - fPacking = 1 (no data averaging)
 * - fStartTimeBin = 0 (first data point)
 * - fEndTimeBin = 0 (no range)
 * - fHandleTag = kEmpty (uninitialized)
 * - fRawRunData = nullptr (no data loaded)
 *
 * This constructor is needed for creating vectors of PRunNonMusr objects.
 * The resulting object cannot be used until properly initialized via the
 * main constructor.
 *
 * \see PRunNonMusr(PMsrHandler*, PRunDataHandler*, UInt_t, EPMusrHandleTag, Bool_t)
 */
PRunNonMusr::PRunNonMusr() : PRunBase()
{
  fNoOfFitBins  = 0;
  fPacking = 1;
  fStartTimeBin = 0;
  fEndTimeBin = 0;

  fHandleTag = kEmpty;

  fRawRunData = nullptr;
}

//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
 * \brief Main constructor initializing a non-μSR run from MSR file and x-y data.
 *
 * Performs comprehensive initialization for general curve fitting:
 *
 * 1. <b>Base Class Initialization:</b>
 *    - Calls PRunBase constructor with MSR/data handlers
 *    - Initializes theory engine and parameter mappings
 *    - Sets up FUNCTIONS block evaluation
 *
 * 2. <b>Raw Data Loading:</b>
 *    - Retrieves raw x-y data using run name from MSR file
 *    - Calls fRawData->GetRunData(*runName)
 *    - Validates data was successfully loaded
 *    - If loading fails → marks run invalid, prints error
 *
 * 3. <b>Data Preparation:</b>
 *    - Calls PrepareData() to process x-y data
 *    - Extracts x-y columns based on MSR file specification
 *    - Applies packing to average data points
 *    - Sets up fit range boundaries
 *    - If preparation fails → marks run invalid
 *
 * The object is marked as invalid (fValid=false) if:
 * - Raw data file cannot be loaded (file not found, wrong format, etc.)
 * - PrepareData() fails (invalid column indices, missing packing, etc.)
 *
 * Key features for non-μSR:
 * - No histogram processing (data is already x-y pairs)
 * - No time-zero determination (not time-differential)
 * - No background subtraction (included in y-data)
 * - No asymmetry calculation (direct y vs. x fitting)
 *
 * \param msrInfo Pointer to MSR file handler (must remain valid for object lifetime)
 * \param rawData Pointer to raw data handler for loading data files
 * \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 x-values, false = high-resolution (minimal effect for non-μSR)
 *
 * \warning Always check IsValid() after construction before using for fitting
 *
 * \see PrepareData() for data processing details
 * \see PRunBase constructor for base class initialization
 */
PRunNonMusr::PRunNonMusr(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData) :
  PRunBase(msrInfo, rawData, runNo, tag), fTheoAsData(theoAsData)
{
  // get the proper run
  fRawRunData = fRawData->GetRunData(*(fRunInfo->GetRunName()));
  if (!fRawRunData) { // couldn't get run
    std::cerr << std::endl << "PRunNonMusr::PRunNonMusr(): **ERROR** Couldn't get raw run  data!";
    std::cerr << std::endl;
    fValid = false;
  }

  // calculate fData
  if (!PrepareData())
    fValid = false;
}

//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
 * \brief Destructor (no cleanup needed for non-μSR).
 *
 * The fRawRunData pointer is not owned by this class and is not deleted here.
 * It is managed by the PRunDataHandler and will be cleaned up externally.
 *
 * Base class destructor (PRunBase) handles cleanup of:
 * - Theory engine objects
 * - Parameter mapping structures
 * - Function value caches
 */
PRunNonMusr::~PRunNonMusr()
{
}

//--------------------------------------------------------------------------
// CalcChiSquare
//--------------------------------------------------------------------------
/**
 * \brief Calculates χ² between non-μSR x-y data and theory.
 *
 * Computes the chi-squared statistic for general x-y data:
 * \f[ \chi^2 = \sum_{i={\rm start}}^{\rm end} \frac{(y_i - f(x_i))^2}{\sigma_i^2} \f]
 *
 * where:
 * - x_i is the independent variable (arbitrary units)
 * - y_i is the measured dependent variable
 * - f(x_i) is the theory function evaluated at x_i
 * - σ_i is the error on y_i (from data file)
 *
 * Algorithm:
 * 1. <b>Evaluate FUNCTIONS block:</b>
 *    - Loop over user-defined functions
 *    - Compute function values using current parameters
 *    - Store in fFuncValues for use by theory
 *
 * 2. <b>Calculate χ² sum:</b>
 *    - Loop from fStartTimeBin to fEndTimeBin (inclusive)
 *    - For each data point i:
 *      a. Get x-value: x = fData.GetX()->at(i)
 *      b. Evaluate theory at x: theo = fTheory->Func(x, par, fFuncValues)
 *      c. Get data and error: y = fData.GetValue()->at(i), σ = fData.GetError()->at(i)
 *      d. Compute contribution: Δχ² = (y - theo)² / σ²
 *    - Sum all contributions
 *
 * Key differences from μSR fits:
 * - x-axis is arbitrary (not necessarily time in μs)
 * - Theory evaluated on-demand (not pre-calculated on grid)
 * - Loop end is INCLUSIVE (i <= fEndTimeBin)
 * - No OpenMP parallelization (simpler, smaller data sets)
 *
 * \param par Parameter vector from MINUIT with current parameter values
 * \return Chi-squared value (minimize during fitting)
 *
 * \see CalcMaxLikelihood() - NOT IMPLEMENTED for non-μSR
 * \see PrepareData() for x-y data loading and fit range setup
 */
Double_t PRunNonMusr::CalcChiSquare(const std::vector<Double_t>& par)
{
  Double_t chisq = 0.0;
  Double_t diff = 0.0;

  // calculate functions
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
  }

  // calculate chi square
  Double_t x(1.0);
  for (UInt_t i=fStartTimeBin; i<=fEndTimeBin; i++) {
    x = fData.GetX()->at(i);
    diff = fData.GetValue()->at(i) - fTheory->Func(x, par, fFuncValues);
    chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
  }

  return chisq;
}

//--------------------------------------------------------------------------
// CalcChiSquareExpected (public)
//--------------------------------------------------------------------------
/**
 * \brief Calculates expected χ² (NOT IMPLEMENTED for non-μSR).
 *
 * This method is not implemented because the concept of "expected χ²" requires
 * knowledge of the underlying statistical distribution of the data, which is
 * not well-defined for general x-y data.
 *
 * For μSR histogram data, the expected χ² can be calculated based on Poisson
 * statistics. For arbitrary non-μSR data, the appropriate statistical model
 * depends on the data source:
 * - Counting experiments → Poisson
 * - Averaged measurements → Gaussian
 * - Other experiments → problem-specific distributions
 *
 * Without knowing the data's statistical nature, a meaningful implementation
 * cannot be provided.
 *
 * \param par Parameter vector from MINUIT (unused)
 * \return Always returns 0.0
 *
 * \note Prints "not implemented yet" message to stdout
 */
Double_t PRunNonMusr::CalcChiSquareExpected(const std::vector<Double_t>& par)
{
  std::cout << std::endl << "PRunNonMusr::CalcChiSquareExpected(): not implemented yet ..." << std::endl;

  return 0.0;
}

//--------------------------------------------------------------------------
// CalcMaxLikelihood
//--------------------------------------------------------------------------
/**
 * \brief Calculates maximum likelihood (NOT IMPLEMENTED for non-μSR).
 *
 * Maximum likelihood fitting is not implemented for general x-y data because:
 *
 * 1. <b>Distribution-dependent:</b> The likelihood function depends on the
 *    underlying statistical distribution of the data points, which varies:
 *    - Poisson: For counting experiments (like μSR histograms)
 *    - Gaussian: For averaged measurements with known errors
 *    - Other: Problem-specific distributions (binomial, exponential, etc.)
 *
 * 2. <b>Gaussian case equivalence:</b> If errors are Gaussian, maximum likelihood
 *    is mathematically equivalent to χ² minimization (already implemented).
 *
 * 3. <b>Implementation complexity:</b> Supporting arbitrary likelihood functions
 *    would require users to specify the distribution, adding complexity without
 *    clear benefit for most non-μSR applications.
 *
 * For general x-y data with Gaussian errors, use χ² fitting via CalcChiSquare().
 * For non-Gaussian statistics, users should implement custom fitting outside musrfit.
 *
 * \param par Parameter vector from MINUIT (unused)
 * \return Always returns 1.0
 *
 * \note Prints "not implemented yet" message to stdout
 *
 * \see CalcChiSquare() for standard least-squares fitting (recommended for non-μSR)
 */
Double_t PRunNonMusr::CalcMaxLikelihood(const std::vector<Double_t>& par)
{
  std::cout << std::endl << "PRunNonMusr::CalcMaxLikelihood(): not implemented yet ..." << std::endl;

  return 1.0;
}

//--------------------------------------------------------------------------
// CalcTheory
//--------------------------------------------------------------------------
/**
 * \brief Evaluates theory function (empty implementation for non-μSR).
 *
 * For non-μSR data, theory calculation is performed on-demand within
 * CalcChiSquare() rather than pre-calculating and storing theory values.
 *
 * This design choice is made because:
 * - Theory is evaluated at arbitrary x-values from the data file
 * - No need for high-resolution theory grid (data defines x-points)
 * - Simpler and more efficient to evaluate during χ² loop
 * - Avoids storing redundant theory array
 *
 * This empty method exists solely to satisfy the PRunBase interface requirement.
 *
 * \see CalcChiSquare() for on-demand theory evaluation at each data point
 */
void PRunNonMusr::CalcTheory()
{
}

//--------------------------------------------------------------------------
// GetNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
 * \brief Calculates and returns the number of x-y points within the fit range.
 *
 * Counts data points where fFitStartTime ≤ x ≤ fFitEndTime (both inclusive).
 * This count is essential for:
 * - Degrees of freedom: ν = N_points - N_params
 * - Reduced χ²: χ²_red = χ² / ν
 * - Statistical quality assessment: χ²/ν ≈ 1 indicates good fit
 *
 * Algorithm:
 * 1. Reset counter: fNoOfFitBins = 0
 * 2. Loop through all x-values in fData.GetX()
 * 3. For each x-value:
 *    - If fFitStartTime ≤ x ≤ fFitEndTime → increment counter
 * 4. Return final count
 *
 * The fit range (fFitStartTime, fFitEndTime) is specified in the MSR file
 * RUN or GLOBAL block "fit" entry, in the same units as the x-axis data.
 *
 * \return Number of data points within fit range
 *
 * \note This method recalculates the count each time it's called (not cached),
 *       allowing for dynamic fit range changes.
 */
UInt_t PRunNonMusr::GetNoOfFitBins()
{
  fNoOfFitBins=0;
  Double_t x;
  for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
    x = fData.GetX()->at(i);
    if ((x >= fFitStartTime) && (x <= fFitEndTime))
      fNoOfFitBins++;
  }

  return fNoOfFitBins;
}

//--------------------------------------------------------------------------
// PrepareData
//--------------------------------------------------------------------------
/**
 * <p>Prepare data for fitting or viewing.
 *
 * <b>return:</b>
 * - true if everthing went smooth
 * - false, otherwise.
 */
Bool_t PRunNonMusr::PrepareData()
{
  Bool_t success = true;

  if (!fValid)
    return false;

  if (fRunInfo->GetRunNameSize() > 1) { // ADDRUN present which is not supported for NonMusr
    std::cerr << std::endl << ">> PRunNonMusr::PrepareData(): **WARNING** ADDRUN NOT SUPPORTED FOR THIS FIT TYPE, WILL IGNORE IT." << std::endl;
  }

  // get packing info
  fPacking = fRunInfo->GetPacking();
  if (fPacking == -1) { // packing not present in the RUN block, will try the GLOBAL block
    fPacking = fMsrInfo->GetMsrGlobal()->GetPacking();
  }
  if (fPacking == -1) { // packing NOT present, in neither the RUN block, nor in the GLOBAL block
    std::cerr << std::endl << ">> PRunNonMusr::PrepareData(): **ERROR** couldn't find any packing information." << std::endl;
    return false;
  }

  // get fit start/end time
  fFitStartTime = fRunInfo->GetFitRange(0);
  fFitEndTime   = fRunInfo->GetFitRange(1);
  if (fFitStartTime == PMUSR_UNDEFINED) { // not present in the RUN block, will try GLOBAL block
    fFitStartTime = fMsrInfo->GetMsrGlobal()->GetFitRange(0);
    fFitEndTime   = fMsrInfo->GetMsrGlobal()->GetFitRange(1);
  }

  if (fHandleTag == kFit)
    success = PrepareFitData();
  else if (fHandleTag == kView)
    success = PrepareViewData();
  else
    success = false;

  return success;
}

//--------------------------------------------------------------------------
// PrepareFitData
//--------------------------------------------------------------------------
/**
 * <p>Prepare data for fitting.
 *
 * <b>return:</b>
 * - true if everthing went smooth
 * - false, otherwise.
 */
Bool_t PRunNonMusr::PrepareFitData()
{
  Bool_t success = true;

  // get x-, y-index
  UInt_t xIndex = GetXIndex();
  UInt_t yIndex = GetYIndex();

  // pack the raw data
  Double_t value  = 0.0;
  Double_t err = 0.0;
  for (UInt_t i=0; i<fRawRunData->fDataNonMusr.GetData()->at(xIndex).size(); i++) {
    if (fPacking == 1) {
      fData.AppendXValue(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i));
      fData.AppendValue(fRawRunData->fDataNonMusr.GetData()->at(yIndex).at(i));
      fData.AppendErrorValue(fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i));
    } else { // packed data, i.e. fPacking > 1
      if ((i % fPacking == 0) && (i != 0)) { // fill data
        fData.AppendXValue(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i)-(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i)-fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i-fPacking))/2.0);
        fData.AppendValue(value);
        fData.AppendErrorValue(TMath::Sqrt(err));
        value = 0.0;
        err = 0.0;
      }
      // sum raw data values
      value += fRawRunData->fDataNonMusr.GetData()->at(yIndex).at(i);
      err += fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i)*fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i);
    }
  }

  // count the number of bins to be fitted
  fNoOfFitBins=0;
  Double_t x;
  for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
    x = fData.GetX()->at(i);
    if ((x >= fFitStartTime) && (x <= fFitEndTime))
      fNoOfFitBins++;
  }

  // get start/end bin
  const PDoubleVector *xx = fData.GetX();
  fStartTimeBin = 0;
  fEndTimeBin = xx->size()-1;
  for (UInt_t i=0; i<xx->size(); i++) {
    if (xx->at(i) < fFitStartTime)
      fStartTimeBin = i;
    if (xx->at(i) < fFitEndTime)
      fEndTimeBin = i;
  }

  return success;
}

//--------------------------------------------------------------------------
// PrepareViewData
//--------------------------------------------------------------------------
/**
 * <p>Prepare data for viewing.
 *
 * <b>return:</b>
 * - true if everthing went smooth
 * - false, otherwise.
 */
Bool_t PRunNonMusr::PrepareViewData()
{
  Bool_t success = true;

  // get x-, y-index
  UInt_t xIndex = GetXIndex();
  UInt_t yIndex = GetYIndex();

  // fill data histo
  // pack the raw data
  Double_t value  = 0.0;
  Double_t err = 0.0;
  for (UInt_t i=0; i<fRawRunData->fDataNonMusr.GetData()->at(xIndex).size(); i++) {
    if (fPacking == 1) {
      fData.AppendXValue(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i));
      fData.AppendValue(fRawRunData->fDataNonMusr.GetData()->at(yIndex).at(i));
      fData.AppendErrorValue(fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i));
    } else { // packed data, i.e. fPacking > 1
      if ((i % fPacking == 0) && (i != 0)) { // fill data
        fData.AppendXValue(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i)-(fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i)-fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i-fPacking))/2.0);
        fData.AppendValue(value);
        fData.AppendErrorValue(TMath::Sqrt(err));
        value = 0.0;
        err = 0.0;
      }
      // sum raw data values
      value += fRawRunData->fDataNonMusr.GetData()->at(yIndex).at(i);
      err += fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i)*fRawRunData->fDataNonMusr.GetErrData()->at(yIndex).at(i);
    }
  }

  // count the number of bins to be fitted
  fNoOfFitBins = fData.GetValue()->size();

  // fill theory histo
  // feed the parameter vector
  std::vector<Double_t> par;
  PMsrParamList *paramList = fMsrInfo->GetMsrParamList();
  for (UInt_t i=0; i<paramList->size(); i++)
    par.push_back((*paramList)[i].fValue);
  // calculate functions
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
  }

  // get plot range
  PMsrPlotList *plotList;
  PMsrPlotStructure plotBlock;
  plotList = fMsrInfo->GetMsrPlotList();
  // find the proper plot block
  // Here a small complication to be handled: there are potentially multiple
  // run blocks and the run might be present in various of these run blocks. In
  // order to get a nice resolution on the theory the following procedure will be
  // followed: the smallest x-interval found will be used to for the fXTheory resolution
  // which is 1000 function points. The function will be calculated from the smallest
  // xmin found up to the largest xmax found.
  Double_t xMin = 0.0, xMax = 0.0;

  // init xMin/xMax, xAbsMin/xAbsMax
  plotBlock = plotList->at(0);
  if (plotBlock.fTmin.size() == 0) { // check if no range information is present
    PMsrRunList *runList = fMsrInfo->GetMsrRunList();
    xMin = runList->at(0).GetFitRange(0);
    xMax = runList->at(0).GetFitRange(1);
    for (UInt_t i=1; i<runList->size(); i++) {
      if (runList->at(i).GetFitRange(0) < xMin)
        xMin = runList->at(i).GetFitRange(0);
      if (runList->at(i).GetFitRange(1) > xMax)
        xMax = runList->at(i).GetFitRange(1);
    }
  } else if (plotBlock.fTmin.size() == 1) { // check if 'range' information is present
    xMin = plotBlock.fTmin[0];
    xMax = plotBlock.fTmax[0];
  } else if (plotBlock.fTmin.size() > 1) { // check if 'sub_ranges' information is present
    xMin = plotBlock.fTmin[0];
    xMax = plotBlock.fTmax[0];
    for (UInt_t i=1; i<plotBlock.fTmin.size(); i++) {
      if (plotBlock.fTmin[i] < xMin)
        xMin = plotBlock.fTmin[i];
      if (plotBlock.fTmax[i] > xMax)
        xMax = plotBlock.fTmax[i];
    }
  }

  if (plotBlock.fUseFitRanges) { // check if 'use_fit_ranges' information is present
    PMsrRunList *runList = fMsrInfo->GetMsrRunList();
    xMin = runList->at(0).GetFitRange(0);
    xMax = runList->at(0).GetFitRange(1);
    for (UInt_t i=1; i<runList->size(); i++) {
      if (runList->at(i).GetFitRange(0) < xMin)
        xMin = runList->at(i).GetFitRange(0);
      if (runList->at(i).GetFitRange(1) > xMax)
        xMax = runList->at(i).GetFitRange(1);
    }
  }

  for (UInt_t i=1; i<plotList->size(); i++) { // go through all the plot blocks
    plotBlock = plotList->at(i);

    if (plotBlock.fTmin.size() == 0) { // check if no range information is present
      PMsrRunList *runList = fMsrInfo->GetMsrRunList();
      for (UInt_t i=0; i<runList->size(); i++) {
        if (runList->at(i).GetFitRange(0) < xMin)
          xMin = runList->at(i).GetFitRange(0);
        if (runList->at(i).GetFitRange(1) > xMax)
          xMax = runList->at(i).GetFitRange(1);
      }
    } else if (plotBlock.fTmin.size() == 1) { // check if 'range' information is present
      if (plotBlock.fTmin[0] < xMin)
        xMin = plotBlock.fTmin[0];
      if (plotBlock.fTmax[0] > xMax)
        xMax = plotBlock.fTmax[0];
    } else if (plotBlock.fTmin.size() > 1) { // check if 'sub_ranges' information is present
      for (UInt_t i=0; i<plotBlock.fTmin.size(); i++) {
        if (plotBlock.fTmin[i] < xMin)
          xMin = plotBlock.fTmin[i];
        if (plotBlock.fTmax[i] > xMax)
          xMax = plotBlock.fTmax[i];
      }
    }

    if (plotBlock.fUseFitRanges) { // check if 'use_fit_ranges' information is present
      PMsrRunList *runList = fMsrInfo->GetMsrRunList();
      for (UInt_t i=0; i<runList->size(); i++) {
        if (runList->at(i).GetFitRange(0) < xMin)
          xMin = runList->at(i).GetFitRange(0);
        if (runList->at(i).GetFitRange(1) > xMax)
          xMax = runList->at(i).GetFitRange(1);
      }
    }
  }

  // typically take 1000 points to calculate the theory, except if there are more data points, than take that number
  Double_t xStep;
  if (fData.GetX()->size() > 1000.0)
    xStep = (xMax-xMin)/fData.GetX()->size();
  else
    xStep = (xMax-xMin)/1000.0;

  if (fTheoAsData) {
    Double_t xx;
    for (UInt_t i=0; i<fRawRunData->fDataNonMusr.GetData()->at(xIndex).size(); i++) {
      // fill x-vector
      xx = fRawRunData->fDataNonMusr.GetData()->at(xIndex).at(i);
      fData.AppendXTheoryValue(xx);
      // fill y-vector
      fData.AppendTheoryValue(fTheory->Func(xx, par, fFuncValues));
    }
  } else {
    Double_t xx = xMin;
    do {
      // fill x-vector
      fData.AppendXTheoryValue(xx);
      // fill y-vector
      fData.AppendTheoryValue(fTheory->Func(xx, par, fFuncValues));
      // calculate next xx
      xx += xStep;
    } while (xx < xMax);
  }

  // clean up
  par.clear();

  return success;
}

//--------------------------------------------------------------------------
// GetXIndex
//--------------------------------------------------------------------------
/**
 * <p>Returns the x-axis data index.
 *
 * <b>return:</b>
 * - x-index
 */
UInt_t PRunNonMusr::GetXIndex()
{
  UInt_t index = 0;
  Bool_t found = false;

  if (fRawRunData->fDataNonMusr.FromAscii()) { // ascii-file format
    index = 0;
    found = true;
  } else { // db-file format
    if (fRunInfo->GetXDataIndex() > 0) { // xy-data already indices
      index = fRunInfo->GetXDataIndex()-1; // since xy-data start with 1 ...
      found = true;
    } else { // xy-data data tags which needs to be converted to an index
      for (UInt_t i=0; i<fRawRunData->fDataNonMusr.GetDataTags()->size(); i++) {
        if (fRawRunData->fDataNonMusr.GetDataTags()->at(i).CompareTo(*fRunInfo->GetXDataLabel()) == 0) {
          index = i;
          found = true;
          break;
        }
      }
    }
  }

  if (!found) {
    std::cerr << std::endl << "PRunNonMusr::GetXIndex(): **ERROR** Couldn't obtain x-data index!";
    std::cerr << std::endl;
    assert(0);
  }

  return index;
}

//--------------------------------------------------------------------------
// GetYIndex
//--------------------------------------------------------------------------
/**
 * <p>Returns the y-axis data index.
 *
 * <b>return:</b>
 * - y-index
 */
UInt_t PRunNonMusr::GetYIndex()
{
  UInt_t index = 0;
  Bool_t found = false;

  if (fRawRunData->fDataNonMusr.FromAscii()) { // ascii-file format
    index = 1;
    found = true;
  } else { // db-file format
    if (fRunInfo->GetYDataIndex() > 0) { // xy-data already indices
      index = fRunInfo->GetYDataIndex()-1; // since xy-data start with 1 ...
      found = true;
    } else { // xy-data data tags which needs to be converted to an index
      for (UInt_t i=0; i<fRawRunData->fDataNonMusr.GetDataTags()->size(); i++) {
        if (fRawRunData->fDataNonMusr.GetDataTags()->at(i).CompareTo(*fRunInfo->GetYDataLabel()) == 0) {
          index = i;
          found = true;
          break;
        }
      }
    }
  }

  if (!found) {
    std::cerr << std::endl << "PRunNonMusr::GetYIndex(): **ERROR** Couldn't obtain y-data index!";
    std::cerr << std::endl;
    assert(0);
  }

  return index;
}