//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #101
// 
// Fitting, plotting, toy data generation on one-dimensional p.d.f
//
// pdf = gauss(x,m,s) 
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooUnitTest.h"

using namespace RooFit ;


// Elementary operations on a gaussian PDF
class TestBasic101 : public RooUnitTest
{
public: 
  TestBasic101(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Fitting,plotting & event generation of basic p.d.f",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

    // S e t u p   m o d e l 
    // ---------------------
    
    // Declare variables x,mean,sigma with associated name, title, initial value and allowed range
    RooRealVar x("x","x",-10,10) ;
    RooRealVar mean("mean","mean of gaussian",1,-10,10) ;
    RooRealVar sigma("sigma","width of gaussian",1,0.1,10) ;
    
    // Build gaussian p.d.f in terms of x,mean and sigma
    RooGaussian gauss("gauss","gaussian PDF",x,mean,sigma) ;  
    
    // Construct plot frame in 'x'
    RooPlot* xframe = x.frame(Title("Gaussian p.d.f.")) ;
    
    
    // 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),Name("another")) ;
    
    
    // G e n e r a t e   e v e n t s 
    // -----------------------------
    
    // Generate a dataset of 1000 events in x from gauss
    RooDataSet* data = gauss.generate(x,10000) ;  
    
    // Make a second plot frame in x and draw both the 
    // data and the p.d.f in the frame
    RooPlot* xframe2 = x.frame(Title("Gaussian p.d.f. with data")) ;
    data->plotOn(xframe2) ;
    gauss.plotOn(xframe2) ;
    
    
    // F i t   m o d e l   t o   d a t a
    // -----------------------------
    
    // Fit pdf to data
    gauss.fitTo(*data) ;
        
    
    // --- Post processing for stressRooFit ---
    regPlot(xframe ,"rf101_plot1") ;
    regPlot(xframe2,"rf101_plot2") ;
    
    delete data ;
    
    return kTRUE ;
  }
} ;


//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #102
// 
// Importing data from ROOT TTrees and THx histograms
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////


#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TTree.h"
#include "TH1D.h"
#include "TRandom.h"
using namespace RooFit ;


