RooFit simultaneous fit problem

Roofit and RooStats
1

Dear RooFitters,

I am trying to do a simultaneous fit over two extended likelihoods.

To test my implementation and model, I first generate some pseudo-data, and then do the fit. As you can see, the best fit values shown on the plots attached agree with the initial values I set. However, the “normalisation” in the right plot does not look right.

Could anyone tell me what I’m missing or doing wrong? I’m also confused about why “cat” appears as a parameter on the plot. My full python code is below.

Thanks a lot,
Christian

[code]from ROOT import *

invariant mass

x = RooRealVar(“x”,“x”,70,130)

efficiency

efficiency = RooRealVar(“efficiency”,“efficiency”,0.9,0,1)

breit-wigner params

mu = RooRealVar(“mu”,“mu”,90,70,130)
sigma = RooRealVar(“sigma”,“sigma”,5,0,20)

signal total, passed, failed

nsigtot = RooRealVar(“nsigtot”,“nsigtot”,10000,0,1E8)
nsigp = RooFormulaVar(“nsigp”,“nsigp”,“efficiency*nsigtot”,RooArgList(efficiency,nsigtot))
nsigf = RooFormulaVar(“nsigf”,“nsigf”,"(1-efficiency)*nsigtot",RooArgList(efficiency,nsigtot))

background

nbkgp = RooRealVar(“nbkgp”,“nbkgp”,100,0,1E5)
nbkgf = RooRealVar(“nbkgf”,“nbkgf”,1000,0,1E5)

exponential params

alphap = RooRealVar(“alphap”,“alphap”,-0.01,-10,10)
alphaf = RooRealVar(“alphaf”,“alphaf”,-0.08,-10,10)

signal pdf

signal = RooBreitWigner(“signal”,“signal”,x,mu,sigma)

background pdf

backgroundp = RooExponential(“backgroundp”,“backgroundp”,x,alphap)
backgroundf = RooExponential(“backgroundf”,“backgroundp”,x,alphaf)

extended pdf’s

totalp = RooAddPdf(“totalp”,“totalp”,RooArgList(signal,backgroundp),RooArgList(nsigp,nbkgp))
totalf = RooAddPdf(“totalf”,“totalf”,RooArgList(signal,backgroundf),RooArgList(nsigf,nbkgf))

cat = RooCategory(“cat”,“cat”)
cat.defineType(“passed”)
cat.defineType(“failed”)

joint likelihood

model = RooSimultaneous(“model”,"",cat)
model.addPdf(totalp,“passed”)
model.addPdf(totalf,“failed”)

generate data for each pdf and combine

data1 = totalp.generate(RooArgSet(x),RooFit.Extended(True))
data2 = totalf.generate(RooArgSet(x),RooFit.Extended(True))
data = RooDataSet(“data”,“combined data”,RooArgSet(x),
RooFit.Index(cat),
RooFit.Import(“passed”,data1),RooFit.Import(“failed”,data2))

fit

results = model.fitTo(data,RooFit.Extended(True))

plot

framep = x.frame(RooFit.Title(“Passed”))
data1.plotOn(framep)
model.plotOn(framep,RooFit.Slice(cat,“passed”),RooFit.ProjWData(RooArgSet(cat),data))
model.plotOn(framep,RooFit.Slice(cat,“passed”),RooFit.ProjWData(RooArgSet(cat),data),RooFit.Components(“backgroundp”),RooFit.LineColor(kRed))

framef = x.frame(RooFit.Title(“Failed”))
data2.plotOn(framef)
model.plotOn(framef,RooFit.Slice(cat,“failed”),RooFit.ProjWData(RooArgSet(cat),data))
model.plotOn(framef,RooFit.Slice(cat,“failed”),RooFit.ProjWData(RooArgSet(cat),data),RooFit.Components(“backgroundf”),RooFit.LineColor(kRed))
model.paramOn(framef)

c = TCanvas(“”,””,800,400)
c.Divide(2)
c.cd(1)
framep.Draw()
c.cd(2)
framef.Draw()
raw_input()[/code]

simulfit.png796×374 10.4 KB

2

