#include <iostream>
#include <fstream>
#include <stdlib.h>

#include "pulse.h"

using namespace std;
using namespace TMath;

Pulse::Pulse()
{
    InitMembers();
}

Pulse::Pulse(TString name)
{
    InitMembers();
    mName           = name;
}

Pulse::Pulse(TString name, TH1* histo)
{
    InitMembers();
    mName           = name;
    mHisto          = histo;
    Fit(mName, mOptions, mType);
}

Pulse::Pulse(TString name, TH1* histo, int type)
{
    InitMembers();
    mName           = name;
    mHisto          = histo;
    mType           = type;
    Fit(mName, mOptions, mType);
}

Pulse::Pulse(TString name, TH1* histo, TString options)
{
    InitMembers();
    mName           = name;
    mHisto          = histo;
    mOptions        += options;
    Fit(mName, mOptions, mType);
}

Pulse::Pulse(TString name, TH1* histo, TString options, int type)
{
    InitMembers();
    mName           = name;
    mHisto          = histo;
    mOptions        += options;
    mType           = type;
    Fit(mName, mOptions, mType);
}


Pulse::~Pulse()
{
    delete mModel;
}

// ===========================================================================
// PUBLIC OPERATIONS
// ===========================================================================

void
Pulse::Fit(
        TString fitName,
        TString fitOptions,
        int     type
        )
{
    int fitMin = mHisto->GetXaxis()->GetFirst();
    int fitMax = mHisto->GetXaxis()->GetLast();

    Fit(
        fitName,
        fitOptions,
        type,
        fitMin,
        fitMax
        );
}

void
Pulse::Fit(
        TString fitName,
        TString fitOptions,
        int     type,
        int     fitMin,
        int     fitMax
        )
{
    fitOptions += "R";
    if (type == 0)
        FitSectionWise(
            fitName,
            fitOptions,
            fitMin,
            fitMax
        );
    else if (type == 1)
        FitContious(
            fitName,
            fitOptions,
            fitMin,
            fitMax
        );
    else{
        cout << "not a correct type number --> fitting skipped!" << endl;
        return;
    }

    CalculateParameters();
    return;
}

// ===========================================================================
// PRIVATE OPERATIONS
// ===========================================================================

void
Pulse::InitMembers()
{
    mName           = "pulse";
    mHisto          = NULL;
    mModel          = NULL;
    mOptions        = "S";
    mBsl            = 0;
    mHeight         = 0;
    mT0             = 0;
    mT1             = 0;
    mTau1           = 0;
    mTau2           = 0;
    mIntegral       = 0;
    mAmplitude      = 0;
    mPhE            = 0;
    mType           = 0;
    mFitProb        = 0;
    mFitNCalls      = 0;
    mFitNdf         = 0;
    mChi2           = 0;
}

void
Pulse::FitSectionWise(
        TString fitName,
        TString fitOptions,
        int     fitMin,
        int     fitMax
        )
{
    mFitMin = fitMin;
    mFitMax = fitMax;
    mModel  = new TF1(fitName, shapeFunc, fitMin, fitMax, 6 );

    // ======================================================================
    // Calculate startvaues for fit
    double tau      = 30.0;
    double bsl = 0; //MW der ersten zehn slices

    int first_bin   = mHisto->GetXaxis()->GetFirst();

    for (int i = 0; i < 30; i++){
        bsl += mHisto->GetBinContent(first_bin+0);
    }
    bsl = bsl/30;

    double stop     = mHisto->GetMaximumBin(); //pos of max
    double height   = mHisto->GetBinContent(stop);
    double start    = stop-10; //pos 10 slices before maximum
    // ======================================================================

    double para[] = {bsl, height, start, stop, tau, tau};



    mModel->SetParameters(para);
    mModel->SetParNames("BSL", "A0", "t0", "t1", "Tau1", "Tau2");
    mModel->SetLineColor(kRed);

    mFitResultPtr = mHisto->Fit(mModel, fitOptions);

    mBsl            = mModel->GetParameter(0);
    mHeight         = mModel->GetParameter(1);
    mT0             = mModel->GetParameter(2);
    mT1             = mModel->GetParameter(3);
    mTau1           = mModel->GetParameter(4);
    mTau2           = mModel->GetParameter(5);
}

void
Pulse::FitContious(
        TString fitName,
        TString fitOptions,
        int     fitMin,
        int     fitMax
        )
{
    mFitMin = fitMin;
    mFitMax = fitMax;
    mModel  = new TF1(fitName, shapeFunc2, fitMin, fitMax, 5 );

    // ======================================================================
    // Calculate startvaues for fit
    double tau      = 30.0;
    double bsl = 0; //MW der ersten zehn slices

    int first_bin   = mHisto->GetXaxis()->GetFirst();

    for (int i = 0; i < 30; i++){
        bsl += mHisto->GetBinContent(first_bin+0);
    }
    bsl = bsl/30;

    double stop     = mHisto->GetMaximumBin(); //pos of max
    double amplitude= mHisto->GetBinContent(stop);
    double start    = stop-10; //pos 10 slices before maximum
    // ======================================================================


    double para[] = {bsl, amplitude, start, tau, tau};
    mModel->SetParameters(para);
    mModel->SetParNames("BSL", "A0", "t0", "Tau1", "Tau2");
    mModel->SetLineColor(kBlue);

    mFitResultPtr = mHisto->Fit(mModel, fitOptions);

    mBsl    = mModel->GetParameter(0);
    mHeight = mModel->GetParameter(1);
    mT0     = mModel->GetParameter(2);
    mTau1   = mModel->GetParameter(3);
    mTau2   = mModel->GetParameter(4);
}

