Added mass scan rge

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

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