Notebook

rf303_conditional¶

Multidimensional models: use of tailored pdf as conditional pdfs.s

pdf = gauss(x,f(y),sx | y ) with f(y) = a0 + a1*y

Author: Wouter Verkerke
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Sunday, February 05, 2023 at 11:14 AM.

In [1]:

%%cpp -d
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooPolyVar.h"
#include "RooProdPdf.h"
#include "RooPlot.h"
#include "TRandom.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH1.h"

using namespace RooFit;

RooDataSet *makeFakeDataXY();

Definition of a helper function:

In [2]:

%%cpp -d
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 tmpy = gRandom->Gaus(0, 10);
      double tmpx = gRandom->Gaus(0.5 * tmpy, 1);
      if (fabs(tmpy) < 10 && fabs(tmpx) < 10) {
         x.setVal(tmpx);
         y.setVal(tmpy);
         d->add(coord);
      }
   }

   return d;
}

Setup composed model gauss(x,m(y),s)¶

Create observables

In [3]:

RooRealVar x("x", "x", -10, 10);
RooRealVar y("y", "y", -10, 10);

Create function f(y) = a0 + a1*y

In [4]:

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 gauss(x,f(y),s)

In [5]:

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

In [6]:

RooDataSet *expDataXY = makeFakeDataXY();

Generate data from conditional p.d.f model(x|y)¶

Make subset of experimental data with only y values

In [7]:

RooDataSet *expDataY = (RooDataSet *)expDataXY->reduce(y);

Generate 10000 events in x obtained from conditional model(x|y) with y values taken from experimental data

In [8]:

RooDataSet *data = model.generate(x, ProtoData(*expDataY));
data->Print();

RooDataSet::modelData[x,y] = 6850 entries

input_line_54:2:2: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration
 RooDataSet *data = model.generate(x, ProtoData(*expDataY));
 ^

Fit conditional p.d.f model(x|y) to data¶

In [9]:

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

[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
 **********
 **    1 **SET PRINT           1
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 a0          -5.00000e-01  1.00000e+00   -5.00000e+00  5.00000e+00
     2 a1          -5.00000e-01  2.00000e-01   -1.00000e+00  1.00000e+00
     3 sigma        5.00000e-01  1.90000e-01    1.00000e-01  2.00000e+00
 **********
 **    3 **SET ERR         0.5
 **********
 **********
 **    4 **SET PRINT           1
 **********
 **********
 **    5 **SET STR           1
 **********
 NOW USING STRATEGY  1: TRY TO BALANCE SPEED AGAINST RELIABILITY
 **********
 **    6 **MIGRAD        1500           1
 **********
 FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
 START MIGRAD MINIMIZATION.  STRATEGY  1.  CONVERGENCE WHEN EDM .LT. 1.00e-03
 FCN=421420 FROM MIGRAD    STATUS=INITIATE       12 CALLS          13 TOTAL
                     EDM= unknown      STRATEGY= 1      NO ERROR MATRIX       
  EXT PARAMETER               CURRENT GUESS       STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  a0          -5.00000e-01   1.00000e+00   2.02430e-01  -7.46265e+04
   2  a1          -5.00000e-01   2.00000e-01   2.35352e-01  -6.95347e+05
   3  sigma        5.00000e-01   1.90000e-01   2.52163e-01  -1.29056e+06
                               ERR DEF= 0.5
 MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=9659.64 FROM MIGRAD    STATUS=CONVERGED     101 CALLS         102 TOTAL
                     EDM=8.18253e-05    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  a0           9.03100e-03   1.19768e-02   1.62548e-04  -2.31829e+00
   2  a1           5.02815e-01   2.21631e-03   1.74026e-04  -2.33580e+00
   3  sigma        9.91234e-01   8.46817e-03   6.05957e-04  -4.27913e-01
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  3    ERR DEF=0.5
  1.434e-04  1.880e-07  4.948e-08 
  1.880e-07  4.912e-06  8.995e-09 
  4.948e-08  8.995e-09  7.171e-05 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3
        1  0.00710   1.000  0.007  0.000
        2  0.00710   0.007  1.000  0.000
        3  0.00068   0.000  0.000  1.000
 **********
 **    7 **SET ERR         0.5
 **********
 **********
 **    8 **SET PRINT           1
 **********
 **********
 **    9 **HESSE        1500
 **********
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=9659.64 FROM HESSE     STATUS=OK             16 CALLS         118 TOTAL
                     EDM=8.17764e-05    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                INTERNAL      INTERNAL  
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   
   1  a0           9.03100e-03   1.19768e-02   3.25095e-05   1.80620e-03
   2  a1           5.02815e-01   2.21631e-03   3.48052e-05   2.61474e+00
   3  sigma        9.91234e-01   8.46818e-03   1.21191e-04  -6.18987e-02
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  3    ERR DEF=0.5
  1.434e-04  1.880e-07  2.151e-09 
  1.880e-07  4.912e-06  2.334e-10 
  2.151e-09  2.334e-10  7.171e-05 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3
        1  0.00708   1.000  0.007  0.000
        2  0.00708   0.007  1.000  0.000
        3  0.00002   0.000  0.000  1.000
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization

Project conditional p.d.f on x and y dimensions¶

Plot x distribution of data and projection of model on x = 1/Ndata sum(data(y_i)) model(x;y_i)

In [10]:

RooPlot *xframe = x.frame();
expDataXY->plotOn(xframe);
model.plotOn(xframe, ProjWData(*expDataY));

[#1] INFO:Plotting -- RooAbsReal::plotOn(model) plot on x averages using data variables (y)
[#1] INFO:Plotting -- RooDataWeightedAverage::ctor(modelDataWgtAvg) constructing data weighted average of function model_Norm[x] over 6850 data points of (y) with a total weight of 6850
.........................................................................................................................................................................................................

Speed up (and approximate) projection by using binned clone of data for projection

In [11]:

RooAbsData *binnedDataY = expDataY->binnedClone();
model.plotOn(xframe, ProjWData(*binnedDataY), LineColor(kCyan), LineStyle(kDotted));

[#1] INFO:Plotting -- RooAbsReal::plotOn(model) plot on x averages using data variables (y)
[#1] INFO:Plotting -- RooDataWeightedAverage::ctor(modelDataWgtAvg) constructing data weighted average of function model_Norm[x] over 100 data points of (y) with a total weight of 6850
.........................................................................................................................................................................................................

Show effect of projection with too coarse binning

In [12]:

((RooRealVar *)expDataY->get()->find("y"))->setBins(5);
RooAbsData *binnedDataY2 = expDataY->binnedClone();
model.plotOn(xframe, ProjWData(*binnedDataY2), LineColor(kRed));

[#1] INFO:Plotting -- RooAbsReal::plotOn(model) plot on x averages using data variables (y)
[#1] INFO:Plotting -- RooDataWeightedAverage::ctor(modelDataWgtAvg) constructing data weighted average of function model_Norm[x] over 5 data points of (y) with a total weight of 6850
.................................................................................................................................................................................................................................................

Make canvas and draw RooPlots

In [13]:

new TCanvas("rf303_conditional", "rf303_conditional", 600, 460);
gPad->SetLeftMargin(0.15);
xframe->GetYaxis()->SetTitleOffset(1.2);
xframe->Draw();

Draw all canvases

In [14]:

%jsroot on
gROOT->GetListOfCanvases()->Draw()