//Be sure to add headfile if you want to compile this script

using namespace RooFit;
using namespace std;

void homework_4()
{	
	TH1D* h10 = new TH1D("significance","sig",200,0.4,1);

	Double_t ann_low_cut=0.4-0.6/100; 
	Double_t ann_up_cut=1.0; 
	Double_t mass,ann;
	TFile* file = new TFile("tree_group_01.root","r");
	TTree* tree = (TTree*) file->Get("tree_data");
	tree->SetBranchAddress("mass",&mass);
	tree->SetBranchAddress("ann",&ann);

	TCanvas* c1 = new TCanvas("c1","c1",1900,1080);
	gStyle->SetOptStat(0);	gStyle->SetStatY(0.98);
	gStyle->SetStatX(0.98);
	gStyle->SetStatW(0.3);
	gStyle->SetStatH(0.3);
	c1->SetGridx(1);
	c1->SetGridy(1);
	c1->SetTickx(1);
	c1->SetTicky(1);
	c1->SetFillColor(0);
	c1->SetBorderMode(0);
	c1->SetBorderSize(2);
	c1->SetLeftMargin(0.11);
	c1->SetFrameBorderMode(0);
	c1->SetFrameLineWidth(2);
	c1->Divide(2,1);

	TH1F* h1 = new TH1F("h1","mass",2000,0,200);
	TH2F* h2 = new TH2F("h2","Invariant mass",2000,0,200,5500,0,550);
	c1->cd(1);
	gPad->SetLeftMargin(0.12);
	gPad->SetRightMargin(0.12);
	gPad->SetTopMargin(0.12);
	gPad->SetBottomMargin(0.12);
	gPad->SetFrameLineWidth(0);
	gPad->SetFrameBorderMode(0);
	gPad->SetTickx(1);
	gPad->SetTicky(1);
	gPad->SetLogx(0);
	gPad->SetLogy(0);
	h2->SetLineStyle(1);
	h2->SetLineColor(2);
	h2->SetLineWidth(1);
	h2->GetXaxis()->SetTitle("mass");
	h2->GetXaxis()->CenterTitle(false);
	h2->GetXaxis()->SetLabelFont(42);
	h2->GetXaxis()->SetLabelSize(0.05);
	h2->GetXaxis()->SetTitleSize(0.05);
	h2->GetXaxis()->SetTitleOffset(0.8);  
	h2->GetXaxis()->SetNdivisions(8,5,0,kTRUE); 
	h2->GetYaxis()->SetNdivisions(7,9,0,kTRUE);
	h2->GetYaxis()->SetTitle( "count" );
	h2->GetYaxis()->CenterTitle(true);
	h2->GetYaxis()->SetLabelFont(42);
	h2->GetYaxis()->SetLabelSize(0.05);
	h2->GetYaxis()->SetTitleSize(0.05);
	h2->GetYaxis()->SetTitleOffset(1.3);
	h2->GetYaxis()->SetNdivisions(8,5,0,kTRUE);
	h2->GetYaxis()->SetTitleFont(42);

	TH2F* h3 = new TH2F("h3","Invariant mass",600,90,150,5500,0,550);
	c1->cd(2);
	gPad->SetLeftMargin(0.12);
	gPad->SetRightMargin(0.12);
	gPad->SetTopMargin(0.12);
	gPad->SetBottomMargin(0.12);
	gPad->SetFrameLineWidth(0);
	gPad->SetFrameBorderMode(0);
	gPad->SetTickx(1);
	gPad->SetTicky(1);
	gPad->SetLogx(0);
	gPad->SetLogy(0);
	h3->SetLineStyle(1);
	h3->SetLineColor(2);
	h3->SetLineWidth(1);
	h3->GetXaxis()->SetRangeUser(90,150);
	h3->GetXaxis()->SetTitle("mass");
	h3->GetXaxis()->CenterTitle(false);
	h3->GetXaxis()->SetLabelFont(42);
	h3->GetXaxis()->SetLabelSize(0.05);
	h3->GetXaxis()->SetTitleSize(0.05);
	h3->GetXaxis()->SetTitleOffset(0.8);  
	h3->GetXaxis()->SetNdivisions(8,5,0,kTRUE); 
	h3->GetYaxis()->SetRangeUser(0.,550);
	h3->GetYaxis()->SetNdivisions(7,9,0,kTRUE);
	h3->GetYaxis()->SetTitle( "count" );
	h3->GetYaxis()->CenterTitle(true);
	h3->GetYaxis()->SetLabelFont(42);
	h3->GetYaxis()->SetLabelSize(0.05);
	h3->GetYaxis()->SetTitleSize(0.05);
	h3->GetYaxis()->SetTitleOffset(1.3);
	h3->GetYaxis()->SetNdivisions(8,5,0,kTRUE);
	h3->GetYaxis()->SetTitleFont(42);

	Int_t n = tree->GetEntriesFast();
	cout<<"The entries number is: "<<n<<endl;
	for(Int_t j = 0; j < 200;j++){
		ann_low_cut = ann_low_cut+0.6/200;
		for(Int_t i = 0; i < n; i++){
			tree->GetEntry(i);
			if (ann > ann_low_cut && ann <ann_up_cut ){
				h1->Fill(mass,1);
			}
		}
	
	c1->cd(1);
	h2->Draw();
	h1->Draw("same");
	
	RooRealVar x1("x1","x1",90,150);
	RooRealVar sigmean("sigmean","mean",125,123,130);
	RooRealVar sigwidth("sigwidth","sigma",2,1,2);
	RooGaussian signalModel("signal","signal PDF",x1,sigmean,sigwidth);

	RooRealVar alpha("alpha","alpha",-3,-0.026,-5);
	RooExponential background("background","background PDF",x1,alpha);

	RooRealVar nsig("nsig","#signal events",0,300000);	
	RooRealVar nbkg("nbkg","#background events",0,300000);

	RooAddPdf model("model","g+e",RooArgList(signalModel,background),RooArgList(nsig,nbkg));

	RooDataHist data("data","data",x1,h1);

	RooFitResult*r = model.fitTo(data,Save(),Extended());

	RooPlot* mesframe = x1.frame();
	data.plotOn(mesframe);
	model.plotOn(mesframe);
	model.plotOn(mesframe,Components(background),LineStyle(kDashed),LineColor(kBlue));
	model.plotOn(mesframe,Components(signalModel),LineStyle(kDotted),LineColor(kRed));
	
	c1->cd(2);
	h3->Draw();
	mesframe->Draw("same");
	c1->Print(Form("./results/fit_result_testroofit_%f.png",ann_low_cut));
	h1->Reset();
	h3->Reset();
//	file<<"Gaussian mean value is :"<<sigmean.getVal()<<endl;
//	file<<"Gaussian #sigma value is :"<<sigwidth.getVal()<<endl;
	cout<<"Gaussian mean value is :"<<sigmean.getVal()<<"\t";
	cout<<"Gaussian #sigma value is :"<<sigwidth.getVal()<<"\t";
	r->Print();

	Double_t NSIG_nb = nsig.getVal();
	Double_t NBKG_nb = nbkg.getVal();
	cout<<"Signal significance is : "<<(NSIG_nb/TMath::Sqrt(NBKG_nb))<<endl;
	h10->Fill(ann_low_cut,(NSIG_nb/TMath::Sqrt(NBKG_nb)));
}
	TCanvas* c = new TCanvas();
	h10->Draw("hist");


/*	TLegend* leg = new TLegend(0.1,0.7,0.48,0.9);
	leg->SetHeader("The Legend ti","C");
	leg->AddEntry("model","background + signal","l");
	leg->AddEntry("signalModel","signal","f");
	leg->AddEntry("background","background","l");
	leg->Draw();
*/
}