void
Pulse::CalculateParameters()
{
    mIntegral   = mModel->Integral(mFitMin, mFitMax);
    mAmplitude  = mModel->GetMaximum() - mBsl;
    mFitProb    = mFitResultPtr->Prob();
    mFitNCalls  = mFitResultPtr->NCalls();
    mFitNdf     = mFitResultPtr->Ndf();
    mChi2       = mFitResultPtr->Chi2();
}

// ===========================================================================
// ACCESS
// ===========================================================================

TString Pulse::GetName(){ return mName;}
double  Pulse::GetBsl(){ return mBsl;}
double  Pulse::GetHeight(){ return mHeight;}
double  Pulse::GetT0(){ return mT0;}
double  Pulse::GetT1(){ return mT1;}
double  Pulse::GetTau1(){ return mTau1;}
double  Pulse::GetTau2(){ return mTau2;}
double  Pulse::GetIntegral(){ return mIntegral;}
double  Pulse::GetAmplitude(){ return mAmplitude;}
int     Pulse::GetPhe(){ return mPhE;}
int     Pulse::GetType(){return mType;}
int     Pulse::GetFitMin(){return mFitMin;}
int     Pulse::GetFitMax(){return mFitMax;}
double  Pulse::GetFitProb(){return mFitProb;}
double  Pulse::GetFitNCalls(){return mFitNCalls;}
double  Pulse::GetFitNdf(){return mFitNdf;}
double  Pulse::GetChi2(){return mChi2;}


// ===========================================================================
// NON CLASS MATH FUNCTIONS
// ===========================================================================

//double
int Heaviside(double val)
{
    if( val >= 0)   return 1;
    else return 0;
}

//Pulse::
double shapeFunc(
        double*     t,
        double*     par
        )
{
    double returnval= 0.0;

    // name parameters
    double bsl      = par[0];
    double height   = par[1];
    double start    = par[2];
//    double rising   = par[3];
    double tau1     = par[4];
    double tau2     = par[5];
    double stop     = par[3];

    // helper variables
//    double max      = start + rising * tau1;
    double e1       = (1 - Exp(-(t[0]-start)/tau1 ) );
//    double e2       = (-1) * height * (1 - Exp(-(x[0]-max)/tau2 ) );
    double e2       = (-1) * (1 - Exp(-(t[0]-stop)/tau2 ) );
    // calculate return value
    returnval += Heaviside(t[0]-start)*e1;
    returnval += Heaviside(t[0]-stop)*e2;
    returnval *= height;
    returnval += bsl;
//    if (t[0] > start)
//        returnval += e1;
//    if (t[0] > stop)
//        returnval += e2;

    return returnval;


}



/*
[18:42:13] Dominik: können wir nachher skypen? haste n headset?
[18:47:34] Dominik: hallo?
[18:53:38] Dominik: so also ich hab das file glaubich gefunden
[18:54:04] Dominik: height ist schonmal nicht die amplitude des pulses.
[18:54:52] Dominik: der puls ist ja sozusagen aus 2 funkionen zusammen gesetzt
[18:55:08] Dominik: die erste funktion, in meinem skript damals e1 genannt
[18:55:22] Dominik: ist ja ne umgekehrte nach obenverschobene e-fkt
[18:55:54] Dominik: und so wie sich ne normale abfallende e-fkt der null annähert, so nähert sich meine e1 dem wert 'height' an
[18:56:11] Dominik: allerdings kommt nun noch die abfallende e-fkt e2 dazu,
[18:57:20] Dominik: im Allgemeinen legt e2 schon los, bevor e1 auch nur in die Nähe von 'height' kommt und beginnt die Kurve nach unten zu drücken.
[18:57:37] Dominik: height ist also ein Skalenparameter, der nicht anschaulich ist
[18:59:27] Dominik: aber ich vermute, wenn man die Funktion nimmt und wie inder 9ten klasse die ableitung macht um die Position und Höhe des Maximums zu bestimmen, dann geht 'height' als linearer Faktor in die Höhe ein und spielt für die Position des Maximums keine Rolle.
Wenn das so ist, dann ist 'height' eine gute Schraube für den Fit-Algorithmus.... aber anschaulich wie gesagt nicht.

p0 = [-1., 10., 50., 70, 30., 30.] # Initial guess for the parameters

bsl = p[0]      # just the baseline
    height = p[1]   # this is not the pulse height, but the maximum of the discharging edge
                    # it *would* be the pulse height only in case the
                    # recharging process would set in later...

    start = p[2]    # start of avalanche aka start of discharge of diode capacitance
    stop = p[3]     # stop of avalanche, or start of recharging of diode
    tau1 = p[4]     # time constant of discharging process
    tau2 = p[5]     # time constant of recharging process

*/


double shapeFunc2(
        double*     t,
        double*     par
        )
{
    double returnval= 0.0;

    // name parameters
    double bsl      = par[0];
    double gain     = par[1];
    double t_0      = par[2];
    double tau1     = par[3];
    double tau2     = par[4];

    // helper variables
    double e1       = 1 - Exp(-(t[0]-t_0)/tau1 ) ;
    double e2       = Exp(-(t[0]-t_0)/tau2 ) ;

    // calculate return value
    returnval += bsl;
    returnval += gain*e1*e2*Heaviside(t[0]-t_0);

    return returnval;
}