#include "TH2F.h"
#include "TH1F.h"
#include "TChain.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TStyle.h"
#include "TDirectory.h"
#include "TF1.h"
#include "TEfficiency.h"
#include "TColor.h"
#include "TRegexp.h"
#include "TMath.h"
#include "TROOT.h"

#include <iostream>
#include <vector>

double frsig, frmean, frchi, fremp;

// -----------------------------------------------------------------------

double Median(const TH2F *h,  const TH2F *h2) 
{
	int n1 = h->GetXaxis()->GetNbins(); 
	int n2 = h->GetYaxis()->GetNbins();
	
	int nfilled=0, N = h->GetEntries();
	
	double x[1000];
	double y[1000];
	
	for (int i=1;i<n1+1;++i)
	for (int j=1;j<n2+1;++j)
	{
		double cont = h->GetBinContent(i,j);
		double cont2 = h2->GetBinContent(i,j);
		
		if (cont>(0.001*N) && cont2>0.)
		{
			x[nfilled]   = cont2/cont;
			y[nfilled++] = 1.;
		}
	}
	return TMath::Median(nfilled, &x[0], &y[0]);
}

// -----------------------------------------------------------------------
void delete_histo(TString name)
{
	TH1 *h = (TH1*)gDirectory->Get(name);
	if (h) delete h;
}

// -----------------------------------------------------------------------

void palette2()
{  
/*   Double_t red[4]    = { 0.00, 176/256., 1.00, 1.00};
   Double_t grn[4]  = { 0.00, 53/256., 168/256., 1.00};
   Double_t blu[4]   = { 0.00, 0.00, 0.00, 0.50};
   Double_t stp[4] = { 0.00, 0.35 ,0.65, 1.00 };*/
//	int hex[8] = {0x482223,0x64394A,0x6D5877,0x5E7E9E,0x3EA5B0,0x48C9A9,0x8DE890,0xE7FE76};
// 	int hex[8] = {0x593432,0x6A4754,0x6A6175,0x5D7E8B,0x5A9A8D,0x7AB17F,0xB4C270,0xF8CB74};
// 	Double_t blu[8], red[8], grn[8], stp[8];
// 	
// 	for (int i=0;i<8;++i)
// 	{
// 		blu[7-i] = (hex[i] & 0xff)/255.;
// 		grn[7-i] = ((hex[i]>>8) & 0xff)/255.;
// 		red[7-i] = ((hex[i]>>16) & 0xff)/255.;
// 		stp[7-i] = (double)(7-i)/7.0;
// // 		cout << red[i]<<" "<<grn[i]<<" "<<blu[i]<<" "<<stp[i]<<endl;
// 	}
    Double_t blu[] = {1,1,0};
    Double_t red[] = {0,1,1};
    Double_t grn[] = {0,1,0};
    Double_t stp[] = {0,0.5,1.0};

	TColor::CreateGradientColorTable(3,stp,red,grn,blu,100) ;  
    gStyle->SetNumberContours(100);
/*	Double_t blu[3] = {1,0.5,0};
	Double_t red[3] = {1,0.5,0};
	Double_t grn[3] = {1,0.5,0};
	Double_t stp[3] = {0,0.5,1};
	
	Int_t n=100;
	TColor::CreateGradientColorTable(3,stp,red,grn,blu,n);*/
 }

// -----------------------------------------------------------------------

void cleanup()
{
	delete_histo("h1");
	delete_histo("h2");
	delete_histo("hp");
	delete_histo("heff");
	//delete_histo("hqa");
}

// -----------------------------------------------------------------------

void config_pad(TVirtualPad* p, double b=0.1, double r=0.12, double t=0.08,  double l=0.1)
{
	p->SetTopMargin(t);
	p->SetBottomMargin(b);
	p->SetLeftMargin(l);
	p->SetRightMargin(r);
/*	p->SetGridx();
	p->SetGridy();*/
}

// -----------------------------------------------------------------------

double qadalitz(TH2F* h1, TH2F* h2, TH1F *hpull)
{
	double globeff = h2->GetEntries()/h1->GetEntries();
	hpull->Reset();
	
	int nx=h1->GetNbinsX();
	int ny=h1->GetNbinsX();
	
	double sum=0;
	double coverage1=0, coverage2=0;
	
	for (int i=1;i<=nx;++i)
	for (int j=1;j<=ny;++j)
	{
		if (h1->GetBinContent(i,j)>0) coverage1+=1.;
		if (h2->GetBinContent(i,j)>0) coverage2+=1.;
		
		double c1 = h1->GetBinContent(i,j);
		double c2 = h2->GetBinContent(i,j);
		
		if (c1>0)
		{	
			double eff = c2/c1;
			double efferr = sqrt(2.)*eff/sqrt(c1);
			if (eff>0.0001) hpull->Fill((eff-globeff)/efferr);
		}
	}
	
	TF1 f1("f1","gaus");
	f1.SetParameters(hpull->GetMaximum(), hpull->GetMean(), hpull->GetRMS());
	hpull->Fit("f1","");
	
	double mean  = f1.GetParameter(1);
	double sigma = f1.GetParameter(2);
	double chindf = f1.GetChisquare()/f1.GetNDF(); 
	if (chindf<1.) chindf=1.;
	
	cout <<"m="<<mean<<" sigma="<<sigma<<" chi/ndf="<<chindf<<" emtpyfrac="<<fabs(1.0 - coverage2/coverage1)<<endl;
	return (fabs(sigma-1.) + fabs(mean) + fabs(1.0-chindf) + fabs(1.0 - coverage2/coverage1));
}
// -----------------------------------------------------------------------