class TestBasic102 : public RooUnitTest
{
public: 
  TestBasic102(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Data import methods",refFile,writeRef,verbose) {} ;

  TH1* makeTH1() 
  {
    // Create ROOT TH1 filled with a Gaussian distribution
    
    TH1D* hh = new TH1D("hh","hh",25,-10,10) ;
    for (int i=0 ; i<100 ; i++) {
      hh->Fill(gRandom->Gaus(0,3)) ;
    }
    return hh ;
  }
  
  
  TTree* makeTTree() 
  {
    // Create ROOT TTree filled with a Gaussian distribution in x and a uniform distribution in y
    
    TTree* tree = new TTree("tree","tree") ;
    Double_t* px = new Double_t ;
    Double_t* py = new Double_t ;
    tree->Branch("x",px,"x/D") ;
    tree->Branch("y",py,"y/D") ;
    for (int i=0 ; i<100 ; i++) {
      *px = gRandom->Gaus(0,3) ;
      *py = gRandom->Uniform()*30 - 15 ;
      tree->Fill() ;
    }

    //delete px ;
    //delete py ;

    return tree ;
  }
  
  Bool_t testCode() {
    
    ////////////////////////////////////////////////////////
    // I m p o r t i n g   R O O T   h i s t o g r a m s  //
    ////////////////////////////////////////////////////////
    
    // I m p o r t   T H 1   i n t o   a   R o o D a t a H i s t
    // ---------------------------------------------------------
    
    // Create a ROOT TH1 histogram
    TH1* hh = makeTH1() ;
    
    // Declare observable x
    RooRealVar x("x","x",-10,10) ;
    
    // Create a binned dataset that imports contents of TH1 and associates its contents to observable 'x'
    RooDataHist dh("dh","dh",x,Import(*hh)) ;
    
    
    // P l o t   a n d   f i t   a   R o o D a t a H i s t
    // ---------------------------------------------------
    
    // Make plot of binned dataset showing Poisson error bars (RooFit default)
    RooPlot* frame = x.frame(Title("Imported TH1 with Poisson error bars")) ;
    dh.plotOn(frame) ; 
    
    // Fit a Gaussian p.d.f to the data
    RooRealVar mean("mean","mean",0,-10,10) ;
    RooRealVar sigma("sigma","sigma",3,0.1,10) ;
    RooGaussian gauss("gauss","gauss",x,mean,sigma) ;
    gauss.fitTo(dh) ;
    gauss.plotOn(frame) ;
    
    // P l o t   a n d   f i t   a   R o o D a t a H i s t   w i t h   i n t e r n a l   e r r o r s
    // ---------------------------------------------------------------------------------------------
    
    // If histogram has custom error (i.e. its contents is does not originate from a Poisson process
    // but e.g. is a sum of weighted events) you can data with symmetric 'sum-of-weights' error instead
    // (same error bars as shown by ROOT)
    RooPlot* frame2 = x.frame(Title("Imported TH1 with internal errors")) ;
    dh.plotOn(frame2,DataError(RooAbsData::SumW2)) ; 
    gauss.plotOn(frame2) ;
    
    // Please note that error bars shown (Poisson or SumW2) are for visualization only, the are NOT used
    // in a maximum likelihood fit
    //
    // A (binned) ML fit will ALWAYS assume the Poisson error interpretation of data (the mathematical definition 
    // of likelihood does not take any external definition of errors). Data with non-unit weights can only be correctly
    // fitted with a chi^2 fit (see rf602_chi2fit.C) 
    
    
    ////////////////////////////////////////////////
    // I m p o r t i n g   R O O T  T T r e e s   //
    ////////////////////////////////////////////////
    
    
    // I m p o r t   T T r e e   i n t o   a   R o o D a t a S e t
    // -----------------------------------------------------------
    
    TTree* tree = makeTTree() ;
    
    // Define 2nd observable y
    RooRealVar y("y","y",-10,10) ;
    
    // Construct unbinned dataset importing tree branches x and y matching between branches and RooRealVars 
    // is done by name of the branch/RRV 
    // 
    // Note that ONLY entries for which x,y have values within their allowed ranges as defined in 
    // RooRealVar x and y are imported. Since the y values in the import tree are in the range [-15,15]
    // and RRV y defines a range [-10,10] this means that the RooDataSet below will have less entries than the TTree 'tree'
    
    RooDataSet ds("ds","ds",RooArgSet(x,y),Import(*tree)) ;
    
    
    // P l o t   d a t a s e t   w i t h   m u l t i p l e   b i n n i n g   c h o i c e s
    // ------------------------------------------------------------------------------------
    
    // Print unbinned dataset with default frame binning (100 bins)
    RooPlot* frame3 = y.frame(Title("Unbinned data shown in default frame binning")) ;
    ds.plotOn(frame3) ;
    
    // Print unbinned dataset with custom binning choice (20 bins)
    RooPlot* frame4 = y.frame(Title("Unbinned data shown with custom binning")) ;
    ds.plotOn(frame4,Binning(20)) ;
    
    // Draw all frames on a canvas
    regPlot(frame ,"rf102_plot1") ;
    regPlot(frame2,"rf102_plot2") ;
    regPlot(frame3,"rf102_plot3") ;
    regPlot(frame4,"rf102_plot4") ;
    
    delete hh ;
    delete tree ;

    return kTRUE ;
  }
} ;







//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #103
// 
// Interpreted functions and p.d.f.s
//
// 
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooConstVar.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "RooGenericPdf.h"

using namespace RooFit ;


class TestBasic103 : public RooUnitTest
{
public: 
  TestBasic103(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Interpreted expression p.d.f.",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {
    
    /////////////////////////////////////////////////////////
    // G e n e r i c   i n t e r p r e t e d   p . d . f . //
    /////////////////////////////////////////////////////////
    
    // Declare observable x
    RooRealVar x("x","x",-20,20) ;
    
    // C o n s t r u c t   g e n e r i c   p d f   f r o m   i n t e r p r e t e d   e x p r e s s i o n
    // -------------------------------------------------------------------------------------------------
    
    // To construct a proper p.d.f, the formula expression is explicitly normalized internally by dividing 
    // it by a numeric integral of the expresssion over x in the range [-20,20] 
    //
    RooRealVar alpha("alpha","alpha",5,0.1,10) ;
    RooGenericPdf genpdf("genpdf","genpdf","(1+0.1*abs(x)+sin(sqrt(abs(x*alpha+0.1))))",RooArgSet(x,alpha)) ;
    
    
    // S a m p l e ,   f i t   a n d   p l o t   g e n e r i c   p d f
    // ---------------------------------------------------------------
    
    // Generate a toy dataset from the interpreted p.d.f
    RooDataSet* data = genpdf.generate(x,10000) ;
    
    // Fit the interpreted p.d.f to the generated data
    genpdf.fitTo(*data) ;
    
    // Make a plot of the data and the p.d.f overlaid
    RooPlot* xframe = x.frame(Title("Interpreted expression pdf")) ;
    data->plotOn(xframe) ;
    genpdf.plotOn(xframe) ;  
    
    
    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    // S t a n d a r d   p . d . f   a d j u s t   w i t h   i n t e r p r e t e d   h e l p e r   f u n c t i o n //
    //                                                                                                             //
    // Make a gauss(x,sqrt(mean2),sigma) from a standard RooGaussian                                               //
    //                                                                                                             //
    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    
    
    // C o n s t r u c t   s t a n d a r d   p d f  w i t h   f o r m u l a   r e p l a c i n g   p a r a m e t e r
    // ------------------------------------------------------------------------------------------------------------
    
    // Construct parameter mean2 and sigma
    RooRealVar mean2("mean2","mean^2",10,0,200) ;
    RooRealVar sigma("sigma","sigma",3,0.1,10) ;
    
    // Construct interpreted function mean = sqrt(mean^2)
    RooFormulaVar mean("mean","mean","sqrt(mean2)",mean2) ;
    
    // Construct a gaussian g2(x,sqrt(mean2),sigma) ;
    RooGaussian g2("g2","h2",x,mean,sigma) ;
    
    
    // G e n e r a t e   t o y   d a t a 
    // ---------------------------------
    
    // Construct a separate gaussian g1(x,10,3) to generate a toy Gaussian dataset with mean 10 and width 3
    RooGaussian g1("g1","g1",x,RooConst(10),RooConst(3)) ;
    RooDataSet* data2 = g1.generate(x,1000) ;
    
    
    // F i t   a n d   p l o t   t a i l o r e d   s t a n d a r d   p d f 
    // -------------------------------------------------------------------
    
    // Fit g2 to data from g1
    RooFitResult* r = g2.fitTo(*data2,Save()) ;
    
    // Plot data on frame and overlay projection of g2
    RooPlot* xframe2 = x.frame(Title("Tailored Gaussian pdf")) ;
    data2->plotOn(xframe2) ;
    g2.plotOn(xframe2) ;
    
    regPlot(xframe,"rf103_plot1") ;
    regPlot(xframe2,"rf103_plot2") ;
    regResult(r,"rf103_fit1") ;

    delete data ;
    delete data2 ;

    return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #105
// 
//  Demonstration of binding ROOT Math functions as RooFit functions
//  and pdfs
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TMath.h"
#include "TF1.h"
#include "Math/DistFunc.h"
#include "RooCFunction1Binding.h" 
#include "RooCFunction3Binding.h" 
#include "RooTFnBinding.h" 

using namespace RooFit ;

class TestBasic105 : public RooUnitTest
{
public: 
  TestBasic105(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("C++ function binding operator p.d.f",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

    // B i n d   T M a t h : : E r f   C   f u n c t i o n
    // ---------------------------------------------------
    
    // Bind one-dimensional TMath::Erf function as RooAbsReal function
    RooRealVar x("x","x",-3,3) ;
    RooAbsReal* erf = bindFunction("erf",TMath::Erf,x) ;
    
    // Plot erf on frame 
    RooPlot* frame1 = x.frame(Title("TMath::Erf bound as RooFit function")) ;
    erf->plotOn(frame1) ;
    
    
    // B i n d   R O O T : : M a t h : : b e t a _ p d f   C   f u n c t i o n
    // -----------------------------------------------------------------------
    
    // Bind pdf ROOT::Math::Beta with three variables as RooAbsPdf function
    RooRealVar x2("x2","x2",0,0.999) ;
    RooRealVar a("a","a",5,0,10) ;
    RooRealVar b("b","b",2,0,10) ;
    RooAbsPdf* beta = bindPdf("beta",ROOT::Math::beta_pdf,x2,a,b) ;
    
    // Generate some events and fit
    RooDataSet* data = beta->generate(x2,10000) ;
    beta->fitTo(*data) ;
    
    // Plot data and pdf on frame
    RooPlot* frame2 = x2.frame(Title("ROOT::Math::Beta bound as RooFit pdf")) ;
    data->plotOn(frame2) ;
    beta->plotOn(frame2) ;
    
    
    
    // B i n d   R O O T   T F 1   a s   R o o F i t   f u n c t i o n
    // ---------------------------------------------------------------
    
    // Create a ROOT TF1 function
    TF1 *fa1 = new TF1("fa1","sin(x)/x",0,10);
    
    // Create an observable 
    RooRealVar x3("x3","x3",0.01,20) ;
    
    // Create binding of TF1 object to above observable
    RooAbsReal* rfa1 = bindFunction(fa1,x3) ;
    
    // Make plot frame in observable, plot TF1 binding function
    RooPlot* frame3 = x3.frame(Title("TF1 bound as RooFit function")) ;
    rfa1->plotOn(frame3) ;
      
      
    regPlot(frame1,"rf105_plot1") ;
    regPlot(frame2,"rf105_plot2") ;
    regPlot(frame3,"rf105_plot3") ;
    
    delete erf ;
    delete beta ;
    delete fa1 ;
    delete rfa1 ;
    delete data ;
    
    return kTRUE ;
  }  
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #108
// 
// Plotting unbinned data with alternate and variable binnings
//
// 
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussModel.h"
#include "RooDecay.h"
#include "RooBMixDecay.h"
#include "RooCategory.h"
#include "RooBinning.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;


class TestBasic108 : public RooUnitTest
{
public: 
  TestBasic108(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Non-standard binning in counting and asymmetry plots",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

    // S e t u p   m o d e l 
    // ---------------------
    
    // Build a B decay p.d.f with mixing
    RooRealVar dt("dt","dt",-20,20) ;
    RooRealVar dm("dm","dm",0.472) ;
    RooRealVar tau("tau","tau",1.547) ;
    RooRealVar w("w","mistag rate",0.1) ;
    RooRealVar dw("dw","delta mistag rate",0.) ;
    
    RooCategory mixState("mixState","B0/B0bar mixing state") ;
    mixState.defineType("mixed",-1) ;
    mixState.defineType("unmixed",1) ;
    RooCategory tagFlav("tagFlav","Flavour of the tagged B0") ;
    tagFlav.defineType("B0",1) ;
    tagFlav.defineType("B0bar",-1) ;
    
    // Build a gaussian resolution model
    RooRealVar dterr("dterr","dterr",0.1,1.0) ;
    RooRealVar bias1("bias1","bias1",0) ;
    RooRealVar sigma1("sigma1","sigma1",0.1) ;  
    RooGaussModel gm1("gm1","gauss model 1",dt,bias1,sigma1) ;
    
    // Construct Bdecay (x) gauss
    RooBMixDecay bmix("bmix","decay",dt,mixState,tagFlav,tau,dm,w,dw,gm1,RooBMixDecay::DoubleSided) ;
    
    
    // S a m p l e   d a t a   f r o m   m o d e l
    // --------------------------------------------
    
    // Sample 2000 events in (dt,mixState,tagFlav) from bmix
    RooDataSet *data = bmix.generate(RooArgSet(dt,mixState,tagFlav),2000) ;
    
    
    
    // S h o w   d t   d i s t r i b u t i o n   w i t h   c u s t o m   b i n n i n g
    // -------------------------------------------------------------------------------
    
    // Make plot of dt distribution of data in range (-15,15) with fine binning for dt>0 and coarse binning for dt<0
    
    // Create binning object with range (-15,15)
    RooBinning tbins(-15,15) ;
    
    // Add 60 bins with uniform spacing in range (-15,0)
    tbins.addUniform(60,-15,0) ;
    
    // Add 15 bins with uniform spacing in range (0,15)
    tbins.addUniform(15,0,15) ;
    
    // Make plot with specified binning
    RooPlot* dtframe = dt.frame(Range(-15,15),Title("dt distribution with custom binning")) ;
    data->plotOn(dtframe,Binning(tbins)) ;
    bmix.plotOn(dtframe) ;
    
    // NB: Note that bin density for each bin is adjusted to that of default frame binning as shown
    // in Y axis label (100 bins --> Events/0.4*Xaxis-dim) so that all bins represent a consistent density distribution
    
    
    // S h o w   m i x s t a t e   a s y m m e t r y  w i t h   c u s t o m   b i n n i n g
    // ------------------------------------------------------------------------------------
    
    // Make plot of dt distribution of data asymmetry in 'mixState' with variable binning 
    
    // Create binning object with range (-10,10)
    RooBinning abins(-10,10) ;
    
    // Add boundaries at 0, (-1,1), (-2,2), (-3,3), (-4,4) and (-6,6)
    abins.addBoundary(0) ;
    abins.addBoundaryPair(1) ;
    abins.addBoundaryPair(2) ;
    abins.addBoundaryPair(3) ;
    abins.addBoundaryPair(4) ;
    abins.addBoundaryPair(6) ;
    
    // Create plot frame in dt
    RooPlot* aframe = dt.frame(Range(-10,10),Title("mixState asymmetry distribution with custom binning")) ;
    
    // Plot mixState asymmetry of data with specified customg binning
    data->plotOn(aframe,Asymmetry(mixState),Binning(abins)) ;
    
    // Plot corresponding property of p.d.f
    bmix.plotOn(aframe,Asymmetry(mixState)) ;
    
    // Adjust vertical range of plot to sensible values for an asymmetry
    aframe->SetMinimum(-1.1) ;
    aframe->SetMaximum(1.1) ;
    
    // NB: For asymmetry distributions no density corrects are needed (and are thus not applied)
    

    regPlot(dtframe,"rf108_plot1") ;
    regPlot(aframe,"rf108_plot2") ;

    delete data ;

    return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #109
// 
// Calculating chi^2 from histograms and curves in RooPlots, 
// making histogram of residual and pull distributions
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooHist.h"
using namespace RooFit ;

class TestBasic109 : public RooUnitTest
{
public: 
  TestBasic109(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Calculation of chi^2 and residuals in plots",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {
    
    // S e t u p   m o d e l 
    // ---------------------
    
    // Create observables
    RooRealVar x("x","x",-10,10) ;
    
    // Create Gaussian
    RooRealVar sigma("sigma","sigma",3,0.1,10) ;
    RooRealVar mean("mean","mean",0,-10,10) ;
    RooGaussian gauss("gauss","gauss",x,RooConst(0),sigma) ;
    
    // Generate a sample of 1000 events with sigma=3
    RooDataSet* data = gauss.generate(x,10000) ;
    
    // Change sigma to 3.15
    sigma=3.15 ;
    
    
    // P l o t   d a t a   a n d   s l i g h t l y   d i s t o r t e d   m o d e l
    // ---------------------------------------------------------------------------
    
    // Overlay projection of gauss with sigma=3.15 on data with sigma=3.0
    RooPlot* frame1 = x.frame(Title("Data with distorted Gaussian pdf"),Bins(40)) ;
    data->plotOn(frame1,DataError(RooAbsData::SumW2)) ;
    gauss.plotOn(frame1) ;
    
    
    // C a l c u l a t e   c h i ^ 2 
    // ------------------------------
    
    // Show the chi^2 of the curve w.r.t. the histogram
    // If multiple curves or datasets live in the frame you can specify
    // the name of the relevant curve and/or dataset in chiSquare()
    regValue(frame1->chiSquare(),"rf109_chi2") ;
    
    
    // S h o w   r e s i d u a l   a n d   p u l l   d i s t s
    // -------------------------------------------------------
    
    // Construct a histogram with the residuals of the data w.r.t. the curve
    RooHist* hresid = frame1->residHist() ;
    
    // Construct a histogram with the pulls of the data w.r.t the curve
    RooHist* hpull = frame1->pullHist() ;
    
    // Create a new frame to draw the residual distribution and add the distribution to the frame
    RooPlot* frame2 = x.frame(Title("Residual Distribution")) ;
    frame2->addPlotable(hresid,"P") ;
    
    // Create a new frame to draw the pull distribution and add the distribution to the frame
    RooPlot* frame3 = x.frame(Title("Pull Distribution")) ;
    frame3->addPlotable(hpull,"P") ;
    
    regPlot(frame1,"rf109_plot1") ;
    regPlot(frame2,"rf109_plot2") ;
    regPlot(frame3,"rf109_plot3") ;

    delete data ;
    //delete hresid ;
    //delete hpull ;
    
    return kTRUE ;    
  }
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #110
// 
// Examples on normalization of p.d.f.s,
// integration of p.d.fs, construction
// of cumulative distribution functions from p.d.f.s
// in one dimension
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooAbsReal.h"
#include "RooPlot.h"
#include "TCanvas.h"
using namespace RooFit ;

class TestBasic110 : public RooUnitTest
{
public: 
  TestBasic110(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Normalization of p.d.f.s in 1D",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

    // S e t u p   m o d e l 
    // ---------------------
    
    // Create observables x,y
    RooRealVar x("x","x",-10,10) ;
    
    // Create p.d.f. gaussx(x,-2,3) 
    RooGaussian gx("gx","gx",x,RooConst(-2),RooConst(3)) ;
    
    
    // R e t r i e v e   r a w  &   n o r m a l i z e d   v a l u e s   o f   R o o F i t   p . d . f . s
    // --------------------------------------------------------------------------------------------------
    
    // Return 'raw' unnormalized value of gx
    regValue(gx.getVal(),"rf110_gx") ;
    
    // Return value of gx normalized over x in range [-10,10]
    RooArgSet nset(x) ;

    regValue(gx.getVal(&nset),"rf110_gx_Norm[x]") ;
    
    // Create object representing integral over gx
    // which is used to calculate  gx_Norm[x] == gx / gx_Int[x]
    RooAbsReal* igx = gx.createIntegral(x) ;
    regValue(igx->getVal(),"rf110_gx_Int[x]") ;
    
    
    // I n t e g r a t e   n o r m a l i z e d   p d f   o v e r   s u b r a n g e
    // ----------------------------------------------------------------------------
    
    // Define a range named "signal" in x from -5,5
    x.setRange("signal",-5,5) ;
    
    // Create an integral of gx_Norm[x] over x in range "signal"
    // This is the fraction of of p.d.f. gx_Norm[x] which is in the
    // range named "signal"
    RooAbsReal* igx_sig = gx.createIntegral(x,NormSet(x),Range("signal")) ;
    regValue(igx_sig->getVal(),"rf110_gx_Int[x|signal]_Norm[x]") ;
    
    
    
    // C o n s t r u c t   c u m u l a t i v e   d i s t r i b u t i o n   f u n c t i o n   f r o m   p d f
    // -----------------------------------------------------------------------------------------------------
    
    // Create the cumulative distribution function of gx
    // i.e. calculate Int[-10,x] gx(x') dx'
    RooAbsReal* gx_cdf = gx.createCdf(x) ;
    
    // Plot cdf of gx versus x
    RooPlot* frame = x.frame(Title("c.d.f of Gaussian p.d.f")) ;
    gx_cdf->plotOn(frame) ;
    
    
    regPlot(frame,"rf110_plot1") ;

    delete igx ;
    delete igx_sig ;
    delete gx_cdf ;

    return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #111 
// 
// Configuration and customization of how numeric (partial) integrals
// are executed
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooNumIntConfig.h"
#include "RooLandau.h"
#include "RooArgSet.h"
#include <iomanip>
using namespace RooFit ;

class TestBasic111 : public RooUnitTest
{
public: 
  TestBasic111(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Numeric integration configuration",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

    // A d j u s t   g l o b a l   1 D   i n t e g r a t i o n   p r e c i s i o n 
    // ----------------------------------------------------------------------------
    
    // Example: Change global precision for 1D integrals from 1e-7 to 1e-6
    //
    // The relative epsilon (change as fraction of current best integral estimate) and
    // absolute epsilon (absolute change w.r.t last best integral estimate) can be specified
    // separately. For most p.d.f integrals the relative change criterium is the most important,
    // however for certain non-p.d.f functions that integrate out to zero a separate absolute
    // change criterium is necessary to declare convergence of the integral
    //
    // NB: This change is for illustration only. In general the precision should be at least 1e-7 
    // for normalization integrals for MINUIT to succeed.
    //
    RooAbsReal::defaultIntegratorConfig()->setEpsAbs(1e-6) ;
    RooAbsReal::defaultIntegratorConfig()->setEpsRel(1e-6) ;
    
    
    // N u m e r i c   i n t e g r a t i o n   o f   l a n d a u   p d f 
    // ------------------------------------------------------------------
    
    // Construct p.d.f without support for analytical integrator for demonstration purposes
    RooRealVar x("x","x",-10,10) ;
    RooLandau landau("landau","landau",x,RooConst(0),RooConst(0.1)) ;
    
    
    // Calculate integral over landau with default choice of numeric integrator
    RooAbsReal* intLandau = landau.createIntegral(x) ;
    Double_t val = intLandau->getVal() ;
    regValue(val,"rf111_val1") ;
    

    // S a m e   w i t h   c u s t o m   c o n f i g u r a t i o n
    // -----------------------------------------------------------
    
    
    // Construct a custom configuration which uses the adaptive Gauss-Kronrod technique
    // for closed 1D integrals
    RooNumIntConfig customConfig(*RooAbsReal::defaultIntegratorConfig()) ;
    customConfig.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
    
    
    // Calculate integral over landau with custom integral specification
    RooAbsReal* intLandau2 = landau.createIntegral(x,NumIntConfig(customConfig)) ;
    Double_t val2 = intLandau2->getVal() ;
    regValue(val2,"rf111_val2") ;


    // A d j u s t i n g   d e f a u l t   c o n f i g   f o r   a   s p e c i f i c   p d f 
    // -------------------------------------------------------------------------------------
    
    
    // Another possibility: associate custom numeric integration configuration as default for object 'landau'
    landau.setIntegratorConfig(customConfig) ;
    
    
    // Calculate integral over landau custom numeric integrator specified as object default
    RooAbsReal* intLandau3 = landau.createIntegral(x) ;
    Double_t val3 = intLandau3->getVal() ;
    regValue(val3,"rf111_val3") ;    

    
    delete intLandau ;
    delete intLandau2 ;
    delete intLandau3 ;
    
    return kTRUE ;
  }
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #201
// 
// Composite p.d.f with signal and background component
//
// pdf = f_bkg * bkg(x,a0,a1) + (1-fbkg) * (f_sig1 * sig1(x,m,s1 + (1-f_sig1) * sig2(x,m,s2)))
//
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic201 : public RooUnitTest
{
public: 
  TestBasic201(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Addition operator p.d.f.",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

    // S e t u p   c o m p o n e n t   p d f s 
    // ---------------------------------------
    
    // 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)) ;
    
    
    ////////////////////////////////////////////////////
    // M E T H O D   1 - T w o   R o o A d d P d f s  //
    ////////////////////////////////////////////////////
    
    
    // A d d   s i g n a l   c o m p o n e n t s 
    // ------------------------------------------
    
    // 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) ;
    
    
    // A d d  s i g n a l   a n d   b a c k g r o u n d
    // ------------------------------------------------
    
    // 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) ;
    
    
    // S a m p l e ,   f i t   a n d   p l o t   m o d e l
    // ---------------------------------------------------
    
    // Generate a data sample of 1000 events in x from model
    RooDataSet *data = model.generate(x,1000) ;
    
    // Fit model to data
    model.fitTo(*data) ;
    
    // Plot data and PDF overlaid
    RooPlot* xframe = x.frame(Title("Example of composite pdf=(sig1+sig2)+bkg")) ;
    data->plotOn(xframe) ;
    model.plotOn(xframe) ;
    
    // Overlay the background component of model with a dashed line
    model.plotOn(xframe,Components(bkg),LineStyle(kDashed)) ;
    
    // Overlay the background+sig2 components of model with a dotted line
    model.plotOn(xframe,Components(RooArgSet(bkg,sig2)),LineStyle(kDotted)) ;
    

    ////////////////////////////////////////////////////////////////////////////////////////////////////
    // M E T H O D   2 - O n e   R o o A d d P d f   w i t h   r e c u r s i v e   f r a c t i o n s  //
    ////////////////////////////////////////////////////////////////////////////////////////////////////
    
    // Construct sum of models on one go using recursive fraction interpretations
    //
    //   model2 = bkg + (sig1 + sig2)
    //
    RooAddPdf  model2("model","g1+g2+a",RooArgList(bkg,sig1,sig2),RooArgList(bkgfrac,sig1frac),kTRUE) ;    
    
    // NB: Each coefficient is interpreted as the fraction of the 
    // left-hand component of the i-th recursive sum, i.e.
    //
    //   sum4 = A + ( B + ( C + D)  with fraction fA, fB and fC expands to
    //
    //   sum4 = fA*A + (1-fA)*(fB*B + (1-fB)*(fC*C + (1-fC)*D))
    
    
    // P l o t   r e c u r s i v e   a d d i t i o n   m o d e l
    // ---------------------------------------------------------
    model2.plotOn(xframe,LineColor(kRed),LineStyle(kDashed)) ;
    model2.plotOn(xframe,Components(RooArgSet(bkg,sig2)),LineColor(kRed),LineStyle(kDashed)) ;

  
    regPlot(xframe,"rf201_plot1") ;
    
    delete data ;
    return kTRUE ;

  }

} ;

//////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #202
// 
// Setting up an extended maximum likelihood fit
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooExtendPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic202 : public RooUnitTest
{
public: 
  TestBasic202(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Extended ML fits to addition operator p.d.f.s",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

    // S e t u p   c o m p o n e n t   p d f s 
    // ---------------------------------------
    
    // 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) ;
    
    /////////////////////
    // M E T H O D   1 //
    /////////////////////
    
    
    // C o n s t r u c t   e x t e n d e d   c o m p o s i t e   m o d e l 
    // -------------------------------------------------------------------
    
    // Sum the composite signal and background into an extended pdf nsig*sig+nbkg*bkg
    RooRealVar nsig("nsig","number of signal events",500,0.,10000) ;
    RooRealVar nbkg("nbkg","number of background events",500,0,10000) ;
    RooAddPdf  model("model","(g1+g2)+a",RooArgList(bkg,sig),RooArgList(nbkg,nsig)) ;
    
    
    
    // S a m p l e ,   f i t   a n d   p l o t   e x t e n d e d   m o d e l 
    // ---------------------------------------------------------------------
    
    // Generate a data sample of expected number events in x from model
    // = model.expectedEvents() = nsig+nbkg
    RooDataSet *data = model.generate(x) ;
    
    // Fit model to data, extended ML term automatically included
    model.fitTo(*data) ; 
    
    // Plot data and PDF overlaid, use expected number of events for p.d.f projection normalization
    // rather than observed number of events (==data->numEntries())
    RooPlot* xframe = x.frame(Title("extended ML fit example")) ;
    data->plotOn(xframe) ;
    model.plotOn(xframe,Normalization(1.0,RooAbsReal::RelativeExpected)) ;
    
    // Overlay the background component of model with a dashed line
    model.plotOn(xframe,Components(bkg),LineStyle(kDashed),Normalization(1.0,RooAbsReal::RelativeExpected)) ;
    
    // Overlay the background+sig2 components of model with a dotted line
    model.plotOn(xframe,Components(RooArgSet(bkg,sig2)),LineStyle(kDotted),Normalization(1.0,RooAbsReal::RelativeExpected)) ;
    
    
    /////////////////////
    // M E T H O D   2 //
    /////////////////////
    
    // C o n s t r u c t   e x t e n d e d   c o m p o n e n t s   f i r s t
    // ---------------------------------------------------------------------
    
    // Associated nsig/nbkg as expected number of events with sig/bkg
    RooExtendPdf esig("esig","extended signal p.d.f",sig,nsig) ;
    RooExtendPdf ebkg("ebkg","extended background p.d.f",bkg,nbkg) ;
    
    
    // S u m   e x t e n d e d   c o m p o n e n t s   w i t h o u t   c o e f s 
    // -------------------------------------------------------------------------
    
    // Construct sum of two extended p.d.f. (no coefficients required)
    RooAddPdf  model2("model2","(g1+g2)+a",RooArgList(ebkg,esig)) ;
    
    
    regPlot(xframe,"rf202_plot1") ;

    delete data ;
    return kTRUE ;

  }
  
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #203
// 
// Fitting and plotting in sub ranges
//
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;


class TestBasic203 : public RooUnitTest
{
public: 
  TestBasic203(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Basic fitting and plotting in ranges",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

    // S e t u p   m o d e l 
    // ---------------------
    
    // Construct observables x
    RooRealVar x("x","x",-10,10) ;
    
    // Construct gaussx(x,mx,1)
    RooRealVar mx("mx","mx",0,-10,10) ;
    RooGaussian gx("gx","gx",x,mx,RooConst(1)) ;
    
    // Construct px = 1 (flat in x)
    RooPolynomial px("px","px",x) ;
    
    // Construct model = f*gx + (1-f)px
    RooRealVar f("f","f",0.,1.) ;
    RooAddPdf model("model","model",RooArgList(gx,px),f) ;
    
    // Generated 10000 events in (x,y) from p.d.f. model
    RooDataSet* modelData = model.generate(x,10000) ;
    
    // F i t   f u l l   r a n g e 
    // ---------------------------
    
    // Fit p.d.f to all data
    RooFitResult* r_full = model.fitTo(*modelData,Save(kTRUE)) ;
    
    
    // F i t   p a r t i a l   r a n g e 
    // ----------------------------------
    
    // Define "signal" range in x as [-3,3]
    x.setRange("signal",-3,3) ;  
    
    // Fit p.d.f only to data in "signal" range
    RooFitResult* r_sig = model.fitTo(*modelData,Save(kTRUE),Range("signal")) ;
    
    
    // P l o t   /   p r i n t   r e s u l t s 
    // ---------------------------------------
    
    // Make plot frame in x and add data and fitted model
    RooPlot* frame = x.frame(Title("Fitting a sub range")) ;
    modelData->plotOn(frame) ;
    model.plotOn(frame,Range("Full"),LineStyle(kDashed),LineColor(kRed)) ; // Add shape in full ranged dashed
    model.plotOn(frame) ; // By default only fitted range is shown
    
    regPlot(frame,"rf203_plot") ;
    regResult(r_full,"rf203_r_full") ;
    regResult(r_sig,"rf203_r_sig") ;

    delete modelData ;    
    return kTRUE;
  }  
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #204
// 
// Extended maximum likelihood fit with alternate range definition
// for observed number of events.
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooExtendPdf.h"
#include "RooFitResult.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic204 : public RooUnitTest
{
public: 
  TestBasic204(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Extended ML fit in sub range",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {
    
    // S e t u p   c o m p o n e n t   p d f s 
    // ---------------------------------------
    
    // 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) ;
    
    
    // C o n s t r u c t   e x t e n d e d   c o m p s   wi t h   r a n g e   s p e c
    // ------------------------------------------------------------------------------
    
    // Define signal range in which events counts are to be defined
    x.setRange("signalRange",4,6) ;
    
    // Associated nsig/nbkg as expected number of events with sig/bkg _in_the_range_ "signalRange"
    RooRealVar nsig("nsig","number of signal events in signalRange",500,0.,10000) ;
    RooRealVar nbkg("nbkg","number of background events in signalRange",500,0,10000) ;
    RooExtendPdf esig("esig","extended signal p.d.f",sig,nsig,"signalRange") ;
    RooExtendPdf ebkg("ebkg","extended background p.d.f",bkg,nbkg,"signalRange") ;
    
    
    // S u m   e x t e n d e d   c o m p o n e n t s
    // ---------------------------------------------
    
    // Construct sum of two extended p.d.f. (no coefficients required)
    RooAddPdf  model("model","(g1+g2)+a",RooArgList(ebkg,esig)) ;
    
    
    // S a m p l e   d a t a ,   f i t   m o d e l
    // -------------------------------------------
    
    // Generate 1000 events from model so that nsig,nbkg come out to numbers <<500 in fit
    RooDataSet *data = model.generate(x,1000) ;
    
    
    // Perform unbinned extended ML fit to data
    RooFitResult* r = model.fitTo(*data,Extended(kTRUE),Save()) ;


    regResult(r,"rf204_result") ;

    delete data ;
    return kTRUE ;
  } 
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #205
// 
// Options for plotting components of composite p.d.f.s.
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooAddPdf.h"
#include "RooChebychev.h"
#include "RooExponential.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic205 : public RooUnitTest
{
public: 
  TestBasic205(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Component plotting variations",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

    // S e t u p   c o m p o s i t e    p d f
    // --------------------------------------
    
    // 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) ;  
    
    // 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) ;
    
    // Build Chebychev polynomial p.d.f.  
    RooRealVar a0("a0","a0",0.5,0.,1.) ;
    RooRealVar a1("a1","a1",-0.2,0.,1.) ;
    RooChebychev bkg1("bkg1","Background 1",x,RooArgSet(a0,a1)) ;
    
    // Build expontential pdf
    RooRealVar alpha("alpha","alpha",-1) ;
    RooExponential bkg2("bkg2","Background 2",x,alpha) ;
    
    // Sum the background components into a composite background p.d.f.
    RooRealVar bkg1frac("sig1frac","fraction of component 1 in background",0.2,0.,1.) ;
    RooAddPdf bkg("bkg","Signal",RooArgList(bkg1,bkg2),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) ;
    
    
    
    // S e t u p   b a s i c   p l o t   w i t h   d a t a   a n d   f u l l   p d f 
    // ------------------------------------------------------------------------------
    
    // Generate a data sample of 1000 events in x from model
    RooDataSet *data = model.generate(x,1000) ;
    
    // Plot data and complete PDF overlaid
    RooPlot* xframe  = x.frame(Title("Component plotting of pdf=(sig1+sig2)+(bkg1+bkg2)")) ;
    data->plotOn(xframe) ;
    model.plotOn(xframe) ;
    
    // Clone xframe for use below
    RooPlot* xframe2 = x.frame(Title("Component plotting of pdf=(sig1+sig2)+(bkg1+bkg2)"),Name("xframe2")) ;
    
    
    // M a k e   c o m p o n e n t   b y   o b j e c t   r e f e r e n c e 
    // --------------------------------------------------------------------
    
    // Plot single background component specified by object reference
    model.plotOn(xframe,Components(bkg),LineColor(kRed)) ;
    
    // Plot single background component specified by object reference
    model.plotOn(xframe,Components(bkg2),LineStyle(kDashed),LineColor(kRed)) ;
    
    // Plot multiple background components specified by object reference
    // Note that specified components may occur at any level in object tree
    // (e.g bkg is component of 'model' and 'sig2' is component 'sig')
    model.plotOn(xframe,Components(RooArgSet(bkg,sig2)),LineStyle(kDotted)) ;
    
    
    
    // M a k e   c o m p o n e n t   b y   n a m e  /   r e g e x p  
    // ------------------------------------------------------------
    
    // Plot single background component specified by name
    model.plotOn(xframe2,Components("bkg"),LineColor(kCyan)) ;
    
    // Plot multiple background components specified by name
    model.plotOn(xframe2,Components("bkg1,sig2"),LineStyle(kDotted),LineColor(kCyan)) ;
    
    // Plot multiple background components specified by regular expression on name
    model.plotOn(xframe2,Components("sig*"),LineStyle(kDashed),LineColor(kCyan)) ;
    
    // Plot multiple background components specified by multiple regular expressions on name
    model.plotOn(xframe2,Components("bkg1,sig*"),LineStyle(kDashed),LineColor(kYellow),Invisible()) ;
    
    
    regPlot(xframe,"rf205_plot1") ;
    regPlot(xframe2,"rf205_plot2") ;

    delete data ;
    return kTRUE ;

  }

} ;
/////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #208
// 
// One-dimensional numeric convolution
// (require ROOT to be compiled with --enable-fftw3)
// 
// pdf = landau(t) (x) gauss(t)
// 
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooLandau.h"
#include "RooFFTConvPdf.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TPluginManager.h"
#include "TROOT.h"

using namespace RooFit ;



class TestBasic208 : public RooUnitTest
{
public: 
  TestBasic208(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("FFT Convolution operator p.d.f.",refFile,writeRef,verbose) {} ;

  Bool_t isTestAvailable() { 

    TPluginHandler *h;
    if ((h = gROOT->GetPluginManager()->FindHandler("TVirtualFFT"))) {
      if (h->LoadPlugin() == -1) {
	gROOT->ProcessLine("new TNamed ;") ;
	return kFALSE;
      } else {
	return kTRUE ;
      }
    }
    return kFALSE ;
  }

  Double_t ctol() { return 1e-2 ; } // Account for difficult shape of Landau distribution

  Bool_t testCode() {

    // 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)) ;
    
    regPlot(frame,"rf208_plot1") ;

    delete data ;
    return kTRUE ;

  }
} ;



/////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #209
// 
// Decay function p.d.fs with optional B physics 
// effects (mixing and CP violation) that can be
// analytically convolved with e.g. Gaussian resolution 
// functions
// 
// pdf1 = decay(t,tau) (x) delta(t)
// pdf2 = decay(t,tau) (x) gauss(t,m,s)
// pdf3 = decay(t,tau) (x) (f*gauss1(t,m1,s1) + (1-f)*gauss2(t,m1,s1))
// 
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussModel.h"
#include "RooAddModel.h"
#include "RooTruthModel.h"
#include "RooDecay.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;



class TestBasic209 : public RooUnitTest
{
public: 
  TestBasic209(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Analytical convolution operator",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {
    
    // B - p h y s i c s   p d f   w i t h   t r u t h   r e s o l u t i o n
    // ---------------------------------------------------------------------
    
    // Variables of decay p.d.f.
    RooRealVar dt("dt","dt",-10,10) ;
    RooRealVar tau("tau","tau",1.548) ;
    
    // Build a truth resolution model (delta function)
    RooTruthModel tm("tm","truth model",dt) ;
    
    // Construct decay(t) (x) delta(t)
    RooDecay decay_tm("decay_tm","decay",dt,tau,tm,RooDecay::DoubleSided) ;
    
    // Plot p.d.f. (dashed)
    RooPlot* frame = dt.frame(Title("Bdecay (x) resolution")) ;
    decay_tm.plotOn(frame,LineStyle(kDashed)) ;
    
    
    // B - p h y s i c s   p d f   w i t h   G a u s s i a n   r e s o l u t i o n
    // ----------------------------------------------------------------------------
    
    // 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 decay(t) (x) gauss1(t)
    RooDecay decay_gm1("decay_gm1","decay",dt,tau,gm1,RooDecay::DoubleSided) ;
    
    // Plot p.d.f. 
    decay_gm1.plotOn(frame) ;
    
    
    // B - p h y s i c s   p d f   w i t h   d o u b l e   G a u s s i a n   r e s o l u t i o n
    // ------------------------------------------------------------------------------------------
    
    // 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 f*gm1+(1-f)*gm2
    RooRealVar gm1frac("gm1frac","fraction of gm1",0.5) ;
    RooAddModel gmsum("gmsum","sum of gm1 and gm2",RooArgList(gm1,gm2),gm1frac) ;
    
    // Construct decay(t) (x) (f*gm1 + (1-f)*gm2)
    RooDecay decay_gmsum("decay_gmsum","decay",dt,tau,gmsum,RooDecay::DoubleSided) ;
    
    // Plot p.d.f. (red)
    decay_gmsum.plotOn(frame,LineColor(kRed)) ;
    
    regPlot(frame,"rf209_plot1") ;

    return kTRUE ;

  } 
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #301
// 
// Multi-dimensional p.d.f.s through composition, e.g. substituting a 
// p.d.f parameter with a function that depends on other observables
// 
// pdf = gauss(x,f(y),s) with f(y) = a0 + a1*y
// 
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolyVar.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;



class TestBasic301 : public RooUnitTest
{
public: 
  TestBasic301(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Composition extension of basic p.d.f",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

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


  regPlot(xframe,"rf301_plot1") ;
  regPlot(yframe,"rf302_plot2") ;
  regTH(hh_model,"rf302_model2d") ;

  delete data ;
  return kTRUE ;
  }
} ;


//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #302
// 
//  Utility functions classes available for use in tailoring
//  of composite (multidimensional) pdfs
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooFormulaVar.h"
#include "RooAddition.h"
#include "RooProduct.h"
#include "RooPolyVar.h"
#include "TCanvas.h"
#include "TH1.h"

using namespace RooFit ;


class TestBasic302 : public RooUnitTest
{
public: 
  TestBasic302(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Sum and product utility functions",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   o b s e r v a b l e s ,   p a r a m e t e r s 
  // -----------------------------------------------------------

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

  // Create parameters
  RooRealVar a0("a0","a0",-1.5,-5,5) ;
  RooRealVar a1("a1","a1",-0.5,-1,1) ;
  RooRealVar sigma("sigma","width of gaussian",0.5) ;


  // U s i n g   R o o F o r m u l a V a r   t o   t a i l o r   p d f 
  // -----------------------------------------------------------------------

  // Create interpreted function f(y) = a0 - a1*sqrt(10*abs(y))
  RooFormulaVar fy_1("fy_1","a0-a1*sqrt(10*abs(y))",RooArgSet(y,a0,a1)) ;

  // Create gauss(x,f(y),s)
  RooGaussian model_1("model_1","Gaussian with shifting mean",x,fy_1,sigma) ;  



  // U s i n g   R o o P o l y V a r   t o   t a i l o r   p d f
  // -----------------------------------------------------------------------

  // Create polynomial function f(y) = a0 + a1*y
  RooPolyVar fy_2("fy_2","fy_2",y,RooArgSet(a0,a1)) ;

  // Create gauss(x,f(y),s)
  RooGaussian model_2("model_2","Gaussian with shifting mean",x,fy_2,sigma) ;  



  // U s i n g   R o o A d d i t i o n   t o   t a i l o r   p d f 
  // -----------------------------------------------------------------------

  // Create sum function f(y) = a0 + y
  RooAddition fy_3("fy_3","a0+y",RooArgSet(a0,y)) ;

  // Create gauss(x,f(y),s)
  RooGaussian model_3("model_3","Gaussian with shifting mean",x,fy_3,sigma) ;  



  // U s i n g   R o o P r o d u c t   t o   t a i l o r   p d f 
  // -----------------------------------------------------------------------

  // Create product function f(y) = a1*y
  RooProduct fy_4("fy_4","a1*y",RooArgSet(a1,y)) ;

  // Create gauss(x,f(y),s)
  RooGaussian model_4("model_4","Gaussian with shifting mean",x,fy_4,sigma) ;  



  // P l o t   a l l   p d f s 
  // ----------------------------

  // Make two-dimensional plots in x vs y
  TH1* hh_model_1 = model_1.createHistogram("hh_model_1",x,Binning(50),YVar(y,Binning(50))) ;
  TH1* hh_model_2 = model_2.createHistogram("hh_model_2",x,Binning(50),YVar(y,Binning(50))) ;
  TH1* hh_model_3 = model_3.createHistogram("hh_model_3",x,Binning(50),YVar(y,Binning(50))) ;
  TH1* hh_model_4 = model_4.createHistogram("hh_model_4",x,Binning(50),YVar(y,Binning(50))) ;
  hh_model_1->SetLineColor(kBlue) ;
  hh_model_2->SetLineColor(kBlue) ;
  hh_model_3->SetLineColor(kBlue) ;
  hh_model_4->SetLineColor(kBlue) ;

  regTH(hh_model_1,"rf202_model2d_1") ;
  regTH(hh_model_2,"rf202_model2d_2") ;
  regTH(hh_model_3,"rf202_model2d_3") ;
  regTH(hh_model_4,"rf202_model2d_4") ;

  return kTRUE ;
  }
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #303
// 
// Use of tailored p.d.f as conditional p.d.fs.s
// 
// pdf = gauss(x,f(y),sx | y ) with f(y) = a0 + a1*y
// 
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolyVar.h"
#include "RooProdPdf.h"
#include "RooPlot.h"
#include "TRandom.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;


class TestBasic303 : public RooUnitTest
{
public: 

RooDataSet* makeFakeDataXY() 
{
  RooRealVar x("x","x",-10,10) ;
  RooRealVar y("y","y",-10,10) ;
  RooArgSet coord(x,y) ;

  RooDataSet* d = new RooDataSet("d","d",RooArgSet(x,y)) ;

  for (int i=0 ; i<10000 ; i++) {
    Double_t tmpy = gRandom->Gaus(0,10) ;
    Double_t tmpx = gRandom->Gaus(0.5*tmpy,1) ;
    if (fabs(tmpy)<10 && fabs(tmpx)<10) {
      x = tmpx ;
      y = tmpy ;
      d->add(coord) ;
    }
      
  }

  return d ;
}



  TestBasic303(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Conditional use of F(x|y)",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // 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",-10,10) ;
  RooRealVar y("y","y",-10,10) ;

  // 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,0.1,2.0) ;
  RooGaussian model("model","Gaussian with shifting mean",x,fy,sigma) ;  


  // Obtain fake external experimental dataset with values for x and y
  RooDataSet* expDataXY = makeFakeDataXY() ;



  // G e n e r a t e   d a t a   f r o m   c o n d i t i o n a l   p . d . f   m o d e l ( x | y )  
  // ---------------------------------------------------------------------------------------------

  // Make subset of experimental data with only y values
  RooDataSet* expDataY= (RooDataSet*) expDataXY->reduce(y) ;

  // Generate 10000 events in x obtained from _conditional_ model(x|y) with y values taken from experimental data
  RooDataSet *data = model.generate(x,ProtoData(*expDataY)) ;



  // F i t   c o n d i t i o n a l   p . d . f   m o d e l ( x | y )   t o   d a t a
  // ---------------------------------------------------------------------------------------------

  model.fitTo(*expDataXY,ConditionalObservables(y)) ;
  


  // P r o j e c t   c o n d i t i o n a l   p . d . f   o n   x   a n d   y   d i m e n s i o n s
  // ---------------------------------------------------------------------------------------------

  // Plot x distribution of data and projection of model on x = 1/Ndata sum(data(y_i)) model(x;y_i)
  RooPlot* xframe = x.frame() ;
  expDataXY->plotOn(xframe) ;
  model.plotOn(xframe,ProjWData(*expDataY)) ; 


  // Speed up (and approximate) projection by using binned clone of data for projection
  RooAbsData* binnedDataY = expDataY->binnedClone() ;
  model.plotOn(xframe,ProjWData(*binnedDataY),LineColor(kCyan),LineStyle(kDotted),Name("Alt1")) ;


  // Show effect of projection with too coarse binning
  ((RooRealVar*)expDataY->get()->find("y"))->setBins(5) ;
  RooAbsData* binnedDataY2 = expDataY->binnedClone() ;
  model.plotOn(xframe,ProjWData(*binnedDataY2),LineColor(kRed),Name("Alt2")) ;


  regPlot(xframe,"rf303_plot1") ;

  delete binnedDataY ;
  delete binnedDataY2 ;
  delete expDataXY ;
  delete expDataY ;
  delete data ;

  return kTRUE ;
  }
} ;




/////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #304
// 
// Simple uncorrelated multi-dimensional p.d.f.s
//
// pdf = gauss(x,mx,sx) * gauss(y,my,sy) 
//
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;



class TestBasic304 : public RooUnitTest
{
public: 
  TestBasic304(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Product operator p.d.f. with uncorrelated terms",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

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

  regPlot(xframe,"rf304_plot1") ;
  regPlot(yframe,"rf304_plot2") ;

  delete data ;

  return kTRUE ;
  }
} ;



/////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #305
// 
// Multi-dimensional p.d.f.s with conditional p.d.fs in product
// 
// pdf = gauss(x,f(y),sx | y ) * gauss(y,ms,sx)    with f(y) = a0 + a1*y
// 
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolyVar.h"
#include "RooProdPdf.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;



class TestBasic305 : public RooUnitTest
{
public: 
  TestBasic305(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Product operator p.d.f. with conditional term",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   c o n d i t i o n a l   p d f   g x ( x | y ) 
  // -----------------------------------------------------------

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

  // Create gaussx(x,f(y),sx)
  RooRealVar sigmax("sigma","width of gaussian",0.5) ;
  RooGaussian gaussx("gaussx","Gaussian in x with shifting mean in y",x,fy,sigmax) ;  



  // C r e a t e   p d f   g y ( y ) 
  // -----------------------------------------------------------

  // Create gaussy(y,0,5)
  RooGaussian gaussy("gaussy","Gaussian in y",y,RooConst(0),RooConst(3)) ;



  // C r e a t e   p r o d u c t   g x ( x | y ) * g y ( y )
  // -------------------------------------------------------

  // Create gaussx(x,sx|y) * gaussy(y)
  RooProdPdf model("model","gaussx(x|y)*gaussy(y)",gaussy,Conditional(gaussx,x)) ;



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

  // Generate 1000 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_rf305",x,Binning(50),YVar(y,Binning(50))) ;
  hh_model->SetLineColor(kBlue) ;

  regTH(hh_model,"rf305_model2d") ;
  regPlot(xframe,"rf305_plot1") ;
  regPlot(yframe,"rf305_plot2") ;

  delete data ;

  return kTRUE ;
  
  }
} ;



//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #306
// 
// Complete example with use of conditional p.d.f. with per-event errors
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooGaussModel.h"
#include "RooDecay.h"
#include "RooLandau.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH2D.h"
using namespace RooFit ;



class TestBasic306 : public RooUnitTest
{
public: 
  TestBasic306(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Conditional use of per-event error p.d.f. F(t|dt)",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // B - p h y s i c s   p d f   w i t h   p e r - e v e n t  G a u s s i a n   r e s o l u t i o n
  // ----------------------------------------------------------------------------------------------

  // Observables
  RooRealVar dt("dt","dt",-10,10) ;
  RooRealVar dterr("dterr","per-event error on dt",0.01,10) ;

  // Build a gaussian resolution model scaled by the per-event error = gauss(dt,bias,sigma*dterr)
  RooRealVar bias("bias","bias",0,-10,10) ;
  RooRealVar sigma("sigma","per-event error scale factor",1,0.1,10) ;
  RooGaussModel gm("gm1","gauss model scaled bt per-event error",dt,bias,sigma,dterr) ;

  // Construct decay(dt) (x) gauss1(dt|dterr)
  RooRealVar tau("tau","tau",1.548) ;
  RooDecay decay_gm("decay_gm","decay",dt,tau,gm,RooDecay::DoubleSided) ;



  // C o n s t r u c t   f a k e   ' e x t e r n a l '   d a t a    w i t h   p e r - e v e n t   e r r o r
  // ------------------------------------------------------------------------------------------------------

  // Use landau p.d.f to get somewhat realistic distribution with long tail
  RooLandau pdfDtErr("pdfDtErr","pdfDtErr",dterr,RooConst(1),RooConst(0.25)) ;
  RooDataSet* expDataDterr = pdfDtErr.generate(dterr,10000) ;



  // S a m p l e   d a t a   f r o m   c o n d i t i o n a l   d e c a y _ g m ( d t | d t e r r )
  // ---------------------------------------------------------------------------------------------

  // Specify external dataset with dterr values to use decay_dm as conditional p.d.f.
  RooDataSet* data = decay_gm.generate(dt,ProtoData(*expDataDterr)) ;

  

  // F i t   c o n d i t i o n a l   d e c a y _ d m ( d t | d t e r r )
  // ---------------------------------------------------------------------

  // Specify dterr as conditional observable
  decay_gm.fitTo(*data,ConditionalObservables(dterr)) ;


  
  // P l o t   c o n d i t i o n a l   d e c a y _ d m ( d t | d t e r r )
  // ---------------------------------------------------------------------


  // Make two-dimensional plot of conditional p.d.f in (dt,dterr)
  TH1* hh_decay = decay_gm.createHistogram("hh_decay",dt,Binning(50),YVar(dterr,Binning(50))) ;
  hh_decay->SetLineColor(kBlue) ;


  // Plot decay_gm(dt|dterr) at various values of dterr
  RooPlot* frame = dt.frame(Title("Slices of decay(dt|dterr) at various dterr")) ;
  for (Int_t ibin=0 ; ibin<100 ; ibin+=20) {
    dterr.setBin(ibin) ;
    decay_gm.plotOn(frame,Normalization(5.),Name(Form("curve_slice_%d",ibin))) ;
  }


  // Make projection of data an dt
  RooPlot* frame2 = dt.frame(Title("Projection of decay(dt|dterr) on dt")) ;
  data->plotOn(frame2) ;

  // Make projection of decay(dt|dterr) on dt. 
  //
  // Instead of integrating out dterr, make a weighted average of curves
  // at values dterr_i as given in the external dataset. 
  // (The kTRUE argument bins the data before projection to speed up the process)
  decay_gm.plotOn(frame2,ProjWData(*expDataDterr,kTRUE)) ;


  regTH(hh_decay,"rf306_model2d") ;
  regPlot(frame,"rf306_plot1") ;
  regPlot(frame2,"rf306_plot2") ;

  delete expDataDterr ;
  delete data ;

  return kTRUE ;
  }
} ;

//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #307
// 
// Complete example with use of full p.d.f. with per-event errors
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooGaussModel.h"
#include "RooDecay.h"
#include "RooLandau.h"
#include "RooProdPdf.h"
#include "RooHistPdf.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;


class TestBasic307 : public RooUnitTest
{
public: 
  TestBasic307(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Full per-event error p.d.f. F(t|dt)G(dt)",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // B - p h y s i c s   p d f   w i t h   p e r - e v e n t  G a u s s i a n   r e s o l u t i o n
  // ----------------------------------------------------------------------------------------------

  // Observables
  RooRealVar dt("dt","dt",-10,10) ;
  RooRealVar dterr("dterr","per-event error on dt",0.1,10) ;

  // Build a gaussian resolution model scaled by the per-event error = gauss(dt,bias,sigma*dterr)
  RooRealVar bias("bias","bias",0,-10,10) ;
  RooRealVar sigma("sigma","per-event error scale factor",1,0.1,10) ;
  RooGaussModel gm("gm1","gauss model scaled bt per-event error",dt,bias,sigma,dterr) ;

  // Construct decay(dt) (x) gauss1(dt|dterr)
  RooRealVar tau("tau","tau",1.548) ;
  RooDecay decay_gm("decay_gm","decay",dt,tau,gm,RooDecay::DoubleSided) ;



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

  // Use landau p.d.f to get empirical distribution with long tail
  RooLandau pdfDtErr("pdfDtErr","pdfDtErr",dterr,RooConst(1),RooConst(0.25)) ;
  RooDataSet* expDataDterr = pdfDtErr.generate(dterr,10000) ;

  // Construct a histogram pdf to describe the shape of the dtErr distribution
  RooDataHist* expHistDterr = expDataDterr->binnedClone() ;
  RooHistPdf pdfErr("pdfErr","pdfErr",dterr,*expHistDterr) ;


  // C o n s t r u c t   c o n d i t i o n a l   p r o d u c t   d e c a y _ d m ( d t | d t e r r ) * p d f ( d t e r r )
  // ----------------------------------------------------------------------------------------------------------------------

  // Construct production of conditional decay_dm(dt|dterr) with empirical pdfErr(dterr)
  RooProdPdf model("model","model",pdfErr,Conditional(decay_gm,dt)) ;

  // (Alternatively you could also use the landau shape pdfDtErr)
  //RooProdPdf model("model","model",pdfDtErr,Conditional(decay_gm,dt)) ;

  

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

  // Specify external dataset with dterr values to use model_dm as conditional p.d.f.
  RooDataSet* data = model.generate(RooArgSet(dt,dterr),10000) ;

  

  // F i t   c o n d i t i o n a l   d e c a y _ d m ( d t | d t e r r )
  // ---------------------------------------------------------------------

  // Specify dterr as conditional observable
  model.fitTo(*data) ;


  
  // P l o t   c o n d i t i o n a l   d e c a y _ d m ( d t | d t e r r )
  // ---------------------------------------------------------------------


  // Make projection of data an dt
  RooPlot* frame = dt.frame(Title("Projection of model(dt|dterr) on dt")) ;
  data->plotOn(frame) ;
  model.plotOn(frame) ;


  regPlot(frame,"rf307_plot1") ;

  delete expDataDterr ;
  delete expHistDterr ;
  delete data ;
  
  return kTRUE ;

  }
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #308
// 
// Examples on normalization of p.d.f.s,
// integration of p.d.fs, construction
// of cumulative distribution functions from p.d.f.s
// in two dimensions
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "RooAbsReal.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;


class TestBasic308 : public RooUnitTest
{
public: 
  TestBasic308(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Normalization of p.d.f.s in 2D",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

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

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

  // Create p.d.f. gaussx(x,-2,3), gaussy(y,2,2) 
  RooGaussian gx("gx","gx",x,RooConst(-2),RooConst(3)) ;
  RooGaussian gy("gy","gy",y,RooConst(+2),RooConst(2)) ;

  // Create gxy = gx(x)*gy(y)
  RooProdPdf gxy("gxy","gxy",RooArgSet(gx,gy)) ;



  // R e t r i e v e   r a w  &   n o r m a l i z e d   v a l u e s   o f   R o o F i t   p . d . f . s
  // --------------------------------------------------------------------------------------------------

  // Return 'raw' unnormalized value of gx
  regValue(gxy.getVal(),"rf308_gxy") ;
  
  // Return value of gxy normalized over x _and_ y in range [-10,10]
  RooArgSet nset_xy(x,y) ;
  regValue(gxy.getVal(&nset_xy),"rf308_gx_Norm[x,y]") ;

  // Create object representing integral over gx
  // which is used to calculate  gx_Norm[x,y] == gx / gx_Int[x,y]
  RooAbsReal* igxy = gxy.createIntegral(RooArgSet(x,y)) ;
  regValue(igxy->getVal(),"rf308_gx_Int[x,y]") ;

  // NB: it is also possible to do the following

  // Return value of gxy normalized over x in range [-10,10] (i.e. treating y as parameter)
  RooArgSet nset_x(x) ;
  regValue(gxy.getVal(&nset_x),"rf308_gx_Norm[x]") ;

  // Return value of gxy normalized over y in range [-10,10] (i.e. treating x as parameter)
  RooArgSet nset_y(y) ;
  regValue(gxy.getVal(&nset_y),"rf308_gx_Norm[y]") ;



  // I n t e g r a t e   n o r m a l i z e d   p d f   o v e r   s u b r a n g e
  // ----------------------------------------------------------------------------

  // Define a range named "signal" in x from -5,5
  x.setRange("signal",-5,5) ;
  y.setRange("signal",-3,3) ;
  
  // Create an integral of gxy_Norm[x,y] over x and y in range "signal"
  // This is the fraction of of p.d.f. gxy_Norm[x,y] which is in the
  // range named "signal"
  RooAbsReal* igxy_sig = gxy.createIntegral(RooArgSet(x,y),NormSet(RooArgSet(x,y)),Range("signal")) ;
  regValue(igxy_sig->getVal(),"rf308_gx_Int[x,y|signal]_Norm[x,y]") ;



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

  // Create the cumulative distribution function of gx
  // i.e. calculate Int[-10,x] gx(x') dx'
  RooAbsReal* gxy_cdf = gxy.createCdf(RooArgSet(x,y)) ;
  
  // Plot cdf of gx versus x
  TH1* hh_cdf = gxy_cdf->createHistogram("hh_cdf",x,Binning(40),YVar(y,Binning(40))) ;
  
  regTH(hh_cdf,"rf308_cdf") ;

  delete igxy_sig ;
  delete igxy ;
  delete gxy_cdf ;

  return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #309
// 
// Projecting p.d.f and data slices in discrete observables 
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussModel.h"
#include "RooDecay.h"
#include "RooBMixDecay.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic310 : public RooUnitTest
{
public: 
  TestBasic310(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Data and p.d.f projection in category slice",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   B   d e c a y   p d f   w it h   m i x i n g 
  // ----------------------------------------------------------

  // Decay time observables
  RooRealVar dt("dt","dt",-20,20) ;

  // Discrete observables mixState (B0tag==B0reco?) and tagFlav (B0tag==B0(bar)?)
  RooCategory mixState("mixState","B0/B0bar mixing state") ;
  RooCategory tagFlav("tagFlav","Flavour of the tagged B0") ;

  // Define state labels of discrete observables
  mixState.defineType("mixed",-1) ;
  mixState.defineType("unmixed",1) ;
  tagFlav.defineType("B0",1) ;
  tagFlav.defineType("B0bar",-1) ;

  // Model parameters
  RooRealVar dm("dm","delta m(B)",0.472,0.,1.0) ;
  RooRealVar tau("tau","B0 decay time",1.547,1.0,2.0) ;
  RooRealVar w("w","Flavor Mistag rate",0.03,0.0,1.0) ;
  RooRealVar dw("dw","Flavor Mistag rate difference between B0 and B0bar",0.01) ;

  // Build a gaussian resolution model
  RooRealVar bias1("bias1","bias1",0) ;
  RooRealVar sigma1("sigma1","sigma1",0.01) ;  
  RooGaussModel gm1("gm1","gauss model 1",dt,bias1,sigma1) ;

  // Construct a decay pdf, smeared with single gaussian resolution model
  RooBMixDecay bmix_gm1("bmix","decay",dt,mixState,tagFlav,tau,dm,w,dw,gm1,RooBMixDecay::DoubleSided) ;
  
  // Generate BMixing data with above set of event errors
  RooDataSet *data = bmix_gm1.generate(RooArgSet(dt,tagFlav,mixState),20000) ;



  // P l o t   f u l l   d e c a y   d i s t r i b u t i o n 
  // ----------------------------------------------------------

  // Create frame, plot data and pdf projection (integrated over tagFlav and mixState)
  RooPlot* frame = dt.frame(Title("Inclusive decay distribution")) ;
  data->plotOn(frame) ;
  bmix_gm1.plotOn(frame) ;



  // P l o t   d e c a y   d i s t r .   f o r   m i x e d   a n d   u n m i x e d   s l i c e   o f   m i x S t a t e
  // ------------------------------------------------------------------------------------------------------------------

  // Create frame, plot data (mixed only)
  RooPlot* frame2 = dt.frame(Title("Decay distribution of mixed events")) ;
  data->plotOn(frame2,Cut("mixState==mixState::mixed")) ;

  // Position slice in mixState at "mixed" and plot slice of pdf in mixstate over data (integrated over tagFlav)
  bmix_gm1.plotOn(frame2,Slice(mixState,"mixed")) ;

  // Create frame, plot data (unmixed only)
  RooPlot* frame3 = dt.frame(Title("Decay distribution of unmixed events")) ;
  data->plotOn(frame3,Cut("mixState==mixState::unmixed")) ;

  // Position slice in mixState at "unmixed" and plot slice of pdf in mixstate over data (integrated over tagFlav)
  bmix_gm1.plotOn(frame3,Slice(mixState,"unmixed")) ;


  regPlot(frame,"rf310_plot1") ;
  regPlot(frame2,"rf310_plot2") ;
  regPlot(frame3,"rf310_plot3") ;

  delete data ;

  return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #310
// 
// Projecting p.d.f and data ranges in continuous observables
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "RooAddPdf.h"
#include "RooPolynomial.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic311 : public RooUnitTest
{
public: 
  TestBasic311(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Data and p.d.f projection in sub range",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   3 D   p d f   a n d   d a t a 
  // -------------------------------------------

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

  // Create signal pdf gauss(x)*gauss(y)*gauss(z) 
  RooGaussian gx("gx","gx",x,RooConst(0),RooConst(1)) ;
  RooGaussian gy("gy","gy",y,RooConst(0),RooConst(1)) ;
  RooGaussian gz("gz","gz",z,RooConst(0),RooConst(1)) ;
  RooProdPdf sig("sig","sig",RooArgSet(gx,gy,gz)) ;

  // Create background pdf poly(x)*poly(y)*poly(z) 
  RooPolynomial px("px","px",x,RooArgSet(RooConst(-0.1),RooConst(0.004))) ;
  RooPolynomial py("py","py",y,RooArgSet(RooConst(0.1),RooConst(-0.004))) ;
  RooPolynomial pz("pz","pz",z) ;
  RooProdPdf bkg("bkg","bkg",RooArgSet(px,py,pz)) ;

  // Create composite pdf sig+bkg
  RooRealVar fsig("fsig","signal fraction",0.1,0.,1.) ;
  RooAddPdf model("model","model",RooArgList(sig,bkg),fsig) ;

  RooDataSet* data = model.generate(RooArgSet(x,y,z),20000) ;



  // P r o j e c t   p d f   a n d   d a t a   o n   x
  // -------------------------------------------------

  // Make plain projection of data and pdf on x observable
  RooPlot* frame = x.frame(Title("Projection of 3D data and pdf on X"),Bins(40)) ;
  data->plotOn(frame) ;
  model.plotOn(frame) ;
  


  // P r o j e c t   p d f   a n d   d a t a   o n   x   i n   s i g n a l   r a n g e
  // ----------------------------------------------------------------------------------
  
  // Define signal region in y and z observables
  y.setRange("sigRegion",-1,1) ;
  z.setRange("sigRegion",-1,1) ;

  // Make plot frame
  RooPlot* frame2 = x.frame(Title("Same projection on X in signal range of (Y,Z)"),Bins(40)) ;

  // Plot subset of data in which all observables are inside "sigRegion"
  // For observables that do not have an explicit "sigRegion" range defined (e.g. observable)
  // an implicit definition is used that is identical to the full range (i.e. [-5,5] for x)
  data->plotOn(frame2,CutRange("sigRegion")) ;

  // Project model on x, integrating projected observables (y,z) only in "sigRegion"
  model.plotOn(frame2,ProjectionRange("sigRegion")) ;


  regPlot(frame,"rf311_plot1") ;
  regPlot(frame2,"rf312_plot2") ;

  delete data ;

  return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #312
// 
// Performing fits in multiple (disjoint) ranges in one or more dimensions
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "RooAddPdf.h"
#include "RooPolynomial.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooFitResult.h"
using namespace RooFit ;


class TestBasic312 : public RooUnitTest
{
public: 
  TestBasic312(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Fit in multiple rectangular ranges",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   2 D   p d f   a n d   d a t a 
  // -------------------------------------------

  // Define observables x,y
  RooRealVar x("x","x",-10,10) ;
  RooRealVar y("y","y",-10,10) ;

  // Construct the signal pdf gauss(x)*gauss(y)
  RooRealVar mx("mx","mx",1,-10,10) ;
  RooRealVar my("my","my",1,-10,10) ;

  RooGaussian gx("gx","gx",x,mx,RooConst(1)) ;
  RooGaussian gy("gy","gy",y,my,RooConst(1)) ;

  RooProdPdf sig("sig","sig",gx,gy) ;

  // Construct the background pdf (flat in x,y)
  RooPolynomial px("px","px",x) ;
  RooPolynomial py("py","py",y) ;
  RooProdPdf bkg("bkg","bkg",px,py) ;

  // Construct the composite model sig+bkg
  RooRealVar f("f","f",0.,1.) ;
  RooAddPdf model("model","model",RooArgList(sig,bkg),f) ;

  // Sample 10000 events in (x,y) from the model
  RooDataSet* modelData = model.generate(RooArgSet(x,y),10000) ;



  // D e f i n e   s i g n a l   a n d   s i d e b a n d   r e g i o n s
  // -------------------------------------------------------------------

  // Construct the SideBand1,SideBand2,Signal regions
  //
  //                    |
  //      +-------------+-----------+                 
  //      |             |           |             
  //      |    Side     |   Sig     |        
  //      |    Band1    |   nal     |             
  //      |             |           |            
  //    --+-------------+-----------+--   
  //      |                         |       
  //      |           Side          |        
  //      |           Band2         |            
  //      |                         |          
  //      +-------------+-----------+            
  //                    |                       

  x.setRange("SB1",-10,+10) ;
  y.setRange("SB1",-10,0) ;

  x.setRange("SB2",-10,0) ;
  y.setRange("SB2",0,+10) ;

  x.setRange("SIG",0,+10) ;
  y.setRange("SIG",0,+10) ;

  x.setRange("FULL",-10,+10) ;
  y.setRange("FULL",-10,+10) ;


  // P e r f o r m   f i t s   i n   i n d i v i d u a l   s i d e b a n d   r e g i o n s
  // -------------------------------------------------------------------------------------

  // Perform fit in SideBand1 region (RooAddPdf coefficients will be interpreted in full range)
  RooFitResult* r_sb1 = model.fitTo(*modelData,Range("SB1"),Save()) ;

  // Perform fit in SideBand2 region (RooAddPdf coefficients will be interpreted in full range)
  RooFitResult* r_sb2 = model.fitTo(*modelData,Range("SB2"),Save()) ;



  // P e r f o r m   f i t s   i n   j o i n t    s i d e b a n d   r e g i o n s
  // -----------------------------------------------------------------------------

  // Now perform fit to joint 'L-shaped' sideband region 'SB1|SB2' 
  // (RooAddPdf coefficients will be interpreted in full range)
  RooFitResult* r_sb12 = model.fitTo(*modelData,Range("SB1,SB2"),Save()) ;


  regResult(r_sb1,"rf312_fit_sb1") ;
  regResult(r_sb2,"rf312_fit_sb2") ;
  regResult(r_sb12,"rf312_fit_sb12") ;

  delete modelData ;

  return kTRUE ;

  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #313
// 
// Working with parameterized ranges to define non-rectangular regions
// for fitting and integration
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooProdPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic313 : public RooUnitTest
{
public: 
  TestBasic313(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Integration over non-rectangular regions",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   3 D   p d f 
  // -------------------------

  // Define observable (x,y,z)
  RooRealVar x("x","x",0,10) ;
  RooRealVar y("y","y",0,10) ;
  RooRealVar z("z","z",0,10) ;

  // Define 3 dimensional pdf
  RooRealVar z0("z0","z0",-0.1,1) ;
  RooPolynomial px("px","px",x,RooConst(0)) ;
  RooPolynomial py("py","py",y,RooConst(0)) ;
  RooPolynomial pz("pz","pz",z,z0) ;
  RooProdPdf pxyz("pxyz","pxyz",RooArgSet(px,py,pz)) ;



  // D e f i n e d   n o n - r e c t a n g u l a r   r e g i o n   R   i n   ( x , y , z ) 
  // -------------------------------------------------------------------------------------

  //
  // R = Z[0 - 0.1*Y^2] * Y[0.1*X - 0.9*X] * X[0 - 10]
  //

  // Construct range parameterized in "R" in y [ 0.1*x, 0.9*x ]
  RooFormulaVar ylo("ylo","0.1*x",x) ;
  RooFormulaVar yhi("yhi","0.9*x",x) ;
  y.setRange("R",ylo,yhi) ;

  // Construct parameterized ranged "R" in z [ 0, 0.1*y^2 ]
  RooFormulaVar zlo("zlo","0.0*y",y) ;
  RooFormulaVar zhi("zhi","0.1*y*y",y) ;
  z.setRange("R",zlo,zhi) ;



  // C a l c u l a t e   i n t e g r a l   o f   n o r m a l i z e d   p d f   i n   R 
  // ----------------------------------------------------------------------------------

  // Create integral over normalized pdf model over x,y,z in "R" region
  RooAbsReal* intPdf = pxyz.createIntegral(RooArgSet(x,y,z),RooArgSet(x,y,z),"R") ;

  // Plot value of integral as function of pdf parameter z0
  RooPlot* frame = z0.frame(Title("Integral of pxyz over x,y,z in region R")) ;
  intPdf->plotOn(frame) ;


  regPlot(frame,"rf313_plot1") ;

  delete intPdf ;

  return kTRUE;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #314
// 
// Working with parameterized ranges in a fit. This an example of a
// fit with an acceptance that changes per-event 
//
//  pdf = exp(-t/tau) with t[tmin,5]
//
//  where t and tmin are both observables in the dataset
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooExponential.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooFitResult.h"

using namespace RooFit ;


class TestBasic314 : public RooUnitTest
{
public: 
  TestBasic314(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Fit with non-rectangular observable boundaries",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // 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,Name("dacc")) ;

  // Print fit results to demonstrate absence of bias
  regResult(r,"rf314_fit") ;
  regPlot(frame,"rf314_plot1") ;

  delete tmp ;
  delete dacc ;
  delete dall ;

  return kTRUE;
  }
} ;


//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #315
// 
// Marginizalization of multi-dimensional p.d.f.s through integration
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "RooPolyVar.h"
#include "TH1.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooNumIntConfig.h"
#include "RooConstVar.h"
using namespace RooFit ;



class TestBasic315 : public RooUnitTest
{
public: 
  TestBasic315(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("P.d.f. marginalization through integration",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

    // C r e a t e   p d f   m ( x , y )  =  g x ( x | y ) * g ( y )
    // --------------------------------------------------------------
    
    // Increase default precision of numeric integration
    // as this exercise has high sensitivity to numeric integration precision
    RooAbsPdf::defaultIntegratorConfig()->setEpsRel(1e-8) ;
    RooAbsPdf::defaultIntegratorConfig()->setEpsAbs(1e-8) ;
    
    // Create observables
    RooRealVar x("x","x",-5,5) ;
    RooRealVar y("y","y",-2,2) ;
    
    // Create function f(y) = a0 + a1*y
    RooRealVar a0("a0","a0",0) ;
    RooRealVar a1("a1","a1",-1.5,-3,1) ;
    RooPolyVar fy("fy","fy",y,RooArgSet(a0,a1)) ;
    
    // Create gaussx(x,f(y),sx)
    RooRealVar sigmax("sigmax","width of gaussian",0.5) ;
    RooGaussian gaussx("gaussx","Gaussian in x with shifting mean in y",x,fy,sigmax) ;  
    
    // Create gaussy(y,0,2)
    RooGaussian gaussy("gaussy","Gaussian in y",y,RooConst(0),RooConst(2)) ;
    
    // Create gaussx(x,sx|y) * gaussy(y)
    RooProdPdf model("model","gaussx(x|y)*gaussy(y)",gaussy,Conditional(gaussx,x)) ;
    
    
    
    // M a r g i n a l i z e   m ( x , y )   t o   m ( x ) 
    // ----------------------------------------------------
    
    // modelx(x) = Int model(x,y) dy
    RooAbsPdf* modelx = model.createProjection(y) ;
    
    
    
    // U s e   m a r g i n a l i z e d   p . d . f .   a s   r e g u l a r   1 - D   p . d . f .
    // ------------------------------------------------------------------------------------------
    
    // Sample 1000 events from modelx
    RooAbsData* data = modelx->generateBinned(x,1000) ;
    
    // Fit modelx to toy data
    modelx->fitTo(*data) ;
    
    // Plot modelx over data
    RooPlot* frame = x.frame(40) ;
    data->plotOn(frame) ;
    modelx->plotOn(frame) ;
    
    regPlot(frame,"rf315_frame") ;

    delete data ;
    delete modelx ;

    return kTRUE ;
  }
} ;





//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #316
// 
// Using the likelihood ratio techique to construct a signal enhanced
// one-dimensional projection of a multi-dimensional p.d.f.
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic316 : public RooUnitTest
{
public: 
  TestBasic316(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Likelihood ratio projection plot",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   3 D   p d f   a n d   d a t a 
  // -------------------------------------------

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

  // Create signal pdf gauss(x)*gauss(y)*gauss(z) 
  RooGaussian gx("gx","gx",x,RooConst(0),RooConst(1)) ;
  RooGaussian gy("gy","gy",y,RooConst(0),RooConst(1)) ;
  RooGaussian gz("gz","gz",z,RooConst(0),RooConst(1)) ;
  RooProdPdf sig("sig","sig",RooArgSet(gx,gy,gz)) ;

  // Create background pdf poly(x)*poly(y)*poly(z) 
  RooPolynomial px("px","px",x,RooArgSet(RooConst(-0.1),RooConst(0.004))) ;
  RooPolynomial py("py","py",y,RooArgSet(RooConst(0.1),RooConst(-0.004))) ;
  RooPolynomial pz("pz","pz",z) ;
  RooProdPdf bkg("bkg","bkg",RooArgSet(px,py,pz)) ;

  // Create composite pdf sig+bkg
  RooRealVar fsig("fsig","signal fraction",0.1,0.,1.) ;
  RooAddPdf model("model","model",RooArgList(sig,bkg),fsig) ;

  RooDataSet* data = model.generate(RooArgSet(x,y,z),20000) ;



  // P r o j e c t   p d f   a n d   d a t a   o n   x
  // -------------------------------------------------

  // Make plain projection of data and pdf on x observable
  RooPlot* frame = x.frame(Title("Projection of 3D data and pdf on X"),Bins(40)) ;
  data->plotOn(frame) ;
  model.plotOn(frame) ;
  


  // D e f i n e   p r o j e c t e d   s i g n a l   l i k e l i h o o d   r a t i o
  // ----------------------------------------------------------------------------------

  // Calculate projection of signal and total likelihood on (y,z) observables
  // i.e. integrate signal and composite model over x
  RooAbsPdf* sigyz = sig.createProjection(x) ;
  RooAbsPdf* totyz = model.createProjection(x) ;

  // Construct the log of the signal / signal+background probability 
  RooFormulaVar llratio_func("llratio","log10(@0)-log10(@1)",RooArgList(*sigyz,*totyz)) ;



  // P l o t   d a t a   w i t h   a   L L r a t i o   c u t 
  // -------------------------------------------------------

  // Calculate the llratio value for each event in the dataset
  data->addColumn(llratio_func) ;

  // Extract the subset of data with large signal likelihood
  RooDataSet* dataSel = (RooDataSet*) data->reduce(Cut("llratio>0.7")) ;

  // Make plot frame
  RooPlot* frame2 = x.frame(Title("Same projection on X with LLratio(y,z)>0.7"),Bins(40)) ;

  // Plot select data on frame
  dataSel->plotOn(frame2) ;



  // M a k e   M C   p r o j e c t i o n   o f   p d f   w i t h   s a m e   L L r a t i o   c u t 
  // ---------------------------------------------------------------------------------------------

  // Generate large number of events for MC integration of pdf projection
  RooDataSet* mcprojData = model.generate(RooArgSet(x,y,z),10000) ;

  // Calculate LL ratio for each generated event and select MC events with llratio)0.7
  mcprojData->addColumn(llratio_func) ;
  RooDataSet* mcprojDataSel = (RooDataSet*) mcprojData->reduce(Cut("llratio>0.7")) ;
    
  // Project model on x, integrating projected observables (y,z) with Monte Carlo technique
  // on set of events with the same llratio cut as was applied to data
  model.plotOn(frame2,ProjWData(*mcprojDataSel)) ;


  regPlot(frame,"rf316_plot1") ;
  regPlot(frame2,"rf316_plot2") ;
  
  delete data ;
  delete dataSel ;
  delete mcprojData ;
  delete mcprojDataSel ;
  delete sigyz ;
  delete totyz ;

  return kTRUE ;
  }  
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #402
// 
// Tools for manipulation of (un)binned datasets
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TFile.h"
using namespace RooFit ;


class TestBasic402 : public RooUnitTest
{
public: 
  TestBasic402(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Basic operations on datasets",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // Binned (RooDataHist) and unbinned datasets (RooDataSet) share
  // many properties and inherit from a common abstract base class
  // (RooAbsData), that provides an interface for all operations
  // that can be performed regardless of the data format

  RooRealVar  x("x","x",-10,10) ;
  RooRealVar  y("y","y", 0, 40) ;
  RooCategory c("c","c") ;
  c.defineType("Plus",+1) ;
  c.defineType("Minus",-1) ;



  // B a s i c   O p e r a t i o n s   o n   u n b i n n e d   d a t a s e t s 
  // --------------------------------------------------------------

  // RooDataSet is an unbinned dataset (a collection of points in N-dimensional space)
  RooDataSet d("d","d",RooArgSet(x,y,c)) ;

  // Unlike RooAbsArgs (RooAbsPdf,RooFormulaVar,....) datasets are not attached to 
  // the variables they are constructed from. Instead they are attached to an internal 
  // clone of the supplied set of arguments

  // Fill d with dummy values
  Int_t i ;
  for (i=0 ; i<1000 ; i++) {
    x = i/50 - 10 ;
    y = sqrt(1.0*i) ;
    c.setLabel((i%2)?"Plus":"Minus") ;

    // We must explicitly refer to x,y,c here to pass the values because
    // d is not linked to them (as explained above)
    d.add(RooArgSet(x,y,c)) ;
  }


  // R e d u c i n g ,   A p p e n d i n g   a n d   M e r g i n g
  // -------------------------------------------------------------

  // The reduce() function returns a new dataset which is a subset of the original
  RooDataSet* d1 = (RooDataSet*) d.reduce(RooArgSet(x,c)) ; 
  RooDataSet* d2 = (RooDataSet*) d.reduce(RooArgSet(y)) ;   
  RooDataSet* d3 = (RooDataSet*) d.reduce("y>5.17") ; 
  RooDataSet* d4 = (RooDataSet*) d.reduce(RooArgSet(x,c),"y>5.17") ; 

  regValue(d3->numEntries(),"rf403_nd3") ;
  regValue(d4->numEntries(),"rf403_nd4") ;

  // The merge() function adds two data set column-wise
  d1->merge(d2) ;

  // The append() function addes two datasets row-wise
  d1->append(*d3) ;

  regValue(d1->numEntries(),"rf403_nd1") ;

  


  // O p e r a t i o n s   o n   b i n n e d   d a t a s e t s 
  // ---------------------------------------------------------

  // A binned dataset can be constructed empty, from an unbinned dataset, or
  // from a ROOT native histogram (TH1,2,3)

  // The binning of real variables (like x,y) is done using their fit range
  // 'get/setRange()' and number of specified fit bins 'get/setBins()'.
  // Category dimensions of binned datasets get one bin per defined category state
  x.setBins(10) ;
  y.setBins(10) ;
  RooDataHist dh("dh","binned version of d",RooArgSet(x,y),d) ;

  RooPlot* yframe = y.frame(Bins(10),Title("Operations on binned datasets")) ;
  dh.plotOn(yframe) ; // plot projection of 2D binned data on y

  // Reduce the 2-dimensional binned dataset to a 1-dimensional binned dataset
  //
  // All reduce() methods are interfaced in RooAbsData. All reduction techniques
  // demonstrated on unbinned datasets can be applied to binned datasets as well.
  RooDataHist* dh2 = (RooDataHist*) dh.reduce(y,"x>0") ;

  // Add dh2 to yframe and redraw
  dh2->plotOn(yframe,LineColor(kRed),MarkerColor(kRed),Name("dh2")) ;

  regPlot(yframe,"rf402_plot1") ;

  delete d1 ;
  delete d2 ;
  delete d3 ;
  delete d4 ;
  delete dh2 ;
  return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #403
// 
// Using weights in unbinned datasets
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooFormulaVar.h"
#include "RooGenericPdf.h"
#include "RooPolynomial.h"
#include "RooChi2Var.h"
#include "RooMinuit.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooFitResult.h"
using namespace RooFit ;


class TestBasic403 : public RooUnitTest
{
public: 
  TestBasic403(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Fits with weighted datasets",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   o b s e r v a b l e   a n d   u n w e i g h t e d   d a t a s e t 
  // -------------------------------------------------------------------------------

  // Declare observable
  RooRealVar x("x","x",-10,10) ;
  x.setBins(40) ;

  // Construction a uniform pdf
  RooPolynomial p0("px","px",x) ;

  // Sample 1000 events from pdf
  RooDataSet* data = p0.generate(x,1000) ;

 

  // C a l c u l a t e   w e i g h t   a n d   m a k e   d a t a s e t   w e i g h t e d 
  // -----------------------------------------------------------------------------------

  // Construct formula to calculate (fake) weight for events
  RooFormulaVar wFunc("w","event weight","(x*x+10)",x) ;

  // Add column with variable w to previously generated dataset
  RooRealVar* w = (RooRealVar*) data->addColumn(wFunc) ;

  // Instruct dataset d in interpret w as event weight rather than as observable
  RooDataSet dataw(data->GetName(),data->GetTitle(),data,*data->get(),0,w->GetName()) ;
  //data->setWeightVar(*w) ;


  // U n b i n n e d   M L   f i t   t o   w e i g h t e d   d a t a 
  // ---------------------------------------------------------------

  // Construction quadratic polynomial pdf for fitting
  RooRealVar a0("a0","a0",1) ;
  RooRealVar a1("a1","a1",0,-1,1) ;
  RooRealVar a2("a2","a2",1,0,10) ;
  RooPolynomial p2("p2","p2",x,RooArgList(a0,a1,a2),0) ;

  // Fit quadratic polynomial to weighted data

  // NOTE: Maximum likelihood fit to weighted data does in general 
  //       NOT result in correct error estimates, unless individual
  //       event weights represent Poisson statistics themselves.
  //       In general, parameter error reflect precision of SumOfWeights
  //       events rather than NumEvents events. See comparisons below
  RooFitResult* r_ml_wgt = p2.fitTo(dataw,Save()) ;



  // P l o t   w e i g h e d   d a t a   a n d   f i t   r e s u l t 
  // ---------------------------------------------------------------

  // Construct plot frame
  RooPlot* frame = x.frame(Title("Unbinned ML fit, binned chi^2 fit to weighted data")) ;

  // Plot data using sum-of-weights-squared error rather than Poisson errors
  dataw.plotOn(frame,DataError(RooAbsData::SumW2)) ;

  // Overlay result of 2nd order polynomial fit to weighted data
  p2.plotOn(frame) ;



  // M L  F i t   o f   p d f   t o   e q u i v a l e n t  u n w e i g h t e d   d a t a s e t
  // -----------------------------------------------------------------------------------------
  
  // Construct a pdf with the same shape as p0 after weighting
  RooGenericPdf genPdf("genPdf","x*x+10",x) ;

  // Sample a dataset with the same number of events as data
  RooDataSet* data2 = genPdf.generate(x,1000) ;

  // Sample a dataset with the same number of weights as data
  RooDataSet* data3 = genPdf.generate(x,43000) ;

  // Fit the 2nd order polynomial to both unweighted datasets and save the results for comparison
  RooFitResult* r_ml_unw10 = p2.fitTo(*data2,Save()) ;
  RooFitResult* r_ml_unw43 = p2.fitTo(*data3,Save()) ;


  // C h i 2   f i t   o f   p d f   t o   b i n n e d   w e i g h t e d   d a t a s e t
  // ------------------------------------------------------------------------------------

  // Construct binned clone of unbinned weighted dataset
  RooDataHist* binnedData = dataw.binnedClone() ;

  // Perform chi2 fit to binned weighted dataset using sum-of-weights errors
  // 
  // NB: Within the usual approximations of a chi2 fit, a chi2 fit to weighted
  // data using sum-of-weights-squared errors does give correct error
  // estimates
  RooChi2Var chi2("chi2","chi2",p2,*binnedData) ;
  RooMinuit m(chi2) ;
  m.migrad() ;
  m.hesse() ;

  // Plot chi^2 fit result on frame as well
  RooFitResult* r_chi2_wgt = m.save() ;
  p2.plotOn(frame,LineStyle(kDashed),LineColor(kRed),Name("p2_alt")) ;



  // C o m p a r e   f i t   r e s u l t s   o f   c h i 2 , M L   f i t s   t o   ( u n ) w e i g h t e d   d a t a 
  // ---------------------------------------------------------------------------------------------------------------

  // Note that ML fit on 1Kevt of weighted data is closer to result of ML fit on 43Kevt of unweighted data 
  // than to 1Kevt of unweighted data, whereas the reference chi^2 fit with SumW2 error gives a result closer to
  // that of an unbinned ML fit to 1Kevt of unweighted data. 

  regResult(r_ml_unw10,"rf403_ml_unw10") ;
  regResult(r_ml_unw43,"rf403_ml_unw43") ;
  regResult(r_ml_wgt  ,"rf403_ml_wgt") ;
  regResult(r_chi2_wgt ,"rf403_ml_chi2") ;
  regPlot(frame,"rf403_plot1") ;
  
  delete binnedData ;
  delete data ;
  delete data2 ;
  delete data3 ;

  return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #404
// 
// Working with RooCategory objects to describe discrete variables
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooPolynomial.h"
#include "RooCategory.h"
#include "Roo1DTable.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic404 : public RooUnitTest
{
public: 
  TestBasic404(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Categories basic functionality",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C o n s t r u c t    a   c a t e g o r y   w i t h   l a b e l s 
  // ----------------------------------------------------------------

  // Define a category with labels only
  RooCategory tagCat("tagCat","Tagging category") ;
  tagCat.defineType("Lepton") ;
  tagCat.defineType("Kaon") ;
  tagCat.defineType("NetTagger-1") ;
  tagCat.defineType("NetTagger-2") ;



  // C o n s t r u c t    a   c a t e g o r y   w i t h   l a b e l s   a n d   i n d e c e s
  // ----------------------------------------------------------------------------------------

  // Define a category with explicitly numbered states
  RooCategory b0flav("b0flav","B0 flavour eigenstate") ;
  b0flav.defineType("B0",-1) ;
  b0flav.defineType("B0bar",1) ;



  // G e n e r a t e   d u m m y   d a t a  f o r   t a b u l a t i o n   d e m o 
  // ----------------------------------------------------------------------------
  
  // Generate a dummy dataset 
  RooRealVar x("x","x",0,10) ;
  RooDataSet *data = RooPolynomial("p","p",x).generate(RooArgSet(x,b0flav,tagCat),10000) ;



  // P r i n t   t a b l e s   o f   c a t e g o r y   c o n t e n t s   o f   d a t a s e t s
  // ------------------------------------------------------------------------------------------

  // Tables are equivalent of plots for categories
  Roo1DTable* btable = data->table(b0flav) ;

  // Create table for subset of events matching cut expression
  Roo1DTable* ttable = data->table(tagCat,"x>8.23") ;

  // Create table for all (tagCat x b0flav) state combinations
  Roo1DTable* bttable = data->table(RooArgSet(tagCat,b0flav)) ;

  // Retrieve number of events from table
  // Number can be non-integer if source dataset has weighed events
  Double_t nb0 = btable->get("B0") ;
  regValue(nb0,"rf404_nb0") ;

  // Retrieve fraction of events with "Lepton" tag
  Double_t fracLep = ttable->getFrac("Lepton") ;
  regValue(fracLep,"rf404_fracLep") ;  


  // D e f i n i n g   r a n g e s   f o r   p l o t t i n g ,   f i t t i n g   o n   c a t e g o r i e s
  // ------------------------------------------------------------------------------------------------------

  // Define named range as comma separated list of labels
  tagCat.setRange("good","Lepton,Kaon") ;

  // Or add state names one by one
  tagCat.addToRange("soso","NetTagger-1") ;
  tagCat.addToRange("soso","NetTagger-2") ;
  
  // Use category range in dataset reduction specification
  RooDataSet* goodData = (RooDataSet*) data->reduce(CutRange("good")) ;
  Roo1DTable* gtable = goodData->table(tagCat) ;
  

  regTable(btable,"rf404_btable") ;
  regTable(ttable,"rf404_ttable") ;
  regTable(bttable,"rf404_bttable") ;
  regTable(gtable,"rf404_gtable") ;


  delete goodData ;
  delete data ;
  return kTRUE ;

  }

} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #405
// 
// Demonstration of real-->discrete mapping functions
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooCategory.h"
#include "RooThresholdCategory.h"
#include "RooBinningCategory.h"
#include "Roo1DTable.h"
#include "RooArgusBG.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooRealConstant.h"
using namespace RooFit ;


class TestBasic405 : public RooUnitTest
{
public: 
  TestBasic405(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Real-to-category functions",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {


  // D e f i n e   p d f   i n   x ,   s a m p l e   d a t a s e t   i n   x 
  // ------------------------------------------------------------------------


  // Define a dummy PDF in x 
  RooRealVar x("x","x",0,10) ;
  RooArgusBG a("a","argus(x)",x,RooRealConstant::value(10),RooRealConstant::value(-1)) ;

  // Generate a dummy dataset 
  RooDataSet *data = a.generate(x,10000) ;



  // C r e a t e   a   t h r e s h o l d   r e a l - > c a t   f u n c t i o n
  // --------------------------------------------------------------------------

  // A RooThresholdCategory is a category function that maps regions in a real-valued 
  // input observable observables to state names. At construction time a 'default'
  // state name must be specified to which all values of x are mapped that are not
  // otherwise assigned
  RooThresholdCategory xRegion("xRegion","region of x",x,"Background") ;

  // Specify thresholds and state assignments one-by-one. 
  // Each statement specifies that all values _below_ the given value 
  // (and above any lower specified threshold) are mapped to the
  // category state with the given name
  //
  // Background | SideBand | Signal | SideBand | Background
  //           4.23       5.23     8.23       9.23 
  xRegion.addThreshold(4.23,"Background") ;
  xRegion.addThreshold(5.23,"SideBand") ;
  xRegion.addThreshold(8.23,"Signal") ;
  xRegion.addThreshold(9.23,"SideBand") ; 



  // U s e   t h r e s h o l d   f u n c t i o n   t o   p l o t   d a t a   r e g i o n s
  // -------------------------------------------------------------------------------------

  // Add values of threshold function to dataset so that it can be used as observable
  data->addColumn(xRegion) ;

  // Make plot of data in x
  RooPlot* xframe = x.frame(Title("Demo of threshold and binning mapping functions")) ;
  data->plotOn(xframe) ;

  // Use calculated category to select sideband data
  data->plotOn(xframe,Cut("xRegion==xRegion::SideBand"),MarkerColor(kRed),LineColor(kRed),Name("data_cut")) ;



  // C r e a t e   a   b i n n i n g    r e a l - > c a t   f u n c t i o n
  // ----------------------------------------------------------------------

  // A RooBinningCategory is a category function that maps bins of a (named) binning definition 
  // in a real-valued input observable observables to state names. The state names are automatically
  // constructed from the variable name, the binning name and the bin number. If no binning name
  // is specified the default binning is mapped

  x.setBins(10,"coarse") ;
  RooBinningCategory xBins("xBins","coarse bins in x",x,"coarse") ;



  // U s e   b i n n i n g   f u n c t i o n   f o r   t a b u l a t i o n   a n d   p l o t t i n g
  // -----------------------------------------------------------------------------------------------

  // Print table of xBins state multiplicity. Note that xBins does not need to be an observable in data
  // it can be a function of observables in data as well
  Roo1DTable* xbtable = data->table(xBins) ;

  // Add values of xBins function to dataset so that it can be used as observable
  RooCategory* xb = (RooCategory*) data->addColumn(xBins) ;

  // Define range "alt" as including bins 1,3,5,7,9 
  xb->setRange("alt","x_coarse_bin1,x_coarse_bin3,x_coarse_bin5,x_coarse_bin7,x_coarse_bin9") ;
  
  // Construct subset of data matching range "alt" but only for the first 5000 events and plot it on the fram
  RooDataSet* dataSel = (RooDataSet*) data->reduce(CutRange("alt"),EventRange(0,5000)) ;
//   dataSel->plotOn(xframe,MarkerColor(kGreen),LineColor(kGreen),Name("data_sel")) ;


  regTable(xbtable,"rf405_xbtable") ;
  regPlot(xframe,"rf405_plot1") ;
  
  delete data ;
  delete dataSel ;

  return kTRUE ;

  }

} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #406
// 
// Demonstration of discrete-->discrete (invertable) functions
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooPolynomial.h"
#include "RooCategory.h"
#include "RooMappedCategory.h"
#include "RooMultiCategory.h"
#include "RooSuperCategory.h"
#include "Roo1DTable.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic406 : public RooUnitTest
{
public: 
  TestBasic406(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Category-to-category functions",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C o n s t r u c t  t w o   c a t e g o r i e s
  // ----------------------------------------------

  // Define a category with labels only
  RooCategory tagCat("tagCat","Tagging category") ;
  tagCat.defineType("Lepton") ;
  tagCat.defineType("Kaon") ;
  tagCat.defineType("NetTagger-1") ;
  tagCat.defineType("NetTagger-2") ;

  // Define a category with explicitly numbered states
  RooCategory b0flav("b0flav","B0 flavour eigenstate") ;
  b0flav.defineType("B0",-1) ;
  b0flav.defineType("B0bar",1) ;

  // Construct a dummy dataset with random values of tagCat and b0flav
  RooRealVar x("x","x",0,10) ;
  RooPolynomial p("p","p",x) ;
  RooDataSet* data = p.generate(RooArgSet(x,b0flav,tagCat),10000) ;



  // C r e a t e   a   c a t - > c a t   m  a p p i n g   c a t e g o r y 
  // ---------------------------------------------------------------------

  // A RooMappedCategory is category->category mapping function based on string expression
  // The constructor takes an input category an a default state name to which unassigned
  // states are mapped
  RooMappedCategory tcatType("tcatType","tagCat type",tagCat,"Cut based") ;

  // Enter fully specified state mappings
  tcatType.map("Lepton","Cut based") ;
  tcatType.map("Kaon","Cut based") ;

  // Enter a wilcard expression mapping
  tcatType.map("NetTagger*","Neural Network") ;

  // Make a table of the mapped category state multiplicit in data
  Roo1DTable* mtable = data->table(tcatType) ;



  // C r e a t e   a   c a t   X   c a t   p r o d u c t   c a t e g o r y 
  // ----------------------------------------------------------------------

  // A SUPER-category is 'product' of _lvalue_ categories. The state names of a super
  // category is a composite of the state labels of the input categories
  RooSuperCategory b0Xtcat("b0Xtcat","b0flav X tagCat",RooArgSet(b0flav,tagCat)) ;

  // Make a table of the product category state multiplicity in data
  Roo1DTable* stable = data->table(b0Xtcat) ;

  // Since the super category is an lvalue, assignment is explicitly possible
  b0Xtcat.setLabel("{B0bar;Lepton}") ;



  // A MULTI-category is a 'product' of any category (function). The state names of a super
  // category is a composite of the state labels of the input categories
  RooMultiCategory b0Xttype("b0Xttype","b0flav X tagType",RooArgSet(b0flav,tcatType)) ;
  
  // Make a table of the product category state multiplicity in data
  Roo1DTable* xtable = data->table(b0Xttype) ;

  regTable(mtable,"rf406_mtable") ;
  regTable(stable,"rf406_stable") ;
  regTable(xtable,"rf406_xtable") ;

  delete data ;

  return kTRUE ;
  }

} ;
//////////////////////////////////////////////////////////////////////////
//
// '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 "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic501 : public RooUnitTest
{
public: 
  TestBasic501(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Simultaneous p.d.f. operator",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

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

  // Create observables
  RooRealVar x("x","x",-8,8) ;

  // Construct signal pdf
  RooRealVar mean("mean","mean",0,-8,8) ;
  RooRealVar sigma("sigma","sigma",0.3,0.1,10) ;
  RooGaussian gx("gx","gx",x,mean,sigma) ;

  // Construct background pdf
  RooRealVar a0("a0","a0",-0.1,-1,1) ;
  RooRealVar a1("a1","a1",0.004,-1,1) ;
  RooChebychev px("px","px",x,RooArgSet(a0,a1)) ;

  // Construct composite pdf
  RooRealVar f("f","f",0.2,0.,1.) ;
  RooAddPdf model("model","model",RooArgList(gx,px),f) ;



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

  // Construct signal pdf. 
  // NOTE that sigma is shared with the signal sample model
  RooRealVar mean_ctl("mean_ctl","mean_ctl",-3,-8,8) ;
  RooGaussian gx_ctl("gx_ctl","gx_ctl",x,mean_ctl,sigma) ;

  // Construct the background pdf
  RooRealVar a0_ctl("a0_ctl","a0_ctl",-0.1,-1,1) ;
  RooRealVar a1_ctl("a1_ctl","a1_ctl",0.5,-0.1,1) ;
  RooChebychev px_ctl("px_ctl","px_ctl",x,RooArgSet(a0_ctl,a1_ctl)) ;

  // Construct the composite model
  RooRealVar f_ctl("f_ctl","f_ctl",0.5,0.,1.) ;
  RooAddPdf model_ctl("model_ctl","model_ctl",RooArgList(gx_ctl,px_ctl),f_ctl) ;



  // G e n e r a t e   e v e n t s   f o r   b o t h   s a m p l e s 
  // ---------------------------------------------------------------

  // Generate 1000 events in x and y from model
  RooDataSet *data = model.generate(RooArgSet(x),100) ;
  RooDataSet *data_ctl = model_ctl.generate(RooArgSet(x),2000) ;



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

  // Define category to distinguish physics and control samples events
  RooCategory sample("sample","sample") ;
  sample.defineType("physics") ;
  sample.defineType("control") ;

  // Construct combined dataset in (x,sample)
  RooDataSet combData("combData","combined data",x,Index(sample),Import("physics",*data),Import("control",*data_ctl)) ;



  // C o n s t r u c t   a   s i m u l t a n e o u s   p d f   i n   ( x , s a m p l e )
  // -----------------------------------------------------------------------------------

  // Construct a simultaneous pdf using category sample as index
  RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;

  // Associate model with the physics state and model_ctl with the control state
  simPdf.addPdf(model,"physics") ;
  simPdf.addPdf(model_ctl,"control") ;



  // P e r f o r m   a   s i m u l t a n e o u s   f i t
  // ---------------------------------------------------

  // Perform simultaneous fit of model to data and model_ctl to data_ctl
  simPdf.fitTo(combData) ;



  // P l o t   m o d e l   s l i c e s   o n   d a t a    s l i c e s 
  // ----------------------------------------------------------------

  // Make a frame for the physics sample
  RooPlot* frame1 = x.frame(Bins(30),Title("Physics sample")) ;

  // Plot all data tagged as physics sample
  combData.plotOn(frame1,Cut("sample==sample::physics")) ;

  // Plot "physics" slice of simultaneous pdf. 
  // NBL You _must_ project the sample index category with data using ProjWData 
  // as a RooSimultaneous makes no prediction on the shape in the index category 
  // and can thus not be integrated
  simPdf.plotOn(frame1,Slice(sample,"physics"),ProjWData(sample,combData)) ;
  simPdf.plotOn(frame1,Slice(sample,"physics"),Components("px"),ProjWData(sample,combData),LineStyle(kDashed)) ;

  // The same plot for the control sample slice
  RooPlot* frame2 = x.frame(Bins(30),Title("Control sample")) ;
  combData.plotOn(frame2,Cut("sample==sample::control")) ;
  simPdf.plotOn(frame2,Slice(sample,"control"),ProjWData(sample,combData)) ;
  simPdf.plotOn(frame2,Slice(sample,"control"),Components("px_ctl"),ProjWData(sample,combData),LineStyle(kDashed)) ;


  regPlot(frame1,"rf501_plot1") ;
  regPlot(frame2,"rf501_plot2") ;

  delete data ;
  delete data_ctl ;

  return kTRUE ;

  }

} ;
//////////////////////////////////////////////////////////////////////////
//
// '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 RooUnitTest
{
public: 
  TestBasic599(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("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 ;
  }
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #601
// 
// Interactive minimization with MINUIT
//
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "RooAddPdf.h"
#include "RooMinuit.h"
#include "RooNLLVar.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;


class TestBasic601 : public RooUnitTest
{
public: 
  TestBasic601(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Interactive Minuit",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // S e t u p   p d f   a n d   l i k e l i h o o d 
  // -----------------------------------------------

  // Observable
  RooRealVar x("x","x",-20,20) ;

  // Model (intentional strong correlations)
  RooRealVar mean("mean","mean of g1 and g2",0) ;
  RooRealVar sigma_g1("sigma_g1","width of g1",3) ;
  RooGaussian g1("g1","g1",x,mean,sigma_g1) ;

  RooRealVar sigma_g2("sigma_g2","width of g2",4,3.0,6.0) ;
  RooGaussian g2("g2","g2",x,mean,sigma_g2) ;

  RooRealVar frac("frac","frac",0.5,0.0,1.0) ;
  RooAddPdf model("model","model",RooArgList(g1,g2),frac) ;

  // Generate 1000 events
  RooDataSet* data = model.generate(x,1000) ;
  
  // Construct unbinned likelihood
  RooNLLVar nll("nll","nll",model,*data) ;


  // I n t e r a c t i v e   m i n i m i z a t i o n ,   e r r o r   a n a l y s i s
  // -------------------------------------------------------------------------------

  // Create MINUIT interface object
  RooMinuit m(nll) ;

  // Call MIGRAD to minimize the likelihood
  m.migrad() ;

  // Run HESSE to calculate errors from d2L/dp2
  m.hesse() ;

  // Run MINOS on sigma_g2 parameter only
  m.minos(sigma_g2) ;


  // S a v i n g   r e s u l t s ,   c o n t o u r   p l o t s 
  // ---------------------------------------------------------

  // Save a snapshot of the fit result. This object contains the initial
  // fit parameters, the final fit parameters, the complete correlation
  // matrix, the EDM, the minimized FCN , the last MINUIT status code and
  // the number of times the RooFit function object has indicated evaluation
  // problems (e.g. zero probabilities during likelihood evaluation)
  RooFitResult* r = m.save() ;


  // C h a n g e   p a r a m e t e r   v a l u e s ,   f l o a t i n g
  // -----------------------------------------------------------------

  // At any moment you can manually change the value of a (constant)
  // parameter
  mean = 0.3 ;
  
  // Rerun MIGRAD,HESSE
  m.migrad() ;
  m.hesse() ;

  // Now fix sigma_g2
  sigma_g2.setConstant(kTRUE) ;

  // Rerun MIGRAD,HESSE
  m.migrad() ;
  m.hesse() ;

  RooFitResult* r2 = m.save() ;

  regResult(r,"rf601_r") ;
  regResult(r2,"rf601_r2") ;

  delete data ;

  return kTRUE ;
  } 
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #602
// 
// Setting up a binning chi^2 fit
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooChi2Var.h"
#include "RooMinuit.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic602 : public RooUnitTest
{
public: 
  TestBasic602(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Chi2 minimization",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // 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,
  // which is the input probability density scaled
  // by the number of events in the dataset
  RooChi2Var chi2("chi2","chi2",model,*dh) ;

  // Use RooMinuit interface to minimize chi^2
  RooMinuit m(chi2) ;
  m.migrad() ;
  m.hesse() ;

  RooFitResult* r = m.save() ;

  regResult(r,"rf602_r") ;
  
  delete d ;
  delete dh ;

  return kTRUE ;
  }
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #604
// 
// Fitting with constraints
//
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;


class TestBasic604 : public RooUnitTest
{
public: 
  TestBasic604(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Auxiliary observable constraints",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   m o d e l  a n d   d a t a s e t 
  // ----------------------------------------------

  // Construct a Gaussian p.d.f
  RooRealVar x("x","x",-10,10) ;

  RooRealVar m("m","m",0,-10,10) ;
  RooRealVar s("s","s",2,0.1,10) ;
  RooGaussian gauss("gauss","gauss(x,m,s)",x,m,s) ;

  // Construct a flat p.d.f (polynomial of 0th order)
  RooPolynomial poly("poly","poly(x)",x) ;

  // Construct model = f*gauss + (1-f)*poly
  RooRealVar f("f","f",0.5,0.,1.) ;
  RooAddPdf model("model","model",RooArgSet(gauss,poly),f) ;

  // Generate small dataset for use in fitting below
  RooDataSet* d = model.generate(x,50) ;



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

  // Construct Gaussian constraint p.d.f on parameter f at 0.8 with resolution of 0.1
  RooGaussian fconstraint("fconstraint","fconstraint",f,RooConst(0.8),RooConst(0.1)) ;



  // M E T H O D   1   -   A d d   i n t e r n a l   c o n s t r a i n t   t o   m o d e l 
  // -------------------------------------------------------------------------------------

  // Multiply constraint term with regular p.d.f using RooProdPdf
  // Specify in fitTo() that internal constraints on parameter f should be used

  // Multiply constraint with p.d.f
  RooProdPdf modelc("modelc","model with constraint",RooArgSet(model,fconstraint)) ;
  
  // Fit modelc without use of constraint term
  RooFitResult* r1 = modelc.fitTo(*d,Save()) ;

  // Fit modelc with constraint term on parameter f
  RooFitResult* r2 = modelc.fitTo(*d,Constrain(f),Save()) ;



  // M E T H O D   2   -     S p e c i f y   e x t e r n a l   c o n s t r a i n t   w h e n   f i t t i n g
  // -------------------------------------------------------------------------------------------------------

  // Construct another Gaussian constraint p.d.f on parameter f at 0.8 with resolution of 0.1
  RooGaussian fconstext("fconstext","fconstext",f,RooConst(0.2),RooConst(0.1)) ;

  // Fit with external constraint
  RooFitResult* r3 = model.fitTo(*d,ExternalConstraints(fconstext),Save()) ;


  regResult(r1,"rf604_r1") ;
  regResult(r2,"rf604_r2") ;
  regResult(r3,"rf604_r3") ;

  delete d ;
  
  return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #605
// 
// Working with the profile likelihood estimator
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooAddPdf.h"
#include "RooNLLVar.h"
#include "RooProfileLL.h"
#include "RooMinuit.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic605 : public RooUnitTest
{
public: 
  TestBasic605(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Profile Likelihood operator",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   m o d e l   a n d   d a t a s e t 
  // -----------------------------------------------

  // Observable
  RooRealVar x("x","x",-20,20) ;

  // Model (intentional strong correlations)
  RooRealVar mean("mean","mean of g1 and g2",0,-10,10) ;
  RooRealVar sigma_g1("sigma_g1","width of g1",3) ;
  RooGaussian g1("g1","g1",x,mean,sigma_g1) ;

  RooRealVar sigma_g2("sigma_g2","width of g2",4,3.0,6.0) ;
  RooGaussian g2("g2","g2",x,mean,sigma_g2) ;

  RooRealVar frac("frac","frac",0.5,0.0,1.0) ;
  RooAddPdf model("model","model",RooArgList(g1,g2),frac) ;

  // Generate 1000 events
  RooDataSet* data = model.generate(x,1000) ;
  


  // C o n s t r u c t   p l a i n   l i k e l i h o o d
  // ---------------------------------------------------

  // Construct unbinned likelihood
  RooNLLVar nll("nll","nll",model,*data) ;

  // Minimize likelihood w.r.t all parameters before making plots
  RooMinuit(nll).migrad() ;

  // Plot likelihood scan frac 
  RooPlot* frame1 = frac.frame(Bins(10),Range(0.01,0.95),Title("LL and profileLL in frac")) ;
  nll.plotOn(frame1,ShiftToZero()) ;

  // Plot likelihood scan in sigma_g2
  RooPlot* frame2 = sigma_g2.frame(Bins(10),Range(3.3,5.0),Title("LL and profileLL in sigma_g2")) ;
  nll.plotOn(frame2,ShiftToZero()) ;



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

  // The profile likelihood estimator on nll for frac will minimize nll w.r.t
  // all floating parameters except frac for each evaluation
  RooProfileLL pll_frac("pll_frac","pll_frac",nll,frac) ;

  // Plot the profile likelihood in frac
  pll_frac.plotOn(frame1,LineColor(kRed)) ;

  // Adjust frame maximum for visual clarity
  frame1->SetMinimum(0) ;
  frame1->SetMaximum(3) ;



  // C o n s t r u c t   p r o f i l e   l i k e l i h o o d   i n   s i g m a _ g 2 
  // -------------------------------------------------------------------------------

  // The profile likelihood estimator on nll for sigma_g2 will minimize nll
  // w.r.t all floating parameters except sigma_g2 for each evaluation
  RooProfileLL pll_sigmag2("pll_sigmag2","pll_sigmag2",nll,sigma_g2) ;

  // Plot the profile likelihood in sigma_g2
  pll_sigmag2.plotOn(frame2,LineColor(kRed)) ;

  // Adjust frame maximum for visual clarity
  frame2->SetMinimum(0) ;
  frame2->SetMaximum(3) ;


  regPlot(frame1,"rf605_plot1") ;
  regPlot(frame2,"rf605_plot2") ;

  delete data ;

  return kTRUE ;
  }
} ;

//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #606
// 
// Understanding and customizing error handling in likelihood evaluations
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooArgusBG.h"
#include "RooNLLVar.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


class TestBasic606 : public RooUnitTest
{
public: 
  TestBasic606(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("NLL error handling",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   m o d e l  a n d   d a t a s e t 
  // ----------------------------------------------

  // Observable
  RooRealVar m("m","m",5.20,5.30) ;

  // Parameters
  RooRealVar m0("m0","m0",5.291,5.20,5.30) ;
  RooRealVar k("k","k",-30,-50,-10) ;

  // Pdf
  RooArgusBG argus("argus","argus",m,m0,k) ;

  // Sample 1000 events in m from argus
  RooDataSet* data = argus.generate(m,1000) ;



  // P l o t   m o d e l   a n d   d a t a 
  // --------------------------------------

  RooPlot* frame1 = m.frame(Bins(40),Title("Argus model and data")) ;
  data->plotOn(frame1) ;
  argus.plotOn(frame1) ;



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

  argus.fitTo(*data,PrintEvalErrors(10),Warnings(kFALSE)) ;    
  m0.setError(0.1) ;
  argus.fitTo(*data,PrintEvalErrors(0),EvalErrorWall(kFALSE),Warnings(kFALSE)) ;    



  // P l o t   l i k e l i h o o d   a s   f u n c t i o n   o f   m 0 
  // ------------------------------------------------------------------

  // Construct likelihood function of model and data
  RooNLLVar nll("nll","nll",argus,*data) ;

  // Plot likelihood in m0 in range that includes problematic values
  // In this configuration the number of errors per likelihood point 
  // evaluated for the curve is shown. A positive number in PrintEvalErrors(N)
  // will show details for up to N events. By default the values for likelihood
  // evaluations with errors are shown normally (unlike fitting), but the shape
  // of the curve can be erratic in these regions.

  RooPlot* frame2 = m0.frame(Range(5.288,5.293),Title("-log(L) scan vs m0")) ;
  nll.plotOn(frame2,PrintEvalErrors(0),ShiftToZero(),LineColor(kRed),Precision(1e-4)) ; 


  // Plot likelihood in m0 in range that includes problematic values
  // In this configuration no messages are printed for likelihood evaluation errors,
  // but if an likelihood value evaluates with error, the corresponding value
  // on the curve will be set to the value given in EvalErrorValue().

  RooPlot* frame3 = m0.frame(Range(5.288,5.293),Title("-log(L) scan vs m0, problematic regions masked")) ;
  nll.plotOn(frame3,PrintEvalErrors(-1),ShiftToZero(),EvalErrorValue(nll.getVal()+10),LineColor(kRed)) ; 


  regPlot(frame1,"rf606_plot1") ;
  regPlot(frame2,"rf606_plot2") ;
  regPlot(frame3,"rf606_plot3") ;

  delete data ;
  return kTRUE ;

  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #607
// 
// Demonstration of options of the RooFitResult class
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooAddPdf.h"
#include "RooChebychev.h"
#include "RooFitResult.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TFile.h"
#include "TStyle.h"
#include "TH2.h"

using namespace RooFit ;


class TestBasic607 : public RooUnitTest
{
public: 
  TestBasic607(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Fit Result functionality",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   p d f ,   d a t a
  // --------------------------------

  // 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,-10,10) ;
  RooRealVar sigma1("sigma1","width of gaussians",0.5,0.1,10) ;
  RooRealVar sigma2("sigma2","width of gaussians",1,0.1,10) ;

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

  // Generate 1000 events
  RooDataSet* data = model.generate(x,1000) ;



  // F i t   p d f   t o   d a t a ,   s a v e   f i t r e s u l t 
  // -------------------------------------------------------------

  // Perform fit and save result
  RooFitResult* r = model.fitTo(*data,Save()) ;


  // V i s u a l i z e   c o r r e l a t i o n   m a t r i x
  // -------------------------------------------------------

  // Construct 2D color plot of correlation matrix
  gStyle->SetOptStat(0) ;
  gStyle->SetPalette(1) ;
  TH2* hcorr = r->correlationHist() ;


  // Sample dataset with parameter values according to distribution
  // of covariance matrix of fit result
  RooDataSet randPars("randPars","randPars",r->floatParsFinal()) ;
  for (Int_t i=0 ; i<10000 ; i++) {
    randPars.add(r->randomizePars()) ;    
  }

  // make histogram of 2D distribution in sigma1 vs sig1frac
  TH1* hhrand = randPars.createHistogram("hhrand",sigma1,Binning(35,0.25,0.65),YVar(sig1frac,Binning(35,0.3,1.1))) ;

  regTH(hcorr,"rf607_hcorr") ;
  regTH(hhrand,"rf607_hhand") ;

  delete data ;
  delete r ;

  return kTRUE ;

  }
} ;


//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #609
// 
// Setting up a chi^2 fit to an unbinned dataset with X,Y,err(Y)
// values (and optionally err(X) values)
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////


#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooPolyVar.h"
#include "RooChi2Var.h"
#include "RooMinuit.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TRandom.h"

using namespace RooFit ;

class TestBasic609 : public RooUnitTest
{
public: 
  TestBasic609(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Chi^2 fit to X-Y dataset",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   d a t a s e t   w i t h   X   a n d   Y   v a l u e s
  // -------------------------------------------------------------------

  // Make weighted XY dataset with asymmetric errors stored
  // The StoreError() argument is essential as it makes
  // the dataset store the error in addition to the values
  // of the observables. If errors on one or more observables
  // are asymmetric, one can store the asymmetric error
  // using the StoreAsymError() argument

  RooRealVar x("x","x",-11,11) ;
  RooRealVar y("y","y",-10,200) ;
  RooDataSet dxy("dxy","dxy",RooArgSet(x,y),StoreError(RooArgSet(x,y))) ;

  // Fill an example dataset with X,err(X),Y,err(Y) values
  for (int i=0 ; i<=10 ; i++) {

    // Set X value and error
    x = -10 + 2*i;
    x.setError( i<5 ? 0.5/1. : 1.0/1. ) ;
    
    // Set Y value and error 
    y = x.getVal() * x.getVal() + 4*fabs(RooRandom::randomGenerator()->Gaus()) ;
    y.setError(sqrt(y.getVal())) ;

    dxy.add(RooArgSet(x,y)) ;    
  }



  // P e r f o r m   c h i 2   f i t   t o   X + / - d x   a n d   Y + / - d Y   v a l u e s
  // ---------------------------------------------------------------------------------------

  // Make fit function
  RooRealVar a("a","a",0.0,-10,10) ;
  RooRealVar b("b","b",0.0,-100,100) ;
  RooPolyVar f("f","f",x,RooArgList(b,a,RooConst(1))) ;

  // Plot dataset in X-Y interpretation
  RooPlot* frame = x.frame(Title("Chi^2 fit of function set of (X#pmdX,Y#pmdY) values")) ;
  dxy.plotOnXY(frame,YVar(y)) ;

  // Fit chi^2 using X and Y errors 
  f.chi2FitTo(dxy,YVar(y)) ;

  // Overlay fitted function
  f.plotOn(frame) ;
  
  // Alternative: fit chi^2 integrating f(x) over ranges defined by X errors, rather
  // than taking point at center of bin
  f.chi2FitTo(dxy,YVar(y),Integrate(kTRUE)) ;

  // Overlay alternate fit result
  f.plotOn(frame,LineStyle(kDashed),LineColor(kRed),Name("alternate")) ;  

  regPlot(frame,"rf609_frame") ;  

  return kTRUE ;
  }
} ;



//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #701
// 
// Unbinned maximum likelihood fit of an efficiency eff(x) function to 
// a dataset D(x,cut), where cut is a category encoding a selection, of which
// the efficiency as function of x should be described by eff(x)
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooFormulaVar.h"
#include "RooProdPdf.h"
#include "RooEfficiency.h"
#include "RooPolynomial.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


// Elementary operations on a gaussian PDF
class TestBasic701 : public RooUnitTest
{
public: 
  TestBasic701(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Efficiency operator p.d.f. 1D",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C o n s t r u c t   e f f i c i e n c y   f u n c t i o n   e ( x ) 
  // -------------------------------------------------------------------

  // Declare variables x,mean,sigma with associated name, title, initial value and allowed range
  RooRealVar x("x","x",-10,10) ;

  // Efficiency function eff(x;a,b) 
  RooRealVar a("a","a",0.4,0,1) ;
  RooRealVar b("b","b",5) ;
  RooRealVar c("c","c",-1,-10,10) ;
  RooFormulaVar effFunc("effFunc","(1-a)+a*cos((x-c)/b)",RooArgList(a,b,c,x)) ;



  // C o n s t r u c t   c o n d i t i o n a l    e f f i c i e n c y   p d f   E ( c u t | x ) 
  // ------------------------------------------------------------------------------------------

  // Acceptance state cut (1 or 0)
  RooCategory cut("cut","cutr") ;
  cut.defineType("accept",1) ;
  cut.defineType("reject",0) ;

  // Construct efficiency p.d.f eff(cut|x)
  RooEfficiency effPdf("effPdf","effPdf",effFunc,cut,"accept") ;



  // G e n e r a t e   d a t a   ( x ,   c u t )   f r o m   a   t o y   m o d e l 
  // -----------------------------------------------------------------------------

  // Construct global shape p.d.f shape(x) and product model(x,cut) = eff(cut|x)*shape(x) 
  // (These are _only_ needed to generate some toy MC here to be used later)
  RooPolynomial shapePdf("shapePdf","shapePdf",x,RooConst(-0.095)) ;
  RooProdPdf model("model","model",shapePdf,Conditional(effPdf,cut)) ;

  // Generate some toy data from model
  RooDataSet* data = model.generate(RooArgSet(x,cut),10000) ;



  // F i t   c o n d i t i o n a l   e f f i c i e n c y   p d f   t o   d a t a 
  // --------------------------------------------------------------------------

  // Fit conditional efficiency p.d.f to data
  effPdf.fitTo(*data,ConditionalObservables(x)) ;



  // P l o t   f i t t e d ,   d a t a   e f f i c i e n c y  
  // --------------------------------------------------------

  // Plot distribution of all events and accepted fraction of events on frame
  RooPlot* frame1 = x.frame(Bins(20),Title("Data (all, accepted)")) ;
  data->plotOn(frame1) ;
  data->plotOn(frame1,Cut("cut==cut::accept"),MarkerColor(kRed),LineColor(kRed)) ;

  // Plot accept/reject efficiency on data overlay fitted efficiency curve
  RooPlot* frame2 = x.frame(Bins(20),Title("Fitted efficiency")) ;
  data->plotOn(frame2,Efficiency(cut)) ; // needs ROOT version >= 5.21
  effFunc.plotOn(frame2,LineColor(kRed)) ;


  regPlot(frame1,"rf701_plot1") ;
  regPlot(frame2,"rf701_plot2") ;

  
  delete data ;

  return kTRUE ;  
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #702
// 
// Unbinned maximum likelihood fit of an efficiency eff(x) function to 
// a dataset D(x,cut), where cut is a category encoding a selection whose 
// efficiency as function of x should be described by eff(x)
//
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooCategory.h"
#include "RooEfficiency.h"
#include "RooPolynomial.h"
#include "RooProdPdf.h"
#include "RooFormulaVar.h"
#include "TCanvas.h"
#include "TH1.h"
#include "RooPlot.h"
using namespace RooFit ;


// Elementary operations on a gaussian PDF
class TestBasic702 : public RooUnitTest
{
public: 
  TestBasic702(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Efficiency operator p.d.f. 2D",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  Bool_t flat=kFALSE ;

  // C o n s t r u c t   e f f i c i e n c y   f u n c t i o n   e ( x , y ) 
  // -----------------------------------------------------------------------

  // Declare variables x,mean,sigma with associated name, title, initial value and allowed range
  RooRealVar x("x","x",-10,10) ;
  RooRealVar y("y","y",-10,10) ;

  // Efficiency function eff(x;a,b) 
  RooRealVar ax("ax","ay",0.6,0,1) ;
  RooRealVar bx("bx","by",5) ;
  RooRealVar cx("cx","cy",-1,-10,10) ;

  RooRealVar ay("ay","ay",0.2,0,1) ;
  RooRealVar by("by","by",5) ;
  RooRealVar cy("cy","cy",-1,-10,10) ;

  RooFormulaVar effFunc("effFunc","((1-ax)+ax*cos((x-cx)/bx))*((1-ay)+ay*cos((y-cy)/by))",RooArgList(ax,bx,cx,x,ay,by,cy,y)) ; 

  // Acceptance state cut (1 or 0)
  RooCategory cut("cut","cutr") ;
  cut.defineType("accept",1) ;
  cut.defineType("reject",0) ;



  // C o n s t r u c t   c o n d i t i o n a l    e f f i c i e n c y   p d f   E ( c u t | x , y ) 
  // ---------------------------------------------------------------------------------------------

  // Construct efficiency p.d.f eff(cut|x)
  RooEfficiency effPdf("effPdf","effPdf",effFunc,cut,"accept") ;



  // G e n e r a t e   d a t a   ( x , y , c u t )   f r o m   a   t o y   m o d e l 
  // -------------------------------------------------------------------------------

  // Construct global shape p.d.f shape(x) and product model(x,cut) = eff(cut|x)*shape(x) 
  // (These are _only_ needed to generate some toy MC here to be used later)
  RooPolynomial shapePdfX("shapePdfX","shapePdfX",x,RooConst(flat?0:-0.095)) ;
  RooPolynomial shapePdfY("shapePdfY","shapePdfY",y,RooConst(flat?0:+0.095)) ;
  RooProdPdf shapePdf("shapePdf","shapePdf",RooArgSet(shapePdfX,shapePdfY)) ;
  RooProdPdf model("model","model",shapePdf,Conditional(effPdf,cut)) ;

  // Generate some toy data from model
  RooDataSet* data = model.generate(RooArgSet(x,y,cut),10000) ;



  // F i t   c o n d i t i o n a l   e f f i c i e n c y   p d f   t o   d a t a 
  // --------------------------------------------------------------------------

  // Fit conditional efficiency p.d.f to data
  effPdf.fitTo(*data,ConditionalObservables(RooArgSet(x,y))) ;



  // P l o t   f i t t e d ,   d a t a   e f f i c i e n c y  
  // --------------------------------------------------------

  // Make 2D histograms of all data, selected data and efficiency function
  TH1* hh_data_all = data->createHistogram("hh_data_all",x,Binning(8),YVar(y,Binning(8))) ;
  TH1* hh_data_sel = data->createHistogram("hh_data_sel",x,Binning(8),YVar(y,Binning(8)),Cut("cut==cut::accept")) ;
  TH1* hh_eff      = effFunc.createHistogram("hh_eff",x,Binning(50),YVar(y,Binning(50))) ;

  // Some adjustsment for good visualization
  hh_data_all->SetMinimum(0) ;
  hh_data_sel->SetMinimum(0) ;
  hh_eff->SetMinimum(0) ;
  hh_eff->SetLineColor(kBlue) ;

  regTH(hh_data_all,"rf702_hh_data_all") ;
  regTH(hh_data_sel,"rf702_hh_data_sel") ;
  regTH(hh_eff,"rf702_hh_eff") ;
  
  delete data ;

  return kTRUE;

  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #703
// 
// Using a product of an (acceptance) efficiency and a p.d.f as p.d.f.
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooExponential.h"
#include "RooEffProd.h"
#include "RooFormulaVar.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


// Elementary operations on a gaussian PDF
class TestBasic703 : public RooUnitTest
{
public: 
  TestBasic703(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Efficiency product operator p.d.f",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

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

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



  // D e f i n e   e f f i c i e n c y   f u n c t i o n
  // ---------------------------------------------------
  
  // Use error function to simulate turn-on slope
  RooFormulaVar eff("eff","0.5*(TMath::Erf((t-1)/0.5)+1)",t) ;



  // D e f i n e   d e c a y   p d f   w i t h   e f f i c i e n c y 
  // ---------------------------------------------------------------

  // Multiply pdf(t) with efficiency in t
  RooEffProd modelEff("modelEff","model with efficiency",model,eff) ;

  

  // P l o t   e f f i c i e n c y ,   p d f  
  // ----------------------------------------

  RooPlot* frame1 = t.frame(Title("Efficiency")) ;
  eff.plotOn(frame1,LineColor(kRed)) ;

  RooPlot* frame2 = t.frame(Title("Pdf with and without efficiency")) ;

  model.plotOn(frame2,LineStyle(kDashed)) ;
  modelEff.plotOn(frame2) ;



  // G e n e r a t e   t o y   d a t a ,   f i t   m o d e l E f f   t o   d a t a
  // ------------------------------------------------------------------------------

  // Generate events. If the input pdf has an internal generator, the internal generator
  // is used and an accept/reject sampling on the efficiency is applied. 
  RooDataSet* data = modelEff.generate(t,10000) ;

  // Fit pdf. The normalization integral is calculated numerically. 
  modelEff.fitTo(*data) ;

  // Plot generated data and overlay fitted pdf
  RooPlot* frame3 = t.frame(Title("Fitted pdf with efficiency")) ;
  data->plotOn(frame3) ;
  modelEff.plotOn(frame3) ;

  
  regPlot(frame1,"rf703_plot1") ;
  regPlot(frame2,"rf703_plot2") ;
  regPlot(frame3,"rf703_plot3") ;


  delete data ;
  return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #704
// 
// Using a p.d.f defined by a sum of real-valued amplitude components
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooTruthModel.h"
#include "RooFormulaVar.h"
#include "RooRealSumPdf.h"
#include "RooPolyVar.h"
#include "RooProduct.h"
#include "TH1.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


// Elementary operations on a gaussian PDF
class TestBasic704 : public RooUnitTest
{
public: 
  TestBasic704(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Amplitude sum operator p.d.f",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // S e t u p   2 D   a m p l i t u d e   f u n c t i o n s
  // -------------------------------------------------------

  // Observables
  RooRealVar t("t","time",-1.,15.);
  RooRealVar cosa("cosa","cos(alpha)",-1.,1.);

  // Use RooTruthModel to obtain compiled implementation of sinh/cosh modulated decay functions
  RooRealVar tau("tau","#tau",1.5);  
  RooRealVar deltaGamma("deltaGamma","deltaGamma", 0.3);
  RooTruthModel tm("tm","tm",t) ;
  RooFormulaVar coshGBasis("coshGBasis","exp(-@0/ @1)*cosh(@0*@2/2)",RooArgList(t,tau,deltaGamma));
  RooFormulaVar sinhGBasis("sinhGBasis","exp(-@0/ @1)*sinh(@0*@2/2)",RooArgList(t,tau,deltaGamma));
  RooAbsReal* coshGConv = tm.convolution(&coshGBasis,&t);
  RooAbsReal* sinhGConv = tm.convolution(&sinhGBasis,&t);
    
  // Construct polynomial amplitudes in cos(a) 
  RooPolyVar poly1("poly1","poly1",cosa,RooArgList(RooConst(0.5),RooConst(0.2),RooConst(0.2)),0);
  RooPolyVar poly2("poly2","poly2",cosa,RooArgList(RooConst(1),RooConst(-0.2),RooConst(3)),0);

  // Construct 2D amplitude as uncorrelated product of amp(t)*amp(cosa)
  RooProduct  ampl1("ampl1","amplitude 1",RooArgSet(poly1,*coshGConv));
  RooProduct  ampl2("ampl2","amplitude 2",RooArgSet(poly2,*sinhGConv));



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

  // Amplitude strengths
  RooRealVar f1("f1","f1",1,0,2) ;
  RooRealVar f2("f2","f2",0.5,0,2) ;
  
  // Construct pdf
  RooRealSumPdf pdf("pdf","pdf",RooArgList(ampl1,ampl2),RooArgList(f1,f2)) ;

  // Generate some toy data from pdf
  RooDataSet* data = pdf.generate(RooArgSet(t,cosa),10000);

  // Fit pdf to toy data with only amplitude strength floating
  pdf.fitTo(*data) ;



  // P l o t   a m p l i t u d e   s u m   p d f 
  // -------------------------------------------

  // Make 2D plots of amplitudes
  TH1* hh_cos = ampl1.createHistogram("hh_cos",t,Binning(50),YVar(cosa,Binning(50))) ;
  TH1* hh_sin = ampl2.createHistogram("hh_sin",t,Binning(50),YVar(cosa,Binning(50))) ;
  hh_cos->SetLineColor(kBlue) ;
  hh_sin->SetLineColor(kBlue) ;

  
  // Make projection on t, plot data, pdf and its components
  // Note component projections may be larger than sum because amplitudes can be negative
  RooPlot* frame1 = t.frame();
  data->plotOn(frame1);
  pdf.plotOn(frame1);
  pdf.plotOn(frame1,Components(ampl1),LineStyle(kDashed));
  pdf.plotOn(frame1,Components(ampl2),LineStyle(kDashed),LineColor(kRed));
  
  // Make projection on cosa, plot data, pdf and its components
  // Note that components projection may be larger than sum because amplitudes can be negative
  RooPlot* frame2 = cosa.frame();
  data->plotOn(frame2);
  pdf.plotOn(frame2);
  pdf.plotOn(frame2,Components(ampl1),LineStyle(kDashed));
  pdf.plotOn(frame2,Components(ampl2),LineStyle(kDashed),LineColor(kRed));
  

  regPlot(frame1,"rf704_plot1") ;
  regPlot(frame2,"rf704_plot2") ;
  regTH(hh_cos,"rf704_hh_cos") ;
  regTH(hh_sin,"rf704_hh_sin") ;

  delete data ;
  delete coshGConv ;
  delete sinhGConv ;

  return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #705
// 
// Linear interpolation between p.d.f shapes using the 'Alex Read' algorithm
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooIntegralMorph.h"
#include "RooNLLVar.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TH1.h"
using namespace RooFit ;


// Elementary operations on a gaussian PDF
class TestBasic705 : public RooUnitTest
{
public: 

  Double_t ctol() { return 5e-2 ; } // very conservative, this is a numerically difficult test

  TestBasic705(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Linear morph operator p.d.f.",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // 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),Name("alt1")) ;
  alpha.setVal(0.25) ;
  lmorph.plotOn(frame1,LineColor(kRed),Name("alt2")) ;
  alpha.setVal(0.375) ;
  lmorph.plotOn(frame1,LineColor(kRed),Name("alt3")) ;
  alpha.setVal(0.50) ;
  lmorph.plotOn(frame1,LineColor(kRed),Name("alt4")) ;
  alpha.setVal(0.625) ;
  lmorph.plotOn(frame1,LineColor(kRed),Name("alt5")) ;
  alpha.setVal(0.75) ;
  lmorph.plotOn(frame1,LineColor(kRed),Name("alt6")) ;
  alpha.setVal(0.875) ;
  lmorph.plotOn(frame1,LineColor(kRed),Name("alt7")) ;
  alpha.setVal(0.95) ;
  lmorph.plotOn(frame1,LineColor(kRed),Name("alt8")) ;



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

  // 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.5,0.9)) ;
  
  // Make 2D pdf of histogram  
  RooNLLVar nll("nll","nll",lmorph,*data) ;  
  nll.plotOn(frame3,ShiftToZero()) ;    

  lmorph.setCacheAlpha(kFALSE) ;


  regPlot(frame1,"rf705_plot1") ;
  regPlot(frame2,"rf705_plot2") ;
  regPlot(frame3,"rf705_plot3") ;
  regTH(hh,"rf705_hh") ;

  delete data ;
  
  return kTRUE;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #706
// 
// Histogram based p.d.f.s and functions
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooHistPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;


// Elementary operations on a gaussian PDF
class TestBasic706 : public RooUnitTest
{
public: 
  TestBasic706(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Histogram based p.d.f.s",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   p d f   f o r   s a m p l i n g 
  // ---------------------------------------------

  RooRealVar x("x","x",0,20) ;
  RooPolynomial p("p","p",x,RooArgList(RooConst(0.01),RooConst(-0.01),RooConst(0.0004))) ;



  // C r e a t e   l o w   s t a t s   h i s t o g r a m
  // ---------------------------------------------------

  // Sample 500 events from p
  x.setBins(20) ;
  RooDataSet* data1 = p.generate(x,500) ;
  
  // Create a binned dataset with 20 bins and 500 events
  RooDataHist* hist1 = data1->binnedClone() ;

  // Represent data in dh as pdf in x
  RooHistPdf histpdf1("histpdf1","histpdf1",x,*hist1,0) ;

  // Plot unbinned data and histogram pdf overlaid
  RooPlot* frame1 = x.frame(Title("Low statistics histogram pdf"),Bins(100)) ;
  data1->plotOn(frame1) ;
  histpdf1.plotOn(frame1) ;  
  

  // C r e a t e   h i g h   s t a t s   h i s t o g r a m
  // -----------------------------------------------------

  // Sample 100000 events from p
  x.setBins(10) ;
  RooDataSet* data2 = p.generate(x,100000) ;

  // Create a binned dataset with 10 bins and 100K events
  RooDataHist* hist2 = data2->binnedClone() ;

  // Represent data in dh as pdf in x, apply 2nd order interpolation  
  RooHistPdf histpdf2("histpdf2","histpdf2",x,*hist2,2) ;

  // Plot unbinned data and histogram pdf overlaid
  RooPlot* frame2 = x.frame(Title("High stats histogram pdf with interpolation"),Bins(100)) ;
  data2->plotOn(frame2) ;
  histpdf2.plotOn(frame2) ;  


  regPlot(frame1,"rf607_plot1") ;
  regPlot(frame2,"rf607_plot2") ;

  delete data1 ;
  delete hist1 ;
  delete data2 ;
  delete hist2 ;

  return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #707
// 
// Using non-parametric (multi-dimensional) kernel estimation p.d.f.s
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooKeysPdf.h"
#include "RooNDKeysPdf.h"
#include "RooProdPdf.h"
#include "TCanvas.h"
#include "TH1.h"
#include "RooPlot.h"
using namespace RooFit ;


// Elementary operations on a gaussian PDF
class TestBasic707 : public RooUnitTest
{
public: 
  TestBasic707(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Kernel estimation p.d.f.s",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   l o w   s t a t s   1 - D   d a t a s e t 
  // -------------------------------------------------------

  // Create a toy pdf for sampling
  RooRealVar x("x","x",0,20) ;
  RooPolynomial p("p","p",x,RooArgList(RooConst(0.01),RooConst(-0.01),RooConst(0.0004))) ;

  // Sample 500 events from p
  RooDataSet* data1 = p.generate(x,200) ;



  // C r e a t e   1 - D   k e r n e l   e s t i m a t i o n   p d f
  // ---------------------------------------------------------------

  // Create adaptive kernel estimation pdf. In this configuration the input data
  // is mirrored over the boundaries to minimize edge effects in distribution
  // that do not fall to zero towards the edges
  RooKeysPdf kest1("kest1","kest1",x,*data1,RooKeysPdf::MirrorBoth) ;

  // An adaptive kernel estimation pdf on the same data without mirroring option
  // for comparison
  RooKeysPdf kest2("kest2","kest2",x,*data1,RooKeysPdf::NoMirror) ;


  // Adaptive kernel estimation pdf with increased bandwidth scale factor
  // (promotes smoothness over detail preservation)
  RooKeysPdf kest3("kest3","kest3",x,*data1,RooKeysPdf::MirrorBoth,2) ;


  // Plot kernel estimation pdfs with and without mirroring over data
  RooPlot* frame = x.frame(Title("Adaptive kernel estimation pdf with and w/o mirroring"),Bins(20)) ;
  data1->plotOn(frame) ;
  kest1.plotOn(frame) ;    
  kest2.plotOn(frame,LineStyle(kDashed),LineColor(kRed)) ;    


  // Plot kernel estimation pdfs with regular and increased bandwidth
  RooPlot* frame2 = x.frame(Title("Adaptive kernel estimation pdf with regular, increased bandwidth")) ;
  kest1.plotOn(frame2) ;    
  kest3.plotOn(frame2,LineColor(kMagenta)) ;    



  // C r e a t e   l o w   s t a t s   2 - D   d a t a s e t 
  // -------------------------------------------------------

  // Construct a 2D toy pdf for sampleing
  RooRealVar y("y","y",0,20) ;
  RooPolynomial py("py","py",y,RooArgList(RooConst(0.01),RooConst(0.01),RooConst(-0.0004))) ;
  RooProdPdf pxy("pxy","pxy",RooArgSet(p,py)) ;
  RooDataSet* data2 = pxy.generate(RooArgSet(x,y),1000) ;



  // C r e a t e   2 - D   k e r n e l   e s t i m a t i o n   p d f
  // ---------------------------------------------------------------

  // Create 2D adaptive kernel estimation pdf with mirroring 
  RooNDKeysPdf kest4("kest4","kest4",RooArgSet(x,y),*data2,"am") ;

  // Create 2D adaptive kernel estimation pdf with mirroring and double bandwidth
  RooNDKeysPdf kest5("kest5","kest5",RooArgSet(x,y),*data2,"am",2) ;

  // Create a histogram of the data
  TH1* hh_data = data2->createHistogram("hh_data",x,Binning(10),YVar(y,Binning(10))) ;

  // Create histogram of the 2d kernel estimation pdfs
  TH1* hh_pdf = kest4.createHistogram("hh_pdf",x,Binning(25),YVar(y,Binning(25))) ;
  TH1* hh_pdf2 = kest5.createHistogram("hh_pdf2",x,Binning(25),YVar(y,Binning(25))) ;
  hh_pdf->SetLineColor(kBlue) ;
  hh_pdf2->SetLineColor(kMagenta) ;
    
  regPlot(frame,"rf707_plot1") ;
  regPlot(frame2,"rf707_plot2") ;
  regTH(hh_data,"rf707_hhdata") ;
  regTH(hh_pdf,"rf707_hhpdf") ;
  regTH(hh_pdf2,"rf707_hhpdf2") ;
  
  delete data1 ;
  delete data2 ;

  return kTRUE ;
  }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #708
// 
// Special decay pdf for B physics with mixing and/or CP violation
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooCategory.h"
#include "RooBMixDecay.h"
#include "RooBCPEffDecay.h"
#include "RooBDecay.h"
#include "RooFormulaVar.h"
#include "RooTruthModel.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;

// Elementary operations on a gaussian PDF
class TestBasic708 : public RooUnitTest
{
public: 
  TestBasic708(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("B Physics p.d.f.s",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  ////////////////////////////////////////////////////
  // B - D e c a y   w i t h   m i x i n g          //
  ////////////////////////////////////////////////////

  // C o n s t r u c t   p d f 
  // -------------------------
  
  // Observable
  RooRealVar dt("dt","dt",-10,10) ;
  dt.setBins(40) ;

  // Parameters
  RooRealVar dm("dm","delta m(B0)",0.472) ;
  RooRealVar tau("tau","tau (B0)",1.547) ;
  RooRealVar w("w","flavour mistag rate",0.1) ;
  RooRealVar dw("dw","delta mistag rate for B0/B0bar",0.1) ;

  RooCategory mixState("mixState","B0/B0bar mixing state") ;
  mixState.defineType("mixed",-1) ;
  mixState.defineType("unmixed",1) ;

  RooCategory tagFlav("tagFlav","Flavour of the tagged B0") ;
  tagFlav.defineType("B0",1) ;
  tagFlav.defineType("B0bar",-1) ;

  // Use delta function resolution model
  RooTruthModel tm("tm","truth model",dt) ;

  // Construct Bdecay with mixing
  RooBMixDecay bmix("bmix","decay",dt,mixState,tagFlav,tau,dm,w,dw,tm,RooBMixDecay::DoubleSided) ;



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

  // Generate some data
  RooDataSet* data = bmix.generate(RooArgSet(dt,mixState,tagFlav),10000) ;

  // Plot B0 and B0bar tagged data separately
  // For all plots below B0 and B0 tagged data will look somewhat differently
  // if the flavor tagging mistag rate for B0 and B0 is different (i.e. dw!=0)
  RooPlot* frame1 = dt.frame(Title("B decay distribution with mixing (B0/B0bar)")) ;

  data->plotOn(frame1,Cut("tagFlav==tagFlav::B0")) ;
  bmix.plotOn(frame1,Slice(tagFlav,"B0")) ;

  data->plotOn(frame1,Cut("tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
  bmix.plotOn(frame1,Slice(tagFlav,"B0bar"),LineColor(kCyan),Name("alt")) ;


  // Plot mixed slice for B0 and B0bar tagged data separately
  RooPlot* frame2 = dt.frame(Title("B decay distribution of mixed events (B0/B0bar)")) ;

  data->plotOn(frame2,Cut("mixState==mixState::mixed&&tagFlav==tagFlav::B0")) ;
  bmix.plotOn(frame2,Slice(tagFlav,"B0"),Slice(mixState,"mixed")) ;

  data->plotOn(frame2,Cut("mixState==mixState::mixed&&tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
  bmix.plotOn(frame2,Slice(tagFlav,"B0bar"),Slice(mixState,"mixed"),LineColor(kCyan),Name("alt")) ;


  // Plot unmixed slice for B0 and B0bar tagged data separately
  RooPlot* frame3 = dt.frame(Title("B decay distribution of unmixed events (B0/B0bar)")) ;

  data->plotOn(frame3,Cut("mixState==mixState::unmixed&&tagFlav==tagFlav::B0")) ;
  bmix.plotOn(frame3,Slice(tagFlav,"B0"),Slice(mixState,"unmixed")) ;

  data->plotOn(frame3,Cut("mixState==mixState::unmixed&&tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
  bmix.plotOn(frame3,Slice(tagFlav,"B0bar"),Slice(mixState,"unmixed"),LineColor(kCyan),Name("alt")) ;





  ///////////////////////////////////////////////////////
  // B - D e c a y   w i t h   C P   v i o l a t i o n //
  ///////////////////////////////////////////////////////

  // C o n s t r u c t   p d f 
  // -------------------------
  
  // Additional parameters needed for B decay with CPV
  RooRealVar CPeigen("CPeigen","CP eigen value",-1) ;
  RooRealVar absLambda("absLambda","|lambda|",1,0,2) ;
  RooRealVar argLambda("absLambda","|lambda|",0.7,-1,1) ;
  RooRealVar effR("effR","B0/B0bar reco efficiency ratio",1) ;

  // Construct Bdecay with CP violation
  RooBCPEffDecay bcp("bcp","bcp", dt, tagFlav, tau, dm, w, CPeigen, absLambda, argLambda, effR, dw, tm, RooBCPEffDecay::DoubleSided) ;
      


  // P l o t   s c e n a r i o   1   -   s i n ( 2 b )   =   0 . 7 ,   | l | = 1 
  // ---------------------------------------------------------------------------

  // Generate some data
  RooDataSet* data2 = bcp.generate(RooArgSet(dt,tagFlav),10000) ;

  // Plot B0 and B0bar tagged data separately
  RooPlot* frame4 = dt.frame(Title("B decay distribution with CPV(|l|=1,Im(l)=0.7) (B0/B0bar)")) ;

  data2->plotOn(frame4,Cut("tagFlav==tagFlav::B0")) ;
  bcp.plotOn(frame4,Slice(tagFlav,"B0")) ;

  data2->plotOn(frame4,Cut("tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
  bcp.plotOn(frame4,Slice(tagFlav,"B0bar"),LineColor(kCyan),Name("alt")) ;



  // P l o t   s c e n a r i o   2   -   s i n ( 2 b )   =   0 . 7 ,   | l | = 0 . 7 
  // -------------------------------------------------------------------------------

  absLambda=0.7 ;

  // Generate some data
  RooDataSet* data3 = bcp.generate(RooArgSet(dt,tagFlav),10000) ;

  // Plot B0 and B0bar tagged data separately (sin2b = 0.7 plus direct CPV |l|=0.5)
  RooPlot* frame5 = dt.frame(Title("B decay distribution with CPV(|l|=0.7,Im(l)=0.7) (B0/B0bar)")) ;  

  data3->plotOn(frame5,Cut("tagFlav==tagFlav::B0")) ;
  bcp.plotOn(frame5,Slice(tagFlav,"B0")) ;

  data3->plotOn(frame5,Cut("tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
  bcp.plotOn(frame5,Slice(tagFlav,"B0bar"),LineColor(kCyan),Name("alt")) ;



  //////////////////////////////////////////////////////////////////////////////////
  // G e n e r i c   B   d e c a y  w i t h    u s e r   c o e f f i c i e n t s  //
  //////////////////////////////////////////////////////////////////////////////////

  // C o n s t r u c t   p d f 
  // -------------------------
  
  // Model parameters
  RooRealVar DGbG("DGbG","DGamma/GammaAvg",0.5,-1,1);
  RooRealVar Adir("Adir","-[1-abs(l)**2]/[1+abs(l)**2]",0);
  RooRealVar Amix("Amix","2Im(l)/[1+abs(l)**2]",0.7);
  RooRealVar Adel("Adel","2Re(l)/[1+abs(l)**2]",0.7);
  
  // Derived input parameters for pdf
  RooFormulaVar DG("DG","Delta Gamma","@1/@0",RooArgList(tau,DGbG));
  
  // Construct coefficient functions for sin,cos,sinh modulations of decay distribution
  RooFormulaVar fsin("fsin","fsin","@0*@1*(1-2*@2)",RooArgList(Amix,tagFlav,w));
  RooFormulaVar fcos("fcos","fcos","@0*@1*(1-2*@2)",RooArgList(Adir,tagFlav,w));
  RooFormulaVar fsinh("fsinh","fsinh","@0",RooArgList(Adel));
  
  // Construct generic B decay pdf using above user coefficients
  RooBDecay bcpg("bcpg","bcpg",dt,tau,DG,RooConst(1),fsinh,fcos,fsin,dm,tm,RooBDecay::DoubleSided);
  
  
  
  // P l o t   -   I m ( l ) = 0 . 7 ,   R e ( l ) = 0 . 7   | l | = 1,   d G / G = 0 . 5 
  // -------------------------------------------------------------------------------------
  
  // Generate some data
  RooDataSet* data4 = bcpg.generate(RooArgSet(dt,tagFlav),10000) ;
  
  // Plot B0 and B0bar tagged data separately 
  RooPlot* frame6 = dt.frame(Title("B decay distribution with CPV(Im(l)=0.7,Re(l)=0.7,|l|=1,dG/G=0.5) (B0/B0bar)")) ;  
  
  data4->plotOn(frame6,Cut("tagFlav==tagFlav::B0")) ;
  bcpg.plotOn(frame6,Slice(tagFlav,"B0")) ;
  
  data4->plotOn(frame6,Cut("tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
  bcpg.plotOn(frame6,Slice(tagFlav,"B0bar"),LineColor(kCyan),Name("alt")) ;
  
 
  regPlot(frame1,"rf708_plot1") ;
  regPlot(frame2,"rf708_plot2") ;
  regPlot(frame3,"rf708_plot3") ;
  regPlot(frame4,"rf708_plot4") ;
  regPlot(frame5,"rf708_plot5") ;
  regPlot(frame6,"rf708_plot6") ;

  delete data ;
  delete data2 ;
  delete data3 ;
  delete data4 ;

  return kTRUE ;
  }  
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #801
// 
// A Toy Monte Carlo study that perform cycles of
// event generation and fittting
//
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooMCStudy.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH2.h"
#include "RooFitResult.h"
#include "TStyle.h"
#include "TDirectory.h"

using namespace RooFit ;


// Elementary operations on a gaussian PDF
class TestBasic801 : public RooUnitTest
{
public: 
  TestBasic801(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Automated MC studies",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   m o d e l
  // -----------------------

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

  // Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
  RooRealVar mean("mean","mean of gaussians",5,0,10) ;
  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,-1,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 nbkg("nbkg","number of background events,",150,0,1000) ;
  RooRealVar nsig("nsig","number of signal events",150,0,1000) ;
  RooAddPdf  model("model","g1+g2+a",RooArgList(bkg,sig),RooArgList(nbkg,nsig)) ;



  // C r e a t e   m a n a g e r
  // ---------------------------

  // Instantiate RooMCStudy manager on model with x as observable and given choice of fit options
  //
  // The Silence() option kills all messages below the PROGRESS level, leaving only a single message
  // per sample executed, and any error message that occur during fitting
  //
  // The Extended() option has two effects: 
  //    1) The extended ML term is included in the likelihood and 
  //    2) A poisson fluctuation is introduced on the number of generated events 
  //
  // The FitOptions() given here are passed to the fitting stage of each toy experiment.
  // If Save() is specified, the fit result of each experiment is saved by the manager  
  //
  // A Binned() option is added in this example to bin the data between generation and fitting
  // to speed up the study at the expemse of some precision

  RooMCStudy* mcstudy = new RooMCStudy(model,x,Binned(kTRUE),Silence(),Extended(),
				       FitOptions(Save(kTRUE),PrintEvalErrors(0))) ;
  

  // G e n e r a t e   a n d   f i t   e v e n t s
  // ---------------------------------------------

  // Generate and fit 100 samples of Poisson(nExpected) events
  mcstudy->generateAndFit(100) ;



  // E x p l o r e   r e s u l t s   o f   s t u d y 
  // ------------------------------------------------

  // Make plots of the distributions of mean, the error on mean and the pull of mean
  RooPlot* frame1 = mcstudy->plotParam(mean,Bins(40)) ;
  RooPlot* frame2 = mcstudy->plotError(mean,Bins(40)) ;
  RooPlot* frame3 = mcstudy->plotPull(mean,Bins(40),FitGauss(kTRUE)) ;

  // Plot distribution of minimized likelihood
  RooPlot* frame4 = mcstudy->plotNLL(Bins(40)) ;

  regPlot(frame1,"rf801_plot1") ;
  regPlot(frame2,"rf801_plot2") ;
  regPlot(frame3,"rf801_plot3") ;
  regPlot(frame4,"rf801_plot4") ;

  delete mcstudy ;

  return kTRUE ;
  }
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #802
// 
// RooMCStudy: using separate fit and generator models, using the chi^2 calculator model 
//
// 
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooMCStudy.h"
#include "RooChi2MCSModule.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TDirectory.h"

using namespace RooFit ;


// Elementary operations on a gaussian PDF
class TestBasic802 : public RooUnitTest
{
public: 
  TestBasic802(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("MC Study with chi^2 calculator",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   m o d e l 
  // -----------------------

  // Observables, parameters
  RooRealVar x("x","x",-10,10) ;
  x.setBins(10) ;
  RooRealVar mean("mean","mean of gaussian",0) ;
  RooRealVar sigma("sigma","width of gaussian",5,1,10) ;

  // Create Gaussian pdf
  RooGaussian gauss("gauss","gaussian PDF",x,mean,sigma) ;  



  // C r e a t e   m a n a g e r  w i t h   c h i ^ 2   a d d - o n   m o d u l e
  // ----------------------------------------------------------------------------

  // Create study manager for binned likelihood fits of a Gaussian pdf in 10 bins
  RooMCStudy* mcs = new RooMCStudy(gauss,x,Silence(),Binned()) ;

  // Add chi^2 calculator module to mcs
  RooChi2MCSModule chi2mod ;
  mcs->addModule(chi2mod) ;

  // Generate 200 samples of 1000 events
  mcs->generateAndFit(200,1000) ;
  
  // Fill histograms with distributions chi2 and prob(chi2,ndf) that
  // are calculated by RooChiMCSModule

  RooRealVar* chi2 = (RooRealVar*) mcs->fitParDataSet().get()->find("chi2") ;
  RooRealVar* prob = (RooRealVar*) mcs->fitParDataSet().get()->find("prob") ;

  TH1* h_chi2  = new TH1F("h_chi2","",40,0,20) ;
  TH1* h_prob  = new TH1F("h_prob","",40,0,1) ;

  mcs->fitParDataSet().fillHistogram(h_chi2,*chi2) ; 
  mcs->fitParDataSet().fillHistogram(h_prob,*prob) ;   



  // C r e a t e   m a n a g e r  w i t h   s e p a r a t e   f i t   m o d e l 
  // ----------------------------------------------------------------------------

  // Create alternate pdf with shifted mean
  RooRealVar mean2("mean2","mean of gaussian 2",0.5) ;
  RooGaussian gauss2("gauss2","gaussian PDF2",x,mean2,sigma) ;  

  // Create study manager with separate generation and fit model. This configuration
  // is set up to generate bad fits as the fit and generator model have different means
  // and the mean parameter is not floating in the fit
  RooMCStudy* mcs2 = new RooMCStudy(gauss2,x,FitModel(gauss),Silence(),Binned()) ;

  // Add chi^2 calculator module to mcs
  RooChi2MCSModule chi2mod2 ;
  mcs2->addModule(chi2mod2) ;

  // Generate 200 samples of 1000 events
  mcs2->generateAndFit(200,1000) ;
  
  // Fill histograms with distributions chi2 and prob(chi2,ndf) that
  // are calculated by RooChiMCSModule

  TH1* h2_chi2  = new TH1F("h2_chi2","",40,0,20) ;
  TH1* h2_prob  = new TH1F("h2_prob","",40,0,1) ;
  
  mcs2->fitParDataSet().fillHistogram(h2_chi2,*chi2) ; 
  mcs2->fitParDataSet().fillHistogram(h2_prob,*prob) ;   

  h_chi2->SetLineColor(kRed) ;
  h_prob->SetLineColor(kRed) ;

  regTH(h_chi2,"rf802_hist_chi2") ;
  regTH(h2_chi2,"rf802_hist2_chi2") ;
  regTH(h_prob,"rf802_hist_prob") ;
  regTH(h2_prob,"rf802_hist2_prob") ;

  delete mcs ;
  delete mcs2 ;

  return kTRUE ;  
  }
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #803
// 
// RooMCStudy: Using the randomizer and profile likelihood add-on models
//
// 
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooMCStudy.h"
#include "RooRandomizeParamMCSModule.h"
#include "RooDLLSignificanceMCSModule.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TDirectory.h"

using namespace RooFit ;


class TestBasic803 : public RooUnitTest
{
public: 
  TestBasic803(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("MC Study with param rand. and Z calc",refFile,writeRef,verbose) {} ;
  Bool_t testCode() {

  // C r e a t e   m o d e l 
  // -----------------------

  // Simulation of signal and background of top quark decaying into
  // 3 jets with background

  // Observable
  RooRealVar mjjj("mjjj","m(3jet) (GeV)",100,85.,350.) ;

  // Signal component (Gaussian)
  RooRealVar mtop("mtop","m(top)",162) ;
  RooRealVar wtop("wtop","m(top) resolution",15.2) ;
  RooGaussian sig("sig","top signal",mjjj,mtop,wtop) ;

  // Background component (Chebychev)
  RooRealVar c0("c0","Chebychev coefficient 0",-0.846,-1.,1.) ;
  RooRealVar c1("c1","Chebychev coefficient 1", 0.112, 0.,1.) ;
  RooRealVar c2("c2","Chebychev coefficient 2", 0.076, 0.,1.) ;
  RooChebychev bkg("bkg","combinatorial background",mjjj,RooArgList(c0,c1,c2)) ;

  // Composite model
  RooRealVar nsig("nsig","number of signal events",53,0,1e3) ;
  RooRealVar nbkg("nbkg","number of background events",103,0,5e3) ;
  RooAddPdf model("model","model",RooArgList(sig,bkg),RooArgList(nsig,nbkg)) ;



  // C r e a t e   m a n a g e r
  // ---------------------------

  // Configure manager to perform binned extended likelihood fits (Binned(),Extended()) on data generated
  // with a Poisson fluctuation on Nobs (Extended())
  RooMCStudy* mcs = new RooMCStudy(model,mjjj,Binned(),Silence(),Extended(kTRUE),
				   FitOptions(Extended(kTRUE),PrintEvalErrors(-1))) ;



  // C u s t o m i z e   m a n a g e r
  // ---------------------------------

  // Add module that randomizes the summed value of nsig+nbkg 
  // sampling from a uniform distribution between 0 and 1000
  //
  // In general one can randomize a single parameter, or a 
  // sum of N parameters, using either a uniform or a Gaussian
  // distribution. Multiple randomization can be executed
  // by a single randomizer module
  
  RooRandomizeParamMCSModule randModule ;
  randModule.sampleSumUniform(RooArgSet(nsig,nbkg),50,500) ;
  mcs->addModule(randModule) ;  


  // Add profile likelihood calculation of significance. Redo each
  // fit while keeping parameter nsig fixed to zero. For each toy,
  // the difference in -log(L) of both fits is stored, as well
  // a simple significance interpretation of the delta(-logL)
  // using Dnll = 0.5 sigma^2

  RooDLLSignificanceMCSModule sigModule(nsig,0) ;
  mcs->addModule(sigModule) ;



  // R u n   m a n a g e r ,   m a k e   p l o t s
  // ---------------------------------------------

  mcs->generateAndFit(50) ;

  // Make some plots  
  RooRealVar* ngen    = (RooRealVar*) mcs->fitParDataSet().get()->find("ngen") ;
  RooRealVar* dll     = (RooRealVar*) mcs->fitParDataSet().get()->find("dll_nullhypo_nsig") ;
  RooRealVar* z       = (RooRealVar*) mcs->fitParDataSet().get()->find("significance_nullhypo_nsig") ;
  RooRealVar* nsigerr = (RooRealVar*) mcs->fitParDataSet().get()->find("nsigerr") ;

  TH1* dll_vs_ngen     = new TH2F("h_dll_vs_ngen"    ,"",40,0,500,40,0,50) ;
  TH1* z_vs_ngen       = new TH2F("h_z_vs_ngen"      ,"",40,0,500,40,0,10) ;
  TH1* errnsig_vs_ngen = new TH2F("h_nsigerr_vs_ngen","",40,0,500,40,0,30) ;
  TH1* errnsig_vs_nsig = new TH2F("h_nsigerr_vs_nsig","",40,0,200,40,0,30) ;

  mcs->fitParDataSet().fillHistogram(dll_vs_ngen,RooArgList(*ngen,*dll)) ;
  mcs->fitParDataSet().fillHistogram(z_vs_ngen,RooArgList(*ngen,*z)) ;
  mcs->fitParDataSet().fillHistogram(errnsig_vs_ngen,RooArgList(*ngen,*nsigerr)) ;
  mcs->fitParDataSet().fillHistogram(errnsig_vs_nsig,RooArgList(nsig,*nsigerr)) ;

  regTH(dll_vs_ngen,"rf803_dll_vs_ngen") ;
  regTH(z_vs_ngen,"rf803_z_vs_ngen") ;
  regTH(errnsig_vs_ngen,"rf803_errnsig_vs_ngen") ;
  regTH(errnsig_vs_nsig,"rf803_errnsig_vs_nsig") ;

  delete mcs ;

  return kTRUE ;

  }
} ;




/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #804
// 
// Using RooMCStudy on models with constrains
//
// 
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooMCStudy.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"

using namespace RooFit ;


class TestBasic804 : public RooUnitTest
{
public: 
  TestBasic804(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("MC Studies with aux. obs. constraints",refFile,writeRef,verbose) {} ;

  Double_t htol() { return 0.1 ; } // numerically very difficult test

  Bool_t testCode() {

  // C r e a t e   m o d e l   w i t h   p a r a m e t e r   c o n s t r a i n t
  // ---------------------------------------------------------------------------

  // Observable
  RooRealVar x("x","x",-10,10) ;

  // Signal component
  RooRealVar m("m","m",0,-10,10) ;
  RooRealVar s("s","s",2,0.1,10) ;
  RooGaussian g("g","g",x,m,s) ;

  // Background component
  RooPolynomial p("p","p",x) ;

  // Composite model
  RooRealVar f("f","f",0.4,0.,1.) ;
  RooAddPdf sum("sum","sum",RooArgSet(g,p),f) ;

  // Construct constraint on parameter f
  RooGaussian fconstraint("fconstraint","fconstraint",f,RooConst(0.7),RooConst(0.1)) ;

  // Multiply constraint with p.d.f
  RooProdPdf sumc("sumc","sum with constraint",RooArgSet(sum,fconstraint)) ;



  // S e t u p   t o y   s t u d y   w i t h   m o d e l
  // ---------------------------------------------------

  // Perform toy study with internal constraint on f
  //RooMCStudy mcs(sumc,x,Constrain(f),Silence(),Binned(),FitOptions(PrintLevel(-1))) ;
  RooMCStudy mcs(sumc,x,Constrain(f),Binned()) ;

  // Run 50 toys of 2000 events.  
  // Before each toy is generated, a value for the f is sampled from the constraint pdf and 
  // that value is used for the generation of that toy.
  mcs.generateAndFit(50,2000) ;

  // Make plot of distribution of generated value of f parameter
  RooRealVar* f_gen = (RooRealVar*) mcs.fitParDataSet().get()->find("f_gen") ;
  TH1* h_f_gen = new TH1F("h_f_gen","",40,0,1) ;
  mcs.fitParDataSet().fillHistogram(h_f_gen,*f_gen) ;

  // Make plot of distribution of fitted value of f parameter
  RooPlot* frame1  = mcs.plotParam(f,Bins(40),Range(0.4,1)) ;
  frame1->SetTitle("Distribution of fitted f values") ;

  // Make plot of pull distribution on f
  RooPlot* frame2 = mcs.plotPull(f,Bins(40),Range(-3,3)) ;
  frame1->SetTitle("Distribution of f pull values") ;

  regTH(h_f_gen,"rf804_h_f_gen") ;
  regPlot(frame1,"rf804_plot1") ;
  regPlot(frame2,"rf804_plot2") ;
  
  return kTRUE ;
  }
} ;