C++ RooDataSet类代码示例

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

}

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

RooDataHist * genHistFromModelPdf(const char * name, RooAbsPdf *model,  RooRealVar *var,    double ScaleLumi,  int range, int rebin, int seed ) {
  double genEvents =  model->expectedEvents(*var);
  TRandom3 *rndm = new TRandom3();
  rndm->SetSeed(seed);
  double nEvt = rndm->PoissonD( genEvents) ;
  int intEvt = ( (nEvt- (int)nEvt) >= 0.5) ? (int)nEvt +1 : int(nEvt);
  RooDataSet * data = model->generate(*var ,   intEvt   );
  cout<< " expected events for " << name << " = "<< genEvents << endl; 
  cout<< " data->numEntries() for name " << name << " == " << data->numEntries()<< endl;
  // cout<< " nEvt from PoissonD for" << name << " == " << nEvt<< endl;
  //cout<< " cast of nEvt  for" << name << " == " << intEvt<< endl; 
  RooAbsData *binned_data = data->binnedClone();
  TH1 * toy_hist = binned_data->createHistogram( name, *var, Binning(range/rebin )  );
  for(int i = 1; i <= toy_hist->GetNbinsX(); ++i) {
    toy_hist->SetBinError( i,  sqrt( toy_hist->GetBinContent(i)) );
    if(toy_hist->GetBinContent(i) == 0.00) {
      cout<< " WARNING: histo " << name << " has 0 enter in bin number " << i << endl;   
    }
    if(toy_hist->GetBinContent(i) < 0.1) {
      toy_hist->SetBinContent(i, 0.0);
      toy_hist->SetBinError(i, 0.0);
      cout<< " WARNING: setting value 0.0 to histo " << name << " for bin number " << i << endl;   
    }  
  }
  RooDataHist * toy_rooHist = new RooDataHist(name, name , RooArgList(*var), toy_hist );
  return toy_rooHist; 
}

void AnalyzeToy::extract_signal()
{
   calculate_yield();
   MakeSpinSPlot splotter(toyData);
   splotter.addSpecies("signal",ws->pdf("model_signal_mass"),signalYield);
   splotter.addSpecies("background",ws->pdf("model_bkg_mass"),backgroundYield);
   splotter.addVariable(ws->var("mass"));
   splotter.calculate();
   RooDataSet *sweights = splotter.getSWeightDataSet();
   sweights->SetName("sweights");

   RooRealVar weight("weight","",-5.,5.);
   RooArgSet event;
   event.add(*mass);
   event.add(*cosT);
   event.add(weight);

   extractedData = new RooDataSet("extractedData","",event,WeightVar("weight"));

   Long64_t nEntries = toyData->numEntries();
   for(int i=0;i<nEntries;i++)
   {
      double weight_double=0;
      weight_double += sweights->get(i)->getRealValue("signal_sw");
//      weight_double += sweights->get(i)->getRealValue("background_sw");

      mass->setVal(toyData->get(i)->getRealValue("mass"));
      cosT->setVal(toyData->get(i)->getRealValue("cosT"));
      extractedData->add(event,weight_double);
   }
   delete toyData;
}

