Software 3 Hands-on:   Modeling distributions in RooFit,  Combination & Reparameterization

(1)

Software 3  Hands-on:  

Modeling distributions in RooFit, 

Combination &

Reparameterization

• Wouter Verkerke, NIKHEF

(2)

Goal of this hands-on session

1.  Learn RooFit model building for distributions

–  Introduction to building unbinned probability density models

–  Introduction to building binned models with systematic uncertainties"

(template morphing)

2.  Leverage the power of workspace for combination and reinterpretation

–  Workspace combination and re-interpretation tools in RooFit"

are very powerful – you can build your own combination in minutes

•  We only have 90 minutes – so tailored approach to that

–  First some introductory slides to familiarize you with the syntax of"

RooFit model building and RooFit model usage

–  10 prepared macros that are fully functional and execute progressively complex tasks

–  You start from a functional working point – goal of your exercise time is to understand what they do and how they do it and work on some extensions and modifications of the macros

(3)

• Wouter Verkerke, UCSB

Unbinned models – Probability density functions

•  Fundamental property of any "

probability density function g(x,p):

–  Easy to construct for 1-dim. PDF"

much more effort for >1 dim.

–  RooFit automatically takes care of this

•  User supplied function need not be normalized"

1 )

, (

max

min

∫ ≡

x

x d p x

g

!

R G dx dy(p~) = 1

R G dx (p~) = 1

(4)

Variables → Parameter or observable?

•  PDF objects have no intrinsic notion of a variable"

begin a parameter or observable

•  But, PDF normalization depends on parameter/observable interpretation of variables

" " "

•  Parameter/observable interpretation is automatic and implicit when a PDF is used together with a dataset

–  All PDF variables that are member of the dataset are observables –  All other PDF variables are parameters

–  Limits are normalization range if variable is observable"

Limits are MINUIT bounds if variable is parameter

• RooGaussian b0sig(“b0sig”,”B0 sig PDF”,mass,mb0,width);

1 )

