C++ RooDataSet::plotOn方法代码示例

void TestNonCentral(){

  RooWorkspace w("w");
  // k <2, must use sum
  w.factory("NonCentralChiSquare::nc(x[0,50],k[1.99,0,5],lambda[5])");
  // kk > 2 can use bessel
  w.factory("NonCentralChiSquare::ncc(x,kk[2.01,0,5],lambda)");
  // kk > 2, force sum
  w.factory("NonCentralChiSquare::nccc(x,kk,lambda)");
  ((RooNonCentralChiSquare*)w.pdf("nccc"))->SetForceSum(true);

  // a normal "central" chi-square for comparision when lambda->0
  w.factory("ChiSquarePdf::cs(x,k)");

  //w.var("kk")->setVal(4.); // test a large kk

  RooDataSet* ncdata = w.pdf("nc")->generate(*w.var("x"),100);
  RooDataSet* csdata = w.pdf("cs")->generate(*w.var("x"),100);
  RooPlot* plot = w.var("x")->frame();
  ncdata->plotOn(plot,MarkerColor(kRed));
  csdata->plotOn(plot,MarkerColor(kBlue));
  w.pdf("nc")->plotOn(plot,LineColor(kRed));
  w.pdf("ncc")->plotOn(plot,LineColor(kGreen));
  w.pdf("nccc")->plotOn(plot,LineColor(kYellow),LineStyle(kDashed));
  w.pdf("cs")->plotOn(plot,LineColor(kBlue),LineStyle(kDotted));
  plot->Draw();

}

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 rf314_paramfitrange()
{

  // D e f i n e   o b s e r v a b l e s   a n d   d e c a y   p d f 
  // ---------------------------------------------------------------

  // Declare observables
  RooRealVar t("t","t",0,5) ;
  RooRealVar tmin("tmin","tmin",0,0,5) ;

  // Make parameterized range in t : [tmin,5]
  t.setRange(tmin,RooConst(t.getMax())) ;

  // Make pdf
  RooRealVar tau("tau","tau",-1.54,-10,-0.1) ;
  RooExponential model("model","model",t,tau) ;



  // C r e a t e   i n p u t   d a t a 
  // ------------------------------------

  // Generate complete dataset without acceptance cuts (for reference)
  RooDataSet* dall = model.generate(t,10000) ;

  // Generate a (fake) prototype dataset for acceptance limit values
  RooDataSet* tmp = RooGaussian("gmin","gmin",tmin,RooConst(0),RooConst(0.5)).generate(tmin,5000) ;

  // Generate dataset with t values that observe (t>tmin)
  RooDataSet* dacc = model.generate(t,ProtoData(*tmp)) ;



  // F i t   p d f   t o   d a t a   i n   a c c e p t a n c e   r e g i o n
  // -----------------------------------------------------------------------

  RooFitResult* r = model.fitTo(*dacc,Save()) ;


 
  // P l o t   f i t t e d   p d f   o n   f u l l   a n d   a c c e p t e d   d a t a 
  // ---------------------------------------------------------------------------------

  // Make plot frame, add datasets and overlay model
  RooPlot* frame = t.frame(Title("Fit to data with per-event acceptance")) ;
  dall->plotOn(frame,MarkerColor(kRed),LineColor(kRed)) ;
  model.plotOn(frame) ;
  dacc->plotOn(frame) ;

  // Print fit results to demonstrate absence of bias
  r->Print("v") ;


  new TCanvas("rf314_paramranges","rf314_paramranges",600,600) ;
  gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.6) ; frame->Draw() ;

  return ;
}

