// --------------------------------------------------------------------
// counts different events for a cut
long int cntEvt(TTree *t1, TString cut, TString var="ev")
{
	float fbev;
	int   ibev;
	
	TBranch *br = t1->GetBranch(var);
	if (br==0) return 0;
	
	TString tit = br->GetTitle(); // should end with '/I' or '/F'
	
	t1->SetEventList(0);
	t1->Draw(">>el",cut);
	TEventList *el=(TEventList*)gROOT->FindObject("el");
	
	t1->SetBranchStatus("*",0);
	t1->SetBranchStatus(var,1);
	
	// int or float branch
	bool isint=tit.EndsWith("/I");
	
	if (isint)	
		t1->SetBranchAddress(var,&ibev);
	else
		t1->SetBranchAddress(var,&fbev);
	
	// last and current event number; counter
	int lev = -1, cev=-1, cnt=0;
	for (int i=0;i<el->GetN();++i)
	{
		t1->GetEvent(el->GetEntry(i));
		
		cev = isint ? ibev : fbev;
		
		if (lev!=cev) cnt++;
		lev = cev;
	}
	
	t1->SetBranchStatus("*",1);
	
	return cnt;
}

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

void confgraph(TGraph *g, TString tit, int col=1, int marker=20)
{
	g->GetHistogram()->SetTitle(tit);
	g->GetHistogram()->SetMinimum(0);
	g->SetLineColor(col);
	g->SetMarkerColor(col);
	g->SetMarkerStyle(marker);
	g->SetMarkerSize(1);
	g->SetLineWidth(2);
}

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

TH1F* createHistoGraph(TGraph *g, TString tit="", double xmin=0, double xmax=0)
{
	static int cnt=0;
	if (tit=="") tit=g->GetTitle();
	
	double ymax  = TMath::MaxElement(g->GetN(), g->GetY());
	double dymax = g->GetErrorY(TMath::LocMax(g->GetN(), g->GetEY()));
	
	if (xmin>=xmax)
	{
		xmin = TMath::MinElement(g->GetN(), g->GetX());
		xmax = TMath::MaxElement(g->GetN(), g->GetX());
	}	
	TH1F *h=new TH1F(Form("h%03d",cnt++),tit,g->GetN(),xmin-0.5*(xmax-xmin)/g->GetN(), xmax+0.5*(xmax-xmin)/g->GetN());
	
	h->GetXaxis()->SetLabelSize(0.05);
	h->GetXaxis()->SetTitleSize(0.05);
	h->GetYaxis()->SetLabelSize(0.05);
	h->GetYaxis()->SetTitleSize(0.05);
	h->GetYaxis()->SetTitleOffset(1.4);
	
	h->SetMaximum((ymax+dymax)*1.05);
	h->SetMinimum(0);
	
	h->SetStats(0);
	
	return h;
}

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

