#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()->GetBinLowEdge( mHisto->GetXaxis()->GetFirst() );
    int fitMax = mHisto->GetXaxis()->GetBinUpEdge(  mHisto->GetXaxis()->GetLast() );

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

void
Pulse::Fit(
        TString fitName,
        TString fitOptions,
        int     type,
        int     fitMin,
        int     fitMax
        )
{
//    mHisto->Sumw2();

    fitOptions += "R";
    if (type == 0)          //use sectionwise function
        FitSectionWise(
            fitName,
            fitOptions,
            fitMin,
            fitMax
        );
    else if (type == 1)     //use continous function
        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.first          = 0;
    mHeight.first       = 0;
    mT0.first           = 0;
    mT1.first           = 0;
    mTau1.first         = 0;
    mTau2.first         = 0;
    mIntegral.first     = 0;
    mAmplitude.first    = 0;
    mPhE.first          = 0;
    mBsl.second         = 0;
    mHeight.second      = 0;
    mT0.second          = 0;
    mT1.second          = 0;
    mTau1.second        = 0;
    mTau2.second        = 0;
    mIntegral.second    = 0;
    mAmplitude.second   = 0;
    mPhE.second         = 0;
    mRisetime.first     = 0;
    mRisetime.second    = 0;
    mOrder              = 0;
    mType               = 0;
    mFitProb            = 0;
    mFitNCalls          = 0;
    mFitNdf             = 0;
    mChi2               = 0;
    mTimeResolution     = 0.5;
}

// Fit pulseshape with a section wise defined funciton
// expressed by a sum of section wise defined exp-functions
//controled by heaviside function
void
Pulse::FitSectionWise(
        TString fitName,
        TString fitOptions,
        int     fitMin,
        int     fitMax
        )
{
    mFitMin = fitMin;
    mFitMax = fitMax;
    mModel  = new TF1(fitName, shapeFunc, fitMin, fitMax, 7 );

    // ======================================================================
    // Calculate startvaues for fit
    double tau = 30.0;
    double bsl = 0.8;     //shall be around 0

    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
    int    phe          = amplitude/10;
    // ======================================================================

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



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

    mFitResultPtr   = mHisto->Fit(mModel, fitOptions);
    mBsl.first      = mModel->GetParameter(0);
    mBsl.second     = mModel->GetParError(0);
    mHeight.first   = mModel->GetParameter(1);
    mHeight.second  = mModel->GetParError(1);
    mT0.first       = mModel->GetParameter(2);
    mT0.second      = mModel->GetParError(2);
    mT1.first       = mModel->GetParameter(3);
    mT1.second      = mModel->GetParameter(3);
    mTau1.first     = mModel->GetParameter(4);
    mTau1.second    = mModel->GetParError(4);
    mTau2.first     = mModel->GetParameter(5);
    mTau2.second    = mModel->GetParError(5);
    mPhE.first      = mModel->GetParameter(6);
    mPhE.second     = mModel->GetParError(6);
}