void plot_toys()
{
  TFile* f = new TFile("/project/atlas/users/tlenz/HWWStatisticsCode_v31/CPmix_kHWW2dot65m_0j_trainv31_v31_allsys_topLumi_toys.root");
  RooWorkspace* w = (RooWorkspace*)f->Get("combined");

  RooDataSet* data = (RooDataSet*)w->data("obsData");
  RooDataSet* toydata = (RooDataSet*)w->data("toyData");
  RooDataSet* toydata2 = (RooDataSet*)w->data("toyData4");
  RooDataSet* thisAsimovData = (RooDataSet*)w->data("asimivData_profiled_for_epsilon_0"); //asimivData_profiled_for_epsilon_0
  
  RooPlot* frame1 = w->var("obs_x_em_signalLike_0j_2012")->frame();
  //  RooPlot* frame2 = w->var("obs_x_em_mainControl_0j_2012")->frame();
  // RooPlot* frame3 = w->var("obs_x_em_zbox_0j_2012")->frame();
  
  data->plotOn(frame1, Name("data"), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_signalLike_0j_2012"));
  //data->plotOn(frame2, Name("data"), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_mainControl_0j_2012"));
  //data->plotOn(frame3, Name("data"), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_zbox_0j_2012"));
  
  toydata->plotOn(frame1, Name("toy1"), LineColor(kRed), MarkerStyle(21), MarkerColor(kRed), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_signalLike_0j_2012"));
  //toydata->plotOn(frame2, Name("toy1"), LineColor(kRed), MarkerStyle(21), MarkerColor(kRed), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_mainControl_0j_2012"));
  //toydata->plotOn(frame3, Name("toy1"), LineColor(kRed), MarkerStyle(21), MarkerColor(kRed), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_zbox_0j_2012"));
  
  toydata2->plotOn(frame1, Name("toy4"), LineColor(kGreen), MarkerStyle(21), MarkerColor(kGreen), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_signalLike_0j_2012"));
  //toydata2->plotOn(frame2, Name("toy4"), LineColor(kGreen), MarkerStyle(21), MarkerColor(kGreen), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_mainControl_0j_2012"));
  //toydata2->plotOn(frame3, Name("toy4"), LineColor(kGreen), MarkerStyle(21), MarkerColor(kGreen), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_zbox_0j_2012"));
  
  thisAsimovData->plotOn(frame1, Name("asimov"), LineColor(kBlue), MarkerStyle(22), MarkerColor(kBlue), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_signalLike_0j_2012"));
  //thisAsimovData->plotOn(frame2, Name("asimov"), LineColor(kBlue), MarkerStyle(22), MarkerColor(kBlue), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_mainControl_0j_2012"));
  //thisAsimovData->plotOn(frame3, Name("asimov"), LineColor(kBlue), MarkerStyle(22), MarkerColor(kBlue), DataError(RooAbsData::Poisson), Cut("merged_cat==merged_cat::em_zbox_0j_2012"));
  
  TCanvas* c = new TCanvas("c","c",1800,600);
  //c->Divide(3,1);
  
  /*c->cd(1);*/ /*gPad->SetLogy();*/ frame1->Draw();
  
  
  //c->cd(2); frame2->Draw();
  //draw legend
  TLegend *leg = new TLegend(0.48, 0.60, 0.84, 0.80, NULL,"brNDC");
  leg->SetBorderSize(0);
  leg->SetTextFont(62);
  leg->SetTextSize(0.04);
  leg->SetLineColor(1);
  leg->SetLineStyle(1);
  leg->SetLineWidth(1);
  leg->SetFillColor(0);
  leg->SetFillStyle(1001);
  leg->SetNColumns(1);
  leg->AddEntry(frame1->findObject("data"),"data","pl");
  leg->AddEntry(frame1->findObject("toy1"),"toy1","pl");
  leg->AddEntry(frame1->findObject("toy4"),"toy4","pl");
  leg->AddEntry(frame1->findObject("asimov"),"asimov","pl");
	leg->Draw();
	//c->cd(3); frame3->Draw();
	c->Print("toys_linear.pdf");
}

void rf304_uncorrprod()
{

   // C r e a t e   c o m p o n e n t   p d f s   i n   x   a n d   y 
   // ----------------------------------------------------------------

   // Create two p.d.f.s gaussx(x,meanx,sigmax) gaussy(y,meany,sigmay) and its variables
   RooRealVar x("x","x",-5,5) ;
   RooRealVar y("y","y",-5,5) ;

   RooRealVar meanx("mean1","mean of gaussian x",2) ;
   RooRealVar meany("mean2","mean of gaussian y",-2) ;
   RooRealVar sigmax("sigmax","width of gaussian x",1) ;
   RooRealVar sigmay("sigmay","width of gaussian y",5) ;

   RooGaussian gaussx("gaussx","gaussian PDF",x,meanx,sigmax) ;  
   RooGaussian gaussy("gaussy","gaussian PDF",y,meany,sigmay) ;  



   // C o n s t r u c t   u n c o r r e l a t e d   p r o d u c t   p d f
   // -------------------------------------------------------------------

   // Multiply gaussx and gaussy into a two-dimensional p.d.f. gaussxy
   RooProdPdf  gaussxy("gaussxy","gaussx*gaussy",RooArgList(gaussx,gaussy)) ;



   // S a m p l e   p d f ,   p l o t   p r o j e c t i o n   o n   x   a n d   y  
   // ---------------------------------------------------------------------------

   // Generate 10000 events in x and y from gaussxy
   RooDataSet *data = gaussxy.generate(RooArgSet(x,y),10000) ;

   // Plot x distribution of data and projection of gaussxy on x = Int(dy) gaussxy(x,y)
   RooPlot* xframe = x.frame(Title("X projection of gauss(x)*gauss(y)")) ;
   data->plotOn(xframe) ;
   gaussxy.plotOn(xframe) ; 

   // Plot x distribution of data and projection of gaussxy on y = Int(dx) gaussxy(x,y)
   RooPlot* yframe = y.frame(Title("Y projection of gauss(x)*gauss(y)")) ;
   data->plotOn(yframe) ;
   gaussxy.plotOn(yframe) ; 



   // Make canvas and draw RooPlots
   TCanvas *c = new TCanvas("rf304_uncorrprod","rf304_uncorrprod",800, 400);
   c->Divide(2);
   c->cd(1) ; gPad->SetLeftMargin(0.15) ; xframe->GetYaxis()->SetTitleOffset(1.4) ; xframe->Draw() ;
   c->cd(2) ; gPad->SetLeftMargin(0.15) ; yframe->GetYaxis()->SetTitleOffset(1.4) ; yframe->Draw() ;

}

void rf301_composition()
{
  // S e t u p   c o m p o s e d   m o d e l   g a u s s ( x , m ( y ) , s )
  // -----------------------------------------------------------------------

  // Create observables
  RooRealVar x("x","x",-5,5) ;
  RooRealVar y("y","y",-5,5) ;

  // Create function f(y) = a0 + a1*y
  RooRealVar a0("a0","a0",-0.5,-5,5) ;
  RooRealVar a1("a1","a1",-0.5,-1,1) ;
  RooPolyVar fy("fy","fy",y,RooArgSet(a0,a1)) ;

  // Creat gauss(x,f(y),s)
  RooRealVar sigma("sigma","width of gaussian",0.5) ;
  RooGaussian model("model","Gaussian with shifting mean",x,fy,sigma) ;  


  // S a m p l e   d a t a ,   p l o t   d a t a   a n d   p d f   o n   x   a n d   y 
  // ---------------------------------------------------------------------------------

  // Generate 10000 events in x and y from model
  RooDataSet *data = model.generate(RooArgSet(x,y),10000) ;

  // Plot x distribution of data and projection of model on x = Int(dy) model(x,y)
  RooPlot* xframe = x.frame() ;
  data->plotOn(xframe) ;
  model.plotOn(xframe) ; 

  // Plot x distribution of data and projection of model on y = Int(dx) model(x,y)
  RooPlot* yframe = y.frame() ;
  data->plotOn(yframe) ;
  model.plotOn(yframe) ; 

  // Make two-dimensional plot in x vs y
  TH1* hh_model = model.createHistogram("hh_model",x,Binning(50),YVar(y,Binning(50))) ;
  hh_model->SetLineColor(kBlue) ;



  // Make canvas and draw RooPlots
  TCanvas *c = new TCanvas("rf301_composition","rf301_composition",1200, 400);
  c->Divide(3);
  c->cd(1) ; gPad->SetLeftMargin(0.15) ; xframe->GetYaxis()->SetTitleOffset(1.4) ; xframe->Draw() ;
  c->cd(2) ; gPad->SetLeftMargin(0.15) ; yframe->GetYaxis()->SetTitleOffset(1.4) ; yframe->Draw() ;
  c->cd(3) ; gPad->SetLeftMargin(0.20) ; hh_model->GetZaxis()->SetTitleOffset(2.5) ; hh_model->Draw("surf") ;

}

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 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 MakeSpinPlots::DrawSpinBackground(TString tag, TString mcName,bool signal){
  bool drawSM = (smName!="" && smName!=mcName);

  TCanvas cv;
  double thisN  = ws->data(mcName+"_Combined")->reduce(TString("evtcat==evtcat::")+tag)->sumEntries();
  float norm = thisN; //607*lumi/12.*thisN/(totEB+totEE);
  cout << norm <<endl;
  if(signal) norm = ws->data(Form("Data_%s_%s_sigWeight",tag.Data(),mcName.Data()))->sumEntries();
  RooPlot *frame = ws->var("cosT")->frame(0,1,5);

  RooDataSet* bkgWeight = (RooDataSet*)ws->data(Form("Data_%s_%s_bkgWeight",tag.Data(),mcName.Data()));
  RooDataSet* tmp = (RooDataSet*)ws->data("Data_Combined")->reduce(TString("((mass>115 && mass<120) || (mass>130 && mass<135)) && evtcat==evtcat::")+tag);
  tmp->plotOn(frame,RooFit::Rescale(norm/tmp->sumEntries()));
  cout << "b" <<endl;
  ws->pdf(Form("%s_FIT_%s_cosTpdf",mcName.Data(),tag.Data()))->plotOn(frame,RooFit::LineColor(kGreen),RooFit::Normalization(norm/tmp->sumEntries()));
  if(drawSM) ws->pdf(Form("%s_FIT_%s_cosTpdf",smName.Data(),tag.Data()))->plotOn(frame,RooFit::LineColor(kRed),RooFit::Normalization(norm/tmp->sumEntries()));
  cout << "c   " <<bkgWeight <<endl;

  bkgWeight->plotOn(frame,RooFit::Rescale(norm/bkgWeight->sumEntries()),RooFit::MarkerColor(kBlue) );  
  if(signal){
    cout << "d" <<endl;
      
    ws->data(Form("Data_%s_%s_sigWeight",tag.Data(),mcName.Data()))->plotOn(frame,RooFit::MarkerStyle(4));
  }
  cout << "d" <<endl;
  
  frame->SetMaximum(frame->GetMaximum()*(signal?0.8:0.4)*norm/tmp->sumEntries());
  frame->SetMinimum(-1*frame->GetMaximum());
  TLegend l(0.6,0.2,0.95,0.45);
  l.SetFillColor(0);
  l.SetBorderSize(0);
  l.SetHeader(tag);
  l.AddEntry(frame->getObject(0),"Data m#in [115,120]#cup[130,135]","p");
  l.AddEntry(frame->getObject(1),mcName,"l");
  if(drawSM) l.AddEntry(frame->getObject(2),"SM Higgs","l");
  l.AddEntry(frame->getObject(2+drawSM),"background weighted Data","p");
  if(signal) l.AddEntry(frame->getObject(3+drawSM),"signal weighted Data","p");

  cout << "e" <<endl;

  frame->Draw();
  l.Draw("SAME");
  cv.SaveAs( basePath+Form("/cosThetaPlots/CosThetaDist_%s%s_%s_%s.png",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data(),tag.Data()) );
  cv.SaveAs( basePath+Form("/cosThetaPlots/C/CosThetaDist_%s%s_%s_%s.C",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data(),tag.Data()) );
  cv.SaveAs( basePath+Form("/cosThetaPlots/CosThetaDist_%s%s_%s_%s.pdf",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data(),tag.Data()) );
}

void test01()
{
  // S e t u p   m o d e l 
  // ---------------------

	TFile *openFile = new TFile("/tmp/kyolee/dimuonTree_upsiMiniTree_pA5tev_14nb_Run210498-210909_trigBit1_allTriggers0_pt4.root");
	TTree* tree = (TTree*)openFile->Get("UpsilonTree");	

  // Declare variables with associated name, title, initial value and allowed range
  RooRealVar mass("invariantMass","M_{#mu#mu} [GeV/c]",9.4,7,14) ;
  RooRealVar muPlusPt("muPlusPt","muPlusPt",0,40) ;
  RooRealVar muMinusPt("muMinusPt","muMinusPt",0,40) ;
  RooRealVar mean("mean","mean of gaussian",9.4,8,11) ;
  RooRealVar sigma("sigma","width of gaussian",1,0.1,10) ;

	RooDataSet* data = new RooDataSet("data", "data", RooArgSet(mass,muPlusPt,muMinusPt), Import(*tree), Cut("muPlusPt>4 && muMinusPt>4"));

 	data->Print();

  // Build gaussian p.d.f in terms of x,mean and sigma
  RooGaussian gauss("gauss","gaussian PDF",mass,mean,sigma) ;  

  // Construct plot frame in 'x'
  RooPlot* xframe = mass.frame(Title("Gaussian p.d.f.")) ;

	data->plotOn(xframe);

	/*
  // P l o t   m o d e l   a n d   c h a n g e   p a r a m e t e r   v a l u e s
  // ---------------------------------------------------------------------------

  // Plot gauss in frame (i.e. in x) 
  gauss.plotOn(xframe) ;

  // Change the value of sigma to 3
  sigma.setVal(3) ;

  // Plot gauss in frame (i.e. in x) and draw frame on canvas
  gauss.plotOn(xframe,LineColor(kRed)) ;
*/  

  // F i t   m o d e l   t o   d a t a
  // -----------------------------

  // Fit pdf to data
  gauss.fitTo(*data,Range(9,10)) ;
  gauss.plotOn(xframe) ;

  // Print values of mean and sigma (that now reflect fitted values and errors)
  mean.Print() ;
  sigma.Print() ;


  // Draw all frames on a canvas
  TCanvas* c = new TCanvas("rf101_basics","rf101_basics",800,400) ;
  c->cd(1) ; gPad->SetLeftMargin(0.15) ; xframe->GetYaxis()->SetTitleOffset(1.6) ; xframe->Draw() ;
 
}

void rf903_numintcache(Int_t mode=0)
{
   // Mode = 0 : Run plain fit (slow)
   // Mode = 1 : Generate workspace with pre-calculated integral and store it on file (prepare for accelerated running)
   // Mode = 2 : Run fit from previously stored workspace including cached integrals (fast, requires run in mode=1 first)

   // C r e a t e ,   s a v e   o r   l o a d   w o r k s p a c e   w i t h   p . d . f .
   // -----------------------------------------------------------------------------------

   // Make/load workspace, exit here in mode 1
   RooWorkspace* w1 = getWorkspace(mode) ;
   if (mode==1) {

      // Show workspace that was created
      w1->Print() ;

      // Show plot of cached integral values
      RooDataHist* hhcache = (RooDataHist*) w1->expensiveObjectCache().getObj(1) ;
      if (hhcache) {

         new TCanvas("rf903_numintcache","rf903_numintcache",600,600) ;
         hhcache->createHistogram("a")->Draw() ;

      }
      else {
         Error("rf903_numintcache","Cached histogram is not existing in workspace");
      }
         return ;
   }

   // U s e   p . d . f .   f r o m   w o r k s p a c e   f o r   g e n e r a t i o n   a n d   f i t t i n g
   // -----------------------------------------------------------------------------------

   // This is always slow (need to find maximum function value empirically in 3D space)
   RooDataSet* d = w1->pdf("model")->generate(RooArgSet(*w1->var("x"),*w1->var("y"),*w1->var("z")),1000) ;

   // This is slow in mode 0, but fast in mode 1
   w1->pdf("model")->fitTo(*d,Verbose(kTRUE),Timer(kTRUE)) ;

   // Projection on x (always slow as 2D integral over Y,Z at fitted value of a is not cached)
   RooPlot* framex = w1->var("x")->frame(Title("Projection of 3D model on X")) ;
   d->plotOn(framex) ;
   w1->pdf("model")->plotOn(framex) ;

   // Draw x projection on canvas
   new TCanvas("rf903_numintcache","rf903_numintcache",600,600) ;
   framex->Draw() ;

   // Make workspace available on command line after macro finishes
   gDirectory->Add(w1) ;

   return ;

}

void rf208_convolution()
{
  // S e t u p   c o m p o n e n t   p d f s 
  // ---------------------------------------

  // Construct observable
  RooRealVar t("t","t",-10,30) ;

  // Construct landau(t,ml,sl) ;
  RooRealVar ml("ml","mean landau",5.,-20,20) ;
  RooRealVar sl("sl","sigma landau",1,0.1,10) ;
  RooLandau landau("lx","lx",t,ml,sl) ;
  
  // Construct gauss(t,mg,sg)
  RooRealVar mg("mg","mg",0) ;
  RooRealVar sg("sg","sg",2,0.1,10) ;
  RooGaussian gauss("gauss","gauss",t,mg,sg) ;


  // C o n s t r u c t   c o n v o l u t i o n   p d f 
  // ---------------------------------------

  // Set #bins to be used for FFT sampling to 10000
  t.setBins(10000,"cache") ; 

  // Construct landau (x) gauss
  RooFFTConvPdf lxg("lxg","landau (X) gauss",t,landau,gauss) ;



  // S a m p l e ,   f i t   a n d   p l o t   c o n v o l u t e d   p d f 
  // ----------------------------------------------------------------------

  // Sample 1000 events in x from gxlx
  RooDataSet* data = lxg.generate(t,10000) ;

  // Fit gxlx to data
  lxg.fitTo(*data) ;

  // Plot data, landau pdf, landau (X) gauss pdf
  RooPlot* frame = t.frame(Title("landau (x) gauss convolution")) ;
  data->plotOn(frame) ;
  lxg.plotOn(frame) ;
  landau.plotOn(frame,LineStyle(kDashed)) ;


  // Draw frame on canvas
  new TCanvas("rf208_convolution","rf208_convolution",600,600) ;
  gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4) ; frame->Draw() ;

}

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 rs602_HLFactoryCombinationexample() {

   using namespace RooStats;
   using namespace RooFit;

   // create a card
   TString card_name("HLFavtoryCombinationexample.rs");
   ofstream ofile(card_name);
   ofile << "// The simplest card for combination\n\n"
         << "gauss1 = Gaussian(x[0,100],mean1[50,0,100],4);\n"
         << "flat1 = Polynomial(x,0);\n"
         << "sb_model1 = SUM(nsig1[120,0,300]*gauss1 , nbkg1[100,0,1000]*flat1);\n"
         << "gauss2 = Gaussian(x,mean2[80,0,100],5);\n"
         << "flat2 = Polynomial(x,0);\n"
         << "sb_model2 = SUM(nsig2[90,0,400]*gauss2 , nbkg2[80,0,1000]*flat2);\n";

   ofile.close();

   HLFactory hlf("HLFavtoryCombinationexample",
               card_name,
               false);

   hlf.AddChannel("model1","sb_model1","flat1");
   hlf.AddChannel("model2","sb_model2","flat2");
   RooAbsPdf* pdf=hlf.GetTotSigBkgPdf();
   RooCategory* thecat = hlf.GetTotCategory();
   RooRealVar* x= static_cast<RooRealVar*>(hlf.GetWs()->arg("x"));

   RooDataSet* data = pdf->generate(RooArgSet(*x,*thecat),Extended());

   // --- Perform extended ML fit of composite PDF to toy data ---
   pdf->fitTo(*data) ;

   // --- Plot toy data and composite PDF overlaid ---
   RooPlot* xframe = x->frame() ;

   data->plotOn(xframe);
   thecat->setIndex(0);
   pdf->plotOn(xframe,Slice(*thecat),ProjWData(*thecat,*data)) ;

   thecat->setIndex(1);
   pdf->plotOn(xframe,Slice(*thecat),ProjWData(*thecat,*data)) ;

   gROOT->SetStyle("Plain");
   xframe->Draw();
}

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();
}