#include "TTree.h"

#include "TH1.h"



void optimal_bins()

{

	TChain myTree("myTree");  
	myTree.Add("/scratch_menkar/CMS_Trees/OniaTrees_2013_5TeV02_pPb/pPb_MC/OniaTree_MC_Ups1S_PA_5TeV02_WithFSR_tuneD6T.root");

    double nmin = 4.7e5;

	myTree.Draw("SumET_HF>>hpt(250,0,250)");

	TH1F *hpt = (TH1F*)gDirectory->Get("hpt");

	int nentries = hpt->Integral();
//int y = (int)round(x);
	int nbins = (int)TMath::Nint(nentries / nmin);
	cout<<"nbins: "<<nbins<<endl;

	double *binedges = new double[nbins + 1];

	double *bincontent = new double[nbins];

	cout << "Total events: " << nentries << ", " << nmin << " events per bin." << endl;

	int cnt = 0;

	int tot = 0;

	int ibin = 0;

	bool isfirst = true;

	for (int i = 1; i < hpt->GetNbinsX() + 1; i++)

	{

		cnt += hpt->GetBinContent(i);

		tot += hpt->GetBinContent(i);

		if (isfirst&&cnt > 0)

		{

			isfirst = false;

			binedges[0] = hpt->GetBinLowEdge(i);

		}

		if (cnt >= nmin)

		{

			bincontent[ibin] = cnt;

			cnt = 0;

			ibin++;

			binedges[ibin] = hpt->GetBinLowEdge(i + 1);

			// cout << ibin << " " << binedges[ibin] << " " << bincontent[ibin] << endl;

		}

	}



	if (ibin < nbins)

	{

		bincontent[ibin] = cnt;

		cnt = 0;

		ibin++;

		binedges[ibin] = 200.;

		// cout << ibin << "..." << endl;

	}



	for (ibin = 0; ibin < nbins; ibin++)

	{

		cout << ibin << ": [" << binedges[ibin] << "," << binedges[ibin + 1] << "]: " << bincontent[ibin] << endl;

		if (binedges[ibin + 1] == 200.) break;

	}

	for (int ibin = 0; ibin < nbins; ibin++)

	{

		if (binedges[ibin] == 200.) break;

		cout << binedges[ibin] << ", ";

	}

	cout << 200 << endl;
	cout<<nentries-tot<<endl;

}