Repository: sibirrer/hierArc (language: python)

Issue:
Minor PR
Here are propositions of small adjustments that I find to be useful when including hierArc in our pipeline.


Discussion / hints:
Thanks!

Fix patch:
diff --git a/hierarc/LensPosterior/kin_constraints_composite.py b/hierarc/LensPosterior/kin_constraints_composite.py
index f7b2e77..ff88859 100644
--- a/hierarc/LensPosterior/kin_constraints_composite.py
+++ b/hierarc/LensPosterior/kin_constraints_composite.py
@@ -15,9 +15,9 @@ def __init__(
         z_lens,
         z_source,
         gamma_in_array,
-        m2l_array,
-        rho0_array,
-        r_s_array,
+        log_m2l_array,
+        kappa_s_array,
+        r_s_angle_array,
         theta_E,
         theta_E_error,
         gamma,
@@ -39,17 +39,20 @@ def __init__(
         num_psf_sampling=100,
         num_kin_sampling=1000,
         multi_observations=False,
+        rho0_array=None,
+        r_s_array=None,
+        is_m2l_population_level=True,
     ):
         """
 
         :param z_lens: lens redshift
         :param z_source: source redshift
         :param gamma_in_array: array of power-law slopes of the mass model
-        :param m2l_array: array of mass-to-light ratios of the stellar component,
-            needs to be in the unit/scaling such that m2l * amp in the
+        :param log_m2l_array: array of log10(mass-to-light ratios) of the stellar component,
+            needs to be in the unit/scaling such that m2l / sigma_crit * amp in the
             kwargs_lens_light provides the convergence amplitude of the stars
-        :param rho0_array: array of halo mass normalizations in M_sun / Mpc^3
-        :param r_s_array: array of halo scale radii in arc seconds
+        :param kappa_s_array: array of generalized NFW profile's convergence normalization at the scale radius
+        :param r_s_angle_array: array of halo scale radii in arcsecond
         :param theta_E: Einstein radius (in arc seconds)
         :param theta_E_error: 1-sigma error on Einstein radius
         :param gamma: power-law slope (2 = isothermal) estimated from imaging data
@@ -80,22 +83,31 @@ def __init__(
             routine
         :param multi_observations: bool, if True, interprets kwargs_aperture and
             kwargs_seeing as lists of multiple observations
+        :param rho0_array: array of halo mass normalizations in M_sun / Mpc^3
+        :param r_s_array: array of halo scale radii in Mpc
         """
-        self._light_profile_analysis = LightProfileAnalysis(
-            light_model=LightModel(light_model_list=lens_light_model_list)
-        )
 
-        (
-            amps,
-            sigmas,
-            center_x,
-            center_y,
-        ) = self._light_profile_analysis.multi_gaussian_decomposition(
-            kwargs_lens_light, r_h=r_eff, **kwargs_mge_light
-        )
+        if (
+            len(lens_light_model_list) == 1
+            and lens_light_model_list[0] == "MULTI_GAUSSIAN"
+        ):
+            kwargs_lens_light = kwargs_lens_light
+        else:
+            self._light_profile_analysis = LightProfileAnalysis(
+                light_model=LightModel(light_model_list=lens_light_model_list)
+            )
+
+            (
+                amps,
+                sigmas,
+                center_x,
+                center_y,
+            ) = self._light_profile_analysis.multi_gaussian_decomposition(
+                kwargs_lens_light, r_h=r_eff, **kwargs_mge_light
+            )
 
-        lens_light_model_list = ["MULTI_GAUSSIAN"]
-        kwargs_lens_light = [{"amp": amps, "sigma": sigmas}]
+            lens_light_model_list = ["MULTI_GAUSSIAN"]
+            kwargs_lens_light = [{"amp": amps, "sigma": sigmas}]
 
         lens_model_list = ["GNFW", "MULTI_GAUSSIAN_KAPPA"]
 
@@ -128,16 +140,39 @@ def __init__(
             multi_observations=multi_observations,
         )
 
-        if len(rho0_array) != len(r_s_array):
-            raise ValueError("rho0 and r_s must have the same length!")
+        if self._check_arrays(kappa_s_array, r_s_angle_array):
+            self._kappa_s_array = kappa_s_array
+            self._r_scale_angle_array = r_s_angle_array
+        elif self._check_arrays(rho0_array, r_s_array):
+            self._kappa_s_array, self._r_scale_angle_array = self.get_kappa_s_r_s_angle(
+                rho0_array, r_s_array
+            )
+        else:
+            raise ValueError(
+                "Both kappa_s_array and r_s_angle_array, or rho0_array and r_s_array must be arrays of the same length!"
+            )
 
-        self._rho_s_array = rho0_array
-        self._r_scale_array = r_s_array
-        self._kappa_s_array, self._r_scale_angle_array = self.get_kappa_s_r_s_angle(
-            rho0_array, r_s_array
-        )
         self.gamma_in_array = gamma_in_array
-        self.m2l_array = m2l_array
+        self.log_m2l_array = log_m2l_array
+        self._is_m2l_population_level = is_m2l_population_level
+
+        if not is_m2l_population_level and not self._check_arrays(
+            self._kappa_s_array, log_m2l_array
+        ):
+            raise ValueError(
+                "log_m2l_array must have the same length as rho0_array or kappa_s_array!"
+            )
+
+    @staticmethod
+    def _check_arrays(array0, array1):
+        """Checks if two arrays are valid sample distributions for pairs of parameters.
+
+        :param array0, array1: arrays representing samples from respective distributions
+        :return: bool, True if arrays are the same length and > 0; False otherwise
+        """
+        if array0 is None or array1 is None:
+            return False
+        return len(array0) == len(array1) and len(array0) > 0
 
     def get_kappa_s_r_s_angle(self, rho_s, r_scale):
         """Computes the surface mass density of the NFW halo at the scale radius.
@@ -158,14 +193,23 @@ def draw_lens(self, no_error=False):
             the mean values instead
         """
         if no_error is True:
-            return (
-                np.mean(self._rho_s_array),
-                np.mean(self._r_scale_array),
-                self._r_eff,
-                1,
-            )
+            if self._is_m2l_population_level:
+                return (
+                    np.mean(self._kappa_s_array),
+                    np.mean(self._r_scale_angle_array),
+                    self._r_eff,
+                    1,
+                )
+            else:
+                return (
+                    np.mean(self._kappa_s_array),
+                    np.mean(self._r_scale_angle_array),
+                    np.mean(self.log_m2l_array),
+                    self._r_eff,
+                    1,
+                )
 
-        random_index = np.random.randint(low=0, high=len(self._rho_s_array))
+        random_index = np.random.randint(low=0, high=len(self._kappa_s_array))
         kappa_s_draw = self._kappa_s_array[random_index]
         r_scale_angle_draw = self._r_scale_angle_array[random_index]
 
@@ -175,7 +219,17 @@ def draw_lens(self, no_error=False):
         )
         r_eff_draw = delta_r_eff * self._r_eff
 
-        return kappa_s_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff
+        if self._is_m2l_population_level:
+            return kappa_s_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff
+        else:
+            log_m2l_draw = self.log_m2l_array[random_index]
+            return (
+                kappa_s_draw,
+                r_scale_angle_draw,
+                log_m2l_draw,
+                r_eff_draw,
+                delta_r_eff,
+            )
 
     def model_marginalization(self, num_sample_model=20):
         """
@@ -190,24 +244,33 @@ def model_marginalization(self, num_sample_model=20):
             (num_sample_model, num_data)
         )  # matrix that contains the sampled J() distribution
         for i in range(num_sample_model):
-            j_kin = self.j_kin_draw_composite(
-                self.kwargs_anisotropy_base,
-                np.mean(self.gamma_in_array),
-                np.mean(self.m2l_array),
-                no_error=False,
-            )
+            if self._is_m2l_population_level:
+                j_kin = self.j_kin_draw_composite(
+                    self.kwargs_anisotropy_base,
+                    np.mean(self.gamma_in_array),
+                    np.mean(self.log_m2l_array),
+                    no_error=False,
+                )
+            else:
+                j_kin = self.j_kin_draw_composite_m2l(
+                    self.kwargs_anisotropy_base,
+                    np.mean(self.gamma_in_array),
+                    no_error=False,
+                )
             j_kin_matrix[i, :] = j_kin
 
         error_cov_j_sqrt = np.cov(np.sqrt(j_kin_matrix.T))
-        j_model_list = np.mean(j_kin_matrix, axis=0)
+        j_model_list = np.nanmean(j_kin_matrix, axis=0)
         return j_model_list, error_cov_j_sqrt
 
-    def j_kin_draw_composite(self, kwargs_anisotropy, gamma_in, m2l, no_error=False):
+    def j_kin_draw_composite(
+        self, kwargs_anisotropy, gamma_in, log_m2l, no_error=False
+    ):
         """One simple sampling realization of the dimensionless kinematics of the model.
 
         :param kwargs_anisotropy: keyword argument of anisotropy setting
         :param gamma_in: power-law slope of the mass model
-        :param m2l: mass-to-light ratio of the stellar component
+        :param log_m2l: log10(mass-to-light ratio) of the stellar component
         :param no_error: bool, if True, does not render from the uncertainty but uses
             the mean values instead
         :return: dimensionless kinematic component J() Birrer et al. 2016, 2019
@@ -218,7 +281,58 @@ def j_kin_draw_composite(self, kwargs_anisotropy, gamma_in, m2l, no_error=False)
 
         kwargs_lens_stars = copy.deepcopy(self._kwargs_lens_light[0])
 
