#ifndef MARS_MExtralgoSpline
#define MARS_MExtralgoSpline

#ifndef ROOT_TROOT
#include <TROOT.h>
#endif

#include <iostream>
class TComplex;

class MExtralgoSpline
{
public:  
    enum ExtractionType_t { kAmplitude, kIntegral };    //! Possible time and charge extraction types

private:
    ExtractionType_t fExtractionType;

private:
    //Bool_t fIsOwner; // Owner of derivatives....

    // Input
    Float_t *fVal;
    Int_t    fNum;

    Float_t *fDer1;
    Float_t *fDer2;

    Float_t fRiseTime;
    Float_t fFallTime;

    Float_t fResolution;

    // Result
    Float_t fTime;
    Float_t fTimeDev;
    Float_t fSignal;
    Float_t fSignalDev;

    Double_t ReMul(const TComplex &c1, const TComplex &th) const;

    inline Float_t Eval(Float_t val, Float_t a, Float_t deriv) const
    {
        return a*val + (a*a*a-a)*deriv;
    }

    // Evaluate value of spline in the interval i with x=[0;1[
    inline Float_t Eval(const Int_t i, const Float_t x) const
    {
        // Eval(i,x) =  (fDer2[i+1]-fDer2[i])*x*x*x + 3*fDer2[i]*x*x +
        //              (fVal[i+1]-fVal[i] -2*fDer2[i]-fDer2[i+1])*x + fVal[i];

        // x := [0; 1[
        return Eval(fVal[i], 1-x, fDer2[i]) + Eval(fVal[i+1], x, fDer2[i+1]);
    }

    /*
    inline Float_t EvalAt(const Float_t x) const
    {
        Int_t i = TMath::FloorNint(x);

        // handle under- and overflow of the array-range by extrapolation
        if (i<0)
            i=0;
        if (i>fNum-2)
            i = fNum-2;

        return Eval(i, x-i);
    }
    */

    // Evaluate first derivative of spline in the interval i with x=[0;1[
    inline Double_t EvalDeriv1(const Float_t x, const Int_t i) const
    {
        // x := [0; 1[
        const Double_t difval = fVal[i+1]-fVal[i];
        const Double_t difder = fDer2[i+1]-fDer2[i];

        return 3*difder*x*x + 6*fDer2[i]*x - 2*fDer2[i] - fDer2[i+1] + difval;
    }

    // Evaluate second derivative of spline in the interval i with x=[0;1[
    inline Double_t EvalDeriv2(const Float_t x, const Int_t i) const
    {
        // x := [0; 1[
        return 6*(fDer2[i+1]*x + fDer2[i]*(1-x));
    }

    Double_t FindY(Int_t i, Double_t y=0, Double_t min=0, Double_t max=1) const;
    Double_t SearchY(Float_t maxpos, Float_t y) const;
/*
    // Evaluate first solution for a possible maximum (x|first deriv==0)
    inline Double_t EvalDerivEq0S1(const Int_t i) const
    {
        // return the x value [0;1[ at which the derivative is zero (solution1)

        Double_t sumder = fDer2[i]+fDer2[i+1];
        Double_t difder = fDer2[i]-fDer2[i+1];

        Double_t sqt1 = sumder*sumder - fDer2[i]*fDer2[i+1];
        Double_t sqt2 = difder*(fVal[i+1]-fVal[i]);

        Double_t x = 3*fDer2[i] - sqrt(3*sqt1 + 3*sqt2);

        Double_t denom = 3*(fDer2[i+1]-fDer2[i]);

        return -x/denom;
    }

    // Evaluate second solution for a possible maximum (x|first deriv==0)
    inline Double_t EvalDerivEq0S2(const Int_t i) const
    {
        // return the x value [0;1[ at which the derivative is zero (solution2)

        Double_t sumder = fDer2[i]+fDer2[i+1];
        Double_t difder = fDer2[i]-fDer2[i+1];

        Double_t sqt1 = sumder*sumder - fDer2[i]*fDer2[i+1];
        Double_t sqt2 = difder*(fVal[i+1]-fVal[i]);

        Double_t x = 3*fDer2[i] + sqrt(3*sqt1 + 3*sqt2);

        Double_t denom = 3*(fDer2[i+1]-fDer2[i]);

        return -x/denom;
    }
    */

    inline void EvalDerivEq0(const Int_t i, Float_t &rc1, Float_t &rc2) const
    {
        Double_t sumder = fDer2[i]+fDer2[i+1];
        Double_t difder = fDer2[i]-fDer2[i+1];

        Double_t sqt1  = sumder*sumder - fDer2[i]*fDer2[i+1];
        Double_t sqt2  = difder*(fVal[i+1]-fVal[i]);
        Double_t sqt3  = sqrt(3*sqt1 + 3*sqt2);
        Double_t denom = 3*(fDer2[i+1]-fDer2[i]);

        rc1 = -(3*fDer2[i] + sqt3)/denom;
        rc2 = -(3*fDer2[i] - sqt3)/denom;
    }

