I have the following example in PyROOT which involves constructing simple background and signal models, constructing a fit model from them, and then using RooMCStudy::generateAndFit to perform some accept-reject MC and see if the fit model is able to recognize the injected signal by checking the pull distribution. This is all fine. However, I then want to access the fit results of each of the toys. But when this fails
[#0] ERROR:InputArguments -- RooMCStudy::fitResult: ERROR, invalid sample number: 1
<ROOT.RooFitResult object at 0x(nil)>
Since I am getting a null pointer here, this means that the fits aren’t being saved even though I am using the 'r' option. Can someone explain to me why, and how to fix this?
PyROOT script (also source is here: RooMCStudy_example.py):
import ROOT
from ROOT import RooFit
from ROOT import gROOT
from ROOT import TRandom3, RooRandom
from ROOT import RooRealVar, RooArgSet, RooArgList
from ROOT import RooExponential, RooGaussian, RooAddPdf, RooMCStudy
def n_bins(range, bin_width):
return int((range[1] - range[0]) / bin_width)
def quick_save(hist, name, path, log_y=False, legend_info=None):
canvas = ROOT.TCanvas()
canvas.cd()
hist.Draw()
if legend_info:
legend = ROOT.TLegend(0.5, 0.6, 0.89, 0.89)
legend.SetBorderSize(0) # no border
legend.SetFillStyle(0) # make transparent
for info in legend_info:
legend.AddEntry(None, info, '')
legend.Draw()
canvas.Draw()
extensions = ['.pdf', '.png', '.eps']
if path is None:
path = ''
for extension in extensions:
canvas.SaveAs(path + name + extension)
if log_y:
canvas.SetLogy()
for extension in extensions:
canvas.SaveAs(path + name + '_log' + extension)
if __name__ == '__main__':
gROOT.SetBatch(ROOT.kTRUE)
# Set error ignore level above kInfo
ROOT.gErrorIgnoreLevel = ROOT.kWarning
output_dir = '<PATH>'
fit_range = [60, 200]
x = RooRealVar('x', 'Large-R Jet Mass',
fit_range[0], fit_range[1], 'GeV')
# Set binning of jet mass to be 5 GeV
x.setBins(n_bins(fit_range, 5))
# background model
decay_rate = RooRealVar('decay_rate', 'exponential decay constant',
-1.5e-02, -1e-01, 1e-01)
background_model = RooExponential(
'background_model', 'background model', x, decay_rate)
# signal model
mean = RooRealVar('mean', 'Gaussian mean', 90, 50, 120)
width = RooRealVar('width', 'Gaussian width', 10, 1, 40)
signal_model = RooGaussian(
'signal_model', 'signal model', x, mean, width)
RNG = TRandom3(0)
expected_signal_events = RNG.Poisson(500)
n_signal = RooRealVar('n_signal', 'Number of signal events',
expected_signal_events,
-1e1 * expected_signal_events,
1e3 * expected_signal_events)
expected_background_events = 1e4
n_background = RooRealVar('n_background', 'Number of QCD background events',
expected_background_events,
0.4 * expected_background_events,
1e3 * expected_background_events)
# fit model
fit_model = RooAddPdf('fit_model', 'signal + background model',
RooArgList(signal_model, background_model),
RooArgList(n_signal, n_background))
# RooMCStudy
RooRandom.randomGenerator().SetSeed(0)
n_toys = 1000
mc_study = RooMCStudy(fit_model, RooArgSet(x),
RooFit.FitOptions('shr'),
RooFit.Minos(ROOT.kTRUE),
RooFit.Binned(ROOT.kTRUE),
RooFit.Silence(ROOT.kTRUE))
ROOT.SetOwnership(mc_study, ROOT.kFALSE)
expected_total_events = int(
expected_signal_events + expected_background_events)
bool_keep_data = ROOT.kTRUE
mc_study.generateAndFit(n_toys, expected_total_events, bool_keep_data)
pull_signal_yield = mc_study.plotPull(n_signal,
RooFit.FitGauss(ROOT.kTRUE),
RooFit.FrameRange(-4., 4.))
quick_save(pull_signal_yield,
'pull_distribution',
path=output_dir,
legend_info=None)
# Why is this a null pointer?
fit_result = mc_study.fitResult(1)
print(fit_result)