void roofitQuestion(){

  gSystem->Load("libRooFit") ;
  using namespace RooFit;

  TFile* fInput  = new TFile("Input.root");

  TH1D *hdata = (TH1D*) fInput->Get("hdata")->Clone();

  TH1D *hs1 = (TH1D*) fInput->Get("hs1")->Clone();
  TH1D *hs2 = (TH1D*) fInput->Get("hs2")->Clone();
  TH1D *hs3 = (TH1D*) fInput->Get("hs3")->Clone();

  TH1D *hb1 = (TH1D*) fInput->Get("hb1")->Clone();
  TH1D *hb2 = (TH1D*) fInput->Get("hb2")->Clone();

  float n_hs  = hs1->Integral() + hs2->Integral() + hs3->Integral();
  float n_hb1 = hb1->Integral();

  double ratio_hs1 = hs1->Integral()/n_hs; 

  // Variable
  RooRealVar *x = new RooRealVar("x", "x", hs1->GetXaxis()->GetXmin(), hs1->GetXaxis()->GetXmax()); 
  
  // Fitted ratios
  RooRealVar  *fit_ratio_hs1   = new RooRealVar("fit_ratio_hs1",   "Fitted ratio 1",   ratio_hs1, 0.0, 0.5);   
  RooFormulaVar *fit_ratio_hs2 = new RooFormulaVar("fit_ratio_hs2", "FITTED ratio contrained", 
						   "fit_ratio_hs1*(2.01052)", 
						   RooArgList(*fit_ratio_hs1));
  
  // Other parameters 
  RooRealVar *n_hs_var  = new RooRealVar("n_hs_var",   "number of hs events",  n_hs);
  RooRealVar *n_hb1_var = new RooRealVar("n_hb1_var",  "number of hb1 events", hb1->Integral());
  RooRealVar *n_hb2_var = new RooRealVar("n_hb2_var",  "number of hb2 events", hb2->Integral());
  
   
  ///////////////////////////
  // WorkSpace: RooHistPdf //
  ///////////////////////////

  RooWorkspace *WS_PDF = new RooWorkspace("Fit using RooHistPdf");

  WS_PDF->import(*x);

  WS_PDF->import(*fit_ratio_hs1);
  WS_PDF->import(*fit_ratio_hs2);

  WS_PDF->import(*n_hs_var);
  WS_PDF->import(*n_hb1_var);
  WS_PDF->import(*n_hb2_var);

  // ----------------------------------------------------------------------------------
  RooDataHist *DH_hs1 = new RooDataHist("DH_hs1",   "hs1 Histogram", *x, hs1);
  RooDataHist *DH_hs2 = new RooDataHist("DH_hs2",   "hs2 Histogram", *x, hs2);
  RooDataHist *DH_hs3 = new RooDataHist("DH_hs3",   "hs3 Histogram", *x, hs3);

  RooDataHist *DH_hb1 = new RooDataHist("DH_hb1",   "hb1 Histogram", *x, hb1);
  RooDataHist *DH_hb2 = new RooDataHist("DH_hb2",   "hb2 Histogram", *x, hb2);

  RooDataHist *DH_hdata = new RooDataHist("DH_hdata",  "data Histogram", *x, hdata);
  // ----------------------------------------------------------------------------------
  RooHistPdf *PDF_hs1  = new RooHistPdf("PDF_hs1",   "PDF for DH_hs1", RooArgSet(*x), *DH_hs1);
  RooHistPdf *PDF_hs2  = new RooHistPdf("PDF_hs2",   "PDF for DH_hs2", RooArgSet(*x), *DH_hs2);
  RooHistPdf *PDF_hs3  = new RooHistPdf("PDF_hs3",   "PDF for DH_hs3", RooArgSet(*x), *DH_hs3);

  RooHistPdf *PDF_hb1  = new RooHistPdf("PDF_hb1",   "PDF for DH_hb1", RooArgSet(*x), *DH_hb1);
  RooHistPdf *PDF_hb2  = new RooHistPdf("PDF_hb2",   "PDF for DH_hb2", RooArgSet(*x), *DH_hb2);

  WS_PDF->import(*PDF_hs1);
  WS_PDF->import(*PDF_hs2);
  WS_PDF->import(*PDF_hs3);

  WS_PDF->import(*PDF_hb1);
  WS_PDF->import(*PDF_hb2);

  // Normalization constant K
  WS_PDF->factory("prod::kn_hs_var(k[1.0,0.0,2.0],n_hs_var)");
  
  // Models
  WS_PDF->factory("SUM::m_hs(fit_ratio_hs1*PDF_hs1,fit_ratio_hs2*PDF_hs2,PDF_hs3)");
  WS_PDF->factory("SUM::model(kn_hs_var*m_hs,n_hb1_var*PDF_hb1,n_hb2_var*PDF_hb2)");

  // Fit
  WS_PDF->pdf("model")->fitTo(*DH_hdata);


  //////////////////////////
  // WorkSpace: Function //
  /////////////////////////

  RooWorkspace *WS_FUN = new RooWorkspace("Fit using RooFunction");

  WS_FUN->import(*x);

  WS_FUN->import(*fit_ratio_hs1);
  WS_FUN->import(*fit_ratio_hs2);

  WS_FUN->import(*n_hs_var);
  WS_FUN->import(*n_hb1_var);
  WS_FUN->import(*n_hb2_var);

  TH1D * hs1_Norm   = (TH1D*)hs1->Clone();
  TH1D * hs2_Norm   = (TH1D*)hs2->Clone();
  TH1D * hs3_Norm   = (TH1D*)hs3->Clone();

  TH1D * hb1_Norm   = (TH1D*)hb1->Clone();
  TH1D * hb2_Norm   = (TH1D*)hb2->Clone();

  // Normalize Histograms
  hs1_Norm->Scale(1.0/hs1_Norm->Integral());
  hs2_Norm->Scale(1.0/hs2_Norm->Integral());
  hs3_Norm->Scale(1.0/hs3_Norm->Integral());

  hb1_Norm->Scale(1.0/hb1_Norm->Integral());
  hb2_Norm->Scale(1.0/hb2_Norm->Integral());

  RooDataHist *DH_hs1_Norm = new RooDataHist("DH_hs1_Norm",   "hs1 Histogram", *x, hs1_Norm);
  RooDataHist *DH_hs2_Norm = new RooDataHist("DH_hs2_Norm",   "hs2 Histogram", *x, hs2_Norm);
  RooDataHist *DH_hs3_Norm = new RooDataHist("DH_hs3_Norm",   "hs3 Histogram", *x, hs3_Norm);

  RooDataHist *DH_hb1_Norm = new RooDataHist("DH_hb1_Norm",   "hb1 Histogram", *x, hb1_Norm);
  RooDataHist *DH_hb2_Norm = new RooDataHist("DH_hb2_Norm",   "hb2 Histogram", *x, hb2_Norm);


  WS_FUN->import(*DH_hs1_Norm);
  WS_FUN->import(*DH_hs2_Norm);
  WS_FUN->import(*DH_hs3_Norm);

  WS_FUN->import(*DH_hb1_Norm);
  WS_FUN->import(*DH_hb2_Norm);

  // Functions
  WS_FUN->factory("HistFunc::f_hs1(x,DH_hs1_Norm)") ;
  WS_FUN->factory("HistFunc::f_hs2(x,DH_hs2_Norm)") ;
  WS_FUN->factory("HistFunc::f_hs3(x,DH_hs3_Norm)") ;

  WS_FUN->factory("HistFunc::f_hb1(x,DH_hb1_Norm)") ;
  WS_FUN->factory("HistFunc::f_hb2(x,DH_hb2_Norm)") ;


  // Normalization constant K
  WS_FUN->factory("prod::kn_hs_var(k[1.0,0.0,2.0],n_hs_var)");
  
  // Models
  WS_FUN->factory("ASUM::m_fhs(fit_ratio_hs1*f_hs1,fit_ratio_hs2*f_hs2,f_hs3)");
  WS_FUN->factory("ASUM::fmodel(kn_hs_var*m_fhs,n_hb1_var*f_hb1,n_hb2_var*f_hb2)");

  // Fit
  WS_FUN->pdf("fmodel")->fitTo(*DH_hdata);



  
}