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

using namespace RooFit;
using namespace std;

void homework_3()
{
	Double_t ann_low_cut=0.7732; 
	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,3500,0,350);
	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,3500,0,350);
	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 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",1.5,1,1.7);
	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("fit_result_testroofit_0.7732.eps");
	c1->Print("fit_result_testroofit_0.7732.bmp");

	cout<<"Gaussian mean value is :"<<sigmean.getVal()<<endl;
	cout<<"Gaussian #sigma value is :"<<sigwidth.getVal()<<endl;
	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;

	TLegend* leg = new TLegend(0.1,0.7,0.48,0.9);
	leg->SetHeader("The Legend ti","C");
	leg->AddEntry(mesframe->findObject("signalModel"),"signalModel","L");
	leg->AddEntry(mesframe->findObject("background"),"background","l");
	leg->AddEntry("background","background","l");
	mesframe->addObject(leg);
	leg->Draw();

}