diff --git a/alea/blueice_extended_model.py b/alea/blueice_extended_model.py index d1af0cdd..03bf7413 100644 --- a/alea/blueice_extended_model.py +++ b/alea/blueice_extended_model.py @@ -1,4 +1,6 @@ from pydoc import locate # to lookup likelihood class +from typing import List + from alea.statistical_model import StatisticalModel from alea.simulators import BlueiceDataGenerator from alea.utils import adapt_likelihood_config_for_blueice @@ -76,7 +78,12 @@ def get_expectation_values(self, **kwargs) -> dict: """ ret = dict() for ll in self._likelihood.likelihood_list[:-1]: # ancillary likelihood does not contribute - mus = ll(full_output=True, **kwargs)[1] + + ll_pars = list(ll.rate_parameters.keys()) + list(ll.shape_parameters.keys()) + ll_pars += ["livetime_days"] + call_args = {k: i for k, i in kwargs.items() if k in ll_pars} + + mus = ll(full_output=True, **call_args)[1] for n, mu in zip(ll.source_name_list, mus): ret[n] = ret.get(n, 0) + mu return ret @@ -105,13 +112,27 @@ def _build_ll_from_config(self, likelihood_config: dict) -> "LogLikelihoodSum": config, template_folder_list) blueice_config["livetime_days"] = self.parameters[ blueice_config["livetime_parameter"]].nominal_value + for p in self.parameters: + blueice_config[p.name] = blueice_config.get(p.name, p.nominal_value) + + # add all parameters to extra_dont_hash for each source unless it is used: + for i, source in enumerate(config["sources"]): + parameters_to_ignore: List[str] = [p.name for p in self.parameters if (p.ptype == "shape") + & (p.name not in source["parameters"])] + # no efficiency affects PDF: + parameters_to_ignore += [p.name for p in self.parameters if (p.ptype == "efficiency")] + parameters_to_ignore += source.get("extra_dont_hash_settings", []) + + # ignore all shape parameters known to this model not named specifically in the source: + blueice_config["sources"][i]["extra_dont_hash_settings"] = parameters_to_ignore + ll = likelihood_object(blueice_config) for source in config["sources"]: # Set rate parameters rate_parameters = [ - p for p in source["parameters"] if self.parameters[p].type == "rate"] + p for p in source["parameters"] if self.parameters[p].ptype == "rate"] if len(rate_parameters) != 1: raise ValueError( f"Source {source['name']} must have exactly one rate parameter.") @@ -125,11 +146,15 @@ def _build_ll_from_config(self, likelihood_config: dict) -> "LogLikelihoodSum": # Set shape parameters shape_parameters = [ - p for p in source["parameters"] if self.parameters[p].type == "shape"] + p for p in source["parameters"] if self.parameters[p].ptype == "shape"] if shape_parameters: # TODO: Implement setting shape parameters raise NotImplementedError("Shape parameters are not yet supported.") + # Set efficiency parameters + if source.get("apply_efficiency", False): + self._set_efficiency(source, ll) + ll.prepare() lls.append(ll) @@ -179,6 +204,21 @@ def _generate_ancillary_measurements(self, **generate_values) -> dict: return ancillary_measurements + def _set_efficiency(self, source, ll): + assert "efficiency_name" in source, "Unspecified efficiency_name for source {:s}".format(source["name"]) + efficiency_name = source["efficiency_name"] + assert efficiency_name in source[ + "parameters"], "The efficiency_name for source {:s} is not in its parameter list".format(source["name"]) + efficiency_parameter = self.parameters[efficiency_name] + assert efficiency_parameter.ptype == "efficiency", "The parameter {:s} must" \ + " be an efficiency".format(efficiency_name) + limits = efficiency_parameter.fit_limits + assert 0 <= limits[0], 'Efficiency parameters including {:s} must be' \ + ' constrained to be nonnegative'.format(efficiency_name) + assert np.isfinite(limits[1]), 'Efficiency parameters including {:s} must be' \ + ' constrained to be finite'.format(efficiency_name) + ll.add_shape_parameter(efficiency_name, anchors=(limits[0], limits[1])) + class CustomAncillaryLikelihood(LogAncillaryLikelihood): """ diff --git a/alea/examples/unbinned_wimp_statistical_model.yaml b/alea/examples/unbinned_wimp_statistical_model.yaml index f482b7d1..58b38c2b 100644 --- a/alea/examples/unbinned_wimp_statistical_model.yaml +++ b/alea/examples/unbinned_wimp_statistical_model.yaml @@ -33,6 +33,18 @@ parameter_definition: - null fit_guess: 1.0 + signal_efficiency: + nominal_value: 1.0 + ptype: efficiency + uncertainty: 0.1 + relative_uncertainty: true + fittable: true + fit_limits: + - 0 + - 10. + fit_guess: 1.0 + description: Parameter to account for the uncertain signal expectation given a certain cross-section + # er_band_shift: # nominal_value: 0 # ptype: shape @@ -76,9 +88,10 @@ likelihood_config: parameters: - wimp_rate_multiplier - wimp_mass + - signal_efficiency template_filename: wimp50gev_template.h5 - apply_efficiency: False - efficiency_name: 'signal_eff' # TODO: Check + apply_efficiency: True + efficiency_name: signal_efficiency # SR1 - name: sr1 @@ -104,6 +117,7 @@ likelihood_config: parameters: - wimp_rate_multiplier - wimp_mass + - signal_efficiency template_filename: wimp50gev_template.h5 - apply_efficiency: False - efficiency_name: 'wimp_eff' # TODO: Check \ No newline at end of file + apply_efficiency: True + efficiency_name: signal_efficiency \ No newline at end of file diff --git a/alea/parameters.py b/alea/parameters.py index 2367c889..17448ffa 100644 --- a/alea/parameters.py +++ b/alea/parameters.py @@ -9,7 +9,7 @@ class Parameter: name (str): The name of the parameter. nominal_value (float, optional): The nominal value of the parameter. fittable (bool, optional): Indicates if the parameter is fittable or always fixed. - ptype (str, optional): The type of the parameter. + ptype (str, optional): The ptype of the parameter. uncertainty (float or str, optional): The uncertainty of the parameter. If a string, it can be evaluated as a numpy or scipy function to define non-gaussian constraints. @@ -39,7 +39,7 @@ def __init__( self.name = name self.nominal_value = nominal_value self.fittable = fittable - self.type = ptype + self.ptype = ptype self._uncertainty = uncertainty self.relative_uncertainty = relative_uncertainty self.blueice_anchors = blueice_anchors diff --git a/alea/statistical_model.py b/alea/statistical_model.py index ce0542a4..9e201b63 100644 --- a/alea/statistical_model.py +++ b/alea/statistical_model.py @@ -163,7 +163,7 @@ def store_data( kw = {'metadata': metadata} if metadata is not None else dict() toydata_to_file(file_name, data_list, data_name_list, **kw) - def get_expectation_values(self): + def get_expectation_values(self, **parameter_values): return NotImplementedError("get_expectation_values is optional to implement") @property @@ -284,9 +284,9 @@ def _confidence_interval_checks( "You must set parameter_interval_bounds in the parameter config" " or when calling confidence_interval") - if parameter_of_interest.type == "rate": + if parameter_of_interest.ptype == "rate": try: - if parameter_of_interest.type == "rate" and poi_name.endswith("_rate_multiplier"): + if parameter_of_interest.ptype == "rate" and poi_name.endswith("_rate_multiplier"): source_name = poi_name.replace("_rate_multiplier", "") else: source_name = poi_name