    // Calculate the "Stammfunktion" of the Eval-function
    inline Double_t EvalPrimitive(Int_t i, Float_t x) const
    {
        /* TO BE CHECKED!
        if (x==0)
            return 0;

        if (x==1)
            return (fVal[i+1]+fVal[i])/2 - fDer2[i+1]/4;
            */
        Align(i, x);

        const Double_t x2  = x*x;
        const Double_t x4  = x2*x2;
        const Double_t x1  = 1-x;
        const Double_t x14 = x1*x1*x1*x1;

        return x2*fVal[i+1]/2 + (x4/2-x2)*fDer2[i+1]/2 + (x-x2/2)*fVal[i] + (x2/2-x-x14/4)*fDer2[i];
    }

    inline void Align(Int_t &i, Float_t &x) const
    {
        if (i<0)
        {
            x += i;
            i=0;
        }
        if (i>=fNum-1)
        {
            x += i-(fNum-2);
            i=fNum-2;
        }
    }

    // Calculate the intgeral of the Eval-function in
    // bin i from a=[0;1[ to b=[0;1[
    inline Double_t EvalInteg(Int_t i, Float_t a=0, Float_t b=1) const
    {
        // This is to make sure that we never access invalid
        // memory, even if this should never happen.
        // If it happens anyhow we extraolate the spline
        Align(i, a);
        Align(i, b);

        return EvalPrimitive(i, b)-EvalPrimitive(i, a);
    }

    // Calculate the intgeral of the Eval-function betwen x0 and x1
    inline Double_t EvalInteg(Float_t x0, Float_t x1) const
    {
        const Int_t min = TMath::CeilNint(x0);
        const Int_t max = TMath::FloorNint(x1);

        // This happens if x0 and x1 are in the same interval
        if (min>max)
            return EvalInteg(max, x0-max, x1-max);

        // Sum complete intervals
        Double_t sum = 0;
        for (int i=min; i<max; i++)
            sum += EvalInteg(i);

        // Sum the incomplete intervals at the beginning and end
        sum += EvalInteg(min-1, 1-(min-x0), 1);
        sum += EvalInteg(max,   0, x1-max);

        // return result
        return sum;
    }

    // We search for the maximum from x=i-1 to x=i+1
    // (Remeber: i corresponds to the value in bin i, i+1 to the
    //  next bin and i-1 to the last bin)
    inline void GetMaxAroundI(Int_t i, Float_t &xmax, Float_t &ymax) const
    {
        Float_t xmax1, xmax2;
        Float_t ymax1, ymax2;

        Bool_t rc1 = i>0      && GetMax(i-1, xmax1, ymax1);
        Bool_t rc2 = i<fNum-1 && GetMax(i,   xmax2, ymax2);

        // In case the medium bin is the first or last bin
        // take the lower or upper edge of the region into account.
        if (i==0)
        {
            xmax1 = 0;
            ymax1 = fVal[0];
            rc1 = kTRUE;
        }
        if (i>fNum-2)
        {
            xmax2 = fNum-1;
            ymax2 = fVal[fNum-1];
            rc2 = kTRUE;
        }

        // Take a default in case no maximum is found
        xmax=i;
        ymax=fVal[i];

        if (rc1)
        {
            ymax = ymax1;
            xmax = xmax1;
        }
        else
            if (rc2)
            {
                ymax = ymax2;
                xmax = xmax2;
            }

        if (rc2 && ymax2>ymax)
        {
            ymax = ymax2;
            xmax = xmax2;
        }
  /*
        // Search real maximum in [i-0.5, i+1.5]
        Float_t xmax1, xmax2, xmax3;
        Float_t ymax1, ymax2, ymax3;

        Bool_t rc1 = i>0 && GetMax(i-1, xmax1, ymax1, 0.5, 1.0);
        Bool_t rc2 = GetMax(i,   xmax2, ymax2, 0.0, 1.0);
        Bool_t rc3 = i<fNum-1 && GetMax(i+1, xmax3, ymax3, 0.0, 0.5);

        // In case the medium bin is the first or last bin
        // take the lower or upper edge of the region into account.
        if (i==0)
        {
            xmax1 = 0;
            ymax1 = Eval(0, 0);
            rc1 = kTRUE;
        }
        if (i==fNum-1)
        {
            xmax3 = fNum-1e-5;
            ymax3 = Eval(fNum-1, 1);
            rc3 = kTRUE;
        }

        // Take a real default in case no maximum is found
        xmax=i+0.5;
        ymax=Eval(i, 0.5);

        //if (!rc1 && !rc2 && !rc3)
        //    cout << "!!!!!!!!!!!!!!!" << endl;

        if (rc1)
        {
            ymax = ymax1;
            xmax = xmax1;
        }
        else
            if (rc2)
            {
                ymax = ymax2;
                xmax = xmax2;
            }
            else
                if (rc3)
                {
                    ymax = ymax3;
                    xmax = xmax3;
                }

        if (rc2 && ymax2>ymax)
        {
            ymax = ymax2;
            xmax = xmax2;
        }
        if (rc3 && ymax3>ymax)
        {
            ymax = ymax3;
            xmax = xmax3;
        }
*/    }