double qadalitz2(TH2F* h1, TH2F* h2, TH1F *hpull, TH2F* heff)
{
	int nx=h1->GetNbinsX();
	int ny=h1->GetNbinsX();
	
	double N = h1->GetEntries();
	
	double coverage1=0, coverage2=0;
	for (int i=1;i<=nx;++i)
	for (int j=1;j<=ny;++j)
	{
		if (h1->GetBinContent(i,j)>0.0001*N) coverage1+=1.;
		if (h1->GetBinContent(i,j)>0.0001*N && h2->GetBinContent(i,j)>0) coverage2+=1.;
	}
	
	if (coverage2==0 || h1->GetEntries()<10 || h2->GetEntries()<10) return 10;
	
	double globeff = Median(h1,h2);//h2->GetEntries()/h1->GetEntries()*coverage1/coverage2;
	
	cout <<"c1="<<coverage1<<" c2="<<coverage2<<" geff="<<globeff<<endl;
	
	hpull->Reset();
	
	TEfficiency teff(*h2, *h1);
	teff.SetStatisticOption(TEfficiency::kFNormal);

	double perfbins=0.;
	
	for (int i=1;i<=nx;++i)
	for (int j=1;j<=ny;++j)
	{
		double c1 = h1->GetBinContent(i,j);
		double c2 = h2->GetBinContent(i,j);
		
		if (c1>0.0001*N)
		{	
			int gbin = teff.GetGlobalBin(i,j);
			
			double eff = teff.GetEfficiency(gbin);
			double efferr = teff.GetEfficiencyErrorUp(gbin);
			if (eff>0.0001 && efferr>0) 
			{
				hpull->Fill((eff-globeff)/efferr);
				heff->SetBinContent(i,j,(eff-globeff)/efferr);
			}
			if (eff>0.99999) perfbins+=1.;
		}
		else 
		{
			h1->SetBinContent(i,j,0);
			h2->SetBinContent(i,j,0);
		}
	}
	
	if (hpull->GetEntries()<5) return 0.001;
	
	TF1 f1("f1","gaus");
	f1.SetParameters(hpull->GetMaximum(), hpull->GetMean(), hpull->GetRMS());
	f1.SetParLimits(1,-8,8);
	f1.SetParLimits(2,0,20);
	
	hpull->Fit("f1","q");
	
	frmean   = f1.GetParameter(1);
	frsig    = f1.GetParameter(2);
	frchi    = f1.GetChisquare()/f1.GetNDF(); 
	fremp    = fabs(1.0 - coverage2/coverage1);
	double facperf  = fabs(1.0 - perfbins/coverage1);
	
	if (frchi<1.) frchi=1.;
	
	cout <<"m="<<frmean<<" sigma="<<frsig<<" chi/ndf="<<frchi<<" emtpyfrac="<<fremp<<"  not perf="<<facperf<<endl;
	return (fabs(frsig-1.) + fabs(frmean) + fabs(1.0-frchi) + fremp)*facperf;
}

// -----------------------------------------------------------------------
// -----------------------------------------------------------------------

