diff --git a/src/smefit/chi2.py b/src/smefit/chi2.py index 35979898..ee4bbfbe 100644 --- a/src/smefit/chi2.py +++ b/src/smefit/chi2.py @@ -83,24 +83,40 @@ def __init__(self, run_card, n_replica, scan_points=100): self.use_multiplicative_prescription = run_card.get( "use_multiplicative_prescription", False ) + self.mass_scan = False + self.run_card = run_card - rge_dict = run_card.get("rge", None) - operators_to_keep = run_card["coefficients"] - cutoff_scale = run_card.get("cutoff_scale", None) + # set all the coefficients to 0 + self.coefficients = CoefficientManager.from_dict(run_card["coefficients"]) + self.coefficients.set_free_parameters( + np.full(self.coefficients.free_parameters.shape[0], 0) + ) - if rge_dict is not None and run_card.get("datasets") is not None: - _logger.info("Loading RGE matrix") - rgemat, operators_to_keep = load_rge_matrix( - rge_dict, - list(operators_to_keep.keys()), - run_card["datasets"], - run_card.get("theory_path", None), - cutoff_scale, - ) - _logger.info("The operators generated by the RGE are: ") - _logger.info(list(operators_to_keep.keys())) - else: + self.rge_dict = run_card.get("rge", None) + operators_to_keep = run_card["coefficients"] + self.cutoff_scale = run_card.get("cutoff_scale", None) + + # Here we assume that if a mass is present in the coefficients + # this is a 1 parameter mass scan only + if np.any(self.coefficients.is_mass) and self.rge_dict is not None: + _logger.info("Mass scan detected") + _logger.info("RGE evolution will be computed for each mass in the scan") + self.mass_scan = True rgemat = None + else: + if self.rge_dict is not None and run_card.get("datasets") is not None: + _logger.info("Loading RGE matrix") + rgemat, operators_to_keep = load_rge_matrix( + self.rge_dict, + list(operators_to_keep.keys()), + run_card["datasets"], + run_card.get("theory_path", None), + self.cutoff_scale, + ) + _logger.info("The operators generated by the RGE are: ") + _logger.info(list(operators_to_keep.keys())) + else: + rgemat = None self.datasets = load_datasets( run_card["data_path"], @@ -116,13 +132,7 @@ def __init__(self, run_card, n_replica, scan_points=100): run_card.get("uv_couplings", False), run_card.get("external_chi2", False), rgemat=rgemat, - cutoff_scale=cutoff_scale, - ) - - # set all the coefficients to 0 - self.coefficients = CoefficientManager.from_dict(run_card["coefficients"]) - self.coefficients.set_free_parameters( - np.full(self.coefficients.free_parameters.shape[0], 0) + cutoff_scale=self.cutoff_scale, ) # build empty dict to store results @@ -166,6 +176,54 @@ def regularized_chi2_func(self, coeff, xs, use_replica): ) return np.array(chi2_list) + def chi2_mass_scan(self, coeff, xs): + + chi2_list = [] + for x in xs: + coeff.value = x + self.coefficients.set_constraints() + # compute rge matrix + _logger.info("Loading RGE matrix") + # assuming that the mass is in TeV + self.rge_dict.update({"init_scale": x * 1e3}) + rgemat, operators_to_keep = load_rge_matrix( + self.rge_dict, + list(self.run_card["coefficients"].keys()), + self.run_card["datasets"], + self.run_card.get("theory_path", None), + self.cutoff_scale, + ) + _logger.info("The operators generated by the RGE are: ") + _logger.info(list(operators_to_keep.keys())) + + self.datasets = load_datasets( + self.run_card["data_path"], + self.run_card["datasets"], + operators_to_keep, + self.run_card["use_quad"], + self.run_card["use_theory_covmat"], + False, + self.use_multiplicative_prescription, + self.run_card.get("default_order", "LO"), + self.run_card.get("theory_path", None), + self.run_card.get("rot_to_fit_basis", None), + self.run_card.get("uv_couplings", False), + self.run_card.get("external_chi2", False), + rgemat=rgemat, + cutoff_scale=self.cutoff_scale, + ) + + chi2_list.append( + compute_chi2( + self.datasets, + self.coefficients.value, + self.use_quad, + self.use_multiplicative_prescription, + use_replica=False, + ) + ) + return np.array(chi2_list) + def compute_bounds(self): r"""Compute individual bounds solving. @@ -243,9 +301,14 @@ def compute_scan(self): for coeff in self.coefficients: if coeff.name not in self.chi2_dict: continue - self.chi2_dict[coeff.name][rep] = self.regularized_chi2_func( - coeff, self.chi2_dict[coeff.name]["x"], use_replica - ) + if self.mass_scan: + self.chi2_dict[coeff.name][rep] = self.chi2_mass_scan( + coeff, self.chi2_dict[coeff.name]["x"] + ) + else: + self.chi2_dict[coeff.name][rep] = self.regularized_chi2_func( + coeff, self.chi2_dict[coeff.name]["x"], use_replica + ) # Parameter is scanned, # reset to 0 to not affect the next scan coeff.value = 0.0 diff --git a/src/smefit/coefficients.py b/src/smefit/coefficients.py index 69bcc358..a730ac6c 100644 --- a/src/smefit/coefficients.py +++ b/src/smefit/coefficients.py @@ -24,10 +24,13 @@ class Coefficient: """ - def __init__(self, name, minimum, maximum, value=None, constrain=False): + def __init__( + self, name, minimum, maximum, value=None, constrain=False, is_mass=False + ): self.name = name self.minimum = minimum self.maximum = maximum + self.is_mass = is_mass # determine if the parameter is free self.is_free = False @@ -142,6 +145,7 @@ def __init__(self, input_array): ) self._table.index = np.array([o.name for o in input_array], dtype=str) self.is_free = np.array([o.is_free for o in input_array], dtype=bool) + self.is_mass = np.array([o.is_mass for o in input_array], dtype=bool) # NOTE: this will not be updated. self._objlist = input_array @@ -193,6 +197,7 @@ def from_dict(cls, coefficient_config): property_dict["max"], constrain=constrain, value=property_dict.get("value", None), + is_mass=property_dict.get("is_mass", False), ) ) # make sure elements are sorted by names