#include <iostream>
#include "TFile.h"
#include "TDirectoryFile.h"
#include "TMath.h"
#include "TH1D.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TCut.h"
#include "TStyle.h"
#include "TMath.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TROOT.h"
#include "TLine.h"
#include "TFitResult.h"

const double pi = TMath::Pi();

double continuum(double* var, double* par){ // power law for the continuum
	double x = var[0];
	return par[0]*TMath::Power(x,par[1]);
}

double sumexpo(double* var, double* par){ // sum of 2 gaussian normalized
	double x = var[0];
	double arg1 = (x-par[1])/par[2];
	double arg2 = (x-par[1]-0.589)/par[4];
	return par[0]*0.05/(sqrt(2*pi)*par[2])*TMath::Exp(-0.5*arg1*arg1) + par[3]*0.05/(sqrt(2*pi)*par[4])*TMath::Exp(-0.5*arg2*arg2);
}

double expPp(double* var, double* par){ // first gaussian, Psi' pick
	double x = var[0];
	double arg = (x-par[1]-0.589)/par[2];
	return par[0]*0.05/(sqrt(2*pi)*par[2])*TMath::Exp(-0.5*arg*arg);
}

double expJP(double* var, double* par){ // second gaussian, J/Psi pick
	double x = var[0];
	double arg = (x-par[1])/par[2];
	return par[0]*0.05/(sqrt(2*pi)*par[2])*TMath::Exp(-0.5*arg*arg);
}

double fitfun(double* var, double* par){ // function for the fit, sum of the 2 gaussians and the power law
	return continuum(var,par) + sumexpo(var,&par[2]);
}

