#include <TStyle.h>
#include <TTree.h>
#include <TFile.h>
#include <TClonesArray.h>
#include "TMath.h"

#include <TH1F.h>
#include <TH1I.h>
#include <TH2F.h>
#include <TH2I.h>

#include "TGraph.h"

#include "TStopwatch.h"
#include "TVirtualPad.h"

#include <cassert>
#include <cmath>

#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <unistd.h>
#include <complex>

using namespace std;
using std::vector;

#include "MpdADC64ECalDigit.h"

//------------------------------------------------------------------------------------------------------------------------
Double_t Novosibirsk(Double_t* x, Double_t* p) 
{
    //As Defined in RooNovosibirsk.cxx
    //If tail=eta=0 the Belle distribution becomes gaussia
    double tail = p[1];
    double width = p[2];
    double peak = p[3];
    if (TMath::Abs(tail) < 1.e-7)
        return p[0] * TMath::Exp(-0.5 * TMath::Power(((x[0] - peak) / width), 2));
    Double_t arg = 1.0 - (x[0] - peak) * tail / width;
    if (arg < 1.e-6) return 0.0; //Argument of logaritem negative. Real continuation -> function equals zero

    Double_t log = TMath::Log(arg);
    static const Double_t xi = 2.3548200450309494; // 2 Sqrt( Ln(4) )
    Double_t width_zero = (2.0 / xi) * TMath::ASinH(tail * xi * 0.5);
    Double_t width_zero2 = width_zero * width_zero;
    Double_t exponent = (-0.5 / (width_zero2) * log * log) - (width_zero2 * 0.5);

    return p[0] * TMath::Exp(exponent);
}

void FitOscillograms()
{
    vector<Double_t> data {12, 0, 15, -2, 6, -2, 6, 0, -14, -11, 22, -3, 21, 3, -22, -651, -2131, -3031, -3323, -2923, -2579, -2267, -1899, -1603, -1299, -1075, -875, -719, -563, -459, -347, -299, -219, -195, -147, -143, -103, -91, -48, -50, -23, -10, 3, 2, 16, 10, 32, 14, 17, -10, 6, -3, 4, 3, 17, 8, 7, -1, 13, 3}; 

    for (int j = 0; j < 60000 ; j++)
    {
        if(j%10 == 0) {
            cout << "\nCount# " << j << endl;
        }

        TGraphErrors* graph = new TGraphErrors(60);
        for(int i = 0; i < 60; i++)
        {
            graph->SetPoint (i, i, 0 - data[i]);
            graph->SetPointError(i, 0, 5);
        }

        Double_t MaxDelta = 3323;
        Double_t BinMax = 19;
                    
        TF1* pFit = new TF1("Novosibirsk", Novosibirsk, 0., 60, 4);
        pFit->SetParameters(MaxDelta, 0.05, 8., BinMax);

        TFitResultPtr r = graph->Fit(pFit, "SQR0", "");
        //r->Print(); 
        //TF1* pRes = graph->GetFunction("Novosibirsk");

        //(TVirtualFitter::GetFitter())->SetCache(0, 0);
        (TVirtualFitter::GetFitter())->Clear("s");

        delete (TVirtualFitter::GetFitter());
        delete graph;
        delete pFit;

    }
    
}

void example_fit
()
{
    TStopwatch timer;
    timer.Start();

    FitOscillograms();

    Double_t rtime = timer.RealTime();
    Double_t ctime = timer.CpuTime();
    cout << endl << endl;
    cout << "Macro finished successfully." << endl;      // marker of successful execution for CDASH
    cout << "Real time " << rtime << " s ~ " << rtime/60 << " mins, CPU time " << ctime/60 << " mins" << endl;
    cout << endl;
}