rf599_wspacepersist.cxx

Go to the documentation of this file.

//////////////////////////////////////////////////////////////////////////
//
// 'ORGANIZATION AND SIMULTANEOUS FITS' RooFit tutorial macro #501
// 
// Using simultaneous p.d.f.s to describe simultaneous fits to multiple
// datasets
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooWorkspace.h"
#include "RooProdPdf.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooGaussModel.h"
#include "RooAddModel.h"
#include "RooDecay.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic599 : public RooFitTestUnit
{
public: 
  TestBasic599(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooFitTestUnit("Workspace and p.d.f. persistence",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

    if (_write) {
            
      RooWorkspace *w = new RooWorkspace("TestBasic11_ws") ;

      regWS(w,"Basic11_ws") ;

      // Build Gaussian PDF in X
      RooRealVar x("x","x",-10,10) ;
      RooRealVar meanx("meanx","mean of gaussian",-1) ;
      RooRealVar sigmax("sigmax","width of gaussian",3) ;
      RooGaussian gaussx("gaussx","gaussian PDF",x,meanx,sigmax) ;  

      // Build Gaussian PDF in Y
      RooRealVar y("y","y",-10,10) ;
      RooRealVar meany("meany","mean of gaussian",-1) ;
      RooRealVar sigmay("sigmay","width of gaussian",3) ;
      RooGaussian gaussy("gaussy","gaussian PDF",y,meany,sigmay) ;  

      // Make product of X and Y
      RooProdPdf gaussxy("gaussxy","gaussx*gaussy",RooArgSet(gaussx,gaussy)) ;

      // Make flat bkg in X and Y
      RooPolynomial flatx("flatx","flatx",x) ;
      RooPolynomial flaty("flaty","flaty",x) ;
      RooProdPdf flatxy("flatxy","flatx*flaty",RooArgSet(flatx,flaty)) ;

      // Make sum of gaussxy and flatxy
      RooRealVar frac("frac","frac",0.5,0.,1.) ;
      RooAddPdf sumxy("sumxy","sumxy",RooArgList(gaussxy,flatxy),frac) ;

      // Store p.d.f in workspace
      w->import(gaussx) ;
      w->import(gaussxy,RenameConflictNodes("set2")) ;
      w->import(sumxy,RenameConflictNodes("set3")) ;

      // Make reference plot of GaussX
      RooPlot* frame1 = x.frame() ;
      gaussx.plotOn(frame1) ;
      regPlot(frame1,"Basic11_gaussx_framex") ;
      
      // Make reference plots for GaussXY
      RooPlot* frame2 = x.frame() ;
      gaussxy.plotOn(frame2) ;
      regPlot(frame2,"Basic11_gaussxy_framex") ;
      
      RooPlot* frame3 = y.frame() ;
      gaussxy.plotOn(frame3) ;
      regPlot(frame3,"Basic11_gaussxy_framey") ;
      
      // Make reference plots for SumXY
      RooPlot* frame4 = x.frame() ;
      sumxy.plotOn(frame4) ;
      regPlot(frame4,"Basic11_sumxy_framex") ;
      
      RooPlot* frame5 = y.frame() ;
      sumxy.plotOn(frame5) ;
      regPlot(frame5,"Basic11_sumxy_framey") ;

      // Analytically convolved p.d.f.s

      // Build a simple decay PDF
      RooRealVar dt("dt","dt",-20,20) ;
      RooRealVar tau("tau","tau",1.548) ;
      
      // Build a gaussian resolution model
      RooRealVar bias1("bias1","bias1",0) ;
      RooRealVar sigma1("sigma1","sigma1",1) ;
      RooGaussModel gm1("gm1","gauss model 1",dt,bias1,sigma1) ;
      
      // Construct a decay PDF, smeared with single gaussian resolution model
      RooDecay decay_gm1("decay_gm1","decay",dt,tau,gm1,RooDecay::DoubleSided) ;
          
      // Build another gaussian resolution model
      RooRealVar bias2("bias2","bias2",0) ;
      RooRealVar sigma2("sigma2","sigma2",5) ;
      RooGaussModel gm2("gm2","gauss model 2",dt,bias2,sigma2) ;
      
      // Build a composite resolution model
      RooRealVar gm1frac("gm1frac","fraction of gm1",0.5) ;
      RooAddModel gmsum("gmsum","sum of gm1 and gm2",RooArgList(gm1,gm2),gm1frac) ;
    
      // Construct a decay PDF, smeared with double gaussian resolution model
      RooDecay decay_gmsum("decay_gmsum","decay",dt,tau,gmsum,RooDecay::DoubleSided) ;

      w->import(decay_gm1) ;
      w->import(decay_gmsum,RenameConflictNodes("set3")) ;

      RooPlot* frame6 = dt.frame() ;
      decay_gm1.plotOn(frame6) ;
      regPlot(frame6,"Basic11_decay_gm1_framedt") ;
    
      RooPlot* frame7 = dt.frame() ;
      decay_gmsum.plotOn(frame7) ;
      regPlot(frame7,"Basic11_decay_gmsum_framedt") ;

      // Construct simultaneous p.d.f
      RooCategory cat("cat","cat") ;
      cat.defineType("A") ;
      cat.defineType("B") ;
      RooSimultaneous sim("sim","sim",cat) ;
      sim.addPdf(gaussxy,"A") ;
      sim.addPdf(flatxy,"B") ;

      w->import(sim,RenameConflictNodes("set4")) ;

      // Make plot with dummy dataset for index projection
      RooDataHist dh("dh","dh",cat) ;
      cat.setLabel("A") ;
      dh.add(cat) ;
      cat.setLabel("B") ;
      dh.add(cat) ;
      
      RooPlot* frame8 = x.frame() ;
      sim.plotOn(frame8,ProjWData(cat,dh),Project(cat)) ;
      
      regPlot(frame8,"Basic11_sim_framex") ;
      
    
    } else {

      RooWorkspace* w = getWS("Basic11_ws") ;
      if (!w) return kFALSE ;

      // Retrieve p.d.f from workspace
      RooAbsPdf* gaussx = w->pdf("gaussx") ;

      // Make test plot and offer for comparison against ref plot
      RooPlot* frame1 = w->var("x")->frame() ;
      gaussx->plotOn(frame1) ;
      regPlot(frame1,"Basic11_gaussx_framex") ;
      
      // Retrieve p.d.f from workspace
      RooAbsPdf* gaussxy = w->pdf("gaussxy") ;

      // Make test plot and offer for comparison against ref plot
      RooPlot* frame2 = w->var("x")->frame() ;
      gaussxy->plotOn(frame2) ;
      regPlot(frame2,"Basic11_gaussxy_framex") ;

      // Make test plot and offer for comparison against ref plot
      RooPlot* frame3 = w->var("y")->frame() ;
      gaussxy->plotOn(frame3) ;
      regPlot(frame3,"Basic11_gaussxy_framey") ;

      // Retrieve p.d.f from workspace
      RooAbsPdf* sumxy = w->pdf("sumxy") ;

      // Make test plot and offer for comparison against ref plot
      RooPlot* frame4 = w->var("x")->frame() ;
      sumxy->plotOn(frame4) ;
      regPlot(frame4,"Basic11_sumxy_framex") ;

      // Make test plot and offer for comparison against ref plot
      RooPlot* frame5 = w->var("y")->frame() ;
      sumxy->plotOn(frame5) ;
      regPlot(frame5,"Basic11_sumxy_framey") ;

      // Retrieve p.d.f from workspace
      RooAbsPdf* decay_gm1 = w->pdf("decay_gm1") ;

      // Make test plot and offer for comparison against ref plot
      RooPlot* frame6 = w->var("dt")->frame() ;
      decay_gm1->plotOn(frame6) ;
      regPlot(frame6,"Basic11_decay_gm1_framedt") ;

      // Retrieve p.d.f from workspace
      RooAbsPdf* decay_gmsum = w->pdf("decay_gmsum") ;

      // Make test plot and offer for comparison against ref plot
      RooPlot* frame7 = w->var("dt")->frame() ;
      decay_gmsum->plotOn(frame7) ;
      regPlot(frame7,"Basic11_decay_gmsum_framedt") ;

      // Retrieve p.d.f. from workspace
      RooAbsPdf* sim = w->pdf("sim") ;
      RooCategory* cat = w->cat("cat") ;

      // Make plot with dummy dataset for index projection
      RooPlot* frame8 = w->var("x")->frame() ;

      RooDataHist dh("dh","dh",*cat) ;
      cat->setLabel("A") ;
      dh.add(*cat) ;
      cat->setLabel("B") ;
      dh.add(*cat) ;
      
      sim->plotOn(frame8,ProjWData(*cat,dh),Project(*cat)) ;

      regPlot(frame8,"Basic11_sim_framex") ;
      
    }

    // "Workspace persistence"
    return kTRUE ;
  }
} ;