int checkphsp2_2(TString filepatt, double step=.1, TString cut="")
{
	palette2();
	
	TCanvas *c2=new TCanvas("c2","c2",1000,10,800,600);
	gStyle->SetOptStat(0);
	gStyle->SetTitleH(0.07);
	
	if (cut=="") cut="1";
	
	TRegexp reg("MI[0-9][0-9][0-9]");
	TRegexp regn("ntp[0-9]");
	
	TString mname=filepatt(reg);
	mname = mname(2,3);
	
	TString ntp=filepatt(regn);

	TChain *n=new TChain(ntp);
	if (filepatt.EndsWith(".root")) filepatt.ReplaceAll(".root","");
	n->Add(filepatt+"*.root");

	
	double minE = int(n->GetMinimum("beamm")*10+0.1)/10.;
	double maxE = 5.5;
	int bins = int((maxE-minE+0.05)/step)+1;
	
	cout <<"minE="<<minE<<" maxE="<<maxE<<" bins="<<bins<<endl;

	TH1F *hqa = new TH1F("hqa","Phsp QA for M"+mname, bins, minE-step/2., maxE+step/2.);
	hqa->SetMarkerStyle(8);
	hqa->SetMarkerSize(1.5);
	hqa->SetMarkerColor(2);
	hqa->SetStats(0);
	hqa->SetYTitle("QA [AU]");
	hqa->SetXTitle("E_{cm} [GeV]");
	hqa->GetYaxis()->SetTitleOffset(1.2);
	
	int cnt=1;
	
	c2->cd();
	TH1F *htsten = new TH1F("htsten","tst", bins, minE-step/2., maxE+step/2.);
	n->Draw("beamm>>htsten",cut,"");

	int Npoints=0;
	for (int j=0;j<bins;++j) if (htsten->GetBinContent(j+1)>0) Npoints++;
	
	int wid = 230;
	TCanvas *c1=new TCanvas("c1","c1",10,10,wid*Npoints,wid*4);
	c1->Divide(Npoints,4);

	c2->cd();
	
	for (int j=0;j<bins;++j)
	{
		cleanup();
		if (htsten->GetBinContent(j+1)<1) continue;
		
		double ecm = hqa->GetBinCenter(j+1);
		cout <<mname<<" Ecm:"<<ecm<<endl;
		TString ecut = cut+TString::Format("&&abs(beamm-%f)<0.01",ecm);
		
		// find limits for Dalitz plot
		int nsel = n->Draw("xdal02:xdal12>>hhh",ecut,"goff");
		if (nsel>500000) nsel=500000;
		TH2F *h = (TH2F*)gDirectory->Get("hhh");
		if (h->GetEntries()==0) continue;
		
		double xmin=1000., ymin=1000., xmax=0., ymax=0.;
		double *x02=n->GetV2(), *x12 = n->GetV1();
		for (int i=0;i<nsel;++i)
		{
			if (x02[i]>xmax) xmax=x02[i];
			if (x02[i]<xmin) xmin=x02[i];
			if (x12[i]>ymax) ymax=x12[i];
			if (x12[i]<ymin) ymin=x12[i];
		}

		TH2F *h1 = new TH2F("h1",TString::Format("M%s @ %3.1f GeV",mname.Data(),ecm),20,xmin,xmax,20,ymin,ymax);
		TH2F *h2 = new TH2F("h2",TString::Format("M%s @ %3.1f GeV triggered",mname.Data(),ecm),20,xmin,xmax,20,ymin,ymax);
		TH2F *heff = new TH2F("heff",TString::Format("rel eff. M%s @ %3.1f GeV ",mname.Data(),ecm),20,xmin,xmax,20,ymin,ymax);
		TH1F *hp = new TH1F("hp",TString::Format("M%s @ %3.1f GeV pull",mname.Data(),ecm),50,-8,8);
		
		h1->Sumw2();
		h2->Sumw2();
		
		int N=n->GetEntries();
		n->Project("h1","xdal02:xdal12",ecut);
		n->Project("h2","xdal02:xdal12",ecut+"&&trig");
		h1->SetContour(100);
		h2->SetContour(100);
		
		//cout <<ecut<<"  : "<<h1->GetEntries()<<" "<<h2->GetEntries()<<endl;
		//c1->cd(1);h1->DrawCopy("colz");c1->cd(2);h2->DrawCopy("colz");c1->Update();
		//continue;
		double qa = qadalitz2(h1, h2, hp, heff);
		c1->cd(cnt); config_pad(gPad); h1->DrawCopy("colz");
		//h2->SetMaximum(h1->GetMaximum());
		c1->cd(cnt+Npoints); config_pad(gPad); h2->DrawCopy("colz");
		//h2->Divide(h1);
		//h2->SetMaximum(1.0);
		// convert in plot with sigma deviation
		//h2->SetTitle(TString::Format("rel eff. M%s @ %3.1f GeV ",mname.Data(),ecm));
		heff->SetMinimum(-8);
		heff->SetMaximum(8);
		
		c1->cd(cnt+Npoints*2); config_pad(gPad); heff->DrawCopy("colz");
		c1->cd(cnt+Npoints*3); gPad->SetTopMargin(0.08);//hp->Draw();
		hp->DrawCopy();
		
		TLatex lat;
		lat.DrawLatex(-7.5,hp->GetMaximum()*0.9,TString::Format("QA = %4.2f",qa));
		cout <<"QA = "<<qa<<endl;
		c1->Update();
		
		if (qa!=qa || qa>20) qa=20;
		
		hqa->SetBinContent(j+1,qa);
		
		c2->cd();
		gPad->SetGridx();gPad->SetGridy();
		config_pad(gPad,0.14,0.05,0.08,0.12);
		hqa->SetMaximum(21);
		hqa->Draw("P");
		c2->Update();
		cnt++;
	}
	
	c1->SaveAs("QA_plots_M"+mname+".gif");
	c2->SaveAs("QA_value_M"+mname+".gif");
  return 0;
}