-        kwargs_lens_stars["amp"] *= m2l
+        kwargs_lens_stars["amp"] *= 10**log_m2l / self.lensCosmo.sigma_crit_angle
+
+        kwargs_lens_stars["sigma"] *= delta_r_eff
+
+        kwargs_light = copy.deepcopy(self._kwargs_lens_light)
+        for kwargs in kwargs_light:
+            kwargs["sigma"] *= delta_r_eff
+
+        kwargs_lens = [
+            {
+                "Rs": r_scale_angle_draw,
+                "gamma_in": gamma_in,
+                "kappa_s": kappa_s_draw,
+                "center_x": 0,
+                "center_y": 0,
+            },
+            kwargs_lens_stars,
+        ]
+
+        print(kwargs_lens)
+
+        j_kin = self.velocity_dispersion_map_dimension_less(
+            kwargs_lens=kwargs_lens,
+            kwargs_lens_light=kwargs_light,
+            kwargs_anisotropy=kwargs_anisotropy,
+            r_eff=r_eff_draw,
+            theta_E=self._theta_E,  # send this to avoid unnecessary recomputation
+            gamma=self._gamma,  # send this to avoid unnecessary recomputation
+        )
+        return j_kin
+
+    def j_kin_draw_composite_m2l(self, kwargs_anisotropy, gamma_in, no_error=False):
+        """Similar function to j_kin_draw_composite, but now drawing from log_m2l
+        distribution.
+
+        :param kwargs_anisotropy: keyword argument of anisotropy setting
+        :param gamma_in: power-law slope of the mass model
+        :param no_error: bool, if True, does not render from the uncertainty but uses
+            the mean values instead
+        :return: dimensionless kinematic component J() Birrer et al. 2016, 2019
+        """
+        (
+            kappa_s_draw,
+            r_scale_angle_draw,
+            log_m2l_draw,
+            r_eff_draw,
+            delta_r_eff,
+        ) = self.draw_lens(no_error=no_error)
+
+        kwargs_lens_stars = copy.deepcopy(self._kwargs_lens_light[0])
+
+        kwargs_lens_stars["amp"] *= log_m2l_draw / self.lensCosmo.sigma_crit_angle
 
         kwargs_lens_stars["sigma"] *= delta_r_eff
 
@@ -275,9 +389,12 @@ def hierarchy_configuration(self, num_sample_model=20):
             "error_cov_j_sqrt": error_cov_j_sqrt,
             "ani_param_array": self.ani_param_array,
             "gamma_in_array": self.gamma_in_array,
-            "m2l_array": self.m2l_array,
+            "log_m2l_array": self.log_m2l_array,
             "param_scaling_grid_list": ani_scaling_grid_list,
         }
+
+        if not self._is_m2l_population_level:
+            kwargs_likelihood["log_m2l_array"] = None
         return kwargs_likelihood
 
     def anisotropy_scaling(self):
@@ -285,13 +402,21 @@ def anisotropy_scaling(self):
 
         :return: anisotropy scaling grid along the axes defined by ani_param_array
         """
-        j_ani_0 = self.j_kin_draw_composite(
-            self.kwargs_anisotropy_base,
-            np.mean(self.gamma_in_array),
-            np.mean(self.m2l_array),
-            no_error=True,
-        )
-        return self._anisotropy_scaling_relative(j_ani_0)
+        if self._is_m2l_population_level:
+            j_ani_0 = self.j_kin_draw_composite(
+                self.kwargs_anisotropy_base,
+                np.mean(self.gamma_in_array),
+                np.mean(self.log_m2l_array),
+                no_error=True,
+            )
+            return self._anisotropy_scaling_relative(j_ani_0)
+        else:
+            j_ani_0 = self.j_kin_draw_composite_m2l(
+                self.kwargs_anisotropy_base,
+                np.mean(self.gamma_in_array),
+                no_error=True,
+            )
+            return self._anisotropy_scaling_relative_m2l(j_ani_0)
 
     def _anisotropy_scaling_relative(self, j_ani_0):
         """Anisotropy scaling relative to a default J prediction.
@@ -309,7 +434,7 @@ def _anisotropy_scaling_relative(self, j_ani_0):
                         len(self.ani_param_array[0]),
                         len(self.ani_param_array[1]),
                         len(self.gamma_in_array),
-                        len(self.m2l_array),
+                        len(self.log_m2l_array),
                     )
                 )
                 for _ in range(num_data)
@@ -317,12 +442,12 @@ def _anisotropy_scaling_relative(self, j_ani_0):
             for i, a_ani in enumerate(self.ani_param_array[0]):
                 for j, beta_inf in enumerate(self.ani_param_array[1]):
                     for k, g_in in enumerate(self.gamma_in_array):
-                        for l, m2l in enumerate(self.m2l_array):
+                        for l, log_m2l in enumerate(self.log_m2l_array):
                             kwargs_anisotropy = self.anisotropy_kwargs(
                                 a_ani=a_ani, beta_inf=beta_inf
                             )
                             j_kin_ani = self.j_kin_draw_composite(
-                                kwargs_anisotropy, g_in, m2l, no_error=True
+                                kwargs_anisotropy, g_in, log_m2l, no_error=True
                             )
 
                             for m, j_kin in enumerate(j_kin_ani):
@@ -336,20 +461,76 @@ def _anisotropy_scaling_relative(self, j_ani_0):
                     (
                         len(self.ani_param_array),
                         len(self.gamma_in_array),
-                        len(self.m2l_array),
+                        len(self.log_m2l_array),
                     )
                 )
                 for _ in range(num_data)
             ]
             for i, a_ani in enumerate(self.ani_param_array):
                 for k, g_in in enumerate(self.gamma_in_array):
-                    for l, m2l in enumerate(self.m2l_array):
+                    for l, log_m2l in enumerate(self.log_m2l_array):
                         kwargs_anisotropy = self.anisotropy_kwargs(a_ani)
                         j_kin_ani = self.j_kin_draw_composite(
-                            kwargs_anisotropy, g_in, m2l, no_error=True
+                            kwargs_anisotropy, g_in, log_m2l, no_error=True
                         )
                         for m, j_kin in enumerate(j_kin_ani):
                             ani_scaling_grid_list[m][i, k, l] = j_kin / j_ani_0[m]
         else:
             raise ValueError("anisotropy model %s not valid." % self._anisotropy_model)
         return ani_scaling_grid_list
