rf705_linearmorph.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_roofit
/// \notebook
///
/// Special p.d.f.'s: linear interpolation between p.d.f shapes using the 'Alex Read' algorithm
///
/// \macro_image
/// \macro_output
/// \macro_code
///
/// \date 07/2008
/// \author Wouter Verkerke
 
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooPolynomial.h"
#include "RooIntegralMorph.h"
#include "RooNLLVar.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "TH1.h"
using namespace RooFit;
 
void rf705_linearmorph()
{
   // C r e a t e   e n d   p o i n t   p d f   s h a p e s
   // ------------------------------------------------------
 
   // Observable
   RooRealVar x("x", "x", -20, 20);
 
   // Lower end point shape: a Gaussian
   RooRealVar g1mean("g1mean", "g1mean", -10);
   RooGaussian g1("g1", "g1", x, g1mean, RooConst(2));
 
   // Upper end point shape: a Polynomial
   RooPolynomial g2("g2", "g2", x, RooArgSet(RooConst(-0.03), RooConst(-0.001)));
 
   // C r e a t e   i n t e r p o l a t i n g   p d f
   // -----------------------------------------------
 
   // Create interpolation variable
   RooRealVar alpha("alpha", "alpha", 0, 1.0);
 
   // Specify sampling density on observable and interpolation variable
   x.setBins(1000, "cache");
   alpha.setBins(50, "cache");
 
   // Construct interpolating pdf in (x,a) represent g1(x) at a=a_min
   // and g2(x) at a=a_max
   RooIntegralMorph lmorph("lmorph", "lmorph", g1, g2, x, alpha);
 
   // P l o t   i n t e r p o l a t i n g   p d f   a t   v a r i o u s   a l p h a
   // -----------------------------------------------------------------------------
 
   // Show end points as blue curves
   RooPlot *frame1 = x.frame();
   g1.plotOn(frame1);
   g2.plotOn(frame1);
 
   // Show interpolated shapes in red
   alpha.setVal(0.125);
   lmorph.plotOn(frame1, LineColor(kRed));
   alpha.setVal(0.25);
   lmorph.plotOn(frame1, LineColor(kRed));
   alpha.setVal(0.375);
   lmorph.plotOn(frame1, LineColor(kRed));
   alpha.setVal(0.50);
   lmorph.plotOn(frame1, LineColor(kRed));
   alpha.setVal(0.625);
   lmorph.plotOn(frame1, LineColor(kRed));
   alpha.setVal(0.75);
   lmorph.plotOn(frame1, LineColor(kRed));
   alpha.setVal(0.875);
   lmorph.plotOn(frame1, LineColor(kRed));
   alpha.setVal(0.95);
   lmorph.plotOn(frame1, LineColor(kRed));
 
   // S h o w   2 D   d i s t r i b u t i o n   o f   p d f ( x , a l p h a )
   // -----------------------------------------------------------------------
 
   // Create 2D histogram
   TH1 *hh = lmorph.createHistogram("hh", x, Binning(40), YVar(alpha, Binning(40)));
   hh->SetLineColor(kBlue);
 
   // F i t   p d f   t o   d a t a s e t   w i t h   a l p h a = 0 . 8
   // -----------------------------------------------------------------
 
   // Generate a toy dataset at alpha = 0.8
   alpha = 0.8;
   RooDataSet *data = lmorph.generate(x, 1000);
 
   // Fit pdf to toy data
   lmorph.setCacheAlpha(kTRUE);
   lmorph.fitTo(*data, Verbose(kTRUE));
 
   // Plot fitted pdf and data overlaid
   RooPlot *frame2 = x.frame(Bins(100));
   data->plotOn(frame2);
   lmorph.plotOn(frame2);
 
   // S c a n   - l o g ( L )   v s   a l p h a
   // -----------------------------------------
 
   // Show scan -log(L) of dataset w.r.t alpha
   RooPlot *frame3 = alpha.frame(Bins(100), Range(0.1, 0.9));
 
   // Make 2D pdf of histogram
   RooNLLVar nll("nll", "nll", lmorph, *data);
   nll.plotOn(frame3, ShiftToZero());
 
   lmorph.setCacheAlpha(kFALSE);
 
   TCanvas *c = new TCanvas("rf705_linearmorph", "rf705_linearmorph", 800, 800);
   c->Divide(2, 2);
   c->cd(1);
   gPad->SetLeftMargin(0.15);
   frame1->GetYaxis()->SetTitleOffset(1.6);
   frame1->Draw();
   c->cd(2);
   gPad->SetLeftMargin(0.20);
   hh->GetZaxis()->SetTitleOffset(2.5);
   hh->Draw("surf");
   c->cd(3);
   gPad->SetLeftMargin(0.15);
   frame3->GetYaxis()->SetTitleOffset(1.4);
   frame3->Draw();
   c->cd(4);
   gPad->SetLeftMargin(0.15);
   frame2->GetYaxis()->SetTitleOffset(1.4);
   frame2->Draw();
 
   return;
}