, (

max

min

∫ ≡

x

dx p

x

g • x = observable

• p = parameter

• demo11.cc

(5)

Building unbinned models

•  Complex PDFs be can be trivially composed using operator classes

–  Addition "

"

" " " " " "

–  Multiplication "

+ =

* ⁼

(6)

Building unbinned models

–  Composition (‘plug & play’) ^"

" " " " " " " " " "

–  Convolution

g(x;m,s) m(y;a

₀

,a

₁

)

=

⊗ ⁼

g(x,y;a0,a1,s)

• Possible in any PDF

No explicit support in PDF code needed

(7)

Plotting multi-dimensional unbinned PDFs

**RooPlot* xframe = x.frame() ; data->plotOn(xframe) ;**

prod->plotOn(xframe) ; xframe->Draw() ;

c->cd(2) ;

**RooPlot* yframe = y.frame() ; data->plotOn(yframe) ;**

prod->plotOn(yframe) ; yframe->Draw() ;

∫

= pdf x y dy x

f ( ) ( , )

∫

= pdf x y dx y

f ( ) ( , )

-Plotting a dataset D(x,y) versus x represents a projection over y

-To overlay PDF(x,y),

you must plot Int(dy)PDF(x,y)

- RooFit automatically takes care of this!

• RooPlot remembers dimensions of plotted datasets

(8)

Introduction to slicing

•  With multidimensional p.d.f.s it is also often useful to be able to plot a slice of a p.d.f

•  In RooFit

–  A slice is thin –  A range is thick

•  Slices mostly useful"

in discrete observables

–  A slice in a continuous observable"

has no width and usually no data"

with the corresponding cut "

(e.g. “x=5.234”)"

•  Ranges work for both"

continuous and discrete "

observables

–  Range of discrete observable"

can be list of >=1 state

x = x.getVal()

Slice in x

Range in y

(9)

Plotting a range of a p.d.f and a dataset

RooPlot* xframe = x.frame() ; data->plotOn(xframe) ;

model.plotOn(xframe) ;

y.setRange(“sig”,-1,1) ;

RooPlot* xframe2 = x.frame() ;

data->plotOn(xframe2,CutRange("sig")) ;

model.plotOn(xframe2,ProjectionRange("sig")) ;

model(x,y) = gauss(x)gauss(y) + poly(x)poly(y)

à Works also with >2D projections (just specify projection range on all projected observables) à Works also with multidimensional p.d.fs that have correlations

(10)

Likelihood ratio plots

•  Idea: use information on S/(S+B) ratio in projected observables to define a cut

•  Example: generalize previous toy model "

to 3 dimensions"

•  Express information on S/(S+B) ratio of model in terms of integrals over model components

[ ]

∫

− +

⋅

= ⋅

dx z y x B f z

y x S f

dx z y x S z f

y

LR ( , , ) ( 1 ) ( , , )

) , , ) (

, (

[ ⁽ ^, ^, ⁾ ⁽ ¹ ⁾ ⁽ ^, ^, ⁾ ]

) , , ) (

, ,

( f S x y z f B x y z

z y x S z f

y x

LR ⋅ + −

= ⋅

• Integrate over x

• Plot LR vs (y,z)

(11)

Likelihood ratio plots

•  Decide on s/(s+b) purity"

contour of LR(y,z)

–  Example s/(s+b) > 50%"

•  Plot both data and model "

with corresponding cut.

–  For data: calculate LR(y,z) for each event, plot only event with LR>0.5 –  For model: using Monte Carlo integration technique:

∫ ^≈ ∑

>0.5 ( , )

) , (

) , , 1 (

) , , (

z y D

i i z

y LR

z y x N M

dydz z

y x

M

• Dataset with values of (y,z)

sampled from p.d.f and filtered for events that meet LR(y,z)>0.5

•   All events •   Only LR(y,z)>0.5

(12)

Likelihood ratio plots – Coded example

// Construct likelihood ratio in projection on (y,z) w.factory("expr::LR('fsig*psig/ptot',fsig,

PROJ::psig(sig,x),PROJ::ptot(model,x))") ;

// Generate toy dataset for MC integration over region with LR>68%

RooDataSet* tmpdata = model.generate(RooArgSet(x,y,z),10000) ; tmpdata->addColumn(*w.function(“LR”)) ;

RooDataSet* projdata = (RooDataSet*) tmpdata->reduce(Cut("LR>0.68")) ;

// Add LR to observed data so we can cut on it data->addColumn(*w.function(“LR”)) ;

RooDataSet* seldata = (RooDataSet*) data->reduce(Cut("LR>0.68")) ;

// Make plot for data and pdf

RooPlot* frame3 = x.frame(Title("Projection with LR(y,z)>68%")) ; seldata->plotOn(frame3) ;

model.plotOn(frame3,ProjWData(*projdata)) ;

[ ]

∫

− +

⋅

= ⋅

dx z y x B f z

y x S f

dx z y x S z f

y

LR ( , , ) ( 1 ) ( , , )

) , , ) (

,

(

(13)

Special pdfs – Kernel estimation model

•  Kernel estimation model

–  Construct smooth pdf from unbinned data, "

using kernel estimation technique

•  Example

•  Also available for n-D data

Sample of events Gaussian pdf

for each event Summed pdf for all events

Adaptive Kernel:

width of Gaussian depends on local event density

w.import(myData,Rename(“myData”)) ; w.factory(“KeysPdf::k(x,myData)”) ;

38

(14)

Example 3 : binned L with syst

•  Example of template morphing"

systematic in a binned likelihood

L( 

N | α ,  s

⁻

, 

s

⁰

, 

s

⁺

) = P(N

_i

| s

_i

( α , s

_i⁻

, s

_i⁰

, s

_i⁺

)

bins

∏ ^{)⋅G(0 |} ^α ^,1)

s

_i

( α ,...) = s

_i⁰

+ α ⋅ (s

_i⁺

− s

_i⁰

) ∀ α > 0 s

_i⁰

+ α ⋅ (s

_i⁰

− s

_i⁻

) ∀ α < 0

$

% &

'&

// Import template histograms in workspace w.import(hs_0,hs_p,hs_m) ;

// Construct template models from histograms w.factory(“HistFunc::s_0(x[80,100],hs_0)”) ; w.factory(“HistFunc::s_p(x,hs_p)”) ;

w.factory(“HistFunc::s_m(x,hs_m)”) ; // Construct morphing model

w.factory(“PiecewiseInterpolation::sig(s_0,s_,m,s_p,alpha[-5,5])”) ; // Construct full model

w.factory(“PROD::model(ASUM(sig,bkg,f[0,1]),Gaussian(0,alpha,1))”) ;

(15)

Example 4 – Beeston-Barlow light

•  Beeston-Barlow-(lite) modeling"

of MC statistical uncertainties

L(  N | 

γ ) = P(N

_i

| γ

_i

( s

_i

+ b

_i

))

bins

∏ ^P( ^s

ⁱ

^{+ b}

ⁱ

^| ^γ

ⁱ

⁽ ^s

ⁱ

^{+ b}

ⁱ

bins

∏ ⁾⁾

// Import template histogram in workspace w.import(hs) ;

// Construct parametric template models from histograms // implicitly creates vector of gamma parameters

w.factory(“ParamHistFunc::s(hs)”) ;

// Product of subsidiary measurement w.factory(“HistConstraint::subs(s)”) ; // Construct full model

w.factory(“PROD::model(s,subs)”) ;

(16)

Example 5 – BB-lite + morphing

•  Template morphing model"

with Beeston-Barlow-lite"

MC statistical uncertainties

L(  N | 

s, 

b) = P(N

_i

| γ

_i

⋅[s

_i

( α , s

_i⁻

, s

_i⁰

, s

_i⁺

) + b

_i

])

bins

∏ ^P( ^s

ⁱ

^{+ b}

ⁱ

^| ^γ

ⁱ

^{⋅[ s}

ⁱ

^{+ b}

ⁱ

^]

bins

∏ ^{)G(0 |} ^α ^,1)

s

_i

( α ,...) = s

_i⁰

+ α ⋅ (s

_i⁺

− s

_i⁰

) ∀ α > 0 s

_i⁰

+ α ⋅ (s

_i⁰

− s

_i⁻

) ∀ α < 0

$

% &

'&

// Import template histograms in workspace w.import(hs_0,hs_p,hs_m,hb) ;

// Construct parametric template morphing signal model w.factory(“ParamHistFunc::s_p(hs_p)”) ;

w.factory(“HistFunc::s_m(x,hs_m)”) ;

w.factory(“HistFunc::s_0(x[80,100],hs_0)”) ;

w.factory(“PiecewiseInterpolation::sig(s_0,s_,m,s_p,alpha[-5,5])”) ; // Construct parametric background model (sharing gamma’s with s_p) w.factory(“ParamHistFunc::bkg(hb,s_p)”) ;

// Construct full model with BB-lite MC stats modeling w.factory(“PROD::model(ASUM(sig,bkg,f[0,1]),

HistConstraint({s_0,bkg}),Gaussian(0,alpha,1))”) ;

(17)

The benefits of modularity

•  Technically very straightforward to combine measurements !

! !

RooFit, or RooFit+HistFactory!

RooStats

RooWorkspace! RooWorkspace!

RooWorkspace!

Higgs channel 1 Higgs channel 2

Combiner!

RooStats!

Higgs

Combination Lightweight!

software tool!

using RooFit!

editor tools!

(~500 LOC) Insertion of !

combination !

step does not !

modify workflow !

before/after !

combination step

(18)

More benefits of modularity

•  Technically very straightforward to reparametrize measurements !

! !

RooFit, or RooFit+HistFactory!

RooStats RooWorkspace!

RooWorkspace!

Standard "

Higgs combination

Reparametrize!

RooStats!

Lightweight"

software tool"

using RooFit"

editor tools Reparametrization

step does not modify workflow

BSM"

Higgs combination

(19)

An excursion – Collaborative analyses with workspaces

•  How can you reparametrize existing Higgs likelihoods in practice?

•  Write functions expressions corresponding to new parameterization"

"

•  Import transformation in workspace, edit existing model

w.factory(“expr::mu_gg_func(‘(KF2*Kg2)/

(0.75KF2+0.25KV2)’, KF2,Kg2,KV2) ;

w.import(mu_gg_func) ;

w.factory(“EDIT::newmodel(model,mu_gg=mu_gg_gunc)”) ;

(20)

Exercise 10

•  Macro ex10_explore_unbinned.C

–  Constructs an unbinned probability density model"

consisting of an exponential background and a (small)"

Gaussian signal (comparabale like a Hàγγ signal) –  Sample an unbinned toy dataset from the model –  Perform an unbinned ML fit to the data

–  Visualize data and model, and its components

•  Questions & explorations

–  Visualize the fitted uncertainty on the model by overlaying"

a band on the nominal fit curve that corresponds to a 68% interval"

of the parameter variations as follows: Add to the plotOn() call of the pdf the argument"

"

RooFit::VisualizeError(*r,1)

You can also add this argument to the 2^nd plotOn() call that only visualizes "

the background

–  See how the fit performs with (much) smaller data samples.

–  For further tips & tricks on error visualization look at RooFit tutorial"

macro rf610_visualerror.C

(21)

Exercise 11

•  Macro ex11_build_unbinned_extended.C

–  This macro builds the same model as ex10, but then as an extended likelihood model, i.e the pdf for the distribution in mgg is multiplied"

by a Poisson model that measures the event count (which the pdf explicitly does not).

–  The model then reparametrizes w.r.t fsig and Nexpected into"

Nsig = fsigNexpected and Nbkg=(1-fsig)Nexpected for a more convenient"

interpretation

–  After plotting and fitting, a RooStats ModelConfig is constructed and the workspace with the modelconfig is written to file

•  Questions & explorations

–  Run a selection of the RooStats interval calculators on the model of ex11"

(e.g. PLR and CLS asymptotic)

–  As variation, unfreeze the signal width and/or mean parameters in the model"

and redo the fit & limit calculations (note that you also have to adjust the listed NPs in the ModelConfig!)

(22)

Exercise 12

•  Macro ex12_build_binned.C

–  This macro explores binned likelihood models and build a model"

similar to that of ex11, but now with histogram templates for"

signal and background.

–  The templates are generated on the fly from analytical pdfs"

that reside in a separate simulation workspace

–  The histograms are turned in to functions that predict the yield of signal"

and background in bins of mgg.

–  The sum of the histograms, weighted by scale factors is turned into "

a probability model that can be fit to a binned data distribution in mgg –  The model is persisted together with a ModelConfig.

•  Questions & explorations

–  Study the macro to understand the process. In particular note that"

for histograms functions, the yield corresponding to a bin is"

the bin height times the bin width (thus if a histogram was filled with the"

notion that the height is an event count) a correction factor for the bin"

width is needed when constructing the probability model

–  Run the macro through your favorite RooStats limit calculators

–  Vary the statistics of the template distribution by adjusting the generated event count in the simulation section and observe the effect on the fit. Note that when"

you adjust the scale factor expression in the probability model. To that end,"

introduce an additional real-value constant named ‘L’, with initial value 1, in the workspace and multiply the expressions for S and B with L. Then, as you increase"

the generated event count for the templates by a factor X, you set the value of L to 1/X"

to compensate for this in the interpretation

(23)

Exercise 13

•  Macro ex13_build_binned_morphing.C

–  This macro builds a similar model to that of ex12 but extends it by reintroducing a nuisance parameter that expresses uncertainty on the position of the signal peak –  To this end three signal templates, at peaking 123, 125 and 127 GeV"

and introduced and interpolated with a ‘vertical’ morphing model with"

a nuisance parameter alpha. A unit-Gaussian subsidiary measurement is introduced on alpha that expresses that the known uncertainty on the peak position is 2 GeV.

–  The model is fit to data, visualized and written to a file together with a ModelConfig

•  Questions & explorations

–  Run the macro and study its output. Is the NP alpha that expresses the uncertainty on the signal mass overconstrained?

–  How does the fit behave if you (strongly) increase or decrease the template statistics?

–  Study the effect of turning off the interpolation class’es feature that zeros event yields that come out negative in the extrapolation. This is best done"

by making a separate plot of the signal yield histogram function alone (named ‘sig’) –  Modify the macro so that the signal width is a free parameter.

–  Modify the macro to have both width and mean as free parameters. With multiple!

interpolation directions and parameters the syntax of the interpolation class becomes !

!PiecewiseInterpolation::sig(nom,{loA,loB},{hiA,hiB},{alphaA,alphaB}) ;

(24)

Exercise 16

•  Macro ex16_combined.C

–  This macro constructs a joint likelihood of two previously constructed models by reading their workspace

–  To be able to join workspaces and datasets, a discrete observable"

(like a C++ enum) must be introduced to label the parts of the pdf"

and the dataset that belong to each channel.

–  After a fit to the joint model, the joint data and models are written"

to a new joint workspace file

•  Questions & explorations

–  First construct the input models by running ex06 and ex11 and renaming"

their output files as described in the ex16 macro

–  The default combined workspace has a common parameter mu for the"

signal strength in either channel (because it was a parameter with a common name in both input workspaces). As a result the joint model will now fit for the weighted average of the two individual measurements.

–  Write a ModelConfig file for the joint model and store it in the works[pace, then run!

the RooStats PLR calculator on the joint model

–  Modify the macro so that the mu variables of the input workspaces are explicitly kept separate: in lines of code where the source models are"

imported add the argument RooFit::RenameVariable(“mu”,”mu_exNN”) ;"

where NN=06 or 11 depending on the imported channel

–  Confirm that when you fit the model with separate mu, you find back"

within numeric precision the same values as were fitted from the individual"

channels (NB: Default MINUIT numeric precision on parameter estimates is ca 4%"

of the fitted uncertainty)

(25)

Exercise 17

•  Macro ex17_reparam.C

–  Import the joint model of ex16 with separate mus"

and perform a Higgs coupling-model style reinterpretation"

of the two signal strengths under the assumption that"

one channel measures only vector boson couplings"

and the other channel measures a product of vector boson and fermion"

couplings.

•  Questions & explorations

–  Run the reparameterized fit and observe how correlations between the kappa parameters are different from the correlations between the mu parameters due to the imposed parameter structure

–  Construct a 3

^rd

channel (for example introduce a clone of ex11 with a different signal yield), use the code of ex16 to make a three-way combination, and and your favorite Higgs channels interpretation to the signal strength of that new channel

(26)

Exercise 18

•  Macro ex18_histfactory.C

–  This macro builds a binned likelihood model with 2 bins using the HistFactory framework. The model has a single channel that is composed of signal and two backgrounds, each with separate systematic uncertainties

–  The output of HistFactory is a workspace with a ModelConfig this is very simular to the models explored so far

•  Questions & explorations

–  Build the model and process the resulting model.root file with your favorite RooStats limit calculator

Software 3 Hands-on: Modeling distributions in RooFit, Combination & Reparameterization