void fitmass_root_talk(const char* filein){

	gROOT->SetBatch(kTRUE);
	gStyle->SetOptFit(1);
	gStyle->SetOptStat(11);

	TFile* fin = new TFile(filein); // input file, the name is taken from the function argument
	fin->ls();
	fin->cd("M2_opt"); // there is a directory called M2_opt in this file

  // getting the histograms saved in the root file
	TH1D* hM2_opt_trig[3];
	hM2_opt_trig[0] = (TH1D*)gDirectory->Get("hM2_opt_trig0");
	hM2_opt_trig[1] = (TH1D*)gDirectory->Get("hM2_opt_trig1");
	hM2_opt_trig[2] = (TH1D*)gDirectory->Get("hM2_opt_trig2");

	fin->cd("M2_opt/M2_opt_r"); // idem
	TH1D* hM2_opt_r[4];
	hM2_opt_r[0] = (TH1D*)gDirectory->Get("hM2_opt_r0");
	hM2_opt_r[1] = (TH1D*)gDirectory->Get("hM2_opt_r1");
	hM2_opt_r[2] = (TH1D*)gDirectory->Get("hM2_opt_r2");
	hM2_opt_r[3] = (TH1D*)gDirectory->Get("hM2_opt_r3");

	fin->cd("M2_opt/M2_opt_deltat"); // idem
	TH1D* hM2_opt_delta[4];
	hM2_opt_delta[0] = (TH1D*)gDirectory->Get("hM2_opt_delta0");
	hM2_opt_delta[1] = (TH1D*)gDirectory->Get("hM2_opt_delta1");
	hM2_opt_delta[2] = (TH1D*)gDirectory->Get("hM2_opt_delta2");
	hM2_opt_delta[3] = (TH1D*)gDirectory->Get("hM2_opt_delta3");

	TH1D** hM2_opt[] = {hM2_opt_trig,hM2_opt_r,hM2_opt_delta}; // all the histograms in an array

	TFile* fout = new TFile("Root_talk_fit.root","recreate"); // output file
	cout << "Creato file di output" << endl;
  // inizialising the 3 canvas
  TCanvas* c[3];
	for(int u=0; u<3;u++){
		c[u] = new TCanvas(Form("c%i",u),Form("c%i",u));
		c[u]->Divide(2,2);
	}

  // inizialasing the functions
	TF1* f1 = new TF1("f1",continuum,1.1,8.,2);
 	f1->SetLineStyle(7);
	f1->SetLineWidth(1);
	f1->SetLineColor(kRed+2);
	TF1* f2 = new TF1("f2",expJP,1.1,8.,3);
	f2->SetLineStyle(3);
	f2->SetLineWidth(1);
	f2->SetLineColor(kMagenta+2);
	TF1* f3 = new TF1("f3",expPp,1.1,8.,3);
  f3->SetLineStyle(3);
	f3->SetLineWidth(1);
	f3->SetLineColor(kAzure-4);

	TF1* fall = new TF1("fall",fitfun,2.5,8.,7); // function for all

	TLine *x1fit = new TLine(2.5,0,2.5,4e4); // just two vertical lines at the start adn the end of the fit function
	TLine *x2fit = new TLine(8,0,8,50);
	x1fit->SetLineStyle(1);
	x2fit->SetLineStyle(1);
	x1fit->SetLineColor(kTeal-1);
	x2fit->SetLineColor(kTeal-1);

	for(int g=0;g<3;g++){ // cycle over the elements of the array hM2_opt
		for(int y=0;y<4;y++){ // cycle over the histos in each array in hM2_opt
		if(y==3 && g==0) continue; // the array of the first block of histograms arrives only at 3 histos, the other two blocks have 4 histos
		cout << "--------------------------FIT " << g << "." << y << " -------- " << hM2_opt[g][y]->GetTitle() << " ----------------------" << endl;
		c[g]->cd(y+1);
		c[g]->cd(y+1)->SetLogy();

		hM2_opt[g][y]->Draw(); // draw the histogram

		fall->SetParameter(0,hM2_opt[g][y]->GetBinContent(hM2_opt[g][y]->GetXaxis()->FindBin(1))); // I want this parameter to be the number at x=1 in the power law
    // this second parameter is the exponent of the power law, obtained as a division of log in natural base
    fall->SetParameter(1,TMath::Log((hM2_opt[g][y]->GetBinContent(hM2_opt[g][y]->GetXaxis()->FindBin(5)))/(hM2_opt[g][y]->GetBinContent(hM2_opt[g][y]->GetXaxis()->FindBin(6))))/TMath::Log(5./6));
		fall->SetParameter(2,1e5);
		fall->SetParName(2,"J/#Psi");
		fall->SetParameter(3,3.1);
		fall->SetParName(3,"M_{J/#Psi}");
		fall->SetParameter(4,0.18);
		fall->SetParName(4,"#sigma_{J/#Psi}");
		fall->SetParameter(5,0.5e5);
		fall->SetParName(5,"#Psi^{'}");
		fall->SetParameter(6,0.15);
		fall->SetParName(6,"#sigma_{#Psi^{'}}");

    // saving the fit result in r
		TFitResultPtr r = hM2_opt[g][y]->Fit(fall,"F E M S","",2.5,8.);

		cout << "Chi2: " << r->Chi2() << "\tdof: " << r->Ndf() << "\tchi2/dof: " << r->Chi2()/r->Ndf() << "\tProb: " << r->Prob() << endl;

    // below setting the functions of continuum and the 2 gaussians with the parameter from the fit, in order to plot on the same
    // graph gauss1, gauss2, continuum and the function-all

		//	continuum
		f1->FixParameter(0,r->Parameter(0));
		f1->FixParameter(1,r->Parameter(1));
		//	J/Psi
		f2->FixParameter(0,r->Parameter(2)); // n. J/psi
		f2->FixParameter(1,r->Parameter(3)); // mass j/psi
		f2->FixParameter(2,r->Parameter(4)); // sigma j/psi
		//	Psi'
		f3->FixParameter(0,r->Parameter(5)); // n. Psi'
		f3->FixParameter(1,r->Parameter(3)); // mass J/psi
		f3->FixParameter(2,r->Parameter(6)); // sigma Psi'

		x1fit->Draw("same");
		x2fit->Draw("same");
		f1->Draw("same");
		f2->Draw("same");
		f3->Draw("same");

		hM2_opt[g][y]->Write();
		} // chiuso 1st for

		c[g]->Write();
		c[g]->SaveAs(Form("c_root_talk%i.pdf",g));
	} // close for

	fin->Close();
	fout->Close();
}