#include <iostream>
#include <cmath>
#include <fstream>
#include <vector>

#include "TVector3.h"
#include "TLorentzVector.h"
#include "TF1.h"
#include "TH1D.h"
#include "TH1F.h"
#include "TH2D.h"
#include "TH2I.h"
#include "TTree.h"
#include "TLeaf.h"
#include "TCanvas.h"
#include "TBox.h"
#include "TText.h"
#include "TObjArray.h"
#include "TAxis.h"
#include "TColor.h"
#include "TMath.h"
#include "Math/BrentRootFinder.h"
#include "Math/WrappedTF1.h"
#include "TFile.h"
#include "TString.h"
#include "TLegend.h"

//#include "PndEmcRecoClust.h"

//----- Fit Function -----//

TF1 *gabler_f=0;

//-------------------------------------------------------------------
double _gablerFunction(double *x, double *par) {
	double val,arg1,arg2;
	double ep=par[3];
	arg1 = pow((x[0]-ep),2);
	arg2 = par[1];                                 // This was par[1]*par[1], before JGM, July 2014
	val = TMath::Exp(-4*TMath::Log(2)*arg1/arg2);
	if (x[0] > par[3])
		return par[0]*val;
	else 
		return par[0]*(val+TMath::Exp((x[0]-ep)/par[2])*(1-val));
}
//
double gablerFunction (double *x, double *p){
	return _gablerFunction(x,p) - p[0]/2.0;
}
//-------------------------------------------------------------------
template<class RF>
double findRoot( RF &r ) {
	double absTol = 1E-3; 
	double relTol = 1E-6; 
	r.Solve( 100, absTol, relTol); // status is somehow used by Root, so it needs to be enabled
	double root = r.Root();
	return  root; 
}
//
double gablerSolver(TF1 &funGab) {
	float amplitude,gamma,lambda,epeak,root1,root2;
	float fwhm;
	epeak = funGab.GetParameter(3);
	lambda = funGab.GetParameter(2);
	gamma = std::abs(funGab.GetParameter(1));
	amplitude = funGab.GetParameter(0);
	//
	TF1 *myfunc=new TF1("myfunc",gablerFunction,0,epeak+3*gamma,4);
	myfunc->SetParameter(3,epeak);
	myfunc->SetParameter(2,lambda);
	myfunc->SetParameter(1,gamma);
	myfunc->SetParameter(0,amplitude);
	//
	ROOT::Math::BrentRootFinder rf;
	ROOT::Math::WrappedTF1 F(*myfunc);
//	printf("Find root in %f %f :", 1., epeak);
	rf.SetFunction( F, 0., epeak);
	root1=findRoot(rf);
//	printf("%f\n",root1);
	rf.SetFunction( F, epeak, epeak+3*gamma);
//	printf("Find root in %f %f :", epeak, epeak+3*gamma);
	root2=findRoot(rf);
//	printf("%f\n",root2);
	fwhm=root2-root1;
	//
	delete myfunc;
	//
	return fwhm;
}
//-------------------------------------------------------------------
double Gabler(TH1 *hst, char *opt="", float *Pmean=0, int getmean=0) {
	//hst->GetXaxis()->SetRangeUser(1,4);
	double RMS      = hst->GetRMS();
//	double mean     = hst->GetMean();
	int    maxBin   = hst->GetMaximumBin();
	double maxCounts= hst->GetBinContent(maxBin);
	double maxpos   = hst->GetBinCenter(maxBin);
    
	//cout << "GetPars: mean= " << mean << ", RMS= " << RMS << ", maxBin= " << maxBin << ", maxCounts= " << maxCounts << ", maxpos= " << maxpos << endl;
    //
    if(!gabler_f)
	gabler_f = new TF1("gabler_f",_gablerFunction,0.05,0.15,4); // (char name, function, xmin, xmax, npar)
    //double param[4]={maxCounts,RMS/2,RMS/2,maxpos};
    double param[4]={maxCounts,RMS*RMS/4.,RMS/2.,maxpos};
    gabler_f->SetParameters(param);
    hst->Fit("gabler_f",opt,"",maxpos-RMS,maxpos+RMS); // was 2*RMS
    float fwhm = gablerSolver(*gabler_f);
    //printf("FWHM: %f\n",fwhm);
    //printf("Resolution: %f\n",fwhm/gabler_f->GetParameter(3));
    if(getmean) (*Pmean) = gabler_f->GetParameter(3);
    return fwhm/gabler_f->GetParameter(3);
 }