+
+    def _anisotropy_scaling_relative_m2l(self, j_ani_0):
+        """Similar function to _anisotropy_scaling_relative, but instead drawing log_m2l
+        from a distribution.
+
+        :param j_ani_0: default J() prediction for default anisotropy
+        :return: list of arrays (for the number of measurements) according to anisotropy
+            scaling
+        """
+        num_data = len(self._sigma_v_measured)
+
+        if self._anisotropy_model == "GOM":
+            ani_scaling_grid_list = [
+                np.zeros(
+                    (
+                        len(self.ani_param_array[0]),
+                        len(self.ani_param_array[1]),
+                        len(self.gamma_in_array),
+                    )
+                )
+                for _ in range(num_data)
+            ]
+            for i, a_ani in enumerate(self.ani_param_array[0]):
+                for j, beta_inf in enumerate(self.ani_param_array[1]):
+                    for k, g_in in enumerate(self.gamma_in_array):
+                        kwargs_anisotropy = self.anisotropy_kwargs(
+                            a_ani=a_ani, beta_inf=beta_inf
+                        )
+                        j_kin_ani = self.j_kin_draw_composite_m2l(
+                            kwargs_anisotropy, g_in, no_error=True
+                        )
+
+                        for m, j_kin in enumerate(j_kin_ani):
+                            ani_scaling_grid_list[m][i, j, k] = j_kin / j_ani_0[m]
+                            # perhaps change the order
+        elif self._anisotropy_model in ["OM", "const"]:
+            ani_scaling_grid_list = [
+                np.zeros(
+                    (
+                        len(self.ani_param_array),
+                        len(self.gamma_in_array),
+                    )
+                )
+                for _ in range(num_data)
+            ]
+            for i, a_ani in enumerate(self.ani_param_array):
+                for k, g_in in enumerate(self.gamma_in_array):
+                    kwargs_anisotropy = self.anisotropy_kwargs(a_ani)
+                    j_kin_ani = self.j_kin_draw_composite_m2l(
+                        kwargs_anisotropy, g_in, no_error=True
+                    )
+                    for m, j_kin in enumerate(j_kin_ani):
+                        ani_scaling_grid_list[m][i, k] = j_kin / j_ani_0[m]
+        else:
+            raise ValueError("anisotropy model %s not valid." % self._anisotropy_model)
+        return ani_scaling_grid_list
diff --git a/hierarc/Likelihood/hierarchy_likelihood.py b/hierarc/Likelihood/hierarchy_likelihood.py
index 533f2c3..2e87b26 100644
--- a/hierarc/Likelihood/hierarchy_likelihood.py
+++ b/hierarc/Likelihood/hierarchy_likelihood.py
@@ -21,7 +21,7 @@ def __init__(
         ani_scaling_array=None,
         ani_scaling_array_list=None,
         gamma_in_array=None,
-        m2l_array=None,
+        log_m2l_array=None,
         param_scaling_grid_list=None,
         num_distribution_draws=50,
         kappa_ext_bias=False,
@@ -47,7 +47,7 @@ def __init__(
          provided J() predictions.
         :param param_scaling_grid_list: list of N-dimensional arrays with the
             scalings of J() for each IFU. Needed when simultaneously scaling
-            anisotropy, gamma_in, and m2l. In that case, gamma_in_array and m2l_array need to be provided.
+            anisotropy, gamma_in, and log_m2l. In that case, gamma_in_array and log_m2l_array need to be provided.
         :param num_distribution_draws: int, number of distribution draws from the likelihood that are being averaged
          over
         :param kappa_ext_bias: bool, if True incorporates the global external selection function into the likelihood.
@@ -77,11 +77,22 @@ def __init__(
         #     ani_param_array=ani_param_array,
         #     ani_scaling_array_list=ani_scaling_array_list,
         # )
-        if gamma_in_array is not None and m2l_array is not None:
+        if gamma_in_array is not None and log_m2l_array is not None:
             if isinstance(ani_param_array, list):
-                param_arrays = ani_param_array + [gamma_in_array, m2l_array]
+                param_arrays = ani_param_array + [gamma_in_array, log_m2l_array]
             else:
-                param_arrays = [ani_param_array, gamma_in_array, m2l_array]
+                param_arrays = [ani_param_array, gamma_in_array, log_m2l_array]
+            ParameterScalingIFU.__init__(
+                self,
+                anisotropy_model=anisotropy_model,
+                param_arrays=param_arrays,
+                scaling_grid_list=param_scaling_grid_list,
+            )
+        elif gamma_in_array is not None and log_m2l_array is None:
+            if isinstance(ani_param_array, list):
+                param_arrays = ani_param_array + [gamma_in_array]
+            else:
+                param_arrays = [ani_param_array, gamma_in_array]
             ParameterScalingIFU.__init__(
                 self,
                 anisotropy_model=anisotropy_model,
@@ -119,7 +130,7 @@ def __init__(
         self._lambda_scaling_property = lambda_scaling_property
         self._lambda_scaling_property_beta = lambda_scaling_property_beta
         self._gamma_in_array = gamma_in_array
-        self._m2l_array = m2l_array
+        self._log_m2l_array = log_m2l_array
 
     def lens_log_likelihood(
         self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None
@@ -249,9 +260,12 @@ def log_likelihood_single(
             kappa_ext=kappa_ext,
             mag_source=mag_source,
         )
-        scaling_param_array = self.draw_scaling_params(
-            kwargs_lens=kwargs_lens, **kwargs_kin
-        )
+        try:
+            scaling_param_array = self.draw_scaling_params(
+                kwargs_lens=kwargs_lens, **kwargs_kin
+            )
+        except ValueError:
+            return np.nan_to_num(-np.inf)
         kin_scaling = self.param_scaling(scaling_param_array)
 
         lnlikelihood = self.log_likelihood(
@@ -270,9 +284,12 @@ def draw_scaling_params(self, kwargs_lens=None, **kwargs_kin):
         :return: array of anisotropy parameter scaling
         """
         ani_param = self.draw_anisotropy(**kwargs_kin)
-        if self._gamma_in_array is not None and self._m2l_array is not None:
-            gamma_in, m2l = self.draw_lens_scaling_params(**kwargs_lens)
-            return np.concatenate([ani_param, [gamma_in, m2l]])
+        if self._gamma_in_array is not None and self._log_m2l_array is not None:
+            gamma_in, log_m2l = self.draw_lens_scaling_params(**kwargs_lens)
+            return np.concatenate([ani_param, [gamma_in, log_m2l]])
+        elif self._gamma_in_array is not None and self._log_m2l_array is None:
+            gamma_in = self.draw_lens_scaling_params(**kwargs_lens)
+            return np.concatenate([ani_param, [gamma_in]])
         else:
             return ani_param
 
@@ -290,9 +307,9 @@ def draw_lens_scaling_params(
         gamma_in=1,
         gamma_in_sigma=0,
         alpha_gamma_in=0,
-        m2l=1,
-        m2l_sigma=0,
-        alpha_m2l=0,
+        log_m2l=1,
+        log_m2l_sigma=0,
+        alpha_log_m2l=0,
     ):
         """Draws a realization of the anisotropy parameter scaling from the
         distribution.
@@ -314,22 +331,30 @@ def draw_lens_scaling_params(
         :param gamma_in_sigma: spread in the distribution
         :param alpha_gamma_in: float, linear slope of the gamma_in scaling relation with
             lens quantity self._lambda_scaling_property
-        :param m2l: mass-to-light ratio
-        :param m2l_sigma: spread in the distribution
-        :param alpha_m2l: float, linear slope of the m2l scaling relation with lens
-            quantity self._lambda_scaling_property
+        :param log_m2l: log(mass-to-light ratio)
+        :param log_m2l_sigma: spread in the distribution
+        :param alpha_log_m2l: float, linear slope of the log(m2l) scaling relation with
+            lens quantity self._lambda_scaling_property
         :return: draw from the distributions
         """
-        if self._gamma_in_array is not None and self._m2l_array is not None:
-            gamma_in_draw, m2l_draw = self.draw_lens_parameters(
+        if self._gamma_in_array is not None and self._log_m2l_array is not None:
+            gamma_in_draw, log_m2l_draw = self.draw_lens_parameters(
                 gamma_in + alpha_gamma_in * self._lambda_scaling_property,
                 gamma_in_sigma,
-                m2l + alpha_m2l * self._lambda_scaling_property,
-                m2l_sigma,
+                log_m2l + alpha_log_m2l * self._lambda_scaling_property,
+                log_m2l_sigma,
             )
-            return gamma_in_draw, m2l_draw
+            return gamma_in_draw, log_m2l_draw
+
+        elif self._gamma_in_array is not None and self._log_m2l_array is None:
+            gamma_in_draw = self.draw_lens_parameters(
+                gamma_in + alpha_gamma_in * self._lambda_scaling_property,
+                gamma_in_sigma,
+            )
+            return gamma_in_draw
+
         else:
-            return None, None
+            return None
 
     def angular_diameter_distances(self, cosmo):
         """Time-delay distance Ddt, angular diameter distance to the lens (dd)
@@ -413,9 +438,9 @@ def draw_lens(
         gamma_in=1,
         gamma_in_sigma=0,
         alpha_gamma_in=0,
-        m2l=1,
-        m2l_sigma=0,
-        alpha_m2l=0,
+        log_m2l=1,
+        log_m2l_sigma=0,
+        alpha_log_m2l=0,
     ):
         """Draws a realization of a specific model from the hyper-parameter
         distribution.
@@ -437,10 +462,10 @@ def draw_lens(
         :param gamma_in_sigma: spread in the distribution
         :param alpha_gamma_in: float, linear slope of the gamma_in scaling relation with
             lens quantity self._lambda_scaling_property
-        :param m2l: mass-to-light ratio
-        :param m2l_sigma: spread in the distribution
-        :param alpha_m2l: float, linear slope of the m2l scaling relation with lens
-            quantity self._lambda_scaling_property
+        :param log_m2l: log(mass-to-light ratio)
+        :param log_m2l_sigma: spread in the distribution
+        :param alpha_log_m2l: float, linear slope of the log(m2l) scaling relation with
+            lens quantity self._lambda_scaling_property
         :return: draw from the distributions
         """
         if self._mst_ifu is True:
diff --git a/hierarc/Likelihood/parameter_scaling.py b/hierarc/Likelihood/parameter_scaling.py
index 3a00f82..4739ce5 100644
--- a/hierarc/Likelihood/parameter_scaling.py
+++ b/hierarc/Likelihood/parameter_scaling.py
@@ -64,6 +64,7 @@ def __init__(
         """
         self._anisotropy_model = anisotropy_model
         self._evalute_ani = False
+        self._is_log_m2l_population_level = False
         if (
             param_arrays is not None
             and scaling_grid_list is not None
@@ -94,8 +95,9 @@ def __init__(
                     self._gamma_in_min = np.min(param_arrays[1])
                     self._gamma_in_max = np.max(param_arrays[1])
                 if self._dim_scaling > 2:
-                    self._m2l_min = np.min(param_arrays[2])
-                    self._m2l_max = np.max(param_arrays[2])
+                    self._log_m2l_min = np.min(param_arrays[2])
+                    self._log_m2l_max = np.max(param_arrays[2])
+                    self._is_log_m2l_population_level = True
 
             elif anisotropy_model == "GOM":
                 self._ani_param_min = [min(param_arrays[0]), min(param_arrays[1])]
@@ -105,8 +107,9 @@ def __init__(
                     self._gamma_in_min = np.min(param_arrays[2])
                     self._gamma_in_max = np.max(param_arrays[2])
                 if self._dim_scaling > 3:
-                    self._m2l_min = np.min(param_arrays[3])
-                    self._m2l_max = np.max(param_arrays[3])
+                    self._log_m2l_min = np.min(param_arrays[3])
+                    self._log_m2l_max = np.max(param_arrays[3])
+                    self._is_log_m2l_population_level = True
             else:
                 raise ValueError(
                     f"Anisotropy model {anisotropy_model} is not recognized!"
@@ -175,34 +178,50 @@ def draw_anisotropy(
         return None
 
     def draw_lens_parameters(
-        self, gamma_in=None, gamma_in_sigma=0, m2l=None, m2l_sigma=0
+        self, gamma_in=None, gamma_in_sigma=0, log_m2l=None, log_m2l_sigma=0
     ):
         """Draw Gaussian distribution and re-sample if outside bounds.
 
         :param gamma_in: mean of the distribution
         :param gamma_in_sigma: std of the distribution
-        :param m2l: mean of the distribution
-        :param m2l_sigma: std of the distribution
+        :param log_m2l: mean of the distribution
+        :param log_m2l_sigma: std of the distribution
         :return: random draw from the distribution
         """
-        if gamma_in < self._gamma_in_min or gamma_in > self._gamma_in_max:
-            raise ValueError(
-                "gamma_in parameter is out of bounds of the interpolated range!"
-            )
-        if m2l < self._m2l_min or m2l > self._m2l_max:
-            raise ValueError(
-                "m2l parameter is out of bounds of the interpolated range!"
-            )
-
-        gamma_in_draw = np.random.normal(gamma_in, gamma_in_sigma)
-        m2l_draw = np.random.normal(m2l, m2l_sigma)
+        if self._is_log_m2l_population_level:
+            if gamma_in < self._gamma_in_min or gamma_in > self._gamma_in_max:
+                raise ValueError(
+                    "gamma_in parameter is out of bounds of the interpolated range!"
+                )
+            if log_m2l < self._log_m2l_min or log_m2l > self._log_m2l_max:
+                raise ValueError(
+                    "m2l parameter is out of bounds of the interpolated range!"
+                )
 
-        if (
-            gamma_in_draw < self._gamma_in_min
-            or gamma_in_draw > self._gamma_in_max
-            or m2l_draw < self._m2l_min
-            or m2l_draw > self._m2l_max
-        ):
-            return self.draw_lens_parameters(gamma_in, gamma_in_sigma, m2l, m2l_sigma)
+            gamma_in_draw = np.random.normal(gamma_in, gamma_in_sigma)
+            log_m2l_draw = np.random.normal(log_m2l, log_m2l_sigma)
+
+            if (
+                gamma_in_draw < self._gamma_in_min
+                or gamma_in_draw > self._gamma_in_max
+                or log_m2l_draw < self._log_m2l_min
+                or log_m2l_draw > self._log_m2l_max
+            ):
+                return self.draw_lens_parameters(
+                    gamma_in, gamma_in_sigma, log_m2l, log_m2l_sigma
+                )
+
+            return gamma_in_draw, log_m2l_draw
+
+        else:
+            if gamma_in < self._gamma_in_min or gamma_in > self._gamma_in_max:
+                raise ValueError(
+                    "gamma_in parameter is out of bounds of the interpolated range!"
+                )
+
+            gamma_in_draw = np.random.normal(gamma_in, gamma_in_sigma)
+
+            if gamma_in_draw < self._gamma_in_min or gamma_in_draw > self._gamma_in_max:
+                return self.draw_lens_parameters(gamma_in, gamma_in_sigma)
 
-        return gamma_in_draw, m2l_draw
+            return gamma_in_draw


Test patch:
diff --git a/test/test_LensPosterior/test_kin_constraints_composite.py b/test/test_LensPosterior/test_kin_constraints_composite.py
index f855ee2..0b842b6 100644
--- a/test/test_LensPosterior/test_kin_constraints_composite.py
+++ b/test/test_LensPosterior/test_kin_constraints_composite.py
@@ -70,7 +70,7 @@ def test_likelihoodconfiguration_om(self):
         lens_light_model_list = ["SERSIC"]
 
         gamma_in_array = np.linspace(0.1, 2.9, 5)
-        m2l_array = np.linspace(0.1, 10, 5)
+        log_m2l_array = np.linspace(0.1, 1, 5)
 
         rho0_array = 10 ** np.random.normal(8, 0, 100) / 1e6
         r_s_array = np.random.normal(0.1, 0, 100)
@@ -80,9 +80,9 @@ def test_likelihoodconfiguration_om(self):
             z_lens=z_lens,
             z_source=z_source,
             gamma_in_array=gamma_in_array,
-            m2l_array=m2l_array,
-            rho0_array=rho0_array,
-            r_s_array=r_s_array,
+            log_m2l_array=log_m2l_array,
+            kappa_s_array=rho0_array,
+            r_s_angle_array=r_s_array,
             theta_E=theta_E,
             theta_E_error=0.01,
             gamma=gamma,
@@ -108,6 +108,258 @@ def test_likelihoodconfiguration_om(self):
         kwargs_kin = {"a_ani": 1}
         ln_class.lens_log_likelihood(cosmo, kwargs_lens={}, kwargs_kin=kwargs_kin)
 
+        kwargs_lens_light_test = [{"amp": [1, 1], "sigma": [1, 2]}]
+        lens_light_model_list_test = ["MULTI_GAUSSIAN"]
+
+        kin_constraints_test = KinConstraintsComposite(
+            z_lens=z_lens,
+            z_source=z_source,
+            gamma_in_array=gamma_in_array,
+            log_m2l_array=log_m2l_array,
+            kappa_s_array=rho0_array,
+            r_s_angle_array=r_s_array,
+            theta_E=theta_E,
+            theta_E_error=0.01,
+            gamma=gamma,
+            gamma_error=0.02,
+            r_eff=r_eff,
+            r_eff_error=0.05,
+            sigma_v_measured=[200],
+            sigma_v_error_independent=[10],
+            sigma_v_error_covariant=0,
+            kwargs_aperture=kwargs_aperture,
+            kwargs_seeing=kwargs_seeing,
+            anisotropy_model=anisotropy_model,
+            kwargs_lens_light=kwargs_lens_light_test,
+            lens_light_model_list=lens_light_model_list_test,
+            **kwargs_kin_api_settings
+        )
+
+        kappa_s_array = 10 ** np.random.normal(8, 0, 100) / 1e6
+        r_s_angle_array = np.random.normal(0.1, 0, 100)
+
+        kin_constraints_kappa = KinConstraintsComposite(
+            z_lens=z_lens,
+            z_source=z_source,
+            gamma_in_array=gamma_in_array,
+            log_m2l_array=log_m2l_array,
+            kappa_s_array=[],
+            r_s_angle_array=[],
+            theta_E=theta_E,
+            theta_E_error=0.01,
+            gamma=gamma,
+            gamma_error=0.02,
+            r_eff=r_eff,
+            r_eff_error=0.05,
+            sigma_v_measured=[200],
+            sigma_v_error_independent=[10],
+            sigma_v_error_covariant=0,
+            kwargs_aperture=kwargs_aperture,
+            kwargs_seeing=kwargs_seeing,
+            anisotropy_model=anisotropy_model,
+            kwargs_lens_light=kwargs_lens_light,
+            lens_light_model_list=lens_light_model_list,
+            rho0_array=kappa_s_array,
+            r_s_array=r_s_angle_array,
+            **kwargs_kin_api_settings
+        )
+
+    def test_likelihoodconfiguration_gom(self):
+        anisotropy_model = "GOM"
+        kwargs_aperture = {
+            "aperture_type": "shell",
+            "r_in": 0,
+            "r_out": 3 / 2.0,
+            "center_ra": 0.0,
+            "center_dec": 0,
+        }
+        kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4}
+
+        # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed)
+        kwargs_numerics_galkin = {
+            "interpol_grid_num": 50,  # numerical interpolation, should converge -> infinity
+            "log_integration": True,
+            # log or linear interpolation of surface brightness and mass models
+            "max_integrate": 100,
+            "min_integrate": 0.001,
+        }  # lower/upper bound of numerical integrals
+
+        # redshift
+        z_lens = 0.5
+        z_source = 1.5
+
+        cosmo = FlatLambdaCDM(H0=70, Om0=0.3)
+        theta_E = 1.0
+        r_eff = 1
+        gamma = 2.1
+
+        kwargs_mge_light = {
+            "grid_spacing": 0.1,
+            "grid_num": 10,
+            "n_comp": 20,
+            "center_x": 0,
+            "center_y": 0,
+        }
+
+        kwargs_kin_api_settings = {
+            "multi_observations": False,
+            "kwargs_numerics_galkin": kwargs_numerics_galkin,
+            "kwargs_mge_light": kwargs_mge_light,
+            "sampling_number": 10,
+            "num_kin_sampling": 10,
+            "num_psf_sampling": 10,
+        }
+
+        kwargs_lens_light = [
+            {
+                "R_sersic": 2,
+                "amp": 1,
+                "n_sersic": 2,
+                "center_x": 0,
+                "center_y": 0,
+            }
+        ]
+        lens_light_model_list = ["SERSIC"]
+
+        gamma_in_array = np.linspace(0.1, 2.9, 5)
+        log_m2l_array = np.linspace(0.1, 1, 5)
+
+        rho0_array = 10 ** np.random.normal(8, 0, 100) / 1e6
+        r_s_array = np.random.normal(0.1, 0, 100)
+
+        # compute likelihood
+        kin_constraints = KinConstraintsComposite(
+            z_lens=z_lens,
+            z_source=z_source,
+            gamma_in_array=gamma_in_array,
+            log_m2l_array=log_m2l_array,
+            kappa_s_array=rho0_array,
+            r_s_angle_array=r_s_array,
+            theta_E=theta_E,
+            theta_E_error=0.01,
+            gamma=gamma,
+            gamma_error=0.02,
+            r_eff=r_eff,
+            r_eff_error=0.05,
+            sigma_v_measured=[200],
+            sigma_v_error_independent=[10],
+            sigma_v_error_covariant=0,
+            kwargs_aperture=kwargs_aperture,
+            kwargs_seeing=kwargs_seeing,
+            anisotropy_model=anisotropy_model,
+            kwargs_lens_light=kwargs_lens_light,
+            lens_light_model_list=lens_light_model_list,
+            **kwargs_kin_api_settings
+        )
+
+        kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5)
+        kwargs_likelihood["normalized"] = False
+        ln_class = LensLikelihood(**kwargs_likelihood)
+        kwargs_kin = {"a_ani": 1, "beta_inf": 0.5}
+        ln_class.lens_log_likelihood(cosmo, kwargs_lens={}, kwargs_kin=kwargs_kin)
+
+
+class TestKinConstraintsCompositeM2l(object):
+    def setup(self):
+        pass
+
+    def test_likelihoodconfiguration_om(self):
+        anisotropy_model = "OM"
+        kwargs_aperture = {
+            "aperture_type": "shell",
+            "r_in": 0,
+            "r_out": 3 / 2.0,
+            "center_ra": 0.0,
+            "center_dec": 0,
+        }
+        kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4}
+
+        # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed)
+        kwargs_numerics_galkin = {
+            "interpol_grid_num": 50,  # numerical interpolation, should converge -> infinity
+            "log_integration": True,
+            # log or linear interpolation of surface brightness and mass models
+            "max_integrate": 100,
+            "min_integrate": 0.001,
+        }  # lower/upper bound of numerical integrals
+
+        # redshift
+        z_lens = 0.5
+        z_source = 1.5
+
+        cosmo = FlatLambdaCDM(H0=70, Om0=0.3)
+        theta_E = 1.0
+        r_eff = 1
+        gamma = 2.1
+
+        kwargs_mge_light = {
+            "grid_spacing": 0.1,
+            "grid_num": 10,
+            "n_comp": 20,
+            "center_x": 0,
+            "center_y": 0,
+        }
+
+        kwargs_kin_api_settings = {
+            "multi_observations": False,
+            "kwargs_numerics_galkin": kwargs_numerics_galkin,
+            "kwargs_mge_light": kwargs_mge_light,
+            "sampling_number": 50,
+            "num_kin_sampling": 50,
+            "num_psf_sampling": 50,
+        }
+
+        kwargs_lens_light = [
+            {
+                "R_sersic": 2,
+                "amp": 1,
+                "n_sersic": 2,
+                "center_x": 0,
+                "center_y": 0,
+            }
+        ]
+        lens_light_model_list = ["SERSIC"]
+
+        gamma_in_array = np.linspace(0.1, 2.9, 5)
+
+        log_m2l_array = np.random.uniform(0.1, 1, 100)
+        rho0_array = 10 ** np.random.normal(8, 0, 100) / 1e6
+        r_s_array = np.random.normal(0.1, 0, 100)
+
+        # compute likelihood
+        kin_constraints = KinConstraintsComposite(
+            z_lens=z_lens,
+            z_source=z_source,
+            gamma_in_array=gamma_in_array,
+            log_m2l_array=log_m2l_array,
+            kappa_s_array=rho0_array,
+            r_s_angle_array=r_s_array,
+            theta_E=theta_E,
+            theta_E_error=0.01,
+            gamma=gamma,
+            gamma_error=0.02,
+            r_eff=r_eff,
+            r_eff_error=0.05,
+            sigma_v_measured=[200],
+            sigma_v_error_independent=[10],
+            sigma_v_error_covariant=0,
+            kwargs_aperture=kwargs_aperture,
+            kwargs_seeing=kwargs_seeing,
+            anisotropy_model=anisotropy_model,
+            kwargs_lens_light=kwargs_lens_light,
+            lens_light_model_list=lens_light_model_list,
+            is_m2l_population_level=False,
+            **kwargs_kin_api_settings
+        )
+
+        kin_constraints.draw_lens(no_error=True)
+
+        kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5)
+        kwargs_likelihood["normalized"] = False
+        ln_class = LensLikelihood(**kwargs_likelihood)
+        kwargs_kin = {"a_ani": 1}
+        ln_class.lens_log_likelihood(cosmo, kwargs_lens={}, kwargs_kin=kwargs_kin)
+
     def test_likelihoodconfiguration_gom(self):
         anisotropy_model = "GOM"
         kwargs_aperture = {
@@ -166,8 +418,8 @@ def test_likelihoodconfiguration_gom(self):
         lens_light_model_list = ["SERSIC"]
 
         gamma_in_array = np.linspace(0.1, 2.9, 5)
-        m2l_array = np.linspace(0.1, 10, 5)
 
+        log_m2l_array = np.random.uniform(0.1, 1, 100)
         rho0_array = 10 ** np.random.normal(8, 0, 100) / 1e6
         r_s_array = np.random.normal(0.1, 0, 100)
 
@@ -176,9 +428,9 @@ def test_likelihoodconfiguration_gom(self):
             z_lens=z_lens,
             z_source=z_source,
             gamma_in_array=gamma_in_array,
-            m2l_array=m2l_array,
-            rho0_array=rho0_array,
-            r_s_array=r_s_array,
+            log_m2l_array=log_m2l_array,
+            kappa_s_array=rho0_array,
+            r_s_angle_array=r_s_array,
             theta_E=theta_E,
             theta_E_error=0.01,
             gamma=gamma,
@@ -193,6 +445,7 @@ def test_likelihoodconfiguration_gom(self):
             anisotropy_model=anisotropy_model,
             kwargs_lens_light=kwargs_lens_light,
             lens_light_model_list=lens_light_model_list,
+            is_m2l_population_level=False,
             **kwargs_kin_api_settings
         )
 
@@ -252,7 +505,7 @@ def test_raise(self):
             kwargs_lens_light = [{"Rs": r_eff * 0.551, "amp": 1.0}]
 
             gamma_in_array = np.linspace(0.1, 2.9, 5)
-            m2l_array = np.linspace(0.1, 10, 5)
+            log_m2l_array = np.linspace(0.1, 1, 5)
             rho0_array = 10 ** np.random.normal(8, 0.2, 100) / 1e6
             r_s_array = np.random.normal(0.1, 0.01, 100)
 
@@ -260,9 +513,9 @@ def test_raise(self):
                 z_lens=z_lens,
                 z_source=z_source,
                 gamma_in_array=gamma_in_array,
-                m2l_array=m2l_array,
-                rho0_array=rho0_array,
-                r_s_array=r_s_array,
+                log_m2l_array=log_m2l_array,
+                kappa_s_array=rho0_array,
+                r_s_angle_array=r_s_array,
                 theta_E=theta_E,
                 theta_E_error=0.01,
                 gamma=gamma,
@@ -328,7 +581,7 @@ def test_raise(self):
 
             kwargs_lens_light = [{"Rs": r_eff * 0.551, "amp": 1.0}]
             gamma_in_array = np.linspace(0.1, 2.9, 5)
-            m2l_array = np.linspace(0.1, 10, 5)
+            log_m2l_array = np.linspace(0.1, 1, 5)
 
             rho0_array = 10 ** np.random.normal(8, 0.2, 5) / 1e6
             r_s_array = np.random.normal(0.1, 0.01, 6)
@@ -337,9 +590,238 @@ def test_raise(self):
                 z_lens=z_lens,
                 z_source=z_source,
                 gamma_in_array=gamma_in_array,
-                m2l_array=m2l_array,
-                rho0_array=rho0_array,
-                r_s_array=r_s_array,
+                log_m2l_array=log_m2l_array,
+                kappa_s_array=rho0_array,
+                r_s_angle_array=r_s_array,
+                theta_E=theta_E,
+                theta_E_error=0.01,
+                gamma=gamma,
+                gamma_error=0.02,
+                r_eff=r_eff,
+                r_eff_error=0.05,
+                sigma_v_measured=[200],
+                sigma_v_error_independent=[10],
+                sigma_v_error_covariant=0,
+                kwargs_aperture=kwargs_aperture,
+                kwargs_seeing=kwargs_seeing,
+                anisotropy_model=anisotropy_model,
+                kwargs_lens_light=kwargs_lens_light,
+                **kwargs_kin_api_settings
+            )
+
+        with self.assertRaises(ValueError):
+            anisotropy_model = "OM"
+            kwargs_aperture = {
+                "aperture_type": "shell",
+                "r_in": 0,
+                "r_out": 3 / 2.0,
+                "center_ra": 0.0,
+                "center_dec": 0,
+            }
+            kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4}
+
+            # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed)
+            kwargs_numerics_galkin = {
+                "interpol_grid_num": 50,  # numerical interpolation, should converge -> infinity
+                "log_integration": True,
+                # log or linear interpolation of surface brightness and mass models
+                "max_integrate": 100,
+                "min_integrate": 0.001,
+            }  # lower/upper bound of numerical integrals
+
+            # redshift
+            z_lens = 0.5
+            z_source = 1.5
+            theta_E = 1.0
+            r_eff = 1
+            gamma = 2.1
+
+            kwargs_mge_light = {
+                "grid_spacing": 0.1,
+                "grid_num": 10,
+                "n_comp": 20,
+                "center_x": 0,
+                "center_y": 0,
+            }
+
+            # settings for kinematics calculation with KinematicsAPI of lenstronomy
+            kwargs_kin_api_settings = {
+                "multi_observations": False,
+                "kwargs_numerics_galkin": kwargs_numerics_galkin,
+                "kwargs_mge_light": kwargs_mge_light,
+                "sampling_number": 10,
+                "num_kin_sampling": 10,
+                "num_psf_sampling": 10,
+            }
+            kwargs_lens_light = [{"Rs": r_eff * 0.551, "amp": 1.0}]
+
+            gamma_in_array = np.linspace(0.1, 2.9, 5)
+            log_m2l_array = np.linspace(0.1, 1, 5)
+            kappa_s_array = 10 ** np.random.normal(8, 0.2, 100) / 1e6
+            r_s_angle_array = np.random.normal(0.1, 0.01, 101)
+
+            kin_constraints = KinConstraintsComposite(
+                z_lens=z_lens,
+                z_source=z_source,
+                gamma_in_array=gamma_in_array,
+                log_m2l_array=log_m2l_array,
+                kappa_s_array=[],
+                r_s_angle_array=[],
+                theta_E=theta_E,
+                theta_E_error=0.01,
+                gamma=gamma,
+                gamma_error=0.02,
+                r_eff=r_eff,
+                r_eff_error=0.05,
+                sigma_v_measured=[200],
+                sigma_v_error_independent=[10],
+                sigma_v_error_covariant=0,
+                kwargs_aperture=kwargs_aperture,
+                kwargs_seeing=kwargs_seeing,
+                anisotropy_model=anisotropy_model,
+                kwargs_lens_light=kwargs_lens_light,
+                rho0_array=kappa_s_array,
+                r_s_array=r_s_angle_array,
+                **kwargs_kin_api_settings
+            )
+
+    def test_raise_m2l(self):
+        with self.assertRaises(ValueError):
+            anisotropy_model = "OM"
+            kwargs_aperture = {
+                "aperture_type": "shell",
+                "r_in": 0,
+                "r_out": 3 / 2.0,
+                "center_ra": 0.0,
+                "center_dec": 0,
+            }
+            kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4}
+
+            # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed)
+            kwargs_numerics_galkin = {
+                "interpol_grid_num": 50,  # numerical interpolation, should converge -> infinity
+                "log_integration": True,
+                # log or linear interpolation of surface brightness and mass models
+                "max_integrate": 100,
+                "min_integrate": 0.001,
+            }  # lower/upper bound of numerical integrals
+
+            # redshift
+            z_lens = 0.5
+            z_source = 1.5
+            theta_E = 1.0
+            r_eff = 1
+            gamma = 2.1
+
+            kwargs_mge_light = {
+                "grid_spacing": 0.1,
+                "grid_num": 10,
+                "n_comp": 20,
+                "center_x": 0,
+                "center_y": 0,
+            }
+
+            # settings for kinematics calculation with KinematicsAPI of lenstronomy
+            kwargs_kin_api_settings = {
+                "multi_observations": False,
+                "kwargs_numerics_galkin": kwargs_numerics_galkin,
+                "kwargs_mge_light": kwargs_mge_light,
+                "sampling_number": 10,
+                "num_kin_sampling": 10,
+                "num_psf_sampling": 10,
+            }
+            kwargs_lens_light = [{"Rs": r_eff * 0.551, "amp": 1.0}]
+
+            gamma_in_array = np.linspace(0.1, 2.9, 5)
+            log_m2l_array = np.random.uniform(0.1, 1, 100)
+            rho0_array = 10 ** np.random.normal(8, 0.2, 100) / 1e6
+            r_s_array = np.random.normal(0.1, 0.01, 100)
+
+            kin_constraints = KinConstraintsComposite(
+                z_lens=z_lens,
+                z_source=z_source,
+                gamma_in_array=gamma_in_array,
+                log_m2l_array=log_m2l_array,
+                kappa_s_array=rho0_array,
+                r_s_angle_array=r_s_array,
+                theta_E=theta_E,
+                theta_E_error=0.01,
+                gamma=gamma,
+                gamma_error=0.02,
+                r_eff=r_eff,
+                r_eff_error=0.05,
+                sigma_v_measured=[200],
+                sigma_v_error_independent=[10],
+                sigma_v_error_covariant=0,
+                kwargs_aperture=kwargs_aperture,
+                kwargs_seeing=kwargs_seeing,
+                anisotropy_model=anisotropy_model,
+                kwargs_lens_light=kwargs_lens_light,
+                is_m2l_population_level=False,
+                **kwargs_kin_api_settings
+            )
+            kin_constraints._anisotropy_model = "BAD"
+            kin_constraints._anisotropy_scaling_relative(j_ani_0=1)
+
+        with self.assertRaises(ValueError):
+            anisotropy_model = "OM"
+            kwargs_aperture = {
+                "aperture_type": "shell",
+                "r_in": 0,
+                "r_out": 3 / 2.0,
+                "center_ra": 0.0,
+                "center_dec": 0,
+            }
+            kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4}
+
+            # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed)
+            kwargs_numerics_galkin = {
+                "interpol_grid_num": 100,  # numerical interpolation, should converge -> infinity
+                "log_integration": True,
+                # log or linear interpolation of surface brightness and mass models
+                "max_integrate": 100,
+                "min_integrate": 0.001,
+            }  # lower/upper bound of numerical integrals
+
+            # redshift
+            z_lens = 0.5
+            z_source = 1.5
+            theta_E = 1.0
+            r_eff = 1
+            gamma = 2.1
+
+            kwargs_mge_light = {
+                "grid_spacing": 0.1,
+                "grid_num": 10,
+                "n_comp": 20,
+                "center_x": 0,
+                "center_y": 0,
+            }
+
+            # settings for kinematics calculation with KinematicsAPI of lenstronomy
+            kwargs_kin_api_settings = {
+                "multi_observations": False,
+                "kwargs_numerics_galkin": kwargs_numerics_galkin,
+                "kwargs_mge_light": kwargs_mge_light,
+                "sampling_number": 50,
+                "num_kin_sampling": 50,
+                "num_psf_sampling": 50,
+            }
+
+            kwargs_lens_light = [{"Rs": r_eff * 0.551, "amp": 1.0}]
+            gamma_in_array = np.linspace(0.1, 2.9, 5)
+
+            log_m2l_array = np.random.uniform(0.1, 1, 6)
+            rho0_array = 10 ** np.random.normal(8, 0.2, 5) / 1e6
+            r_s_array = np.random.normal(0.1, 0.01, 5)
+
+            KinConstraintsComposite(
+                z_lens=z_lens,
+                z_source=z_source,
+                gamma_in_array=gamma_in_array,
+                log_m2l_array=log_m2l_array,
+                kappa_s_array=rho0_array,
+                r_s_angle_array=r_s_array,
                 theta_E=theta_E,
                 theta_E_error=0.01,
                 gamma=gamma,
@@ -353,9 +835,87 @@ def test_raise(self):
                 kwargs_seeing=kwargs_seeing,
                 anisotropy_model=anisotropy_model,
                 kwargs_lens_light=kwargs_lens_light,
+                is_m2l_population_level=False,
                 **kwargs_kin_api_settings
             )
 
+        with self.assertRaises(ValueError):
+            anisotropy_model = "OM"
+            kwargs_aperture = {
+                "aperture_type": "shell",
+                "r_in": 0,
+                "r_out": 3 / 2.0,
+                "center_ra": 0.0,
+                "center_dec": 0,
+            }
+            kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4}
+
+            # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed)
+            kwargs_numerics_galkin = {
+                "interpol_grid_num": 50,  # numerical interpolation, should converge -> infinity
+                "log_integration": True,
+                # log or linear interpolation of surface brightness and mass models
+                "max_integrate": 100,
+                "min_integrate": 0.001,
+            }  # lower/upper bound of numerical integrals
+
+            # redshift
+            z_lens = 0.5
+            z_source = 1.5
+            theta_E = 1.0
+            r_eff = 1
+            gamma = 2.1
+
+            kwargs_mge_light = {
+                "grid_spacing": 0.1,
+                "grid_num": 10,
+                "n_comp": 20,
+                "center_x": 0,
+                "center_y": 0,
+            }
+
+            # settings for kinematics calculation with KinematicsAPI of lenstronomy
+            kwargs_kin_api_settings = {
+                "multi_observations": False,
+                "kwargs_numerics_galkin": kwargs_numerics_galkin,
+                "kwargs_mge_light": kwargs_mge_light,
+                "sampling_number": 10,
+                "num_kin_sampling": 10,
+                "num_psf_sampling": 10,
+            }
+            kwargs_lens_light = [{"Rs": r_eff * 0.551, "amp": 1.0}]
+
+            gamma_in_array = np.linspace(0.1, 2.9, 5)
+            log_m2l_array = np.linspace(0.1, 1, 100)
+            rho0_array = 10 ** np.random.normal(8, 0.2, 100) / 1e6
+            r_s_array = np.random.normal(0.1, 0.01, 100)
+
+            kin_constraints = KinConstraintsComposite(
+                z_lens=z_lens,
+                z_source=z_source,
+                gamma_in_array=gamma_in_array,
+                log_m2l_array=log_m2l_array,
+                kappa_s_array=rho0_array,
+                r_s_angle_array=r_s_array,
+                theta_E=theta_E,
+                theta_E_error=0.01,
+                gamma=gamma,
+                gamma_error=0.02,
+                r_eff=r_eff,
+                r_eff_error=0.05,
+                sigma_v_measured=[200],
+                sigma_v_error_independent=[10],
+                sigma_v_error_covariant=0,
+                kwargs_aperture=kwargs_aperture,
+                kwargs_seeing=kwargs_seeing,
+                anisotropy_model=anisotropy_model,
+                kwargs_lens_light=kwargs_lens_light,
+                is_m2l_population_level=False,
+                **kwargs_kin_api_settings
+            )
+            kin_constraints._anisotropy_model = "BAD"
+            kin_constraints._anisotropy_scaling_relative(j_ani_0=1)
+
 
 if __name__ == "__main__":
     pytest.main()
diff --git a/test/test_Likelihood/test_hierarchy_likelihood.py b/test/test_Likelihood/test_hierarchy_likelihood.py
index a0ce60e..119a5ff 100644
--- a/test/test_Likelihood/test_hierarchy_likelihood.py
+++ b/test/test_Likelihood/test_hierarchy_likelihood.py
@@ -1,6 +1,7 @@
 from hierarc.Likelihood.hierarchy_likelihood import LensLikelihood
 from astropy.cosmology import FlatLambdaCDM
 import pytest
+import unittest
 import numpy as np
 import numpy.testing as npt
 
@@ -93,12 +94,12 @@ def setup(self):
         )
 
         gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5)
-        m2l_array = np.linspace(start=1, stop=10, num=10)
+        log_m2l_array = np.linspace(start=0.1, stop=1, num=10)
         param_scaling_array = np.multiply.outer(
             np.ones_like(ani_param_array),
-            np.outer(np.ones_like(gamma_in_array), np.ones_like(m2l_array)),
+            np.outer(np.ones_like(gamma_in_array), np.ones_like(log_m2l_array)),
         )
-        self.likelihood_gamma_in_m2l_list_ani = LensLikelihood(
+        self.likelihood_gamma_in_log_m2l_list_ani = LensLikelihood(
             z_lens,
             z_source,
             name="name",
@@ -108,7 +109,7 @@ def setup(self):
             ani_scaling_array_list=None,
             param_scaling_grid_list=[param_scaling_array],
             gamma_in_array=gamma_in_array,
-            m2l_array=m2l_array,
+            log_m2l_array=log_m2l_array,
             num_distribution_draws=200,
             kappa_ext_bias=False,
             kappa_pdf=None,
@@ -117,7 +118,7 @@ def setup(self):
             **kwargs_likelihood
         )
 
-        self.likelihood_gamma_in_m2l = LensLikelihood(
+        self.likelihood_gamma_in_log_m2l = LensLikelihood(
             z_lens,
             z_source,
             name="name",
@@ -127,7 +128,7 @@ def setup(self):
             ani_scaling_array_list=None,
             param_scaling_grid_list=[param_scaling_array],
             gamma_in_array=gamma_in_array,
-            m2l_array=m2l_array,
+            log_m2l_array=log_m2l_array,
             num_distribution_draws=200,
             kappa_ext_bias=False,
             kappa_pdf=None,
@@ -136,6 +137,26 @@ def setup(self):
             **kwargs_likelihood
         )
 
+        self.likelihood_gamma_in_fail_case = LensLikelihood(
+            z_lens,
+            z_source,
+            name="name",
+            likelihood_type="DdtDdGaussian",
+            anisotropy_model="OM",
+            ani_param_array=ani_param_array,
+            ani_scaling_array_list=None,
+            param_scaling_grid_list=[param_scaling_array],
+            gamma_in_array=gamma_in_array,
+            log_m2l_array=log_m2l_array,
+            num_distribution_draws=200,
+            kappa_ext_bias=False,
+            kappa_pdf=None,
+            kappa_bin_edges=None,
+            mst_ifu=False,
+            lambda_scaling_property=100,
+            **kwargs_likelihood
+        )
+
     def test_lens_log_likelihood(self):
         np.random.seed(42)
         kwargs_lens = {
@@ -189,23 +210,51 @@ def test_lens_log_likelihood(self):
         kwargs_test = self.likelihood._kwargs_init(kwargs=None)
         assert type(kwargs_test) is dict
 
-        gamma_in_draw, m2l_draw = self.likelihood.draw_lens_scaling_params()
+        gamma_in_draw = self.likelihood.draw_lens_scaling_params()
         assert gamma_in_draw is None
-        assert m2l_draw is None
 
         kwargs_lens = {
             "gamma_in": 1,
             "gamma_in_sigma": 0,
             "alpha_gamma_in": 0,
-            "m2l": 1,
-            "m2l_sigma": 0,
-            "alpha_m2l": 0,
+            "log_m2l": 1,
+            "log_m2l_sigma": 0,
+            "alpha_log_m2l": 0,
         }
-        ln_likelihood = self.likelihood_gamma_in_m2l.lens_log_likelihood(
+        ln_likelihood = self.likelihood_gamma_in_log_m2l.lens_log_likelihood(
             self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin
         )
         npt.assert_almost_equal(ln_likelihood, -0.0, decimal=1)
 
+        kwargs_source = self.likelihood_gamma_in_fail_case._kwargs_init(None)
+        z_apparent_m_anchor = kwargs_source.get("z_apparent_m_anchor", 0.1)
+        delta_lum_dist = self.likelihood_gamma_in_fail_case.luminosity_distance_modulus(
+            self.cosmo, z_apparent_m_anchor
+        )
+
+        z_lens = 0.5
+        z_source = 1.5
+
+        kwargs_lens = {
+            "gamma_in": 1,
+            "gamma_in_sigma": 0,
+            "alpha_gamma_in": 0,
+            "log_m2l": 1,
+            "log_m2l_sigma": 0,
+            "alpha_log_m2l": 1000,
+        }
+
+        dd = self.cosmo.angular_diameter_distance(z=z_lens).value
+        ds = self.cosmo.angular_diameter_distance(z=z_source).value
+        dds = self.cosmo.angular_diameter_distance_z1z2(z1=z_lens, z2=z_source).value
+        ddt = (1.0 + z_lens) * dd * ds / dds
+
+        ln_likelihood = self.likelihood_gamma_in_fail_case.log_likelihood_single(
+            ddt, dd, delta_lum_dist, kwargs_lens, kwargs_kin, kwargs_source
+        )
+
+        assert ln_likelihood < -10000000
+
 
 if __name__ == "__main__":
     pytest.main()
diff --git a/test/test_Likelihood/test_parameter_scaling.py b/test/test_Likelihood/test_parameter_scaling.py
index dad4f0c..508fcb6 100644
--- a/test/test_Likelihood/test_parameter_scaling.py
+++ b/test/test_Likelihood/test_parameter_scaling.py
@@ -26,12 +26,12 @@ def setup(self):
 
         ani_param_array = np.linspace(start=0, stop=1, num=10)
         gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5)
-        m2l_array = np.linspace(start=1, stop=10, num=10)
+        log_m2l_array = np.linspace(start=0.1, stop=1, num=10)
 
-        param_arrays = [ani_param_array, gamma_in_array, m2l_array]
+        param_arrays = [ani_param_array, gamma_in_array, log_m2l_array]
         param_scaling_array = np.multiply.outer(
             ani_param_array,
-            np.outer(gamma_in_array, m2l_array),
+            np.outer(gamma_in_array, log_m2l_array),
         )
         self.scaling_nfw = ParameterScalingSingleAperture(
             param_arrays, param_scaling_array
@@ -45,13 +45,13 @@ def setup(self):
             gom_param_array[0],
             gom_param_array[1],
             gamma_in_array,
-            m2l_array,
+            log_m2l_array,
         ]
         param_scaling_array = np.multiply.outer(
             gom_param_array[0],
             np.multiply.outer(
                 gom_param_array[1],
-                np.outer(gamma_in_array, m2l_array),
+                np.outer(gamma_in_array, log_m2l_array),
             ),
         )
 
@@ -69,11 +69,11 @@ def test_param_scaling(self):
         scaling = self.scaling_2d.param_scaling(param_array=[1, 2])
         assert scaling == 2
 
-        scaling = self.scaling_nfw.param_scaling(param_array=[1, 2.9, 10])
-        assert scaling == 1 * 2.9 * 10
+        scaling = self.scaling_nfw.param_scaling(param_array=[1, 2.9, 0.5])
+        assert scaling == 1 * 2.9 * 0.5
 
-        scaling = self.scaling_nfw_2d.param_scaling(param_array=[1, 2, 2.9, 10])
-        assert scaling == 1 * 2 * 2.9 * 10
+        scaling = self.scaling_nfw_2d.param_scaling(param_array=[1, 2, 2.9, 0.5])
+        assert scaling == 1 * 2 * 2.9 * 0.5
 
 
 class TestParameterScalingIFU(object):
@@ -99,12 +99,12 @@ def setup(self):
 
         ani_param_array = np.linspace(start=0, stop=1, num=10)
         gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5)
-        m2l_array = np.linspace(start=1, stop=10, num=10)
+        log_m2l_array = np.linspace(start=0.1, stop=1, num=10)
 
-        param_arrays = [ani_param_array, gamma_in_array, m2l_array]
+        param_arrays = [ani_param_array, gamma_in_array, log_m2l_array]
         param_scaling_array = np.multiply.outer(
             ani_param_array,
-            np.outer(gamma_in_array, m2l_array),
+            np.outer(gamma_in_array, log_m2l_array),
         )
         self.scaling_nfw = ParameterScalingIFU(
             anisotropy_model="OM",
@@ -120,13 +120,13 @@ def setup(self):
             gom_param_array[0],
             gom_param_array[1],
             gamma_in_array,
-            m2l_array,
+            log_m2l_array,
         ]
         param_scaling_array = np.multiply.outer(
             gom_param_array[0],
             np.multiply.outer(
                 gom_param_array[1],
-                np.outer(gamma_in_array, m2l_array),
+                np.outer(gamma_in_array, log_m2l_array),
             ),
         )
         self.scaling_nfw_2d = ParameterScalingIFU(
@@ -135,6 +135,39 @@ def setup(self):
             scaling_grid_list=[param_scaling_array],
         )
 
+        param_arrays = [ani_param_array, gamma_in_array]
+        param_scaling_array = np.multiply.outer(
+            ani_param_array,
+            gamma_in_array,
+        )
+        self.scaling_nfw_no_m2l = ParameterScalingIFU(
+            anisotropy_model="OM",
+            param_arrays=param_arrays,
+            scaling_grid_list=[param_scaling_array],
+        )
+
+        gom_param_array = [
+            np.linspace(start=0, stop=1, num=10),
+            np.linspace(start=1, stop=2, num=5),
+        ]
+        param_arrays = [
+            gom_param_array[0],
+            gom_param_array[1],
+            gamma_in_array,
+        ]
+        param_scaling_array = np.multiply.outer(
+            gom_param_array[0],
+            np.multiply.outer(
+                gom_param_array[1],
+                gamma_in_array,
+            ),
+        )
+        self.scaling_nfw_2d_no_m2l = ParameterScalingIFU(
+            anisotropy_model="GOM",
+            param_arrays=param_arrays,
+            scaling_grid_list=[param_scaling_array],
+        )
+
     def test_param_scaling(self):
         scaling = self.scaling.param_scaling(param_array=None)
         assert scaling[0] == 1
@@ -145,11 +178,17 @@ def test_param_scaling(self):
         scaling = self.scaling_2d.param_scaling(param_array=[1, 2])
         assert scaling[0] == 2
 
-        scaling = self.scaling_nfw.param_scaling(param_array=[1, 2.9, 10])
-        assert scaling[0] == 1 * 2.9 * 10
+        scaling = self.scaling_nfw.param_scaling(param_array=[1, 2.9, 0.5])
+        assert scaling[0] == 1 * 2.9 * 0.5
+
+        scaling = self.scaling_nfw_2d.param_scaling(param_array=[1, 2, 2.9, 0.5])
+        assert scaling[0] == 1 * 2 * 2.9 * 0.5
+
+        scaling = self.scaling_nfw_no_m2l.param_scaling(param_array=[1, 2.9])
+        assert scaling[0] == 1 * 2.9
 
-        scaling = self.scaling_nfw_2d.param_scaling(param_array=[1, 2, 2.9, 10])
-        assert scaling[0] == 1 * 2 * 2.9 * 10
+        scaling = self.scaling_nfw_2d_no_m2l.param_scaling(param_array=[1, 2, 2.9])
+        assert scaling[0] == 1 * 2 * 2.9
 
     def test_draw_anisotropy(self):
         a_ani = 1
@@ -197,6 +236,25 @@ def test_draw_anisotropy(self):
                 a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1
             )
 
+        param_draw = self.scaling_nfw_no_m2l.draw_anisotropy(
+            a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0
+        )
+        assert param_draw[0] == a_ani
+        for i in range(100):
+            param_draw = self.scaling_nfw_no_m2l.draw_anisotropy(
+                a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1
+            )
+
+        param_draw = self.scaling_nfw_2d_no_m2l.draw_anisotropy(
+            a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0
+        )
+        assert param_draw[0] == a_ani
+        assert param_draw[1] == beta_inf
+        for i in range(100):
+            param_draw = self.scaling_nfw_2d_no_m2l.draw_anisotropy(
+                a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1
+            )
+
         scaling = ParameterScalingIFU(anisotropy_model="NONE")
         param_draw = scaling.draw_anisotropy(
             a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0
@@ -205,23 +263,42 @@ def test_draw_anisotropy(self):
 
     def test_draw_lens_parameters(self):
         param_draw = self.scaling_nfw.draw_lens_parameters(
-            gamma_in=1, gamma_in_sigma=0, m2l=5, m2l_sigma=0
+            gamma_in=1, gamma_in_sigma=0, log_m2l=0.5, log_m2l_sigma=0
         )
         assert param_draw[0] == 1
-        assert param_draw[1] == 5
+        assert param_draw[1] == 0.5
         for i in range(100):
             param_draw = self.scaling_nfw.draw_lens_parameters(
-                gamma_in=1, gamma_in_sigma=1, m2l=5, m2l_sigma=3
+                gamma_in=1, gamma_in_sigma=1, log_m2l=0.5, log_m2l_sigma=3
             )
 
         param_draw = self.scaling_nfw_2d.draw_lens_parameters(
-            gamma_in=1, gamma_in_sigma=0, m2l=5, m2l_sigma=0
+            gamma_in=1, gamma_in_sigma=0, log_m2l=0.5, log_m2l_sigma=0
         )
         assert param_draw[0] == 1
-        assert param_draw[1] == 5
+        assert param_draw[1] == 0.5
+
         for i in range(100):
             param_draw = self.scaling_nfw_2d.draw_lens_parameters(
-                gamma_in=1, gamma_in_sigma=1, m2l=5, m2l_sigma=3
+                gamma_in=1, gamma_in_sigma=1, log_m2l=0.5, log_m2l_sigma=3
+            )
+
+        param_draw = self.scaling_nfw_no_m2l.draw_lens_parameters(
+            gamma_in=1, gamma_in_sigma=0, log_m2l=0.5, log_m2l_sigma=0
+        )
+        assert param_draw == 1
+        for i in range(100):
+            param_draw = self.scaling_nfw_no_m2l.draw_lens_parameters(
+                gamma_in=1, gamma_in_sigma=1, log_m2l=0.5, log_m2l_sigma=3
+            )
+
+        param_draw = self.scaling_nfw_2d_no_m2l.draw_lens_parameters(
+            gamma_in=1, gamma_in_sigma=0, log_m2l=0.5, log_m2l_sigma=0
+        )
+        assert param_draw == 1
+        for i in range(100):
+            param_draw = self.scaling_nfw_2d_no_m2l.draw_lens_parameters(
+                gamma_in=1, gamma_in_sigma=1, log_m2l=0.5, log_m2l_sigma=3
             )
 
 
@@ -263,12 +340,12 @@ def test_raise(self):
 
         ani_param_array = np.linspace(start=0, stop=1, num=10)
         gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5)
-        m2l_array = np.linspace(start=1, stop=10, num=10)
+        log_m2l_array = np.linspace(start=0.1, stop=1, num=10)
 
-        param_arrays = [ani_param_array, gamma_in_array, m2l_array]
+        param_arrays = [ani_param_array, gamma_in_array, log_m2l_array]
         param_scaling_array = np.multiply.outer(
             ani_param_array,
-            np.multiply.outer(gamma_in_array, m2l_array),
+            np.multiply.outer(gamma_in_array, log_m2l_array),
         )
         print(param_scaling_array.shape)
         scaling = ParameterScalingIFU(
@@ -278,21 +355,21 @@ def test_raise(self):
         )
         with self.assertRaises(ValueError):
             scaling.draw_lens_parameters(
-                gamma_in=-1, gamma_in_sigma=0.1, m2l=5, m2l_sigma=0.1
+                gamma_in=-1, gamma_in_sigma=0.1, log_m2l=0.5, log_m2l_sigma=0.1
             )
         with self.assertRaises(ValueError):
             scaling.draw_lens_parameters(
-                gamma_in=1, gamma_in_sigma=0.1, m2l=-5, m2l_sigma=0.1
+                gamma_in=1, gamma_in_sigma=0.1, log_m2l=-0.5, log_m2l_sigma=0.1
             )
 
         ani_param_array = np.linspace(start=0, stop=1, num=10)
         gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5)
-        m2l_array = np.linspace(start=1, stop=10, num=10)
+        log_m2l_array = np.linspace(start=0.1, stop=1, num=10)
 
-        param_arrays = [ani_param_array, gamma_in_array, m2l_array]
+        param_arrays = [ani_param_array, gamma_in_array, log_m2l_array]
         param_scaling_array = np.multiply.outer(
             ani_param_array,
-            np.multiply.outer(gamma_in_array, m2l_array),
+            np.multiply.outer(gamma_in_array, log_m2l_array),
         )
 
         scaling = ParameterScalingIFU(
@@ -302,7 +379,37 @@ def test_raise(self):
         )
         with self.assertRaises(ValueError):
             scaling.draw_lens_parameters(
-                gamma_in=1, gamma_in_sigma=0, m2l=0, m2l_sigma=0
+                gamma_in=1, gamma_in_sigma=0, log_m2l=0, log_m2l_sigma=0
+            )
+
+        gom_param_array = [
+            np.linspace(start=0, stop=1, num=10),
+            np.linspace(start=1, stop=2, num=5),
+        ]
+
+        gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5)
+
+        param_arrays = [
+            gom_param_array[0],
+            gom_param_array[1],
+            gamma_in_array,
+        ]
+        param_scaling_array = np.multiply.outer(
+            gom_param_array[0],
+            np.multiply.outer(
+                gom_param_array[1],
+                gamma_in_array,
+            ),
+        )
+        self.scaling_nfw_2d_no_m2l = ParameterScalingIFU(
+            anisotropy_model="GOM",
+            param_arrays=param_arrays,
+            scaling_grid_list=[param_scaling_array],
+        )
+
+        with self.assertRaises(ValueError):
+            param_draw = self.scaling_nfw_2d_no_m2l.draw_lens_parameters(
+                gamma_in=-1, gamma_in_sigma=1, log_m2l=0.5, log_m2l_sigma=3
             )