#include <iostream>
#include <fstream>

#include "TFile.h"
#include "TF1.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TFitResultPtr.h"
#include "TFitResult.h"
#include "Fit/Fitter.h"

using namespace std;

double DoubleTailedStepedGaussian(double*xx,double*pp)
{
	double x   = xx[0];

	Float_t f_tot = 0.;

	Float_t Back_const = pp[0];
	Float_t Back_slope = pp[1];

	f_tot += Back_const + x*Back_slope;

	Float_t Ampli     = pp[2];
	Float_t Mean      = pp[3];
	Float_t Sigma     = pp[4]*1./sqrt(8.*log(2.));
	Float_t Lambda    = pp[5];
	Float_t Rho       = pp[6];
	Float_t S         = pp[7];

	Float_t U         = (x-Mean)/Sigma;
	Float_t f_g       = Ampli*TMath::Exp(-U*U*0.5);
	Float_t f_lambda  = Ampli*TMath::Exp(-0.5*Lambda*(2.*U-Lambda));
	Float_t f_rho     = Ampli*TMath::Exp(-0.5*Rho*(2.*U-Rho));
	Float_t f_S       = Ampli*S*1./((1+TMath::Exp(U))*(1+TMath::Exp(U)));

	if(U<Lambda) f_tot += f_lambda;
	else if(U>Rho) f_tot += f_rho;
	else f_tot += f_g;

	f_tot += f_S;

	return f_tot;
}

double DoubleTailedStepedGaussian2(double*xx,double*pp)
{
    double x   = xx[0];

    int    NSubPeaks = (int)pp[0]; //Number of subpeaks in the peak range

    Float_t f_tot = 0.;

    Int_t Npar = 6;

    Float_t Back_const = pp[1];
    Float_t Back_slope = pp[2];

    f_tot += Back_const + x*Back_slope;

    for(int i=0 ; i<NSubPeaks ; i++) {
	Float_t Ampli     = pp[3+i*Npar+0];
	Float_t Mean      = pp[3+i*Npar+1];
	Float_t Sigma     = pp[3+i*Npar+2]*1./sqrt(8.*log(2.));
	Float_t Lambda    = pp[3+i*Npar+3];
	Float_t Rho       = pp[3+i*Npar+4];
	Float_t S         = pp[3+i*Npar+5];

	Float_t U         = (x-Mean)/Sigma;
	Float_t f_g       = Ampli*TMath::Exp(-U*U*0.5);
	Float_t f_lambda  = Ampli*TMath::Exp(-0.5*Lambda*(2.*U-Lambda));
	Float_t f_rho     = Ampli*TMath::Exp(-0.5*Rho*(2.*U-Rho));
	Float_t f_S       = Ampli*S*1./((1+TMath::Exp(U))*(1+TMath::Exp(U)));

	if(U<Lambda) f_tot += f_lambda;
	else if(U>Rho) f_tot += f_rho;
	else f_tot += f_g;

	f_tot += f_S;
    }

    return f_tot;
}

void FitTest()
{
	ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2","Migrad");
//	ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(100000000);		// Default: 0
//	ROOT::Math::MinimizerOptions::SetDefaultTolerance(1);				// Default: 0.010000
//	ROOT::Math::MinimizerOptions::SetDefaultStrategy(0);				// Default: 1


	TFile::Open("HistTest.root");
	TH1F *fHistogram = (TH1F*)gFile->Get("EGamma");
	fHistogram->Draw("hist");
	fHistogram->GetXaxis()->SetRangeUser(1300,1500);

	Float_t DefFWHM = 5.;
	Float_t fBackgd[2];
	fBackgd[0] = 1380;
	fBackgd[1] = 1430;
	Float_t fEnergies = 1408;

	TF1 *fFitFunction = new TF1("MyFit", DoubleTailedStepedGaussian2, fBackgd[0], fBackgd[1], 9);
	fFitFunction->SetParName(0, "NPeaks");
	fFitFunction->SetParName(1, "BkgConst");
	fFitFunction->SetParName(2, "BkgSlope");
	fFitFunction->SetParName(3, "Height");
	fFitFunction->SetParName(4, "Position");
	fFitFunction->SetParName(5, "FWHM");
	fFitFunction->SetParName(6, "LeftTail");
	fFitFunction->SetParName(7, "RightTail");
	fFitFunction->SetParName(8, "Step");
	fFitFunction->SetNpx(10000);

	fFitFunction->FixParameter(0,1);

	fHistogram->GetXaxis()->SetRangeUser(fBackgd[0],fBackgd[1]);

	//Calc Bckd
	fFitFunction->SetParameter(1, fHistogram->GetBinContent(fHistogram->FindBin(fBackgd[0])));
	fFitFunction->SetParLimits(1, fHistogram->GetMinimum(),fHistogram->GetMaximum());
	fFitFunction->SetParameter(2, 0);
	fFitFunction->SetParLimits(2, -10, 10);

	//Height
	fFitFunction->SetParameter(3, fHistogram->GetBinContent(fHistogram->FindBin(fEnergies)) - (fHistogram->GetBinContent(fHistogram->FindBin(fBackgd[0]))+fHistogram->GetBinContent(fHistogram->FindBin(fBackgd[1])))*0.5 );
	fFitFunction->SetParLimits(3, fHistogram->GetBinContent((fHistogram->FindBin(fEnergies))-(fHistogram->GetBinContent(fHistogram->FindBin(fBackgd[0]))+fHistogram->GetBinContent(fHistogram->FindBin(fBackgd[1])))*0.5)*0.5, fHistogram->GetBinContent(fHistogram->GetMaximumBin()));
	//Position
	fFitFunction->SetParameter(4, fEnergies);
	fFitFunction->SetParLimits(4, fEnergies-DefFWHM, fEnergies+DefFWHM);
	//FWHM
	fFitFunction->SetParameter(5, DefFWHM);
	fFitFunction->SetParLimits(5, 0, 10);
	//LeftTail
	fFitFunction->SetParameter(6, -2);
	fFitFunction->SetParLimits(6, -5, -0.1);
	//RightTail
	fFitFunction->SetParameter(7, 2);
	fFitFunction->SetParLimits(7, 0.1, 5.);
	//AmplitudeStep
	fFitFunction->SetParameter(8, 0.01);
	fFitFunction->SetParLimits(8, -1., 1.);

	fHistogram->GetXaxis()->SetRangeUser(1300,1500);

	TFitResultPtr r = fHistogram->Fit(fFitFunction,"R0SV");

	fFitFunction->Draw("same");
	gPad->SetLogy();

	if(r->Status()==0)cout<<"Successful fit !!!!"<<endl;
	else cout<<"ERROR !!!!! Status = "<< r->Status() << endl;
}