#include <string>
#include "TH1D.h"        //mirar root makefile
#include "RooRealVar.h"
#include "RooConstVar.h"
#include "RooGaussian.h"
#include "RooBreitWigner.h"
//#include "RooRelBreitWigner.h"
#include "RooArgusBG.h"
#include "RooAddPdf.h"
#include "RooAbsPdf.h"
#include "RooDataSet.h"
#include "TTree.h"
#include "TFile.h"
#include "TBranch.h"
#include "RooExponential.h"
#include "RooStats/SPlot.h"
#include "RooDataHist.h"
#include "RooExtendPdf.h"

using namespace RooFit;
using namespace RooStats;

void pruebavalor(){

	//importando el tree

		TFile *file = new TFile("mass.root");
		TTree* tree = (TTree*)file->Get("Tree;1");
		RooRealVar DP_Mass("DP_Mass","DP_Mass",2860,2900);
		//DP_Mass.setBins(200);

		RooDataSet data("data","dataset with x",RooArgSet(DP_Mass),Import(*tree)); //RooDataHist controlar los puntos
		//RooDataHist data("data","dataset with x",RooArgSet(DP_Mass),Import(*tree),200);

		//cout<<tree->GetMaximum("DP_Mass")<<endl;

	// creando BreitWigner

		RooRealVar mass1("mass1","mass1",2875, 2860, 2900);  //se puede binar
		RooRealVar width1("width1","width1",5,1,40);
		RooBreitWigner breit("BreitWigner1","BreitWigner1",DP_Mass,mass1,width1);

	// Ns  numero de eventos

		RooRealVar Ns("Ns", "numero de eventos", 1000, 2262182);
		RooExtendPdf ebreit("ebreit", "extended breit with a range", breit, Ns);

	//background con una exponencial

		RooRealVar a("a","parametro de la exponencial", -4600, 4600);
		RooExponential exponencial("exponencial","exponencial", DP_Mass, a);

	//Nb numero de background

		RooRealVar Nb("Nb","numero de background", 11000, 2262182);
		RooExtendPdf Nexponencial("Nexponencial", "extended exponencial with Nb", exponencial, Nb);

	// "suma" de los pdf

		RooAddPdf sigbac("sigbac", "signal + background", RooArgList(ebreit,Nexponencial));

		std::cout<<"\n -------------- Ajustando la BreitWigner y la exponencial------------- \n"<<std::endl;

		sigbac.fitTo(data);     //signal + background

	// ------------ imprimir la data sin SPlot

		std::cout<<"\n -------- data print antes del SWeight ------- \n"<<std::endl;

		data.Print();

	//Roostats   --Splot----Esta clase calcula los pesos que se utilizan para crear un sPlot

		std::cout<<"\n----------------calculando los Sweigth------------\n"<<std::endl;

		RooStats::SPlot *sData = new RooStats::SPlot("sData", "An SPlot", data, &sigbac, RooArgList(Ns,Nb) );
		std::cout<<" \n --------------  data print --------------- \n "<<std::endl;
		data.Print();

		std::cout << "\n ----------------mirar los SWeights ------------- \n"<<std::endl;
		std::cout << std::endl << "Yield de Ns (signal) " << Ns.getVal() << ". From SWeight it is "
				  << sData->GetYieldFromSWeight("Ns") << std::endl;
		std::cout << " Yield of background (Nb) " << Nb.getVal() << ". From SWeight it is "
				  << sData-> GetYieldFromSWeight("Nb") << std::endl;

	// calculando el peso para cada evento

	/* ----------------------------------------------------
		for (int i=0; i<228825; i++) {

			std::cout << "Ns Weight for event " << i << std::right << std::setw(12) << sData->GetSWeight(i,"Ns") << " Nb Weight "
					  << std::setw(12) << sData->GetSWeight(i, "Nb") << " Total Weight" << std::setw(12) << sData->GetSumOfEventSWeight(i)
					  << std::endl;
		}

		std::cout << std::endl;
		------------------------------------------------------- */
		data.Print();

	// create weighted data set

		RooDataSet * dataw_sig = new RooDataSet(data.GetName(),data.GetTitle(),&data,*data.get(),0,"Ns_sw");
		RooDataSet * dataw_bac = new RooDataSet(data.GetName(),data.GetTitle(),&data,*data.get(),0,"Nb_sw");


	//------------------ meter numero ente 0 y 1 para los eventos-------------  el numero aleatorio se mete en la Rootopla


	// plotear el ajuste de Breit-Wigner y la exponencial---------------------------------------------------------------------------------

		TCanvas *cdata = new TCanvas("sPlot", "sPlot", 2860, 2900);
		cdata->Divide(1);
		cdata->cd(1);

		RooPlot* xframe = DP_Mass.frame();   //se plotea puntos, que corresponden al tree DP_Mass con el ajuste
		data.plotOn(xframe) ;
		sigbac.plotOn(xframe);
		xframe->SetTitle("Fit");
		xframe->Draw();

	// plotear el SWeight
		TCanvas *cdata2 = new TCanvas("sPlot2", "sPlot2", 2860, 2900);
		cdata2->Divide(1);
		cdata2->cd(1);

		RooPlot* frame2 = DP_Mass.frame();
		dataw_sig->plotOn(frame2, RooFit::DataError(RooAbsData::SumW2));
		frame2->SetTitle("pesos signal");
		frame2->Draw();


		TCanvas *cdata3 = new TCanvas("sPlot3", "sPlot3", 2860, 2900);
		cdata3->Divide(1);
		cdata3->cd(1);

		RooPlot* frame3 = DP_Mass.frame();
		dataw_bac->plotOn(frame3, RooFit::DataError(RooAbsData::SumW2));
		frame3->SetTitle("pesos background");
		frame3->Draw();



		// ------------ creando la Rootopla, para el hacerle un RootBranch-----------


		std::cout<<"----------------------creando un arbol con los pesos------------------"<<std::endl;

		TFile f("arbolpruebapesos2.root","recreate");
		TTree tree2("Pesos", "Pesos");
		double S[1];
		double B[1];
		double N[1];
		tree2.Branch("Pesosig1",&S,"S/D");
		tree2.Branch("Pesobac1",&B,"B/D");
		tree2.Branch("numeracion",&N,"N/D");


		for (int i=0;i<228825; i++){
			S[0]= sData->GetSWeight(i,"Ns");
			B[0]=sData->GetSWeight(i, "Nb");
			N[0]= (i+0.0)/228824;
			tree2.Fill();

		}

		tree2.Print();
		//tree2.Scan();
		tree2.Write();
}