diff --git a/bin/rabbit_limit.py b/bin/rabbit_limit.py index 5e77a52..5ae71e0 100755 --- a/bin/rabbit_limit.py +++ b/bin/rabbit_limit.py @@ -24,6 +24,14 @@ def make_parser(): parser = parsing.common_parser() parser.add_argument("--outname", default="limits.hdf5", help="output file name") + parser.add_argument( + "--modes", + nargs="+", + default=["gaussian", "likelihood"], + type=str, + choices=["gaussian", "likelihood"], + help="Different modes of approximations to be used", + ) parser.add_argument( "--asymptoticLimits", nargs="+", @@ -42,10 +50,8 @@ def make_parser(): return parser.parse_args() -def do_asymptotic_limits( - args, fitter, ws, data_values, mapping=None, bkg_only_fit=False -): - if bkg_only_fit: +def do_asymptotic_limits(args, fitter, ws, data_values, mapping=None, fit_data=False): + if fit_data: logger.info("Perform background only fit") # set process to zero and freeze fitter.freeze_params(args.asymptoticLimits) @@ -73,9 +79,10 @@ def do_asymptotic_limits( # Clone fitter to simultaneously minimize on asimov dataset fitter_asimov = copy.deepcopy(fitter) - # unconditional fit to real data - fitter.set_nobs(data_values) - fitter.minimize() + if fit_data: + logger.info("Unconditional fit to real data") + fitter.set_nobs(data_values) + fitter.minimize() # asymptotic limits (CLs) # see: @@ -130,13 +137,15 @@ def do_asymptotic_limits( idx = np.where(fitter.parms.astype(str) == key)[0][0] if key in fitter.poi_model.pois.astype(str): + is_poi = True xbest = fitter_asimov.get_blinded_poi()[idx] xobs = fitter.get_blinded_poi()[idx] elif key in fitter.parms.astype(str): + is_poi = False xbest = fitter_asimov.get_blinded_theta()[ - fitter_asimov.poi_model.npoi + idx + idx - fitter_asimov.poi_model.npoi ] - xobs = fitter.get_blinded_theta()[fitter_asimov.poi_model.npoi + idx] + xobs = fitter.get_blinded_theta()[idx - fitter.poi_model.npoi] xerr = fitter_asimov.cov[idx, idx] ** 0.5 @@ -149,7 +158,7 @@ def do_asymptotic_limits( fitter.cov.assign(cov) xobs_err = fitter.cov[idx, idx] ** 0.5 - if not args.allowNegativePOI: + if is_poi and not args.allowNegativePOI: xerr = 2 * xerr * xbest**0.5 xobs_err = 2 * xobs_err * xobs**0.5 @@ -167,53 +176,58 @@ def do_asymptotic_limits( qmuA = (norm.ppf(cl_b) - norm.ppf(cl_sb)) ** 2 logger.debug(f"Find r with q_(r,A)=-2ln(L)/ln(L0) = {qmuA}") + if "gaussian" in args.modes: + # Gaussian approximation + r = xbest + xerr * qmuA**0.5 + logger.info( + f"Expected (Gaussian) {round((cl_b)*100,1):4.1f}%: {key} < {r}" + ) + limits[i, j, k] = r + + if "likelihood" in args.modes: + # Likelihood approximation + r, _ = fitter_asimov.contour_scan( + key, nllvalreduced_asimov, qmuA, signs=[1], fun=fun + ) + logger.info( + f"Expected (Likelihood) {round((cl_b)*100,1):4.1f}%: {key} < {r}" + ) + limits_nll[i, j, k] = r + + if "gaussian" in args.modes: # Gaussian approximation - r = xbest + xerr * qmuA**0.5 - logger.info( - f"Expected (Gaussian) {round((cl_b)*100,1):4.1f}%: {key} < {r}" + limits_obs[i, j] = asymptotic_limits.compute_gaussian_limit( + key, xobs, xobs_err, xerr, cl_s ) - limits[i, j, k] = r + if "likelihood" in args.modes: # Likelihood approximation - r, _ = fitter_asimov.contour_scan( - key, nllvalreduced_asimov, qmuA, signs=[1], fun=fun - ) - logger.info( - f"Expected (Likelihood) {round((cl_b)*100,1):4.1f}%: {key} < {r}" - ) - limits_nll[i, j, k] = r - - # Gaussian approximation - limits_obs[i, j] = asymptotic_limits.compute_gaussian_limit( - key, xobs, xobs_err, xerr, cl_s - ) + if mapping is not None: + # TODO: implement for channels + pass + else: + # TODO: make it work + # nllvalreduced = fitter.reduced_nll().numpy() + # limits_nll_obs[i, j] = asymptotic_limits.compute_likelihood_limit(fitter, fitter_asimov, nllvalreduced, nllvalreduced_asimov, key, cl_s) + pass - # Likelihood approximation - if mapping is not None: - # TODO: implement for channels - pass - else: - # TODO: make it work - # nllvalreduced = fitter.reduced_nll().numpy() - # limits_nll_obs[i, j] = asymptotic_limits.compute_likelihood_limit(fitter, fitter_asimov, nllvalreduced, nllvalreduced_asimov, key, cl_s) - pass - - ws.add_limits_hist( - limits, - args.asymptoticLimits, - args.levels, - clb_list, - base_name="gaussian_asymptoticLimits_expected", - ) + if "gaussian" in args.modes: + ws.add_limits_hist( + limits, + args.asymptoticLimits, + args.levels, + clb_list, + base_name="gaussian_asymptoticLimits_expected", + ) - ws.add_limits_hist( - limits_obs, - args.asymptoticLimits, - args.levels, - base_name="gaussian_asymptoticLimits_observed", - ) + ws.add_limits_hist( + limits_obs, + args.asymptoticLimits, + args.levels, + base_name="gaussian_asymptoticLimits_observed", + ) - if mapping is not None: + if "likelihood" in args.modes: # TODO: implement for channels ws.add_limits_hist( limits_nll, @@ -348,7 +362,7 @@ def main(): ws, data_values=observations, mapping=mapping, - bkg_only_fit=ifit >= 0, + fit_data=ifit >= 0, ) fit_time.append(time.time()) diff --git a/bin/rabbit_plot_pulls_and_impacts.py b/bin/rabbit_plot_pulls_and_impacts.py index b847c30..c67660d 100755 --- a/bin/rabbit_plot_pulls_and_impacts.py +++ b/bin/rabbit_plot_pulls_and_impacts.py @@ -14,12 +14,14 @@ from rabbit import io_tools -from wums import output_tools, plot_tools # isort: skip +from wums import logging, output_tools, plot_tools # isort: skip # prevent MathJax from bein loaded pio.kaleido.scope.mathjax = None +logger = None + def writeOutput(fig, outfile, extensions=[], postfix=None, args=None, meta_info=None): name, _ = os.path.splitext(outfile) @@ -700,6 +702,17 @@ def parseArgs(): type=str, help="fitresults output hdf5 file from fit", ) + parser.add_argument( + "-v", + "--verbose", + type=int, + default=3, + choices=[0, 1, 2, 3, 4], + help="Set verbosity level with logging, the larger the more verbose", + ) + parser.add_argument( + "--noColorLogger", action="store_true", help="Do not use logging with colors" + ) parser.add_argument( "--result", default=None, @@ -712,14 +725,26 @@ def parseArgs(): default=None, type=str, nargs="+", - help="Print impacts on observables use '-m channel axes' for mapping results.", + help="Plot impacts on observables use '-m channel axes' for mapping results.", + ) + parser.add_argument( + "--channel", + default=None, + type=str, + help="Plot impacts for given channel", ) parser.add_argument( - "--mappingRef", + "--refMapping", default=None, type=str, nargs="+", - help="Print impacts on observables use '-m channel axes' for mapping results for reference.", + help="Plot impacts on observables from mapping results for reference.", + ) + parser.add_argument( + "--refChannel", + default=None, + type=str, + help="Plot impacts for given channel", ) parser.add_argument( "-r", @@ -1228,6 +1253,8 @@ def produce_plots_hist( def main(): args = parseArgs() + global logger + logger = logging.setup_logger(__file__, args.verbose, args.noColorLogger) config = plot_tools.load_config(args.config) @@ -1235,12 +1262,12 @@ def main(): fitresult, meta = io_tools.get_fitresult(args.inputFile, args.result, meta=True) if any( - x is not None for x in [args.referenceFile, args.refResult, args.mappingRef] + x is not None for x in [args.referenceFile, args.refResult, args.refMapping] ): - referenceFile = ( - args.referenceFile if args.referenceFile is not None else args.inputFile - ) - fitresult_ref = io_tools.get_fitresult(referenceFile, args.refResult) + if args.referenceFile is not None: + fitresult_ref = io_tools.get_fitresult(args.referenceFile, args.refResult) + else: + fitresult_ref = fitresult else: fitresult_ref = None @@ -1292,13 +1319,17 @@ def get_mapping_key(result, key): channels = res[key]["channels"] return channels, key else: - keys = [key for key in res.keys() if key.startswith(key)] + keys = [k for k in res.keys() if k.startswith(key)] if len(keys) == 0: raise ValueError( f"Mapping {key} not found, available mappings are: {res.keys()}" ) channels = res[keys[0]]["channels"] + logger.info( + f"Found mapping {keys[0]} with channels {[k for k in channels.keys()]}" + ) + return channels, keys[0] mapping_key = " ".join(args.mapping) @@ -1307,8 +1338,8 @@ def get_mapping_key(result, key): if fitresult_ref: mapping_key_ref = ( - " ".join(args.mappingRef) - if args.mappingRef is not None + " ".join(args.refMapping) + if args.refMapping is not None else mapping_key ) @@ -1317,6 +1348,8 @@ def get_mapping_key(result, key): ) for channel, hists in channels.items(): + if args.channel and channel not in args.channel: + continue modes = ["ungrouped", "group"] if args.mode == "both" else [args.mode] for mode in modes: @@ -1331,7 +1364,9 @@ def get_mapping_key(result, key): hist = hists[key].get() if fitresult_ref: - if channel in channels_ref.keys(): + if args.refChannel: + channel_ref = args.refChannel + elif channel in channels_ref.keys(): channel_ref = channel elif len(channels_ref.keys()) == 1: channel_ref = [v for v in channels_ref.keys()][0] diff --git a/rabbit/poi_models/poi_model.py b/rabbit/poi_models/poi_model.py index f45d997..758e8b8 100644 --- a/rabbit/poi_models/poi_model.py +++ b/rabbit/poi_models/poi_model.py @@ -64,6 +64,7 @@ def __init__(self, indata, **kwargs): self.poidefault = tf.zeros([], dtype=self.indata.dtype) self.allowNegativePOI = False + self.is_linear = True def compute(self, poi): rnorm = tf.ones(self.indata.nproc, dtype=self.indata.dtype) diff --git a/tests/plot_limits.py b/tests/plot_limits.py index cffff1d..2ea9358 100644 --- a/tests/plot_limits.py +++ b/tests/plot_limits.py @@ -5,6 +5,7 @@ import mplhep as hep import numpy as np +from wums import boostHistHelpers as hh from wums import logging from rabbit import io_tools @@ -59,6 +60,13 @@ def parseArgs(): parser.add_argument( "--scale", type=float, default=1.0, help="Scale limits by this value" ) + parser.add_argument( + "--mode", + type=str, + default="gaussian", + choices=["gaussian", "likelihood"], + help="Different modes of approximations to be used", + ) parser.add_argument( "--config", type=str, @@ -127,30 +135,49 @@ def main(): logger = logging.setup_logger(__file__, args.verbose, args.noColorLogger) + key_expected = f"{args.mode}_asymptoticLimits_expected" + key_observed = f"{args.mode}_asymptoticLimits_observed" + cls_list = set() clb_list = set() limits = [] limits_asimov = [] lumi = None + xvals = [] for infile in args.infiles: fitresult_asimov, meta = io_tools.get_fitresult( infile, result="asimov", meta=True ) - h_limits_asimov = fitresult_asimov["asymptoticLimits"].get() + bsm_sample, mixing = meta["meta_info_input"]["meta_info"]["args"][ + "addBSMMixing" + ] + bsm_name, mass = bsm_sample.split("_") + xvals.append(int(mass)) + + scale = float(mixing) * args.scale + + if key_expected in fitresult_asimov.keys(): + h_limits_asimov = fitresult_asimov[key_expected].get() + h_limits_asimov = hh.scaleHist(h_limits_asimov, scale) + else: + raise RuntimeError( + f"Key {key_expected} not found, available keyes are {[key for key in fitresult_asimov.keys()]}" + ) cls_list = cls_list.union(set([cls for cls in h_limits_asimov.axes["cls"]])) clb_list = clb_list.union(set([clb for clb in h_limits_asimov.axes["clb"]])) limits_asimov.append(h_limits_asimov) - fitresult, meta = io_tools.get_fitresult(infile, meta=True) + fitresult = io_tools.get_fitresult(infile) if fitresult != fitresult_asimov and not args.noObs: if ( - "asymptoticLimits" in fitresult.keys() - and "clb" not in fitresult["asymptoticLimits"].get().axes.name + key_observed in fitresult.keys() + and "clb" not in fitresult[key_observed].get().axes.name ): - h_limits = fitresult["asymptoticLimits"].get() + h_limits = fitresult[key_observed].get() + h_limits = hh.scaleHist(h_limits, scale) cls_list = cls_list.union(set([cls for cls in h_limits.axes["cls"]])) limits.append(h_limits) else: @@ -165,33 +192,33 @@ def main(): if args.xvals is not None: x = np.array(args.xvals) else: - x = np.arange(len(limits_asimov)) + x = np.array(xvals) clb_list = list(clb_list) for cls in cls_list: yexp = np.array( - [l[{"cls": cls, "clb": "0.5"}] * args.scale for l in limits_asimov] + [l[{"cls": cls, "clb": "0.5"}] for l in limits_asimov] ).flatten() yexp_m2 = np.array( - [l[{"cls": cls, "clb": "0.025"}] * args.scale for l in limits_asimov] + [l[{"cls": cls, "clb": "0.025"}] for l in limits_asimov] ).flatten() yexp_m1 = np.array( - [l[{"cls": cls, "clb": "0.16"}] * args.scale for l in limits_asimov] + [l[{"cls": cls, "clb": "0.16"}] for l in limits_asimov] ).flatten() yexp_p1 = np.array( - [l[{"cls": cls, "clb": "0.84"}] * args.scale for l in limits_asimov] + [l[{"cls": cls, "clb": "0.84"}] for l in limits_asimov] ).flatten() yexp_p2 = np.array( - [l[{"cls": cls, "clb": "0.975"}] * args.scale for l in limits_asimov] + [l[{"cls": cls, "clb": "0.975"}] for l in limits_asimov] ).flatten() ylist = [yexp, yexp_m1, yexp_m2, yexp_p1, yexp_p2] if len(limits) > 0: - yobs = np.array([l[{"cls": cls}] * args.scale for l in limits]).flatten() + yobs = np.array([l[{"cls": cls}] for l in limits]).flatten() ylist.append(yobs) ylist = np.array(ylist)