// Fit pulseshape with a continous funciton expressed by a product of two exp-functions
void
Pulse::FitContious(
        TString fitName,
        TString fitOptions,
        int     fitMin,
        int     fitMax
        )
{
    mFitMin = fitMin;
    mFitMax = fitMax;
    mModel  = new TF1(fitName, shapeFunc2, 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 amplitude= mHisto->GetBinContent(stop);
    double start    = stop-10; //pos 10 slices before maximum
    int    phe      = amplitude/10;
    // ======================================================================


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

    mFitResultPtr   = mHisto->Fit(mModel, fitOptions);
    mBsl.first      = mModel->GetParameter(0);
    mBsl.second     = mModel->GetParError(0);
    mHeight.first   = mModel->GetParameter(1);
    mHeight.second  = mModel->GetParError(1);
    mT0.first       = mModel->GetParameter(2);
    mT0.second      = mModel->GetParError(2);
    mTau1.first     = mModel->GetParameter(3);
    mTau1.second    = mModel->GetParError(3);
    mTau2.first     = mModel->GetParameter(4);
    mTau2.second    = mModel->GetParError(4);
    mPhE.first      = mModel->GetParameter(5);
    mPhE.second     = mModel->GetParError(5);

//    cout << "test" << mBsl.first << "\t +/-" << mBsl.second << endl;
}

void
Pulse::CalculateParameters()
{
    mIntegral.first     = mModel->Integral(mFitMin, mFitMax);
    //TODO setze konstantes glied auf 0 dann integrieren und fehler bestimmen
    mIntegral.second    = mModel->IntegralError(mFitMin, mFitMax);
    mAmplitude.first    = mModel->GetMaximum() - mBsl.first;
//    mAmplitude.second   = mModel->GetMaximum() - mBsl.first;
    mRisetime.first     = mTimeResolution*CalculateHistEdgeRisetime();
    mFitProb            = mFitResultPtr->Prob();
    mFitNCalls          = mFitResultPtr->NCalls();
    mFitNdf             = mFitResultPtr->Ndf();
    mChi2               = mFitResultPtr->Chi2();
}

double
Pulse::CalculateHistEdgeRisetime() // in bins
{
    //rise-tim for pulse for pulse height from 10% to 90%
    double riseTime = 0;

    // Get Pulse Maximum value and bin
    int maxBin = mHisto->GetMaximumBin();
    double max = mHisto->GetBinContent(maxBin);

    // values of bin with 10% and 90%
    int tenPctBin       = 0;
    int ninetyPctBin    = 0;

    // flag points
    bool finishedLow    = false;
    bool finishedHigh   = false;

    //2nd runvariable
    int j = maxBin;

    //loop over from 20 bins bevore max to maximum
    for (int i = maxBin - 20; i<maxBin; i++){
        j--;
        if (j<0) break;

        if (mHisto->GetBinContent(i) >= 0.1*max && !finishedLow)
        {
            tenPctBin = i;
            finishedLow = true;
        }


        if (mHisto->GetBinContent(j) <= 0.9*max && !finishedHigh)
        {
            ninetyPctBin = j;
            finishedHigh = true;
        }

        if (finishedLow && finishedHigh)
        {
            break;
        }
    }

    riseTime = mHisto->GetBinLowEdge(ninetyPctBin) - mHisto->GetBinLowEdge(tenPctBin);
    return riseTime;
}




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

TString Pulse::GetName(){ return mName;}
double  Pulse::GetBsl(){ return mBsl.first;}
double  Pulse::GetBslErr(){ return mBsl.second;}
double  Pulse::GetHeight(){ return mHeight.first;}
double  Pulse::GetHeightErr(){ return mHeight.second;}
double  Pulse::GetT0(){ return mT0.first;}
double  Pulse::GetT0Err(){ return mT0.second;}
double  Pulse::GetT1(){ return mT1.first;}
double  Pulse::GetT1Err(){ return mT1.second;}
double  Pulse::GetTau1(){ return mTau1.first;}
double  Pulse::GetTau1Err(){ return mTau1.second;}
double  Pulse::GetTau2(){ return mTau2.first;}
double  Pulse::GetTau2Err(){ return mTau2.second;}
double  Pulse::GetIntegral(){ return mIntegral.first;}
double  Pulse::GetIntegralErr(){ return mIntegral.second;}
double  Pulse::GetAmplitude(){ return mAmplitude.first;}
double  Pulse::GetAmplitudeErr(){ return mAmplitude.second;}
int     Pulse::GetPhE(){ return mPhE.first;}
int     Pulse::GetPhEErr(){ return mPhE.second;}
double  Pulse::GetRiseTime(){ return mRisetime.first; }
double  Pulse::GetRiseTimeErr(){ return mRisetime.second;}
int     Pulse::GetOrder(){ return mOrder;}
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;}
double  Pulse::GetTimeResolution(){return mTimeResolution;}
void    Pulse::SetTimeResolution(double timeResolution){mTimeResolution = timeResolution;}


// ===========================================================================
// 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 stop     = par[3];
    double tau1     = par[4];
    double tau2     = par[5];
    int    phe      = par[6];

    // 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 *= phe*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];
    int    phe      = par[5];

    // 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 += phe*gain*e1*e2*Heaviside(t[0]-t_0);

    return returnval;
}