    inline Bool_t GetMax(Int_t i, Float_t &xmax, Float_t &ymax, Float_t min=0, Float_t max=1) const
    {
        // Find analytical maximum in the bin i in the interval [min,max[

        Float_t x1, x2;
        EvalDerivEq0(i, x1, x2);
        // const Float_t x1 = EvalDerivEq0S1(i);
        // const Float_t x2 = EvalDerivEq0S2(i);

        const Bool_t ismax1 = x1>=min && x1<max && EvalDeriv2(x1, i)<0;
        const Bool_t ismax2 = x2>=min && x2<max && EvalDeriv2(x2, i)<0;

        if (!ismax1 && !ismax2)
            return kFALSE;

        if (ismax1 && !ismax2)
        {
            xmax = i+x1;
            ymax = Eval(i, x1);
            return kTRUE;
        }

        if (!ismax1 && ismax2)
        {
            xmax = i+x2;
            ymax = Eval(i, x2);
            return kTRUE;
        }

        // Somehting must be wrong...
        return kFALSE;
        /*
        std::cout << "?????????????" << std::endl;

        const Double_t y1 = Eval(i, x1);
        const Double_t y2 = Eval(i, x2);

        if (y1>y2)
        {
            xmax = i+x1;
            ymax = Eval(i, x1);
            return kTRUE;
        }
        else
        {
            xmax = i+x2;
            ymax = Eval(i, x2);
            return kTRUE;
        }

        return kFALSE;*/
    }
/*
    inline Int_t GetMaxPos(Int_t i, Float_t &xmax, Float_t &ymax) const
    {
        Double_t x[3];

        x[0] = 0;
        // x[1] = 1; // This means we miss a possible maximum at the
                            // upper edge of the last interval...

        x[1] = EvalDerivEq0S1(i);
        x[2] = EvalDerivEq0S2(i);

        //y[0] = Eval(i, x[0]);
        //y[1] = Eval(i, x[1]);
        //y[1] = Eval(i, x[1]);
        //y[2] = Eval(i, x[2]);

        Int_t rc = 0;
        Double_t max = Eval(i, x[0]);

        for (Int_t j=1; j<3; j++)
        {
            if (x[j]<=0 || x[j]>=1)
                continue;

            const Float_t y = Eval(i, x[j]);
            if (y>max)
            {
                max = y;
                rc = j;
            }
        }

        if (max>ymax)
        {
            xmax = x[rc]+i;
            ymax = max;
        }

        return rc;
    }

    inline void GetMaxPos(Int_t min, Int_t max, Float_t &xmax, Float_t &ymax) const
    {
        Float_t xmax=-1;
        Float_t ymax=-FLT_MAX;

        for (int i=min; i<max; i++)
            GetMaxPos(i, xmax, ymax);

        for (int i=min+1; i<max; i++)
        {
            Float_t y = Eval(i, 0);
            if (y>ymax)
            {
                ymax = y;
                xmax = i;
            }
        }

    }*/


    void InitDerivatives() const;
    Float_t CalcIntegral(Float_t start) const;

public:
    MExtralgoSpline(Float_t *val, Int_t n, Float_t *der1, Float_t *der2)
        : fExtractionType(kIntegral), fVal(val), fNum(n), fDer1(der1), fDer2(der2), fTime(0), fTimeDev(-1), fSignal(0), fSignalDev(-1)
    {
        InitDerivatives();
    }

    void SetRiseFallTime(Float_t rise, Float_t fall) { fRiseTime=rise; fFallTime=fall; }
    void SetExtrationType(ExtractionType_t typ)      { fExtractionType = typ; }
    void SetResolution(Float_t res)                  { fResolution=res; }

    Float_t GetTime() const      { return fTime; }
    Float_t GetSignal() const    { return fSignal; }

    Float_t GetTimeDev() const   { return fTimeDev; }
    Float_t GetSignalDev() const { return fSignalDev; }

    void GetSignal(Float_t &sig, Float_t &dsig) const { sig=fSignal; dsig=fSignalDev; }
    void GetTime(Float_t &sig, Float_t &dsig) const   { sig=fTime; dsig=fTimeDev; }

    Float_t ExtractNoise(Int_t iter);
    void Extract(Byte_t sat, Int_t maxpos);
};

#endif