using namespace RooFit;

void bidi_fit(int doTemplates=1)
{
	// declare the generic variable
	const Double_t xLo   = 0;
	const Double_t xHi   = 3.5;
	RooRealVar Mass("dmfit","Lb-Lc",3.1,3.55,"GeV/c^{2}");
  	RooRealVar mrho46_gRandom("mrho46_gRandom","#rho_{46}",0,3.5,"GeV/c^{2}");
  	RooRealVar mrho56_gRandom("mrho56_gRandom","#rho_{56}",0,3.5,"GeV/c^{2}");
  	RooRealVar m3pi("m3pi","m3#pi",0,3.5,"GeV/c^{2}");



	if (doTemplates) {
		// Load Lc->a1 macro and do the template
		gROOT->LoadMacro("bidi_fit_Lca1.C");
		bidi_fit_Lca1(mrho46_gRandom,mrho56_gRandom);
	}
	// load Lca1 template
	TFile *Lca1_file = TFile::Open("Lca1_WS_template2D.root");
	RooWorkspace *wsLca1 = (RooWorkspace *) Lca1_file->Get("Lca1_ws");
	RooAbsPdf *pdf_a1 = wsLca1->pdf("totalPdf_Lca1Pdf");
// 	wsLca1->var("mrho46_gRandom")->setMin(xLo);
// 	wsLca1->var("mrho46_gRandom")->setMax(xHi);
// 	wsLca1->var("mrho56_gRandom")->setMin(xLo);
// 	wsLca1->var("mrho56_gRandom")->setMax(xHi);
	wsLca1->var("dm")->setConstant(false);
	wsLca1->var("k")->setConstant(false);	

	if (doTemplates) {
		// Load Lc->rho pi macro and do the template
		gROOT->LoadMacro("bidi_fit_Lcrho.C");
		bidi_fit_Lcrho(mrho46_gRandom,mrho56_gRandom);
	}
	// load Lca1 template
	TFile *Rhopi_file = TFile::Open("Rhopi_WS_template2D.root");
	RooWorkspace *wsRhopi = (RooWorkspace *) Rhopi_file->Get("Rhopi_ws");
	RooAbsPdf *pdf_rhopi = wsRhopi->pdf("totalPdf_RhopiPdf");
// 	wsRhopi->var("mrho46_gRandom")->setMin(xLo);
// 	wsRhopi->var("mrho46_gRandom")->setMax(xHi);
// 	wsRhopi->var("mrho56_gRandom")->setMin(xLo);
// 	wsRhopi->var("mrho56_gRandom")->setMax(xHi);
	wsRhopi->var("dm")->setConstant(false);
	wsRhopi->var("k")->setConstant(false);	

	if (doTemplates) {
		// Load Lc->nr macro and do the template
		gROOT->LoadMacro("bidi_fit_Lcnr.C");
		bidi_fit_Lcnr(mrho46_gRandom,mrho56_gRandom);
	}
	// load Lca1 template
	TFile *Lcnr_file = TFile::Open("Lcnr_WS_template2D.root");
	RooWorkspace *wsLcnr = (RooWorkspace *) Lcnr_file->Get("Lcnr_ws");
	RooAbsPdf *pdf_nr = wsLcnr->pdf("totalPdf_LcnrPdf");
// 	wsLcnr->var("mrho46_gRandom")->setMin(xLo);
// 	wsLcnr->var("mrho46_gRandom")->setMax(xHi);
// 	wsLcnr->var("mrho56_gRandom")->setMin(xLo);
// 	wsLcnr->var("mrho56_gRandom")->setMax(xHi);
	wsLcnr->var("dm")->setConstant(false);
// //  	wsLcnr->var("k")->setConstant(false);	


//________________rimetti questo_________________________________-
// 	TCanvas *canvas_test = new TCanvas("canvas_test","Test templates",1);
// 	RooPlot *test_frame = mrho46_gRandom.frame(Title("CDF Montecarlo"),Name("fit_templates"));
// 	pdf_a1->plotOn(test_frame,LineColor(kRed));
// 	pdf_rhopi->plotOn(test_frame,LineColor(kGreen));
//  	pdf_nr->plotOn(test_frame,LineColor(kViolet));
// 	test_frame->Draw();


	// load the sPlot rw sample for the final fit

	//______________________pdf per fit 2D_____________________________
	
// 	RooFitResult *fit= totalPdf.fitTo(data,Extended(kFALSE),SumW2Error(kFALSE),Hesse(kTRUE),Save(kTRUE),Minos(kTRUE));
// 	fit->Print("v");
// 
//         TH2D* pdf_tot = totalPdf.createHistogram("pdf_tot",mrho46_gRandom,YVar(mrho56_gRandom)) ; // for pdf
// 	pdf_tot->SetLineColor(kBlue) ;
// // 	TH2* dh2 = data->createHistogram("dg2",mrho46_gRandom,Binning(100),YVar(mrho56_gRandom,Binning(100))); // for data
// 	TH2* dh2 = data.createHistogram("dh2",mrho46_gRandom,Binning(100),YVar(mrho56_gRandom,Binning(100))); // for data ,Cut("m3pi>2")
// 	TCanvas *canvas_plot2D= new TCanvas("cplot2D","Plot",1);
// 	canvas_plot2D->Divide(1,2);
//         canvas_plot2D->cd(1);
// 	pdf_tot->Draw("SURFZ");
// 	canvas_plot2D->cd(2);
// 	dh2->Draw("LEGO");
	//___________________________________________________

// 	TFile *fitres = TFile::Open("gRandom_fitres_agnese.root");
//  	RooDataSet *rwDataset = (RooDataSet*) fitres->Get("dataset_rwLb");
	TFile *data_file = TFile::Open("gRandom_trig_check_all_data.root");
	TTree *data_tree = dynamic_cast<TTree*>(data_file->Get("lb"));
	RooDataSet rwDataset("data","Dataset",data_tree,RooArgSet(mrho46_gRandom,mrho56_gRandom));//
// 	data->reduce(Cut("m3pi>2")) ;
        cout << "=== Data loaded" << endl;
	

	
	RooRealVar Lca1_yield("Lca1_yield","N Lc->a1",200,1,1000);
	RooRealVar Lcrhopi_yield("Lcrhopi_yield","N Lc->rhopi",100,1,1000);
	RooRealVar Lcnr_yield("Lcnr_yield","N Lc->3pi",500,1,1000);
	
// 	RooExtendPdf *sig_a1 = new RooExtendPdf("sig_a1","Extendend signal PDF",*pdf_a1,Lca1_yield);
//   	RooExtendPdf *sig_rho = new RooExtendPdf("sig_rho","Extendend background PDF",*pdf_rhopi,Lcrhopi_yield);
// 	RooExtendPdf *sig_nr = new RooExtendPdf("sig_nr","Extendend background PDF",*pdf_nr,Lcnr_yield);
//   	RooAddPdf * model_2D = new RooAddPdf("Model", "right sign pdf", RooArgList(*sig_a1, *sig_rho, *sig_nr));
// 	model_2D->printCompactTree();
// 	model_2D->fitTo(rwDataset,NumCPU(2),Extended());
// 
//   cout<<" ============= More Accurate FIT ============="<<endl;
//   RooFitResult* res1 = 
//     model_2D->fitTo(rwDataset,NumCPU(2),Extended(),/*Minos(true),*/Strategy(2),Save(true),Verbose(false),Timer(kTRUE));
//   res1->Print("v");

// 	RooFitResult *fit_res = model_2D.fitTo(rwDataset46,Extended(kTRUE),Strategy(2),SumW2Error(kTRUE),Hesse(kTRUE),Save(kTRUE),Minos(kTRUE));
// 	TCanvas *canvas_final46 = new TCanvas("canvas_final","Final fit",1);
// 	RooPlot *fit_frame = Mass.frame(Title("CDF Preliminary 2.4 fb^{-1}"),Name("fit_frame"),Bins(100));//,Bins(50) //Bins(100)
// 	rwDataset46->plotOn(fit_frame,Name("Data"));
// 	model_2D.plotOn(fit_frame,Components("Model"),LineColor(kBlue),LineStyle(1),Name("Model"));
// 	model_2D.plotOn(fit_frame,Components("totalPdf_RhopiPdf"),LineColor(kRed),LineStyle(1),Name("Rhopi"));
//  	model_2D.plotOn(fit_frame,Components("totalPdf_Lca1Pdf"),LineColor(kGreen),LineStyle(1),Name("Lca1"));
// 	model_2D.plotOn(fit_frame,Components("totalPdf_LcnrPdf"),LineColor(kViolet),LineStyle(1),Name("Lcnr"));
// 	model_2D.paramOn(fit_frame) ;
	

	RooAddPdf model_2D("Model","Final Lb->Lc3pi 46",RooArgList(*pdf_a1,*pdf_rhopi,*pdf_nr),
				RooArgList(Lca1_yield,Lcrhopi_yield,Lcnr_yield));


	RooFitResult *fit_res = model_2D->fitTo(rwDataset,Extended(),Strategy(2),SumW2Error(kTRUE),Hesse(kTRUE),Save(kTRUE),Minos(kTRUE));

	TCanvas *canvas_final46 = new TCanvas("canvas_final","Final fit",1);
	RooPlot *fit_frame = Mass.frame(Title("CDF Preliminary 2.4 fb^{-1}"),Name("fit_frame"),Bins(100));//,Bins(50) //Bins(100)
	rwDataset.plotOn(fit_frame,Name("Data"));
	model_2D->plotOn(fit_frame,Components("Model"),LineColor(kBlue),LineStyle(1),Name("Model"));
	model_2D->plotOn(fit_frame,Components("totalPdf_RhopiPdf"),LineColor(kRed),LineStyle(1),Name("Rhopi"));
 	model_2D->plotOn(fit_frame,Components("totalPdf_Lca1Pdf"),LineColor(kGreen),LineStyle(1),Name("Lca1"));
	model_2D->plotOn(fit_frame,Components("totalPdf_LcnrPdf"),LineColor(kViolet),LineStyle(1),Name("Lcnr"));
	model_2D->paramOn(fit_frame) ;
	fit_frame->Draw();
        TLegend *leg = new TLegend(.75,.75,.95,.95);
   	leg->AddEntry(fit_frame->getHist("Data"),"Data","l");
   	leg->AddEntry(fit_frame->getCurve("Model"),"Model","l");
   	leg->AddEntry(fit_frame->getCurve("Lca1"),"Lca1 ","l");
   	leg->AddEntry(fit_frame->getCurve("Rhopi"),"Lcrhopi","l");
	leg->AddEntry(fit_frame->getCurve("Lcnr"),"Lc3pi(nr)","l");
   	
  	leg->Draw();
        
	fit_res->Print("v");

        // C r e a t e   w o r k s p a c e ,   i m p o r t   d a t a   a n d   m o d e l
  // -----------------------------------------------------------------------------

  // Create a new empty workspace
     RooWorkspace *w46 = new RooWorkspace("w46","workspace46") ;

  // Import model and all its components into the workspace
     w46->import(model_46) ;

  // Import data into the workspace
     w46->import(*rwDataset46) ;
     
  // Print workspace contents
     w46->Print() ;
    


  // S a v e   w o r k s p a c e   i n   f i l e
  // -------------------------------------------

  // Save the workspace into a ROOT file
  	w46->writeToFile("workspace46.root");

	TFile *ofile = TFile::Open("Lc3pi_final.root","recreate");
	fit_res->Write("fit_res");
	canvas_final46->Write("canvas_final46");
	canvas_final46->SaveAs("plot_final46.jpg");
        
	ofile->Close();
}