void trig_pt()
{
  TCanvas *myCan=new TCanvas("myCan","myCan");
  myCan->SetGrid();
  TFile *f_MC= new TFile("TnP_WptCutToTrig_MCptminus.root","read");
  RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTree/Tnp_WptCut_to_Mu15_eta2p1_pt/cnt_eff");
  cout<<"ntry: "<<datasetMC->numEntries()<<endl;
  double XMC[Nbin],XMCerrL[Nbin],XMCerrH[Nbin],YMC[Nbin],YMCerrLo[Nbin],YMCerrHi[Nbin];
  for(int i(0); i<datasetMC->numEntries();i++)
  {
    const RooArgSet &pointMC=*datasetMC->get(i);
    RooRealVar &ptMC=pointMC["pt"],&effMC = pointMC["efficiency"];
    XMC[i]=ptMC.getVal();
    XMCerrL[i]=-ptMC.getAsymErrorLo();
    XMCerrH[i]=ptMC.getAsymErrorHi();
    YMC[i]=effMC.getVal();
    YMCerrLo[i]=-effMC.getAsymErrorLo();
    YMCerrHi[i]=effMC.getAsymErrorHi();
  }
  grMC=new TGraphAsymmErrors(11,XMC,YMC,XMCerrL,XMCerrH,YMCerrLo,YMCerrHi);
  grMC->SetLineColor(kRed);
  grMC->SetMarkerColor(kRed);
  grMC->Draw("AP");
  //grMC->Draw("psame");
  myCan->SaveAs("trig_McMinus_pt.png");
  myCan->SaveAs("trig_McMinus_pt.eps");

}

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 glbToId_eta()
{
  TCanvas *myCan=new TCanvas("myCan","myCan");
  myCan->SetGrid();
  TFile *f_MC= new TFile("TnP_GlbToID_MCetaplus_WptTight2012_eta.root","read");
  RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTree/WptTight2012_eta/fit_eff");
  //RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTree/WptTight2012_eta/cnt_eff");
  cout<<"ntry: "<<datasetMC->numEntries()<<endl;
  double XMC[Nbin],XMCerrL[Nbin],XMCerrH[Nbin],YMC[Nbin],YMCerrLo[Nbin],YMCerrHi[Nbin];
  for(int i(0); i<datasetMC->numEntries();i++)
  {
    const RooArgSet &pointMC=*datasetMC->get(i);
    RooRealVar &etaMC=pointMC["eta"],&effMC = pointMC["efficiency"];
    XMC[i]=etaMC.getVal();
    XMCerrL[i]=-etaMC.getAsymErrorLo();
    XMCerrH[i]=etaMC.getAsymErrorHi();
    YMC[i]=effMC.getVal();
    YMCerrLo[i]=-effMC.getAsymErrorLo();
    YMCerrHi[i]=effMC.getAsymErrorHi();
  }
  grMC=new TGraphAsymmErrors(Nbin,XMC,YMC,XMCerrL,XMCerrH,YMCerrLo,YMCerrHi);
  grMC->SetLineColor(kRed);
  grMC->SetMarkerColor(kRed);
  grMC->Draw("AP");
  //grMC->Draw("psame");
  myCan->SaveAs("glbToId_MCplus_eta.png");
  myCan->SaveAs("glbToId_MCplus_eta.eps");

}

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 makeDataset( TString fname, TString tname, TString outfname ) {

  RooWorkspace *w = new RooWorkspace("w","w");

  w->factory( "Dst_M[1950.,2070.]" );
  w->factory( "D0_M[1810.,1920.]" );
  w->factory( "D0_LTIME_ps[0.00,5.]" );

  RooArgSet *observables = new RooArgSet();
  observables->add( *w->var("Dst_M") );
  observables->add( *w->var("D0_M") );
  observables->add( *w->var("D0_LTIME_ps") );
  w->defineSet("observables", *observables);

  w->var("Dst_M")->setBins(240);
  w->var("D0_M")->setBins(220);
  w->var("D0_LTIME_ps")->setBins(200);

  double Dst_M = -999.;
  double D0_M = -999.;
  double D0_LTIME_ps = -999.;

  TFile *tf = TFile::Open(fname);
  TTree *tree = (TTree*)tf->Get(tname);

  tree->SetBranchAddress( "Dst_M", &Dst_M );
  tree->SetBranchAddress( "D0_M" , &D0_M );
  tree->SetBranchAddress( "D0_LTIME_ps", &D0_LTIME_ps );

  RooDataSet *data = new RooDataSet("Data","Data",*observables);
  RooDataHist *dataH = new RooDataHist("DataHist","Data",*observables);

  for ( int ev=0; ev<tree->GetEntries(); ev++) {

    tree->GetEntry(ev);

    if ( ev%10000 == 0 ) cout << ev << " / " << tree->GetEntries() << endl;

    if ( Dst_M < w->var("Dst_M")->getMin() || Dst_M > w->var("Dst_M")->getMax() ) continue;
    if ( D0_M < w->var("D0_M")->getMin() || D0_M > w->var("D0_M")->getMax() ) continue;
    if ( D0_LTIME_ps < w->var("D0_LTIME_ps")->getMin() || D0_LTIME_ps > w->var("D0_LTIME_ps")->getMax() ) continue;

    w->var("Dst_M")->setVal(Dst_M);
    w->var("D0_M")->setVal(D0_M);
    w->var("D0_LTIME_ps")->setVal(D0_LTIME_ps);

    data->add( *observables );
    dataH->add( *observables );

  }
  tf->Close();

  w->import(*data);
  w->import(*dataH);

  w->writeToFile(outfname);

}

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 reweight_diphotonpt_1d(TH1F*weight_diphopt,TH1F*weight_diphopto,TH2F*weight_diphoptn,TH2F*weight_diphopt2o,TH2F*weight_diphoptnr, TH2F* weight_diphopt2,RooDataSet **dset, RooDataSet *dsetdestination){

  if (!(*dset)) return;
  if (!(dsetdestination)) return;
///numerator and denominator
  TH1F *hnum = new TH1F("hnum","hnum",n_diphoptbins_forreweighting,diphoptbins_forreweighting);
  TH1F *hden = new TH1F("hden","hden",n_diphoptbins_forreweighting,diphoptbins_forreweighting);
  hnum->Sumw2();
  hden->Sumw2();

     
// RooAbsData->get*() Load a given row of data

 //getRealValue Get value of a RooAbsReal stored in set with given name. If none is found, value of defVal is returned.
 for (int i=0; i<(*dset)->numEntries(); i++){
      hden->Fill(fabs((*dset)->get(i)->getRealValue("roodiphopt")),(*dset)->store()->weight(i));
  }
  for (int i=0; i<dsetdestination->numEntries(); i++){
    hnum->Fill(fabs(dsetdestination->get(i)->getRealValue("roodiphopt")),dsetdestination->store()->weight(i));
  }
 //normalize to one 
  hnum->Scale(1.0/hnum->Integral());
  hden->Scale(1.0/hden->Integral());
  hnum->Divide(hden);
  TH1F *h = hnum;
  RooDataSet *newdset = new RooDataSet(**dset,Form("%s_diphoptrew",(*dset)->GetName()));
  newdset->reset();
  assert(newdset);
  for (int i=0; i<(*dset)->numEntries(); i++){
    RooArgSet args = *((*dset)->get(i));
    float oldw = (*dset)->store()->weight(i);
    float diphopt = fabs((*dset)->get(i)->getRealValue("roodiphopt"));
    float neww = oldw*h->GetBinContent(h->FindBin(fabs(diphopt)));
    newdset->add(args,neww);
    if(debug){
		weight_diphopt->Fill(neww);
		weight_diphopto->Fill(oldw);
		weight_diphopt2o->Fill((h->FindBin(fabs(diphopt))),oldw);	
		weight_diphoptn->Fill((h->FindBin(fabs(diphopt))),neww);	
		if(oldw!=0)weight_diphoptnr->Fill(h->FindBin(fabs(diphopt)),oldw/neww);
		else {weight_diphoptnr->Fill(-10,1);}//cout << "dipho weight old 0" << endl;}
		if(oldw!=0)weight_diphopt2->Fill(diphopt,oldw/neww);
		else {weight_diphopt2->Fill(-10,1);}//cout << "dipho weight old 0" << endl;}
	 }
		}

  newdset->SetName((*dset)->GetName());
  newdset->SetTitle((*dset)->GetTitle());
  delete hnum; delete hden;
 
  RooDataSet *old_dset = *dset;
  *dset=newdset;
  std::cout << "Diphoton Pt 1d rew: norm from " << old_dset->sumEntries() << " to " << newdset->sumEntries() << std::endl;

  delete old_dset;

};