void analyse_J_slc(TString fname="ntp1_slc_Jee_A.root", TString cut="xd0d0pide>0.5&&xd0d1pide>0.5&&chi24c<50", int sigmode=100, TString tit="")
{
	//data/DPMJee_A_SLC1        ( 366/ 500 files) =  29996980207
	//data/DPMJee_A_SLC2        ( 491/ 616 files) =  27740885200
	//data/DPMJee_A_SLC3        ( 460/ 604 files) =  26471309572
	//data/DPMJee_A_SLC4        ( 491/ 617 files) =  27745096509
	//data/DPMJee_A_SLC5        ( 485/ 624 files) =  27497769591
	//data/DPMJee_A_SLC6        ( 513/ 624 files) =  28645444363
	//data/DPMJee_A_SLC7        ( 564/ 659 files) =  25532064241
	//data/DPMJee_A_SLC8        ( 594/ 662 files) =  24338865274
	//data/DPMJee_A_SLC9        ( 502/ 503 files) =  20568203496
	//--------------------
	//data/DPMJee_B_SLC1        ( 297/ 500 files) =  24341219982
	//data/DPMJee_B_SLC2        ( 659/ 687 files) =  27009017984
	//data/DPMJee_B_SLC3        ( 656/ 682 files) =  26881191552
	//data/DPMJee_B_SLC4        ( 613/ 673 files) =  25119567078
	//data/DPMJee_B_SLC5        ( 612/ 665 files) =  25082906031
	//data/DPMJee_B_SLC6        ( 614/ 616 files) =  25163739980
	//data/DPMJee_B_SLC7        ( 608/ 610 files) =  24914199407
	//data/DPMJee_B_SLC8        ( 511/ 512 files) =  20942409620
	//data/DPMJee_B_SLC9        ( 507/ 507 files) =  20772725573
	//--------------------
	//--------------------
	//data/DPMJmm_A_SLC1        ( 272/ 273 files) =  20255703739
	//data/DPMJmm_A_SLC2        ( 200/ 200 files) =   7446866207
	//data/DPMJmm_A_SLC3        ( 214/ 214 files) =   7970194795
	//data/DPMJmm_A_SLC4        ( 202/ 202 files) =   7520702389
	//data/DPMJmm_A_SLC5        ( 214/ 214 files) =   7965350396
	//data/DPMJmm_A_SLC6        ( 250/ 250 files) =   9310510861
	//data/DPMJmm_A_SLC7        ( 244/ 244 files) =   9084144004
	//data/DPMJmm_A_SLC8        ( 247/ 247 files) =   9197632458
	//data/DPMJmm_A_SLC9        ( 233/ 233 files) =   8676693188
	//--------------------
	//data/DPMJmm_B_SLC1        ( 101/ 101 files) =   7519725318
	//data/DPMJmm_B_SLC2        ( 231/ 231 files) =   8595694989
	//data/DPMJmm_B_SLC3        ( 235/ 235 files) =   8751031806
	//data/DPMJmm_B_SLC4        ( 253/ 253 files) =   9421220844
	//data/DPMJmm_B_SLC5        ( 242/ 242 files) =   9010126829
	//data/DPMJmm_B_SLC6        ( 273/ 273 files) =  10164780274
	//data/DPMJmm_B_SLC7        ( 239/ 239 files) =   8896609302
	//data/DPMJmm_B_SLC8        ( 248/ 248 files) =   9235338848
	//data/DPMJmm_B_SLC9        ( 239/ 239 files) =   8900632658
	
	std::map<long int,long int> evcnts = 
                           { {1100, 29996980207}, {1101, 27740885200}, {1102, 26471309572}, {1103, 27745096509},  {1104, 27497769591}, {1105, 28645444363}, {1106, 25532064241}, {1107, 24338865274}, {1108, 20568203496},    // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
		                     {1120, 24341219982}, {1121, 27009017984}, {1122, 26881191552 }, {1123, 25119567078}, {1124, 25082906031}, {1125, 25163739980}, {1126, 24914199407}, {1127, 20942409620}, {1128, 20772725573},    // Setup B, pbp -> J/psi (-> e+ e-) pi+ pi-  //FIXME modes 1120 ...
		                 
		                     {1200, 20255703739}, {1201,  7446866207}, {1202,  7970194795},  {1203,  7520702389}, {1204,  7965350396}, {1205,  9310510861}, {1206,  8972578036}, {1207,  9197632458}, {1208,  8676693188},    // Setup A, pbp -> J/psi (-> mu+ mu-) pi+ pi- 
		                     {1220,  7519725318}, {1221,  8595694989}, {1222,  8751031806},  {1223,  9421220844}, {1224,  8935730839}, {1225,  9084618361}, {1226,  8896609302}, {1227,  9235338848}, {1228,  8900632658}     // Setup B, pbp -> J/psi (-> mu+ mu-) pi+ pi-  //FIXME modes 1220...
		                     };

		 
		                 
	//for ( auto x:evcnts) cout <<"mode "<<x.first<<" : "<<x.second<<endl;
	//cout <<endl;
		                 
	int bkgmode = sigmode+1000;
	//if (sigmode%100>19) bkgmode+=1;
	
	TFile *f = new TFile(fname);
	TTree *t = (TTree*)f->Get("ntp1");
	
	TFile fana("anaJ_slc.root","UPDATE");
	
	int Nred = 9;
	
	TGraphErrors *g[4];
	g[0] = new TGraphErrors(Nred);
	g[1] = new TGraphErrors(Nred);
	g[2] = new TGraphErrors(Nred); 
	g[3] = new TGraphErrors(Nred); 
	
	confgraph(g[0], "signal to noise");
	confgraph(g[1], "significance");
	confgraph(g[2], "signal efficiency");
	confgraph(g[3], "background efficiency");
	
	TF1 *f1[4];
	
	f1[0]=new TF1("f0","0.5*[0]*(1.0-TMath::Erf((x-[1])/[2]))+[3]",0,100.);	
	f1[0]->SetParameters(10,70, 10, 5);
	
	f1[1]=new TF1("f1","0.5*[0]*(1.0-TMath::Erf((x-[1])/[2]))+[3]",0,100.);	
   //1  p0           9.25779e-05   1.53699e-06   1.20127e-09   1.87008e+02
   //2  p1           5.00805e+01   2.94904e+00   6.21765e-06  -1.66797e-02
   //3  p2          -6.50046e+01   2.02989e+00   2.97643e-04** at limit **
   //4  p3           6.43684e-05   1.96720e-06   6.21319e-10   1.95246e+02
	f1[1]->SetParameters(1e-4, 50, -65, 6e-5);
	f1[1]->SetParLimits(1,20,200);
	f1[1]->SetParLimits(2,20,200);
	
	//f1[2]=new TF1("f2","0.5*[0]*(1.0-TMath::Erf((x-[1])/[2]))+[3]",0,100.);	
	//f1[2]->SetParameters(1e-7, 120, -10, 1e-8);
	//f1[2]->SetParLimits(1,20,200);

	f1[2]=new TF1("f2","expo(0)+pol1(2)",0,100.);	
	f1[2]->SetParameters(-0.5,-0.01,1,1);
	
	f1[3]=new TF1("f3","expo(0)+[2]",0,100.);	
	f1[3]->SetParameters(1.2,-0.02,1.6);  //
	

	//f1->SetParLimits(2,15,100);	
	//TF1 *f1=new TF1("f1","pol2(0)",0,100.);
	//f1->SetParameters(1,1,1);	
	
	// -------------------------------------------------------------------
	// 8x2 phi slices, symmetric
	// -------------------------------------------------------------------

	// missing fraction in [%] of EMC = [ 1.0 - (tht_max - 22°)/118° ]* 100
	double emc_rmv[9] =    {   0.,   12.5,  25.,  37.5,  50.,  62.5, 75., 87.5, 100. };

	//TString lab[4] = { ";supermodules missing;S/B", ";supermodules missing;significance [#sigma]", ";supermodules missing;signal efficiency [%]", ";supermodules missing;background efficiency [%]"};
	TString lab[4] = {  ";EMC missing [%];signal efficiency (#phi) [%]", ";EMC missing (#phi) [%];background efficiency [%]", ";EMC missing (#phi) [%];S/B", ";EMC missing (#phi) [%];significance [#sigma]"};
	
	
	TCanvas *c1 = new TCanvas("c1","c1",1000,800);
	c1->Divide(2,2,0.0001,0.0001);
	
	double sig_S = 50;
	double sig_B = 46e6;
	
	double Lint  = 1170*2;
	double fBR   = 0.05*0.06;
	
	double S_dat = sig_S * fBR * Lint;
	double B_dat = sig_B * Lint;
		
	cout <<"S:B = "<<S_dat/B_dat<<"  mode_S="<<sigmode<<"  mode_B="<<bkgmode<<endl; 
	
	printf (" i  |    S0    |     B0     |     S    |    B    |    fS    |    fB     |    S*   |     B*   |\n");
	printf ("----+----------+------------+----------+---------+----------+-----------+---------+----------+\n");
	for (int i=0;i<Nred;++i)
	{
		long int S0 = 1e5;
		long int B0 = evcnts[sigmode+i+1000];
	
		//cout <<"S0 = "<<S0<<"  B0 = "<<B0<<endl;
		
		double fS = S_dat/S0;
		double fB = B_dat/B0;
		
		TString sigcut = Form("%s && mode==%d", cut.Data(), sigmode+i);
		TString bkgcut = Form("%s && mode==%d", cut.Data(), bkgmode+i);
		
		
		double S  = (double) cntEvt(t, sigcut);
		double B  = (double) cntEvt(t, bkgcut);
		if (B==0) B=1;
		double dS = sqrt(S);
		double dB = sqrt(B);

		printf("%2d  |  %6ld  |  %6.3f G  |  %6d  |  %5d  |  %6.4f  |  %7.4f  |  %5.0f  |  %6.0f  |\n", i, S0, (double)B0/1e9, (int)S, (int)B, fS, fB, S*fS, B*fB);
		//cout <<i<<":  S="<<S<<"  B="<<B<<"  -->  S*"<<fS<<"="<<S*fS<<"  B*"<<fB<<"="<<B*fB<<endl;
		
		double SN  = S*fS/(B*fB);
		double dSN = SN*sqrt(dS*dS/(S*S) + dB*dB/(B*B));
		
		g[2]->SetPoint(i, emc_rmv[i], SN);
		g[2]->SetPointError(i, 0, dSN);
		
		
		double Z  = S*fS/sqrt(S*fS+B*fB);
		double dZ = 0.5 * sqrt( ((fS*fS*S+2*fB*fS*B)*(fS*fS*S+2*fB*fS*B)*S + fS*fS*fB*fB*S*S*B)/pow(fS*S+fB*B,3)) ;//0.5 * sqrt( (pow(2*B*fB + S*fS,2)*S + fS*fS*S*S*B)/pow(fS*S+fB*B,3));
		
		g[3]->SetPoint(i, emc_rmv[i], Z);
		g[3]->SetPointError(i, 0, dZ);
		
		
		double effS  = S/S0;
		double deffS = effS*sqrt( 1./S + 1./S0 ); 
		
		g[0]->SetPoint(i, emc_rmv[i], effS*100.);
		g[0]->SetPointError(i, 0, deffS*100.);
		
		
		double effB  = B/B0;
		double deffB = effB*sqrt( 1./B + 1./B0 ); 
		
		g[1]->SetPoint(i, emc_rmv[i], effB*100.);
		g[1]->SetPointError(i, 0, deffB*100.);
	}	
	
	TH1F *h[4];
	
	gStyle->SetOptFit(0);
	gStyle->SetOptStat(0);
	
	TString plotnam[4]={"effs","effb","sn","sign"};
	
	for (int i=0;i<4;++i) 
	{
		c1->cd(i+1);
		h[i] = createHistoGraph(g[i],tit+lab[i]);
		h[i]->Draw();
		g[i]->Draw("P same");
		
		//f1[i]->SetParameters(TMath::MaxElement(g[i]->GetN(), g[i]->GetY()),70, 10, TMath::MinElement(g[i]->GetN(), g[i]->GetY()));
		//if (i==1)
		//{
			//f1[i]->SetParameters(TMath::MaxElement(g[i]->GetN(), g[i]->GetY()),120, -10, TMath::MinElement(g[i]->GetN(), g[i]->GetY()));
			//f1[i]->SetParLimits(1,20,200);
		//}	
		
		g[i]->Fit(Form("f%d",i),"q");
			
		g[i]->SetName(plotnam[i]+"_"+((TString)fname(9,5)));
		g[i]->Write();
	}
	
	fname.ReplaceAll(".root","");
	cut.ReplaceAll("&&","_AND_");
	cut.ReplaceAll("||","_OR_");
	cut.ReplaceAll("!","_NOT_");
	cut.ReplaceAll(">","_lg_");
	cut.ReplaceAll("<","_sm_");
	cut.ReplaceAll(".","_");
	
	c1->SaveAs(Form("fig/%s__%03d__%s.gif",fname.Data(), sigmode, cut.Data()));
	c1->SaveAs(Form("fig/%s__%03d__%s.C",fname.Data(), sigmode, cut.Data()));
	
	fana.Close();
}


	//data/DPMJee_A_SLC1        ( 366/ 500 files) =  29996980207
	//data/DPMJee_A_SLC2        ( 186/ 500 files) =  15245616034
	//data/DPMJee_A_SLC3        ( 186/ 500 files) =  15245616034
	//data/DPMJee_A_SLC4        ( 186/ 500 files) =  15245616034
	//data/DPMJee_A_SLC5        ( 186/ 500 files) =  15245616034
	//data/DPMJee_A_SLC6        ( 186/ 500 files) =  15245616034
	//data/DPMJee_A_SLC7        (  59/ 500 files) =   4836581051
	//data/DPMJee_A_SLC8        ( 185/ 500 files) =   7580713308
	//data/DPMJee_A_SLC9        ( 211/ 211 files) =   8644837141
	//--------------------
	//data/DPMJee_B_SLC1        ( 297/ 500 files) =  24341219982
	//data/DPMJee_B_SLC2        ( 216/ 500 files) =   8855594650
	//data/DPMJee_B_SLC3        ( 208/ 500 files) =   8524796869
	//data/DPMJee_B_SLC4        ( 182/ 500 files) =   7456922637
	//data/DPMJee_B_SLC5        ( 224/ 500 files) =   9182421996
	//data/DPMJee_B_SLC6        ( 205/ 205 files) =   8400473125
	//data/DPMJee_B_SLC7        ( 207/ 208 files) =   8484522825
	//data/DPMJee_B_SLC8        ( 196/ 196 files) =   8034610317
	//data/DPMJee_B_SLC9        ( 190/ 190 files) =   7785364661
	
	//--------------------
	//--------------------
	//data/DPMJmm_A_SLC1        ( 272/ 273 files) =  20255703739
	//data/DPMJmm_A_SLC2        ( 200/ 200 files) =   7446866207
	//data/DPMJmm_A_SLC3        ( 214/ 214 files) =   7970194795
	//data/DPMJmm_A_SLC4        ( 202/ 202 files) =   7520702389
	//data/DPMJmm_A_SLC5        ( 214/ 214 files) =   7965350396
	//data/DPMJmm_A_SLC6        ( 250/ 250 files) =   9310510861
	//data/DPMJmm_A_SLC7        ( 241/ 241 files) =   8972578036
	//data/DPMJmm_A_SLC8        ( 247/ 247 files) =   9197632458
	//data/DPMJmm_A_SLC9        ( 233/ 233 files) =   8676693188
	//--------------------
	//data/DPMJmm_B_SLC1        ( 101/ 101 files) =   7519725318
	//data/DPMJmm_B_SLC2        ( 231/ 231 files) =   8595694989
	//data/DPMJmm_B_SLC3        ( 235/ 235 files) =   8751031806
	//data/DPMJmm_B_SLC4        ( 253/ 253 files) =   9421220844
	//data/DPMJmm_B_SLC5        ( 240/ 240 files) =   8935730839
	//data/DPMJmm_B_SLC6        ( 244/ 244 files) =   9084618361
	//data/DPMJmm_B_SLC7        ( 239/ 239 files) =   8896609302
	//data/DPMJmm_B_SLC8        ( 248/ 248 files) =   9235338848
	//data/DPMJmm_B_SLC9        ( 239/ 239 files) =   8900632658


	//std::map<long int,long int> evcnts = 
                           //{ {1100, 29996980207}, {1101, 15245616034}, {1102, 15245616034}, {1103, 15245616034},  {1104, 15245616034}, {1105, 15245616034}, {1106, 4836581051}, {1107, 7580713308}, {1108, 8644837141},    // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
		                     //{1120, 24341219982}, {1121,  8855594650}, {1122, 8524796869 }, {1123, 7456922637}, {1124, 9182421996}, {1125, 8400473125}, {1126, 8484522825},  {1127, 8034610317},  {1128, 7785364661},    // Setup B, pbp -> J/psi (-> e+ e-) pi+ pi-  //FIXME modes 1120 ...
		                 
		                     //{1200, 1}, {1201, 1}, {1202, 1},  {1203, 1}, {1204, 1}, {1205, 1}, {1206, 1},  {1207, 1},  {1208, 1},    // Setup A, pbp -> J/psi (-> mu+ mu-) pi+ pi- 
		                     //{1220, 1}, {1221, 1}, {1222, 1},  {1223, 1}, {1224, 1}, {1225, 1}, {1226, 1},  {1227, 1},  {1228, 1}     // Setup B, pbp -> J/psi (-> mu+ mu-) pi+ pi-  //FIXME modes 1220...
		                     //};
		 
	//std::map<long int,long int> evcnts = 
                           //{ {1100, 29996980207}, {1101, 15245616034}, {1102, 15245616034}, {1103, 15245616034},  {1104, 15245616034}, {1105, 15245616034}, {1106, 4836581051}, {1107, 7580713308}, {1108, 8644837141},    // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
		                     //{1120, 24341219982}, {1121,  8855594650}, {1122, 8524796869 }, {1123, 7456922637}, {1124, 9182421996}, {1125, 8400473125}, {1126, 8484522825},  {1127, 8034610317},  {1128, 7785364661},    // Setup B, pbp -> J/psi (-> e+ e-) pi+ pi-  //FIXME modes 1120 ...
		                 
		                     //{1200, 20255703739}, {1201, 7446866207}, {1202, 7970194795},  {1203, 7520702389}, {1204, 7965350396}, {1205, 9310510861}, {1206, 8972578036},  {1207, 9197632458},  {1208, 8676693188},    // Setup A, pbp -> J/psi (-> mu+ mu-) pi+ pi- 
		                     //{1220, 7519725318}, {1221, 8595694989}, {1222, 8751031806},  {1223, 9421220844}, {1224, 8935730839}, {1225, 9084618361}, {1226, 8896609302},  {1227, 9235338848},  {1228, 8900632658}     // Setup B, pbp -> J/psi (-> mu+ mu-) pi+ pi-  //FIXME modes 1220...
		                     //};