Hi,

you are not plotting the combined dataset correctly.
The right way to achieve this is to explicitly select the category, e.g.:

data.plotOn(framep,RooFit.Cut("cat==cat::passed"))

The correct script therefore is:

from ROOT import *

# invariant mass
x = RooRealVar("x","x",70,130)

# efficiency
efficiency = RooRealVar("efficiency","efficiency",0.9,0,1)

# breit-wigner params
mu = RooRealVar("mu","mu",90,70,130)
sigma = RooRealVar("sigma","sigma",5,0,20)

# signal total, passed, failed
nsigtot = RooRealVar("nsigtot","nsigtot",10000,0,1E8)
nsigp = RooFormulaVar("nsigp","nsigp","efficiency*nsigtot",RooArgList(efficiency,nsigtot))
nsigf = RooFormulaVar("nsigf","nsigf","(1-efficiency)*nsigtot",RooArgList(efficiency,nsigtot))

# background
nbkgp = RooRealVar("nbkgp","nbkgp",100,0,1E5)
nbkgf = RooRealVar("nbkgf","nbkgf",1000,0,1E5)

# exponential params
alphap = RooRealVar("alphap","alphap",-0.01,-10,10)
alphaf = RooRealVar("alphaf","alphaf",-0.08,-10,10)

# signal pdf
signal = RooBreitWigner("signal","signal",x,mu,sigma)

# background pdf
backgroundp = RooExponential("backgroundp","backgroundp",x,alphap)
backgroundf = RooExponential("backgroundf","backgroundp",x,alphaf)

# extended pdf's
totalp = RooAddPdf("totalp","totalp",RooArgList(signal,backgroundp),RooArgList(nsigp,nbkgp))
totalf = RooAddPdf("totalf","totalf",RooArgList(signal,backgroundf),RooArgList(nsigf,nbkgf))

cat = RooCategory("cat","cat")
cat.defineType("passed")
cat.defineType("failed")

# joint likelihood
model = RooSimultaneous("model","",cat)
model.addPdf(totalp,"passed")
model.addPdf(totalf,"failed")

# generate data for each pdf and combine
data1 = totalp.generate(RooArgSet(x),RooFit.Extended(True))
data2 = totalf.generate(RooArgSet(x),RooFit.Extended(True))
data = RooDataSet("data","combined data",RooArgSet(x),
                  RooFit.Index(cat),
                  RooFit.Import("passed",data1),RooFit.Import("failed",data2))

# fit
results = model.fitTo(data,RooFit.Extended(True))

# plot
framep = x.frame(RooFit.Title("Passed"))
data.plotOn(framep,RooFit.Cut("cat==cat::passed"))
model.plotOn(framep,RooFit.Slice(cat,"passed"),RooFit.ProjWData(RooArgSet(cat),data))
model.plotOn(framep,RooFit.Slice(cat,"passed"),RooFit.ProjWData(RooArgSet(cat),data),RooFit.Components("backgroundp"),RooFit.LineColor(kRed))

framef = x.frame(RooFit.Title("Failed"))
data.plotOn(framef,RooFit.Cut("cat==cat::failed"))
model.plotOn(framef,RooFit.Slice(cat,"failed"),RooFit.ProjWData(RooArgSet(cat),data))
model.plotOn(framef,RooFit.Slice(cat,"failed"),RooFit.ProjWData(RooArgSet(cat),data),RooFit.Components("backgroundf"),RooFit.LineColor(kRed))
model.paramOn(framef)

c = TCanvas("","",800,400)
c.Divide(2)
c.cd(1)
framep.Draw()
c.cd(2)
framef.Draw()

raw_input()

One very handy utility RooFit provides is the WSFactory: root.cern.ch/doc/master/classRo … f4c6aaf8df
It allows to express complex models with short and expressive strings de facto implementing a powerful domain specific language: the programming model is slimmer - the same result can be achieved with far less lines of code. This tutorial root.cern.ch/doc/master/rf511__ … ic_8C.html also illustrates how to use the factory.

Cheers,
Danilo

3

Hi Danilo, that solved my problem, thanks a lot for your help!