void rf602_chi2fit()
{

   // S e t u p   m o d e l
   // ---------------------

   // Declare observable x
   RooRealVar x("x","x",0,10) ;

   // Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
   RooRealVar mean("mean","mean of gaussians",5) ;
   RooRealVar sigma1("sigma1","width of gaussians",0.5) ;
   RooRealVar sigma2("sigma2","width of gaussians",1) ;

   RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ;  
   RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ;  
   
   // Build Chebychev polynomial p.d.f.  
   RooRealVar a0("a0","a0",0.5,0.,1.) ;
   RooRealVar a1("a1","a1",0.2,0.,1.) ;
   RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ;

   // Sum the signal components into a composite signal p.d.f.
   RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
   RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;

   // Sum the composite signal and background 
   RooRealVar bkgfrac("bkgfrac","fraction of background",0.5,0.,1.) ;
   RooAddPdf  model("model","g1+g2+a",RooArgList(bkg,sig),bkgfrac) ;


   // C r e a t e   b i n n e d   d a t a s e t
   // -----------------------------------------

   RooDataSet* d = model.generate(x,10000) ;
   RooDataHist* dh = d->binnedClone() ;

   // Construct a chi^2 of the data and the model.
   // When a p.d.f. is used in a chi^2 fit, the probability density scaled
   // by the number of events in the dataset to obtain the fit function
   // If model is an extended p.d.f, the expected number events is used
   // instead of the observed number of events.
   model.chi2FitTo(*dh) ;

   // NB: It is also possible to fit a RooAbsReal function to a RooDataHist
   // using chi2FitTo(). 

   // Note that entries with zero bins are _not_ allowed 
   // for a proper chi^2 calculation and will give error
   // messages  
   RooDataSet* dsmall = (RooDataSet*) d->reduce(EventRange(1,100)) ;
   RooDataHist* dhsmall = dsmall->binnedClone() ;
   RooChi2Var chi2_lowstat("chi2_lowstat","chi2",model,*dhsmall) ;
   cout << chi2_lowstat.getVal() << endl ;


}

pair<double,double> datEvents(RooWorkspace *work, int m_hyp, int cat, bool spin=false){
  
  vector<double> result;
  RooDataSet *data = (RooDataSet*)work->data(Form("data_mass_cat%d",cat));
  double evs = data->numEntries();
  double evsPerGev;
  if (!spin) evsPerGev = data->sumEntries(Form("CMS_hgg_mass>=%4.1f && CMS_hgg_mass<%4.1f",double(m_hyp)-0.5,double(m_hyp)+0.5));
  else evsPerGev = data->sumEntries(Form("mass>=%4.1f && mass<%4.1f",double(m_hyp)-0.5,double(m_hyp)+0.5));
  return pair<double,double>(evs,evsPerGev);
}

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 ;

}