///////////////////////////////////////////////////////////////////////
//
// EnergyRThetaFunctional
//
// Uses Segmentor and poisson correction to convert from NHits to an
// energy estimator (H).  H(r, theta) is converted to an effective H(r=0)
// using 3rd order polynomials of H vs r at a given theta (interpolation
// of the two closest H(r) values to find H(r, theta)).  The energy is
// then estimated using a lookup of H(0) vs energy.
//
// Author: Matt Mottram <m.mottram@qmul.ac.uk> -- contact person
//
// REVISION HISTORY:
// - 2015/11/26: Matt Mottram - first instance
// - 2016/05/19: Matt Mottram - revert to lookup H vs E
//
///////////////////////////////////////////////////////////////////////

#ifndef __RAT_Methods_EnergyRThetaFunctional__
#define __RAT_Methods_EnergyRThetaFunctional__

#include <RAT/SeededMethod.hh>

#include <vector>
#include <string>

namespace RAT{

namespace DS
{
  class FitResult;
}

namespace Methods{

  class EnergyRThetaFunctional : public SeededMethod {

  public:

    // Name getters
    std::string GetName() const { return EnergyRThetaFunctional::Name(); }
    static std::string Name() { return std::string("energyRThetaFunctional"); }

    // Initialise the method, index is optional index in the DB
    void Initialise(const std::string& index) { fIndex = index; }

    // BeginOfRun to load constants from the database
    void BeginOfRun(DS::Run& run);

    // Nothing to do in EndOfRun
    void EndOfRun(DS::Run&) {}

    // Set the seed to a default value (origin)
    void DefaultSeed();

    // Calculate and return best fit value
    DS::FitResult GetBestFit();

  private:
    std::string fIndex; // Optional index override for the DB table

    double fRPiecewise; // Radius cutoff to transition from f1 to f2
    double fRCutoff; // Upper radius beyond which fits are invalid

    std::vector<double> fEnergies; // Energies H per MeV is coordinated at
    std::vector<double> fHAtEnergy; // Segmentor H parameter (energy scaling) as a lookup

    std::vector<double> fThetas; // Theta values that index the other vector member variables

    std::vector<double> fF1Pol0; // 0th order polynomial constant for function1
    std::vector<double> fF1Pol1; // 0th order polynomial constant for function1
    std::vector<double> fF1Pol2; // 0th order polynomial constant for function1
    std::vector<double> fF1Pol3; // 0th order polynomial constant for function1

    std::vector<double> fF2Constant; // constant for polynomial x -> (x - [0]) correction
    std::vector<double> fF2Pol0; // 0th order polynomial constant for function2
    std::vector<double> fF2Pol1; // 0th order polynomial constant for function2
    std::vector<double> fF2Pol2; // 0th order polynomial constant for function2
    std::vector<double> fF2Pol3; // 0th order polynomial constant for function2

  };

} // Methods

} // RAT

#endif