//----- Fit Function: END -----//


//-----  Group Clusters  -----//

void emc_reconstruct_fromclusters(const char* inFile) {
cout << endl;
cout << "//----- Start initialisation -----//" << endl;
cout << endl;
cout << "[info] => RUNNING TIMEBASED VERSION" << endl;
cout <<  endl;
cout << "[info] => Description: Function that attempts to reconstruct the invariant mass of the parent system from PndEmcClusters." << endl;

TFile *f = new TFile(inFile); 

TTree *t1 = (TTree*)f->Get("cbmsim"); 	// Get the Tree from the root file

int nEvents = t1->GetEntries();

TH1D *h1 = new TH1D("h1","Energy of clusters",400,0,4);
//TH1I *h7 = new TH1I("h7","Number of clusters per timebunch",40,0,20);
//TH1D *h12 = new TH1D("h12","2-photon invariant mass (selected on clusters with highest energy, timecut)",400,0,4);
TH1D *h13 = new TH1D("h13","N-photon invariant mass (selected on clusters above threshold)",1200,0,12);
TH1D *h14 = new TH1D("h14","h_{c} invariant mass (reco via #pi#pi#eta#gamma)",800,0,8);
TH1D *h15 = new TH1D("h15","#eta_{c} invariant mass",800,0,8);
TH1D *h16 = new TH1D("h16","h_{c} invariant mass (reco via #eta_{c})",800,0,8);
TH1D *h17 = new TH1D("h17","h_{c} invariant mass (reco via #eta_{c}) (angle cut)",800,0,8);
TH1D *h18 = new TH1D("h18","h_{c} invariant mass (reco via #pi#pi#eta#gamma) (angle cut)",800,0,8);
//TH2D *hMassT = new TH2D("hMassT","DeltaT between clusters VS invariant mass",150,0,3,300,0,30);
//TH2D *hE1 = new TH2D("hE1","Energy highest energy clusters VS invariant mass",200,0,2,110,0,2.2);
//TH2D *hE2 = new TH2D("hE2","Energy lowest energy clusters VS invariant mass",200,0,2,110,0,2.2);
//TH1D *hClusT = new TH1D("hClusT","Time between clusters",100,0,100);
//TH1D *hHitT = new TH1D("hHitT","Time between hits in clusters",100,0,100);
TH1I *hClusN = new TH1I("hClusN","Number of high-energy clusters per timebunch (from RECO)",80,0,80);
//TH2D *hNvsEn = new TH2D("hNvsEn","#Clusters/event vs Invariant mass of those clusters",40,0,20,300,2,5);
TH2I *hNparticles = new TH2I("hNparticles","Number of pions vs number of eta's per event",40,0,20,40,0,20);
TH2D *hAngle = new TH2D("hAngle","Angle between particles",360,0,180,800,0,8);
TH1I *hClusPerTB = new TH1I("hClusPerTB","Number of clusters per timebunch",400,0,400);
TH1D *h2fsp = new TH1D("h2fsp","2-photon invariant mass",400,0,4);
TH2I *hRecoPerEvent = new TH2I("hRecoPerEvent","#invariant mass combinations per event",5000,0,5000,20,0,20);
TH1I *hPhotonMultiplicity = new TH1I("hPhotonMultiplicity","Number of high-E clusters per timebunch",20,0,20);

int nProg = 0;
int hit = 0;
int nFinalStateParticles = 2; // tell the reco how many final state particles there should be
int nRecos = 0;
int nTooMuchRecos = 0;

double clusEnergy = 0;
double clusX = 0;
double clusY = 0;
double clusZ = 0;
double clusDir = 0;
double invMass = 0;

TLorentzVector myLorentzVector;
TLorentzVector combiVector;
TLorentzVector* invmassPtr;

std::vector<TLorentzVector> lorentzArray;

//PndEmcRecoClust reco;

int nHits = t1->Draw("EmcClusterTemp.fEnergy>>hist","","goff");

int len = 0;

cout << "[info] => " << nHits << " clusters in " << nEvents << " timebunches" << endl;

if (nHits == 0) {
	cout << "[warning] => No clusters found in this tree!" << endl;
	cout << "\n[info] => Aborting script execution\n" << endl;
	exit(0);
}

//ofstream ofile;
//ofile.open("LogInvMass_ppgg_5k_1GeV(STD).txt");

cout << endl;
cout << "//----- Initialisation successful, start Exec -----//" << endl;
cout << endl;
cout << "           0%              100%"<<endl;
cout << "PROGRESS: [" << std::flush;

for (int i=0; i<nEvents; i++) {
	// get clusters of this event
	TClonesArray* clusterArray=new TClonesArray("PndEmcCluster");
	t1->SetBranchAddress("EmcClusterTemp",&clusterArray);
	t1->GetEntry(i);
    len = clusterArray->GetEntriesFast(); 				// Gets number of clusters
//if (len != 0) cout << "Event: " << i << ", " << len << " hits." << endl;
if (len == 0) continue;
	lorentzArray.clear();

    for (int j=0; j<len; j++) {
		if (hit == nProg*nHits/100) { 	// filling the progress bar
			nProg += 5;
			/*progress += .009; 		// this condition is fulfilled 100 times, so in the end, progress will be increased by 100*0.009=0.9, the length of the progress bar	
			cProg.Draw();
			progText.DrawText(0.04,0.65,"Performing cluster finding (");
			CoverUp.SetFillColor(kWhite);
			CoverUp.Draw();
			progPerc.DrawText(0.52,0.65,Form("%2d%%)...",nProg));
			ProgressBarHolder.SetFillColor(kGray);
			ProgressBarHolder.Draw();
			TBox ProgressBar(0.04,0.3,progress,0.45); // progress bar needs to be redefined each time, or it won't be updated of course
			ProgressBar.SetFillColor(kGreen+2);
			ProgressBar.Draw();
			cProg.Update();*/
			cout << "/" << std::flush;
		}

		hit++; // current hit number; only used for the progress bar atm

		myLorentzVector.Clear();
		PndEmcCluster *myCluster=(PndEmcCluster*)clusterArray->At(j);

	// construct their corresponding Lorentz vectors

		clusEnergy = myCluster->GetEnergy();
/*		clusX = myCluster->x();
		clusY = myCluster->y();
		clusZ = myCluster->z();
		clusDir = TMath::Sqrt(clusX*clusX+clusY*clusY+clusZ*clusZ);

		myLorentzVector.SetPxPyPzE(clusX*clusEnergy/clusDir,clusY*clusEnergy/clusDir,clusZ*clusEnergy/clusDir,clusEnergy);*/

		myLorentzVector = myCluster->GetLorentzVector();

		h1->Fill(clusEnergy);

		if (clusEnergy > 0.6) nRecos++;

		lorentzArray.push_back(myLorentzVector);
	}
	hPhotonMultiplicity->Fill(nRecos);	

	nRecos = 0;
	// and reconstruct
	if (nFinalStateParticles == 2) {
		for (int b=0; b<len-1; b++) { // sum over all combinations of two photons
			if (lorentzArray[b].E()<0.50) continue; // skip low-energy clusters
			for (int d=b+1; d<len; d++) {
				if (lorentzArray[d].E()<0.50) continue;
				combiVector.Clear();
				combiVector = lorentzArray[b] + lorentzArray[d]; // sum Lorentz vectors to get inv mass of this combination
				h2fsp->Fill(combiVector.M();); // fill histo with inv mass of the 2 final state particles
				//ofile << combiVector.M() << endl; // also write values to file
				nRecos++;
			}
		}
	}
	if (nRecos>1) nTooMuchRecos++;
	hRecoPerEvent->Fill(i,nRecos);

	if (nFinalStateParticles > 2) {
		invmassPtr = lorentzArray;
		//invMass = reco.Reconstruct(len, invmassPtr, nFinalStateParticles, h14, hNparticles, h15, h16, h17, h18, hAngle);//, hE1, hE2);
		//h13->Fill(invMass);
		//hClusN->Fill(reco.nClustersPassed);
		invmassPtr = NULL;
	}

}
cout << "]" << endl;
//ofile.close();

//----- Fitting the histograms -----//
cout << "\n[task] => Obtaining resolution and related quantities..." << std::flush;
/*
h1->Fit("gaus","Q");
TF1* fit = h1->GetFunction("gaus");
Double_t Sigma = fit->GetParameter("Sigma");
Double_t Mean = fit->GetParameter("Mean");
Double_t Sum = h1->Integral((Mean-2*Sigma)/0.01,(Mean+2*Sigma)/0.01);
Double_t Res = 0;
Double_t Rest = h1->GetEntries()-Sum;
*/

//cout << " Res..." << std::flush;
double Res = Gabler(h1,"Q"); // Resolution := FWHM/mean
//cout << " Mean..." << std::flush;
double Mean = gabler_f->GetParameter(3); // Position of mean
//cout << " Sigma..." << std::flush;
double Sigma = (Res*Mean)/2.355; 
//cout << " Sum..." << std::flush;
double Sum = h1->Integral((Mean-3*Sigma)/0.01,(Mean+3*Sigma)/0.01); // integral of bin content in the fitted peak (xbin = x/binwidth)
double Rest = h1->GetEntries()-Sum;

cout << " done. " << endl;

cout << "[info] => Integration bounds: [" << Mean-3*Sigma << ", " << Mean+3*Sigma << "]" << endl;

//----- Printing table with fit results -----//

cout << "\nMean\tRes\tSigma\t#events in peak:\t#events outside:" << endl;
cout.precision(4);
cout << Mean << "\t" << Res << "\t" << Sigma << "\t" << Sum << "\t\t\t" << Rest << endl;// " (" << TMath::Floor(Sum/7) << " timebunches with at average 7 clusters)" << endl;

cout << "Multiplicity: " << TMath::Ceil(hPhotonMultiplicity->GetMean()) << endl;

gStyle->SetOptStat(1111111);
TCanvas *c1 = new TCanvas("c1","Cluster energies",1);
c1->cd(1);
h1->GetXaxis()->SetTitle("Cluster energies (GeV/c)");
h1->GetXaxis()->SetRangeUser(0,2);
h1->Draw("");

/*TCanvas *c2 = new TCanvas("c2","Cluster multiplicity",1);
c2->cd(1);
hPhotonMultiplicity->GetXaxis()->SetTitle("Number of clusters per timebunch with E > 600 MeV");
hPhotonMultiplicity->Draw("");*/

// to draw the rest, just create a TCanvas after the macro has finished and draw the histo you want in it
h13->GetXaxis()->SetTitle("M_{N#gamma} (GeV/c^{2})");
h13->GetYaxis()->SetTitle("Counts / 10 MeV/c^{2}");
h14->GetXaxis()->SetTitle("M_{#pi^{0}#pi^{0}#eta#gamma} (GeV/c^{2})");
h15->GetXaxis()->SetTitle("M_{#pi^{0}#pi^{0}#eta} (GeV/c^{2})");
h16->GetXaxis()->SetTitle("M_{#eta_{c}#gamma} (GeV/c^{2})");
h17->GetXaxis()->SetTitle("h_{c} invariant mass (reco via #eta_{c}) (angle cut)");
h18->GetXaxis()->SetTitle("h_{c} invariant mass (reco via #pi#pi#eta#gamma) (angle cut)");
hNparticles->GetXaxis()->SetTitle("##pi^{0}s vs # of clusters");
h2fsp->GetXaxis()->SetTitle("M_{#gamma#gamma} (GeV/c^{2})");
h2fsp->GetXaxis()->SetLabelSize(0.05);
h2fsp->GetYaxis()->SetLabelSize(0.05);
h2fsp->GetXaxis()->SetTitleSize(0.05);
h2fsp->GetYaxis()->SetTitleSize(0.05);
h2fsp->GetYaxis()->SetTitle("Counts per 10 MeV/c^{2}");
/*
TFile fout("histONL.root","recreate");
h2fsp->Write();
fout.Close();*/

cout << "[info] => " << nTooMuchRecos << " timebunches had more than 1 combination in them." << endl;
cout << endl;
cout << "//----- Finished - Script execution successful -----//" << endl;
cout << endl;

}