C++ RooRealVar::frame方法代码示例

void Plot(RooAbsPdf *iGen,RooAbsPdf *iFit,int iN,int iNEvents,RooRealVar &iVar,RooRealVar &iSig,RooRealVar &iMean,
	  RooRealVar &iScale,RooRealVar &iRes,RooDataSet *iData=0) { 
  TRandom1 *lRand = new TRandom1(0xDEADBEEF);
  RooDataSet * lSignal  = iGen->generate(iVar,iNEvents);
  if(iData == 0) {
    iFit->fitTo(*lSignal,Strategy(1));
    if(iMean.getError() < 0.05) iFit->fitTo(*lSignal,Strategy(2));
  } else {
    iFit->fitTo(*iData,Strategy(1));
    if(iMean.getError() < 0.05) iFit->fitTo(*iData,Strategy(2));
  }    
  iVar.setBins(30);
  RooPlot *lFrame1 = iVar.frame(RooFit::Title("XXX")) ;
  if(iData == 0) lSignal->plotOn(lFrame1);
  if(iData != 0) iData->plotOn(lFrame1);
  iFit->plotOn(lFrame1);
  TCanvas *iC =new TCanvas("A","A",800,600);
  iC->cd(); lFrame1->Draw();
  iC->SaveAs("Crap.png");
  if(iData != 0) { 
    RooPlot *lFrame2 = iVar.frame(RooFit::Title("XXX")) ;
    iData->plotOn(lFrame2);
    iGen->plotOn(lFrame2);
    TCanvas *iC1 =new TCanvas("B","B",800,600);
    iC1->cd(); lFrame2->Draw();
    iC1->SaveAs("Crap.png");
  }
}

void FitterUtils::PlotShape2D(RooDataSet& originDataSet, RooDataSet& genDataSet, RooAbsPdf& shape, string plotsfile, string canvName, RooRealVar& B_plus_M, RooRealVar& misPT)
{
   //**************Prepare TFile to save the plots

   TFile f2(plotsfile.c_str(), "UPDATE");

   //**************Plot Signal Zero Gamma

   TH2F* th2fKey = (TH2F*)shape.createHistogram("th2Shape", B_plus_M, Binning(20), YVar(misPT, Binning(20)));
   cout<<genDataSet.sumEntries()<<endl;
   TH2F* th2fGen = (TH2F*)genDataSet.createHistogram("th2fGen", B_plus_M, Binning(20), YVar(misPT, Binning(20)));

   RooPlot* plotM = B_plus_M.frame();
   originDataSet.plotOn(plotM);
   shape.plotOn(plotM);

   RooPlot* plotMisPT = misPT.frame();
   originDataSet.plotOn(plotMisPT);
   shape.plotOn(plotMisPT);

   TCanvas canv(canvName.c_str(), canvName.c_str(), 800, 800);
   canv.Divide(2,2);
   canv.cd(1); th2fGen->Draw("lego");
   canv.cd(2); th2fKey->Draw("surf");
   canv.cd(3); plotM->Draw();
   canv.cd(4); plotMisPT->Draw();

   canv.Write();

   f2.Close();
}

