Dear experts,
I am attempting to use RooStats to set the expected limit on a cross section, based on Monte Carlo studies.
I encounter problems when I try to include the systematic uncertainties on the luminosity.
The model is constituted by a 2nd order polynomial for the background, and a Crystal Ball with fixed parameters for the signal.
I use the factory syntax, in Python environment.
If I run the code without including the systematics, I do get a result that makes sense (this, to say that I am not 100% sure it’s the correct result, but at least is totally reasonable).
When I incorporate the machinery for the systematics, I obtain the same result, and I don’t know if I’m doing something wrong or if that’s the expected behaviour.
The essential parts of the codes are the following.
At the end of the post, I’ll paste the final part of the script, entitled to the actual limit extraction, if that was useful for debugging.
ALP_InvM_sq = w.factory('ALP_InvM_sq['+str(myrange[0])+','+str(myrange[1])+']')
w.factory('lumi_nom[511.0]')
w.factory('alpha_lumi[1., 0.0, 2]')
w.factory('prod::lumi(lumi_nom,alpha_lumi)')
w.factory('Gaussian::constr_lumi(alpha_lumi, glob_lumi[1., 0.0, 2], 0.1)')
w.factory('xsec[0, 0, 10]') #parameter of interest
w.factory('sigeff[0.258]')
w.factory('prod:sig_yield(lumi, xsec, sigeff)')
w.factory('CBShape:sig(ALP_InvM_sq, m0['+str(d["mean"])+'], sigma['+str(d["sigma"])+'], alpha['+str(d["alpha"])+'], n['+str(d["n"])+'])') # d is a dictionary
w.factory('Chebychev:bkg(ALP_InvM_sq, {a0[0, -5, 5], a1[0, -5, 5]})')
# w.factory('bkg_yield_nom[6000, 0, 50000]') # does it make a difference?
w.factory('bkg_yield_nom[6000]')
w.factory('prod::bkg_yield(bkg_yield_nom,alpha_lumi)')
w.factory('SUM::factorymodel_core(sig_yield*sig, bkg_yield*bkg)')
w.factory('PROD::factorymodel(factorymodel_core,constr_lumi)')
pdf = w.pdf('factorymodel')
# ----------------------------------------------------------------------------------------------------------------
# create set of observables (will need it for datasets and ModelConfig later)
pObs = ROOT.RooRealVar(w.var("ALP_InvM_sq"))
obs = ROOT.RooArgSet("observables")
obs.add(pObs)
# create set of global observables (need to be defined as constants?)
w.var("alpha_lumi").setConstant(True)
# w.var("bkg_yield_nom").setConstant(True)
w.var("lumi_nom").setConstant(True)
w.var("sigeff").setConstant(True)
globalObs = ROOT.RooArgSet("global_obs")
globalObs.add( w.var("alpha_lumi") )
globalObs.add( w.var("lumi_nom") )
globalObs.add( w.var("sigeff") )
# globalObs.add( w.var("bkg_yield_nom") )
# create set of nuisance parameters
nuis = ROOT.RooArgSet("nuis")
nuis.add( w.var("glob_lumi") )
nuis.add( w.var("a0") )
nuis.add( w.var("a1") )
nuis.add( w.var("bkg_yield_nom") )
# create set of parameters of interest (POI)
poi = ROOT.RooArgSet("poi")
poi.add( w.var("xsec") )
# ----------------------------------------------------------------------------------------------------------------
# bkg-only model
b_modelNM = ROOT.RooStats.ModelConfig("B_modelNM", w)
b_modelNM.SetPdf(w.pdf("factorymodel"))
b_modelNM.SetObservables( obs )
b_modelNM.SetParametersOfInterest( poi )
w.var("xsec").setVal(0.0)
b_modelNM.SetGlobalObservables( globalObs )
b_modelNM.SetNuisanceParameters( nuis )
b_modelNM.SetSnapshot(ROOT.RooArgSet(w.var("xsec")))
# bkg+signal model
sb_modelNM = ROOT.RooStats.ModelConfig("S+B_modelNM", w)
sb_modelNM.SetPdf(w.pdf("factorymodel"))
sb_modelNM.SetObservables( obs )
sb_modelNM.SetParametersOfInterest( poi )
w.var("xsec").setVal(1) # is any number ok?
sb_modelNM.SetGlobalObservables( globalObs )
sb_modelNM.SetNuisanceParameters( nuis )
sb_modelNM.SetSnapshot(ROOT.RooArgSet(w.var("xsec")))
- Is there any obvious error?
- Are the nuisance parameters and the global observables properly set? In particular alpha_lumi and glob_lumi.
- Is it normal that I extract the same limit (in particular the median limit, see the very last block of code at the end of the post) as without the nuisance parameters?
- I performed some tests, and I saw that if I set the starting value of alpha_lumi to any other value (0.5, 1.5…), then it and glob_lumi get have that same value at the end of fit. Is this correct? Is glob_lumi supposed to converge to the same starting value of alpha_lumi and not being modified?
Thanks for any help!
The rest of the code is here:
# ----------------------------------------------------------------------------------------------------------------
profll = ROOT.RooStats.ProfileLikelihoodTestStat(sb_modelNM.GetPdf())
ac = ROOT.RooStats.AsymptoticCalculator(dh_bg, b_modelNM, sb_modelNM)
ac.SetOneSided(True) # KLUDGE -- should want one sided (true) for limit
ROOT.RooStats.AsymptoticCalculator.SetPrintLevel(-1)
asympResult = ac.GetHypoTest()
asympResult.SetPValueIsRightTail(False) # appears to do nothing?!
asymp_pb = 1. - asympResult.AlternatePValue() # KLUDGE!!! Needs 1 -
asympResult.SetPValueIsRightTail(True)
asymp_psb = asympResult.NullPValue()
print( "Results based on asymptotic formulae:" )
print( "psb = ", asymp_psb )
print( "pb = ", asymp_pb )
asymp_clb = 1. - asymp_pb
asymp_clsb = asymp_psb
asymp_cls = 9999.
if asymp_clb > 0:
asymp_cls = asymp_clsb/asymp_clb
print( "cls = ", asymp_cls )
# create hypotest inverter passing the desired calculator (hc or ac)
calc = ROOT.RooStats.HypoTestInverter(ac)
calc.SetVerbose(False)
calc.SetConfidenceLevel(0.95)
useCLs = True
calc.UseCLs(useCLs)
if useCLs:
profll.SetOneSided(True)
npoints = 200 # number of points to scan
# min and max for scan (better to choose smaller intervals)
poimin = w.var("xsec").getMin()
poimax = w.var("xsec").getMax()
poimin = -10
poimax = 10
print( "Doing a fixed scan in interval : ",poimin," , ",poimax)
calc.SetFixedScan(npoints, poimin, poimax)
r = calc.GetInterval()
upperLimit = r.UpperLimit()
ulError = r.UpperLimitEstimatedError()
print( "The computed upper limit is: ", upperLimit," +/- ",ulError)
# compute expected limit
print("Expected upper limits using b-only model : ")
print("median limit = " , r.GetExpectedUpperLimit(0) )
print("med. limit (-1 sig) " , r.GetExpectedUpperLimit(-1))
print("med. limit (+1 sig) " , r.GetExpectedUpperLimit(1))
print("med. limit (-2 sig) " , r.GetExpectedUpperLimit(-2))
print("med. limit (+2 sig) " , r.GetExpectedUpperLimit(2))