SUSY Limit Plotting

[ekourlit]

This issue describes the way of presenting SUSY limits as it has been agreed between ATLAS and CMS. An example from JHEP 12 (2019) 060 can be seen here:

Description of limit lines

1. Observed limit (thick solid dark-red line): all uncertainties are included in the fit as nuisance parameters, with the exception of the theoretical signal uncertainties (PDF, scales, etc.) on the inclusive cross section.
2. ±1σ lines around observed limit (1) with style "thin dark-red dotted": re-run limit calculation (1) while increasing or decreasing the signal cross section by the theoretical signal uncertainties (PDF, scales, etc.).
3. Expected limit (less thick long-dashed dark-blue line): all uncertainties are included in the fit as nuisance parameters, with the exception of the theoretical signal uncertainties (PDF, scales, etc.) on the inclusive cross section.
4. ±1σ band around expected limit (2) with style "yellow band": the band contours are the ±1σ results of the fit (2).

Implementation

Below an implementation of the calculation of the above CLs value using pyhf and cabinetry is presented.

Workspace

The workspace is created by a cabinetry configuration file. Is should be noted that this workspace includes multiple signal models, i.e. multiple signal Samples with different NormFactors each. An example of such configuration file can be seen here: config_excl_AnalysisChallenge_grid.yml.

Code snippet

# get model and observed data (pre-fit)
model, data = cabinetry.model_utils.model_and_data(ws)
channels = model.config.channels

# workspace and model without signal theory uncertainties. signal theory uncertainties have the string '_Signal_' in their name.
signal_theory_names = [systematic['Name'] for systematic in cabinetry_config['Systematics'] if '_Signal_' in systematic['Name']]
pruned_ws = pyhf.Workspace(ws)._prune_and_rename(prune_modifiers=signal_theory_names)
pruned_model, _ = cabinetry.model_utils.model_and_data(pruned_ws)

# signal normalisation parameter names (having the string 'Signal_norm' in their name).
signal_norms = [param for param in model.config.par_names if 'Signal_norm' in param]
# signal normalisation parameters indices for both full and pruned models
signal_norms_idx = [model.config.par_map[norm]['slice'].start for norm in signal_norms]
pruned_signal_norms_idx = [pruned_model.config.par_map[norm]['slice'].start for norm in signal_norms]
# signal theory modifier indices
signal_theory_idx = [model.config.par_map[modifier]['slice'].start for modifier in signal_theory_names]

def fix_and_calculate_CLs(param):
    '''
    Fix all signal normalisation parameters to zero except the one specified by param
    then calculate expected and observed CLs values
    '''

    # Expected CLs
    print(f"Calculating explected CLs for {param}")

    # get the index of the param, this will be the new POI
    pruned_poi_index = cabinetry.model_utils._poi_index(pruned_model, poi_name=param)
    # set POI to param
    pruned_model.config.set_poi(pruned_model.config.par_names[pruned_poi_index])

    # indices of all the signal normalisation parameters except the poi_idx
    pruned_fixed_signal_norms_idx = [idx for idx in pruned_signal_norms_idx if idx!=pruned_poi_index]

    # set the rest of signal normalisation parameters fixed at zero
    pruned_fix_pars = pruned_model.config.suggested_fixed().copy()
    pruned_init_pars = pruned_model.config.suggested_init().copy()
    for idx in pruned_fixed_signal_norms_idx:
        pruned_fix_pars[idx] = True
        pruned_init_pars[idx] = 0.0
    
    # calculate expected CLs with respect to POI
    try:
        _, CLs_exp_band = pyhf.infer.hypotest(1.0, asimov_dataset, pruned_model, init_pars=pruned_init_pars, fixed_params=pruned_fix_pars, return_expected_set=True)
        # keep only ±1 sigma and convert to simple floats
        CLs_exp_band = list(map(lambda x: float(x), CLs_exp_band[1:4]))
    except:
        CLs_exp_band = None
    
    # Observed CLs
    print(f"Calculating observed CLs for {param}")

    # get the index of the param, this will be the new POI
    poi_index = cabinetry.model_utils._poi_index(model, poi_name=param)
    # set POI to param
    model.config.set_poi(model.config.par_names[poi_index])
    
    # indices of all the signal normalisation parameters except the poi_idx
    fixed_signal_norms_idx = [idx for idx in signal_norms_idx if idx!=poi_index]

    # set the rest of signal normalisation parameters fixed at zero
    fix_pars = model.config.suggested_fixed().copy()
    init_pars = model.config.suggested_init().copy()
    for idx in fixed_signal_norms_idx:
        fix_pars[idx] = True
        init_pars[idx] = 0.0

    CLs_obs_band = []
    # run three fits with different fixed signal theory NPs
    for signal_theory_np in [-1.0, 0.0, 1.0]:
        # fix signal theory NPs
        for idx in signal_theory_idx:
            fix_pars[idx] = True
            init_pars[idx] = signal_theory_np

        # calculate observed CLs with respect to POI
        try:
            CLs_obs = pyhf.infer.hypotest(1.0, data, model, init_pars=init_pars, fixed_params=fix_pars)
            # convert to simple floats
            CLs_obs = float(CLs_obs)
        except:
            CLs_obs = None

        # save results
        CLs_obs_band.append(CLs_obs)

    return [CLs_obs_band, CLs_exp_band]

# results dictionary
post_fit_significance_results = {}

# loop over the signal parameters of the workspace and calculate
for param in signal_norms:
    result = fix_and_calculate_CLs[param]
    post_fit_significance_results[param] = result