//____________________________________
void MakePlots(RooWorkspace* wks) {

  // Make plots of the data and the best fit model in two cases:
  // first the signal+background case
  // second the background-only case.

  // get some things out of workspace
  RooAbsPdf* model = wks->pdf("model");
  RooAbsPdf* sigModel = wks->pdf("sigModel");
  RooAbsPdf* zjjModel = wks->pdf("zjjModel");
  RooAbsPdf* qcdModel = wks->pdf("qcdModel");

  RooRealVar* mu = wks->var("mu");
  RooRealVar* invMass = wks->var("invMass");
  RooAbsData* data = wks->data("data");


  //////////////////////////////////////////////////////////
  // Make plots for the Alternate hypothesis, eg. let mu float

  mu->setConstant(kFALSE);

  model->fitTo(*data,Save(kTRUE),Minos(kFALSE), Hesse(kFALSE),PrintLevel(-1));

  //plot sig candidates, full model, and individual componenets
  new TCanvas();
  RooPlot* frame = invMass->frame() ;
  data->plotOn(frame ) ;
  model->plotOn(frame) ;
  model->plotOn(frame,Components(*sigModel),LineStyle(kDashed), LineColor(kRed)) ;
  model->plotOn(frame,Components(*zjjModel),LineStyle(kDashed),LineColor(kBlack)) ;
  model->plotOn(frame,Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;

  frame->SetTitle("An example fit to the signal + background model");
  frame->Draw() ;
  //  cdata->SaveAs("alternateFit.gif");

  //////////////////////////////////////////////////////////
  // Do Fit to the Null hypothesis.  Eg. fix mu=0

  mu->setVal(0); // set signal fraction to 0
  mu->setConstant(kTRUE); // set constant

  model->fitTo(*data, Save(kTRUE), Minos(kFALSE), Hesse(kFALSE),PrintLevel(-1));

  // plot signal candidates with background model and components
  new TCanvas();
  RooPlot* xframe2 = invMass->frame() ;
  data->plotOn(xframe2, DataError(RooAbsData::SumW2)) ;
  model->plotOn(xframe2) ;
  model->plotOn(xframe2, Components(*zjjModel),LineStyle(kDashed),LineColor(kBlack)) ;
  model->plotOn(xframe2, Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;

  xframe2->SetTitle("An example fit to the background-only model");
  xframe2->Draw() ;
  //  cbkgonly->SaveAs("nullFit.gif");

}

void ComputeUpperLimit(RooAbsData *data, RooStats::ModelConfig *model, float &UpperLimit, float &signif, RooRealVar *mu, RooArgSet *nullParams,RooWorkspace *ws,REGION region,const char* tag) {

  bool StoreEverything=false; // activate if you want to store frames and all
  
  RooStats::ProfileLikelihoodCalculator *plc = new RooStats::ProfileLikelihoodCalculator(*data, *model);
  plc->SetParameters(*mu);
  plc->SetNullParameters(*nullParams);
  plc->SetTestSize(0.05);
  RooStats::LikelihoodInterval *interval = plc->GetInterval();

  bool ComputationSuccessful=false;
  UpperLimit = interval->UpperLimit(*mu,ComputationSuccessful);
  signif = 0.0; // plc->GetHypoTest()->Significance();   // deactivated significance (to make algorithm faster)

  if(!ComputationSuccessful) {
    cout << "There seems to have been a problem. Returned upper limit is " << UpperLimit << " but it will be set to -999" << endl;
    UpperLimit=-999;
    signif=-999;
  }

  if(StoreEverything) {
    // Store it all
    RooRealVar* minv = (RooRealVar*)model->GetObservables()->first();
    minv->setBins(static_cast<int>((minv->getMax()-minv->getMin())/5.));

    RooPlot* frameEE = minv->frame(RooFit::Title("ee sample"));
    frameEE->GetXaxis()->CenterTitle(1);
    frameEE->GetYaxis()->CenterTitle(1);
    
    RooPlot* frameMM = minv->frame(RooFit::Title("mm sample"));
    frameMM->GetXaxis()->CenterTitle(1);
    frameMM->GetYaxis()->CenterTitle(1);
    
    RooPlot* frameOF = minv->frame(RooFit::Title("OF sample"));
    frameOF->GetXaxis()->CenterTitle(1);
    frameOF->GetYaxis()->CenterTitle(1);
    
    data->plotOn(frameMM,RooFit::Cut("catCentral==catCentral::MMCentral"));
    model->GetPdf()->plotOn(frameMM,RooFit::Slice(*ws->cat("catCentral"), "MMCentral"),RooFit::ProjWData(*data));
    
    data->plotOn(frameEE,RooFit::Cut("catCentral==catCentral::EECentral"));
    model->GetPdf()->plotOn(frameEE,RooFit::Slice(*ws->cat("catCentral"), "EECentral"),RooFit::ProjWData(*data));
    
    data->plotOn(frameOF,RooFit::Cut("catCentral==catCentral::OFOSCentral"));
    model->GetPdf()->plotOn(frameOF,RooFit::Slice(*ws->cat("catCentral"), "OFOSCentral"),RooFit::ProjWData(*data));
    
    TFile *fout = new TFile("fout.root","UPDATE");
    frameMM->Write(Concatenate(Concatenate(data->GetName(),"_MM"),tag),TObject::kOverwrite);
    frameEE->Write(Concatenate(Concatenate(data->GetName(),"_EE"),tag),TObject::kOverwrite);
    frameOF->Write(Concatenate(Concatenate(data->GetName(),"_OF"),tag),TObject::kOverwrite);
    fout->Close();
  }

  delete plc;
  plc=0;
}

void rf706_histpdf_modified_makeroohistpdf()
{
  //  gROOT->SetBatch(kTRUE);

  //load data
  TFile *inHistos= new TFile("output/histos.root", "READ");
  TH1F * h100   = (TH1F*)inHistos->Get("cutmassLb7");
  TFile * inHistos2 = new TFile("../lambda/output/histos.root", "READ");
  TH1F * h200 = (TH1F*)inHistos2->Get("cutmassLb7");  

  double massLb = 5500;
  double massLbmin = 5300;
  double massLbmax = 5650;

  RooRealVar *mass = new RooRealVar("mass","m(#Lambda_{b})",massLb,"MeV");
  RooDataHist *hist1 = new RooDataHist("hist1","1D",RooArgList(*mass),h100);
  RooDataHist *hist2 = new RooDataHist("hist2","1D",RooArgList(*mass),h200);
  
  // Represent data in hist1 as pdf in *mass
  RooHistPdf histpdf1 = makeroohistpdf(hist1,mass) ;

  // Plot binned data and histogram pdf overlaid
  RooPlot* frame1 = mass->frame(Title("#Sigma^{0} data with RooHistPdf")) ;
  frame1->SetMaximum(1100);
  hist1->plotOn(frame1) ;
  histpdf1.plotOn(frame1) ;
  histpdf1.paramOn(frame1) ;
  
  TCanvas* c1 = new TCanvas("rf706_histpdf","rf706_histpdf",800,400) ;
  c1->Divide(1,2);
  c1->cd(1);
  frame1->Draw() ;

  //see what happens if I try to fit histpdf1 to lambda data
  RooHistPdf histpdf2 = makeroohistpdf(hist1,mass);
  RooAbsReal * nll = histpdf2.createNLL(*hist2,Extended(kTRUE));
  RooMinuit m(*nll);
  m.setVerbose(kTRUE);
  m.migrad();
  m.minos();
  RooPlot* frame2 = mass->frame(Title("#Lambda^{0} data fit with RooHistPdf from #Sigma^{0} data"));
  frame2->SetMaximum(1100);
  hist2->plotOn(frame2);
  histpdf2.plotOn(frame2);
  histpdf2.paramOn(frame2);
  c1->cd(2);
  frame2->Draw();

  //save histpdf1
  // cout<<"saving histpdf1..."<<endl;
  // histpdf1.RooHistPdf::SavePrimitive(std::cout);

  //  gROOT->SetPrimitive(kFALSE);
}

int main(){
  
  RooMsgService::instance().setGlobalKillBelow(ERROR);
  
  TFile *bkgFile = TFile::Open("comb_svn/hgg.inputbkgdata_8TeV_MVA.root");
  RooWorkspace *bkgWS = (RooWorkspace*)bkgFile->Get("cms_hgg_workspace");
  RooRealVar *mass = (RooRealVar*)bkgWS->var("CMS_hgg_mass");
  RooDataSet *data = (RooDataSet*)bkgWS->data("data_mass_cat0");

  RooRealVar *p1 = new RooRealVar("p1","p1",-2.,-10.,0.);
  RooRealVar *p2 = new RooRealVar("p2","p2",-0.001,-0.01,0.01);
  RooRealVar *p3 = new RooRealVar("p3","p3",-0.0001,-0.01,0.01);

  //RooPowerLawSum *pow1 = new RooPowerLawSum("pow","pow",*mass,RooArgList(*p1));
  //RooPowerLawSum *pow2 = new RooPowerLawSum("pow","pow",*mass,RooArgList(*p1,*p2));
  //RooPowerLawSum *pow3 = new RooPowerLawSum("pow","pow",*mass,RooArgList(*p1,*p2,*p3));

  RooExponentialSum *pow1 = new RooExponentialSum("pow1","pow1",*mass,RooArgList(*p1));
  RooExponentialSum *pow2 = new RooExponentialSum("pow2","pow2",*mass,RooArgList(*p1,*p2));
  RooExponentialSum *pow3 = new RooExponentialSum("pow3","pow3",*mass,RooArgList(*p1,*p2,*p3));

  TCanvas *canv = new TCanvas();
  RooPlot *frame = mass->frame();
  data->plotOn(frame);
  
  cout << "bus" << endl;
  pow1->fitTo(*data,PrintEvalErrors(-1),PrintLevel(-1),Warnings(-1),Verbose(-1));
  cout << "bus" << endl;
  pow1->plotOn(frame);
  pow2->fitTo(*data);
  pow2->plotOn(frame,LineColor(kRed),LineStyle(kDashed));
  pow3->fitTo(*data);
  pow3->plotOn(frame,LineColor(kGreen),LineStyle(7));
  frame->Draw();
  canv->Print("test.pdf");

  TFile *outFile = new TFile("test.root","RECREATE");

  TTree *tree = new TTree("tree","tree");

  vector<double> mu;
  mu.push_back(1.);
  mu.push_back(2.);
  mu.push_back(3.);
  mu.push_back(4.);
  vector<string> label;
  label.push_back("pol");
  label.push_back("pow");
  label.push_back("lau");
  label.push_back("exp");

  tree->Branch("mu",&mu);
  tree->Branch("label",&label);

  tree->Fill();
  outFile->cd();
  tree->Write();
  outFile->Close();
  return 0;
}

void chi2(int xmin = 0, int xmax = 200, TString filename="../DsubMC/met2j0bIso2ewkScale_0_200.root", int nparam = 2){
RooAbsData::ErrorType errorType = RooAbsData::SumW2;

file = new TFile(filename);
RooRealVar* h = new RooRealVar("h","h",xmin,xmax); 

//Get Data
TH1D* hData = file->Get("dataih");
hData->SetName("hData");
//hData->Draw();
RooDataHist* data = new RooDataHist("data","data",*h,hData);

//Get Summed MC
TH1D* hMC = file->Get("hh");
hMC->SetName("hMC");
hMC->Draw();
RooDataHist rdhMC("MC","MC",*h,hMC);
RooHistPdf pdfMC("MCpdf","MCpdf",*h,rdhMC);

//make (plot) the curves
RooPlot* hFrame = h->frame(Name("hFrame"));
data->plotOn(hFrame,RooFit::DataError(errorType));
pdfMC.plotOn(hFrame,ProjWData(*data),Components(pdfMC),Name("h_total"));

//Determine Chi^2
double chi2fit = hFrame->chiSquare("h_total", "h_data", nparam);
cout<<"Chi 2 / dof: "<<chi2fit<<endl;

//Determine K-S
double ks = hMC->KolmogorovTest(hData);
cout<<"Kolmogorov-Smirnov: "<<ks<<endl;
}

pair<double,double> bkgEvPerGeV(RooWorkspace *work, int m_hyp, int cat, int spin=false){
  
  RooRealVar *mass = (RooRealVar*)work->var("CMS_hgg_mass");
  if (spin) mass = (RooRealVar*)work->var("mass");
  mass->setRange(100,180);
  RooAbsPdf *pdf = (RooAbsPdf*)work->pdf(Form("pdf_data_pol_model_8TeV_cat%d",cat));
  RooAbsData *data = (RooDataSet*)work->data(Form("data_mass_cat%d",cat));
  RooPlot *tempFrame = mass->frame();
  data->plotOn(tempFrame,Binning(80));
  pdf->plotOn(tempFrame);
  RooCurve *curve = (RooCurve*)tempFrame->getObject(tempFrame->numItems()-1);
  double nombkg = curve->Eval(double(m_hyp));
 
  RooRealVar *nlim = new RooRealVar(Form("nlim%d",cat),"",0.,0.,1.e5);
  //double lowedge = tempFrame->GetXaxis()->GetBinLowEdge(FindBin(double(m_hyp)));
  //double upedge  = tempFrame->GetXaxis()->GetBinUpEdge(FindBin(double(m_hyp)));
  //double center  = tempFrame->GetXaxis()->GetBinUpCenter(FindBin(double(m_hyp)));

  nlim->setVal(nombkg);
  mass->setRange("errRange",m_hyp-0.5,m_hyp+0.5);
  RooAbsPdf *epdf = 0;
  epdf = new RooExtendPdf("epdf","",*pdf,*nlim,"errRange");
		
  RooAbsReal *nll = epdf->createNLL(*data,Extended(),NumCPU(4));
  RooMinimizer minim(*nll);
  minim.setStrategy(0);
  minim.setPrintLevel(-1);
  minim.migrad();
  minim.minos(*nlim);
  
  double error = (nlim->getErrorLo(),nlim->getErrorHi())/2.;
  data->Print(); 
  return pair<double,double>(nombkg,error); 
}

//#include "/uscms_data/d3/cvernier/DiH_13TeV/CMSSW_7_1_5/src/HiggsAnalysis/CombinedLimit/interface/RooMultiPdf.h"
//#include "HiggsAnalysis/CombinedLimit/interface/RooMultiPdf.h"
void makeRooMultiPdfWorkspace(){

   // Load the combine Library 
   gSystem->Load("libHiggsAnalysisCombinedLimit.so");

   // Open the dummy H->gg workspace 
   TFile *f_hgg = TFile::Open("w_background_Bern.root");
   RooWorkspace *w_hgg = (RooWorkspace*)f_hgg->Get("HbbHbb");
   // The observable (CMS_hgg_mass in the workspace)
   RooRealVar *mass =  w_hgg->var("x");

   // Get three of the functions inside, exponential, linear polynomial, power law
   RooAbsPdf *pdf_exp = w_hgg->pdf("bg_exp");
   RooAbsPdf *pdf_pol = w_hgg->pdf("bg");


   // Fit the functions to the data to set the "prefit" state (note this can and should be redone with combine when doing 
   // bias studies as one typically throws toys from the "best-fit"
   RooAbsData *data = w_hgg->data("data_obs");
   pdf_exp->fitTo(*data);  // index 0
   pdf_pol->fitTo(*data);   // index 2

   // Make a plot (data is a toy dataset)
   RooPlot *plot = mass->frame();   data->plotOn(plot);
   pdf_exp->plotOn(plot,RooFit::LineColor(kBlue));
   pdf_pol->plotOn(plot,RooFit::LineColor(kRed));
   plot->SetTitle("PDF fits to toy data");
   plot->Draw();

   // Make a RooCategory object. This will control which of the pdfs is "active"
   RooCategory cat("pdf_index","Index of Pdf which is active");

   // Make a RooMultiPdf object. The order of the pdfs will be the order of their index, ie for below 
   // 0 == exponential
   // 1 == linear function
   // 2 == powerlaw
   RooArgList mypdfs;
   mypdfs.add(*pdf_exp);
   mypdfs.add(*pdf_pol);
   
   RooMultiPdf multipdf("roomultipdf","All Pdfs",cat,mypdfs);
   
   // As usual make an extended term for the background with _norm for freely floating yield
   RooRealVar norm("roomultipdf_norm","Number of background events",0,10000);
   
   // Save to a new workspace
   TFile *fout = new TFile("background_pdfs.root","RECREATE");
   RooWorkspace wout("backgrounds","backgrounds");
   wout.import(cat);
   wout.import(norm);
   wout.import(multipdf);
   wout.Print();
   wout.Write();

}

void FitterUtils::plot_fit_result(string plotsfile, RooAbsPdf &totPdf, RooDataSet dataGenTot)
{

   //**************Prepare TFile to save the plots

   TFile f2(plotsfile.c_str(), "UPDATE");
   //**************Plot the results of the fit

   RooArgSet *var_set = totPdf.getObservables(dataGenTot);
   TIterator *iter = var_set->createIterator();
   RooRealVar *var;

   std::vector<RooPlot*> plots;
   RooPlot* frame;

   while((var = (RooRealVar*) iter->Next()))
   {

      frame = var->frame();
      dataGenTot.plotOn(frame);
      totPdf.plotOn(frame, Components("histPdfPartReco"), LineColor(kBlue));
      totPdf.plotOn(frame, Components("histPdfSignalZeroGamma"), LineColor(kGreen));
      totPdf.plotOn(frame, Components("histPdfSignalOneGamma"), LineColor(kMagenta));
      totPdf.plotOn(frame, Components("histPdfSignalTwoGamma"), LineColor(kOrange));
      totPdf.plotOn(frame, Components("histPdfJpsiLeak"), LineColor(14));
      totPdf.plotOn(frame, Components("combPDF"), LineColor(kBlack));
      totPdf.plotOn(frame, LineColor(kRed));

      plots.push_back(frame);

   }  

   if (!(plots.size())) return;

   TCanvas cFit("cFit", "cFit", 600, 800);
   cFit.Divide(1,2);
   cFit.cd(1); plots[0]->Draw();
   if (plots.size()>1){ 
      cFit.cd(2); plots[1]->Draw();
   }

   cFit.Write();
   f2.Close();


}

void plotFromWorkspace() {
    TFile* file = new TFile("card_m125_1JetIncl_XX_workspace.root");

    RooWorkspace* w = (RooWorkspace*)file->Get("w");

    RooRealVar* m = (RooRealVar*)w->var("CMS_hmumu_mass");
    RooRealVar* e = (RooRealVar*)w->var("CMS_hmumu_merr");
    RooDataSet* d = (RooDataSet*)w->data("data_pseudo");
    RooAbsPdf * b = w->pdf("bkg_mass_merr_1JetIncl_XX_pdf");
    RooAbsPdf * s = w->pdf("sig_mass_merr_ggH_1JetIncl_XX_pdf");

    RooPlot* frame = m->frame();
    d->plotOn(frame);
    b->plotOn(frame);
    //s->plotOn(frame, RooFit::ProjWData(*e, *d), RooFit::LineColor(kOrange+1));
    //s->plotOn(frame, RooFit::LineColor(kOrange+1));
    frame->Draw();
}

void simul(){

	TFile* in = new TFile( "../bin/Kdata.root", "READ" );

	TH1D* h = in->Get("tof_8");

	RooRealVar * x = new RooRealVar( "x", "x", -10, 10 );

	RooDataHist * rdh = new RooDataHist( "data", "data", RooArgSet( *x ), h);

	RooRealVar m( "mean", "mean", 0, -1, 1 );
	RooRealVar s( "sigma", "sigma", .5, 0, 1 );
	RooGaussian g( "gauss", "gauss", )

	RooPlot * frame = x->frame();
	rdh->plotOn( frame );
	frame->Draw();
	

}

void plot_with_residuals(TCanvas& _c,RooPlot& _frame,RooRealVar& var,Int_t& _nbins,Double_t& r_min,Double_t& r_max) {

    _c.SetName("_c");
    _c.Divide(1,2);
    _c_1->SetPad(0.01,0.16,0.99,0.99);
    _c_2->SetPad(0.01,0.01,0.99,0.16);

    _c.cd(1);
    _frame.Draw();

    _c.cd(2);
    RooHist* hpull = _frame.pullHist();

    RooPlot* frame_pull = var.frame(Title(" "),Range(r_min,r_max),Bins(_nbins));
    frame_pull->addPlotable(hpull,"P");

    // this is if you want residuals in number of events
    /*RooHist* hresid = _frame -> residHist();*/
    /*RooPlot* frame2 = m_Kpipi.frame(1000,2000,25);*/
    /*frame2->addPlotable(hresid,"P");*/

    frame_pull->SetMinimum(-5);
    frame_pull->SetMaximum(+5);
    frame_pull->SetNdivisions(0,"x"); // 510 for having also an x scale here
    frame_pull->SetNdivisions(203,"y");   // axis divisions was 510
    frame_pull->SetXTitle(" ");
    frame_pull->SetLabelSize(0.15,"Y");
    frame_pull->Draw();

    TLine *_line = new TLine(r_min,-3.,r_max,-3);
    _line->SetLineStyle(1);
    _line->SetLineColor(2);
    _line->SetLineWidth(1);
    _line->Draw();
    _line2 = new TLine(r_min,+3.,r_max,+3);
    _line2->SetLineStyle(1);
    _line2->SetLineColor(2);
    _line2->SetLineWidth(1);
    _line2->Draw();

}

void FitterUtilsSimultaneousExpOfPolyTimesX::plot_kemu_fit_result(string plotsfile, RooAbsPdf &totKemuPdf, RooDataSet const& dataGenKemu)
{

   //**************Prepare TFile to save the plots

   TFile f2(plotsfile.c_str(), "UPDATE");
   //**************Plot the results of the fit

   RooArgSet *var_set = totKemuPdf.getObservables(dataGenKemu);
   TIterator *iter = var_set->createIterator();
   RooRealVar *var;

   std::vector<RooPlot*> plots;
   RooPlot* frame;

   while((var = (RooRealVar*) iter->Next()))
   {
      frame = var->frame();
      dataGenKemu.plotOn(frame);
      totKemuPdf.plotOn(frame, LineColor(kRed));

      plots.push_back(frame);
   }

   if (!(plots.size())) return;

   TCanvas cFit("cKemuFit", "cKemuFit", 600, 800);
   cFit.Divide(1,2);
   cFit.cd(1); plots[0]->Draw();
   if (plots.size()>1){ 
      cFit.cd(2); plots[1]->Draw();
   }

   cFit.Write();
   f2.Close();
}

void FitBias(TString CAT,TString CUT,float SIG,float BKG,int NTOYS)
{
  gROOT->ForceStyle();
  
  RooMsgService::instance().setSilentMode(kTRUE);
  RooMsgService::instance().setStreamStatus(0,kFALSE);
  RooMsgService::instance().setStreamStatus(1,kFALSE);
  
  // -----------------------------------------
  TFile *fTemplates = TFile::Open("templates_"+CUT+"_"+CAT+"_workspace.root");
  RooWorkspace *wTemplates = (RooWorkspace*)fTemplates->Get("w");
  RooRealVar *x            = (RooRealVar*)wTemplates->var("mTop");
  RooAbsPdf *pdf_signal    = (RooAbsPdf*)wTemplates->pdf("ttbar_pdf_Nominal");
  RooAbsPdf *pdf_bkg       = (RooAbsPdf*)wTemplates->pdf("qcdCor_pdf"); 
  TRandom *rnd = new TRandom();
  rnd->SetSeed(0);
  x->setBins(250);   
  RooPlot *frame;

  TFile *outf;

  if (NTOYS > 1) { 
    outf = TFile::Open("FitBiasToys_"+CUT+"_"+CAT+".root","RECREATE");
  }

  float nSigInj,nBkgInj,nSigFit,nBkgFit,eSigFit,eBkgFit,nll;

  TTree *tr = new TTree("toys","toys");
  
  tr->Branch("nSigInj",&nSigInj,"nSigInj/F");
  tr->Branch("nSigFit",&nSigFit,"nSigFit/F");
  tr->Branch("nBkgInj",&nBkgInj,"nBkgInj/F");
  tr->Branch("nBkgFit",&nBkgFit,"nBkgFit/F");
  tr->Branch("eSigFit",&eSigFit,"eSigFit/F");
  tr->Branch("eBkgFit",&eBkgFit,"eBkgFit/F");
  tr->Branch("nll"    ,&nll    ,"nll/F");

  for(int itoy=0;itoy<NTOYS;itoy++) {
    // generate pseudodataset
    nSigInj = rnd->Poisson(SIG);
    nBkgInj = rnd->Poisson(BKG);
    RooRealVar *nSig = new RooRealVar("nSig","nSig",nSigInj);
    RooRealVar *nBkg = new RooRealVar("nBkg","nBkg",nBkgInj);
    RooAddPdf *model = new RooAddPdf("model","model",RooArgList(*pdf_signal,*pdf_bkg),RooArgList(*nSig,*nBkg)); 
    RooDataSet *data = model->generate(*x,nSigInj+nBkgInj);
    
    RooDataHist *roohist = new RooDataHist("roohist","roohist",RooArgList(*x),*data);
    // build fit model
    RooRealVar *nFitSig = new RooRealVar("nFitSig","nFitSig",SIG,0,10*SIG);
    RooRealVar *nFitBkg = new RooRealVar("nFitBkg","nFitBkg",BKG,0,10*BKG);
    RooAddPdf *modelFit = new RooAddPdf("modelFit","modelFit",RooArgList(*pdf_signal,*pdf_bkg),RooArgList(*nFitSig,*nFitBkg)); 
    // fit the pseudo dataset
    RooFitResult *res = modelFit->fitTo(*roohist,RooFit::Save(),RooFit::Extended(kTRUE));
    //res->Print();
    nSigFit = nFitSig->getVal();
    nBkgFit = nFitBkg->getVal();
    eSigFit = nFitSig->getError();
    eBkgFit = nFitBkg->getError();
    nll     = res->minNll();
    tr->Fill();
    if (itoy % 100 == 0) {
      cout<<"Toy #"<<itoy<<": injected = "<<nSigInj<<", fitted = "<<nSigFit<<", error = "<<eSigFit<<endl;
    }
    if (NTOYS == 1) {
      frame = x->frame();
      roohist->plotOn(frame); 
      model->plotOn(frame);
    }
  }
  if (NTOYS == 1) {
    TCanvas *can = new TCanvas("Toy","Toy",900,600);
    frame->Draw();
  }  
  else {
    outf->cd();
    tr->Write();
    outf->Close();
    fTemplates->Close();
  }  
}