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/base_config.py b/hierarc/LensPosterior/base_config.py
index 832cca6..9fef09e 100644
--- a/hierarc/LensPosterior/base_config.py
+++ b/hierarc/LensPosterior/base_config.py
@@ -22,6 +22,7 @@ def __init__(
         kwargs_seeing,
         kwargs_numerics_galkin,
         anisotropy_model,
+        lens_model_list=None,
         kwargs_lens_light=None,
         lens_light_model_list=["HERNQUIST"],
         MGE_light=False,
@@ -31,6 +32,7 @@ def __init__(
         num_psf_sampling=100,
         num_kin_sampling=1000,
         multi_observations=False,
+        cosmo_fiducial=None,
     ):
         """
 
@@ -42,27 +44,43 @@ def __init__(
         :param gamma_error: 1-sigma uncertainty on power-law slope
         :param r_eff: half-light radius of the deflector (arc seconds)
         :param r_eff_error: uncertainty on half-light radius
-        :param kwargs_aperture: spectroscopic aperture keyword arguments, see lenstronomy.Galkin.aperture for options
-        :param kwargs_seeing: seeing condition of spectroscopic observation, corresponds to kwargs_psf in the GalKin module specified in lenstronomy.GalKin.psf
-        :param kwargs_numerics_galkin: numerical settings for the integrated line-of-sight velocity dispersion
-        :param anisotropy_model: type of stellar anisotropy model. See details in MamonLokasAnisotropy() class of lenstronomy.GalKin.anisotropy
-        :param multi_observations: bool, if True, interprets kwargs_aperture and kwargs_seeing as lists of multiple
-         observations
+        :param kwargs_aperture: spectroscopic aperture keyword arguments, see
+        lenstronomy.Galkin.aperture for options
+        :param kwargs_seeing: seeing condition of spectroscopic observation, corresponds
+            to kwargs_psf in the GalKin module specified in lenstronomy.GalKin.psf
+        :param kwargs_numerics_galkin: numerical settings for the integrated
+            line-of-sight velocity dispersion
+        :param anisotropy_model: type of stellar anisotropy model. See details in
+            MamonLokasAnisotropy() class of lenstronomy.GalKin.anisotropy
+        :param multi_observations: bool, if True, interprets kwargs_aperture and
+            kwargs_seeing as lists of multiple observations
+        :param lens_model_list: keyword argument list of lens model (optional)
         :param kwargs_lens_light: keyword argument list of lens light model (optional)
-        :param kwargs_mge_light: keyword arguments that go into the MGE decomposition routine
-        :param hernquist_approx: bool, if True, uses the Hernquist approximation for the light profile
+        :param kwargs_mge_light: keyword arguments that go into the MGE decomposition
+            routine
+        :param hernquist_approx: bool, if True, uses the Hernquist approximation for the
+            light profile
+        :param cosmo_fiducial: astropy.cosmology instance, if None,
+            uses astropy's default cosmology
         """
         self._z_lens, self._z_source = z_lens, z_source
-        kwargs_model = {
-            "lens_model_list": ["SPP"],
-            "lens_light_model_list": lens_light_model_list,
-        }
+
+        if lens_model_list is None:
+            kwargs_model = {
+                "lens_model_list": ["SPP"],
+                "lens_light_model_list": lens_light_model_list,
+            }
+        else:
+            kwargs_model = {
+                "lens_model_list": lens_model_list,
+                "lens_light_model_list": lens_light_model_list,
+            }
         TDCosmography.__init__(
             self,
             z_lens,
             z_source,
             kwargs_model,
-            cosmo_fiducial=None,
+            cosmo_fiducial=cosmo_fiducial,
             lens_model_kinematics_bool=None,
             light_model_kinematics_bool=None,
             kwargs_seeing=kwargs_seeing,
diff --git a/hierarc/LensPosterior/kin_constraints.py b/hierarc/LensPosterior/kin_constraints.py
index fe1e28a..74147ff 100644
--- a/hierarc/LensPosterior/kin_constraints.py
+++ b/hierarc/LensPosterior/kin_constraints.py
@@ -27,6 +27,7 @@ def __init__(
         sigma_v_error_cov_matrix=None,
         kwargs_lens_light=None,
         lens_light_model_list=["HERNQUIST"],
+        lens_model_list=None,
         MGE_light=False,
         kwargs_mge_light=None,
         hernquist_approx=True,
@@ -34,6 +35,7 @@ def __init__(
         num_psf_sampling=100,
         num_kin_sampling=1000,
         multi_observations=False,
+        cosmo_fiducial=None,
     ):
         """
 
@@ -45,23 +47,36 @@ def __init__(
         :param gamma_error: 1-sigma uncertainty on power-law slope
         :param r_eff: half-light radius of the deflector (arc seconds)
         :param r_eff_error: uncertainty on half-light radius
-        :param sigma_v_measured: numpy array of IFU velocity dispersion of the main deflector in km/s
-        :param sigma_v_error_independent: numpy array of 1-sigma uncertainty in velocity dispersion of the IFU
+        :param sigma_v_measured: numpy array of IFU velocity dispersion of the main
+            deflector in km/s
+        :param sigma_v_error_independent: numpy array of 1-sigma uncertainty in velocity
+            dispersion of the IFU
          observation independent of each other
-        :param sigma_v_error_covariant: covariant error in the measured kinematics shared among all IFU measurements
-        :param sigma_v_error_cov_matrix: error covariance matrix in the sigma_v measurements (km/s)^2
-        :type sigma_v_error_cov_matrix: nxn matrix with n the length of the sigma_v_measured array
-        :param kwargs_aperture: spectroscopic aperture keyword arguments, see lenstronomy.Galkin.aperture for options
-        :param kwargs_seeing: seeing condition of spectroscopic observation, corresponds to kwargs_psf in the GalKin
+        :param sigma_v_error_covariant: covariant error in the measured kinematics
+            shared among all IFU measurements
+        :param sigma_v_error_cov_matrix: error covariance matrix in the sigma_v
+            measurements (km/s)^2
+        :type sigma_v_error_cov_matrix: nxn matrix with n the length of the
+            sigma_v_measured array
+        :param kwargs_aperture: spectroscopic aperture keyword arguments, see
+            lenstronomy.Galkin.aperture for options
+        :param kwargs_seeing: seeing condition of spectroscopic observation, corresponds
+            to kwargs_psf in the GalKin
          module specified in lenstronomy.GalKin.psf
-        :param kwargs_numerics_galkin: numerical settings for the integrated line-of-sight velocity dispersion
-        :param anisotropy_model: type of stellar anisotropy model. See details in MamonLokasAnisotropy() class of
-         lenstronomy.GalKin.anisotropy
+        :param kwargs_numerics_galkin: numerical settings for the integrated
+            line-of-sight velocity dispersion
+        :param anisotropy_model: type of stellar anisotropy model. See details in
+            MamonLokasAnisotropy() class of lenstronomy.GalKin.anisotropy
+        :param lens_model_list: keyword argument list of lens model (optional)
         :param kwargs_lens_light: keyword argument list of lens light model (optional)
-        :param kwargs_mge_light: keyword arguments that go into the MGE decomposition routine
-        :param hernquist_approx: bool, if True, uses the Hernquist approximation for the light profile
-        :param multi_observations: bool, if True, interprets kwargs_aperture and kwargs_seeing as lists of multiple
-         observations
+        :param kwargs_mge_light: keyword arguments that go into the MGE decomposition
+            routine
+        :param hernquist_approx: bool, if True, uses the Hernquist approximation for the
+            light profile
+        :param multi_observations: bool, if True, interprets kwargs_aperture and
+            kwargs_seeing as lists of multiple observations
+        :param cosmo_fiducial: astropy.cosmology instance, if None,
+            uses astropy's default
         """
         self._sigma_v_measured = np.array(sigma_v_measured)
         self._sigma_v_error_independent = np.array(sigma_v_error_independent)
@@ -70,6 +85,7 @@ def __init__(
 
         self._kwargs_lens_light = kwargs_lens_light
         self._anisotropy_model = anisotropy_model
+
         BaseLensConfig.__init__(
             self,
             z_lens,
@@ -84,6 +100,7 @@ def __init__(
             kwargs_seeing,
             kwargs_numerics_galkin,
             anisotropy_model,
+            lens_model_list=lens_model_list,
             kwargs_lens_light=kwargs_lens_light,
             lens_light_model_list=lens_light_model_list,
             MGE_light=MGE_light,
@@ -93,6 +110,7 @@ def __init__(
             num_psf_sampling=num_psf_sampling,
             num_kin_sampling=num_kin_sampling,
             multi_observations=multi_observations,
+            cosmo_fiducial=cosmo_fiducial,
         )
 
     def j_kin_draw(self, kwargs_anisotropy, no_error=False):
@@ -161,9 +179,10 @@ def hierarchy_configuration(self, num_sample_model=20):
     def model_marginalization(self, num_sample_model=20):
         """
 
-        :param num_sample_model: number of samples drawn from the lens and light model posterior to compute the
-         dimensionless kinematic component J()
-        :return: J() as array for each measurement prediction, covariance matrix in sqrt(J)
+        :param num_sample_model: number of samples drawn from the lens and light model
+            posterior to compute the dimensionless kinematic component J()
+        :return: J() as array for each measurement prediction, covariance matrix in
+            sqrt(J)
         """
         num_data = len(self._sigma_v_measured)
         j_kin_matrix = np.zeros(
diff --git a/hierarc/LensPosterior/kin_constraints_composite.py b/hierarc/LensPosterior/kin_constraints_composite.py
new file mode 100644
index 0000000..f7b2e77
--- /dev/null
+++ b/hierarc/LensPosterior/kin_constraints_composite.py
@@ -0,0 +1,355 @@
+__author__ = "ajshjib"
+
+import copy
+
+import numpy as np
+from hierarc.LensPosterior.kin_constraints import KinConstraints
+from lenstronomy.Util import constants as const
+from lenstronomy.Analysis.light_profile import LightProfileAnalysis
+from lenstronomy.LightModel.light_model import LightModel
+
+
+class KinConstraintsComposite(KinConstraints):
+    def __init__(
+        self,
+        z_lens,
+        z_source,
+        gamma_in_array,
+        m2l_array,
+        rho0_array,
+        r_s_array,
+        theta_E,
+        theta_E_error,
+        gamma,
+        gamma_error,
+        r_eff,
+        r_eff_error,
+        sigma_v_measured,
+        kwargs_aperture,
+        kwargs_seeing,
+        kwargs_numerics_galkin,
+        anisotropy_model,
+        sigma_v_error_independent=None,
+        sigma_v_error_covariant=None,
+        sigma_v_error_cov_matrix=None,
+        kwargs_lens_light=None,
+        lens_light_model_list=["HERNQUIST"],
+        kwargs_mge_light=None,
+        sampling_number=1000,
+        num_psf_sampling=100,
+        num_kin_sampling=1000,
+        multi_observations=False,
+    ):
+        """
+
+        :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
+            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 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
+        :param gamma_error: 1-sigma uncertainty on power-law slope
+        :param r_eff: half-light radius of the deflector (arc seconds)
+        :param r_eff_error: uncertainty on half-light radius
+        :param sigma_v_measured: numpy array of IFU velocity dispersion of the main
+            deflector in km/s
+        :param sigma_v_error_independent: numpy array of 1-sigma uncertainty in velocity
+            dispersion of the IFU
+         observation independent of each other
+        :param sigma_v_error_covariant: covariant error in the measured kinematics
+            shared among all IFU measurements
+        :param sigma_v_error_cov_matrix: error covariance matrix in the sigma_v
+            measurements (km/s)^2
+        :type sigma_v_error_cov_matrix: nxn matrix with n the length of the
+            sigma_v_measured array
+        :param kwargs_aperture: spectroscopic aperture keyword arguments, see
+            lenstronomy.Galkin.aperture for options
+        :param kwargs_seeing: seeing condition of spectroscopic observation, corresponds
+            to kwargs_psf in the GalKin module specified in lenstronomy.GalKin.psf
+        :param kwargs_numerics_galkin: numerical settings for the integrated
+            line-of-sight velocity dispersion
+        :param anisotropy_model: type of stellar anisotropy model. See details in
+            MamonLokasAnisotropy() class of lenstronomy.GalKin.anisotropy
+        :param kwargs_lens_light: keyword argument list of lens light model (optional)
+        :param kwargs_mge_light: keyword arguments that go into the MGE decomposition
+            routine
+        :param multi_observations: bool, if True, interprets kwargs_aperture and
+            kwargs_seeing as lists of multiple observations
+        """
+        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_model_list = ["GNFW", "MULTI_GAUSSIAN_KAPPA"]
+
+        super(KinConstraintsComposite, self).__init__(
+            z_lens,
+            z_source,
+            theta_E,
+            theta_E_error,
+            gamma,
+            gamma_error,
+            r_eff,
+            r_eff_error,
+            sigma_v_measured,
+            kwargs_aperture,
+            kwargs_seeing,
+            kwargs_numerics_galkin,
+            anisotropy_model,
+            sigma_v_error_independent=sigma_v_error_independent,
+            sigma_v_error_covariant=sigma_v_error_covariant,
+            sigma_v_error_cov_matrix=sigma_v_error_cov_matrix,
+            kwargs_lens_light=kwargs_lens_light,
+            lens_light_model_list=lens_light_model_list,
+            lens_model_list=lens_model_list,
+            MGE_light=False,  # set False, as MGE is already done as default
+            kwargs_mge_light=None,
+            hernquist_approx=False,
+            sampling_number=sampling_number,
+            num_psf_sampling=num_psf_sampling,
+            num_kin_sampling=num_kin_sampling,
+            multi_observations=multi_observations,
+        )
+
+        if len(rho0_array) != len(r_s_array):
+            raise ValueError("rho0 and r_s must have 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
+
+    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.
+
+        :param rho_s: halo mass normalization in M_sun / Mpc^3
+        :param r_scale: halo scale radius in arc seconds
+        :return: surface mass density divided by the critical density
+        """
+        r_s_angle = r_scale / self.lensCosmo.dd / const.arcsec  # Rs in arcsec
+        kappa_s = rho_s * r_scale / self.lensCosmo.sigma_crit
+
+        return kappa_s, r_s_angle
+
+    def draw_lens(self, no_error=False):
+        """Draws a lens model from the posterior.
+
+        :param no_error: bool, if True, does not render from the uncertainty but uses
+            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,
+            )
+
+        random_index = np.random.randint(low=0, high=len(self._rho_s_array))
+        kappa_s_draw = self._kappa_s_array[random_index]
+        r_scale_angle_draw = self._r_scale_angle_array[random_index]
+
+        # we make sure no negative r_eff are being sampled
+        delta_r_eff = np.maximum(
+            np.random.normal(loc=1, scale=self._r_eff_error / self._r_eff), 0.001
+        )
+        r_eff_draw = delta_r_eff * self._r_eff
+
+        return kappa_s_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff
+
+    def model_marginalization(self, num_sample_model=20):
+        """
+
+        :param num_sample_model: number of samples drawn from the lens and light model
+            posterior to compute the dimensionless kinematic component J()
+        :return: J() as array for each measurement prediction, covariance matrix in
+            sqrt(J)
+        """
+        num_data = len(self._sigma_v_measured)
+        j_kin_matrix = np.zeros(
+            (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,
+            )
+            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)
+        return j_model_list, error_cov_j_sqrt
+
+    def j_kin_draw_composite(self, kwargs_anisotropy, gamma_in, 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 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, 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"] *= m2l
+
+        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,
+        ]
+
+        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 hierarchy_configuration(self, num_sample_model=20):
+        """Routine to configure the likelihood to be used in the hierarchical sampling.
+        In particular, a default configuration is set to compute the Gaussian
+        approximation of Ds/Dds by sampling the posterior and the estimate of the
+        variance of the sample. The anisotropy scaling is then performed. Different
+        anisotropy models are supported.
+
+        :param num_sample_model: number of samples drawn from the lens and light model
+            posterior to compute the dimensionless kinematic component J()
+        :return: keyword arguments
+        """
+
+        j_model_list, error_cov_j_sqrt = self.model_marginalization(num_sample_model)
+        ani_scaling_grid_list = self.anisotropy_scaling()
+
+        error_cov_measurement = self.error_cov_measurement
+        # configuration keyword arguments for the hierarchical sampling
+        kwargs_likelihood = {
+            "z_lens": self._z_lens,
+            "z_source": self._z_source,
+            "likelihood_type": "IFUKinCov",
+            "sigma_v_measurement": self._sigma_v_measured,
+            "anisotropy_model": self._anisotropy_model,
+            "j_model": j_model_list,
+            "error_cov_measurement": error_cov_measurement,
+            "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,
+            "param_scaling_grid_list": ani_scaling_grid_list,
+        }
+        return kwargs_likelihood
+
+    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)
+
+    def _anisotropy_scaling_relative(self, j_ani_0):
+        """Anisotropy scaling relative to a default J prediction.
+
+        :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),
+                        len(self.m2l_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):
+                        for l, m2l in enumerate(self.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
+                            )
+
+                            for m, j_kin in enumerate(j_kin_ani):
+                                ani_scaling_grid_list[m][i, j, k, l] = (
+                                    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),
+                        len(self.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):
+                        kwargs_anisotropy = self.anisotropy_kwargs(a_ani)
+                        j_kin_ani = self.j_kin_draw_composite(
+                            kwargs_anisotropy, g_in, 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
diff --git a/hierarc/Likelihood/LensLikelihood/base_lens_likelihood.py b/hierarc/Likelihood/LensLikelihood/base_lens_likelihood.py
index 854b7f4..c618311 100644
--- a/hierarc/Likelihood/LensLikelihood/base_lens_likelihood.py
+++ b/hierarc/Likelihood/LensLikelihood/base_lens_likelihood.py
@@ -154,13 +154,13 @@ def num_data(self):
         return self._lens_type.num_data
 
     def log_likelihood(
-        self, ddt, dd, aniso_scaling=None, sigma_v_sys_error=None, mu_intrinsic=None
+        self, ddt, dd, kin_scaling=None, sigma_v_sys_error=None, mu_intrinsic=None
     ):
         """
 
         :param ddt: time-delay distance [physical Mpc]
         :param dd: angular diameter distance to the lens [physical Mpc]
-        :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
+        :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
          dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the
          anisotropy model used to derive the prediction and covariance matrix in the init of this class.
         :param sigma_v_sys_error: unaccounted uncertainty in the velocity dispersion measurement
@@ -175,12 +175,12 @@ def log_likelihood(
         ]:
             return self._lens_type.log_likelihood(ddt, dd)
         elif self.likelihood_type in ["DdtDdKDE", "DdtDdGaussian", "DsDdsGaussian"]:
-            return self._lens_type.log_likelihood(ddt, dd, aniso_scaling=aniso_scaling)
+            return self._lens_type.log_likelihood(ddt, dd, kin_scaling=kin_scaling)
         elif self.likelihood_type in ["DdtHistKin", "IFUKinCov", "DdtGaussKin"]:
             return self._lens_type.log_likelihood(
                 ddt,
                 dd,
-                aniso_scaling=aniso_scaling,
+                kin_scaling=kin_scaling,
                 sigma_v_sys_error=sigma_v_sys_error,
             )
         elif self.likelihood_type in ["Mag"]:
@@ -220,14 +220,14 @@ def sigma_v_measurement(self, sigma_v_sys_error=None):
             )
         return None, None
 
-    def sigma_v_prediction(self, ddt, dd, aniso_scaling=None):
+    def sigma_v_prediction(self, ddt, dd, kin_scaling=None):
         """
 
         :param ddt: ddt in physical Mpc
         :param dd: dd in physical Mpc
-        :param aniso_scaling: anisotropy scaling in J
+        :param kin_scaling: anisotropy scaling in J
         :return: model predicted velocity dispersion (vector) and model covariance matrix thereof
         """
         if self.likelihood_type in ["DdtHistKin", "IFUKinCov", "DdtGaussKin"]:
-            return self._lens_type.sigma_v_prediction(ddt, dd, aniso_scaling)
+            return self._lens_type.sigma_v_prediction(ddt, dd, kin_scaling)
         return None, None
diff --git a/hierarc/Likelihood/LensLikelihood/ddt_dd_gauss_likelihood.py b/hierarc/Likelihood/LensLikelihood/ddt_dd_gauss_likelihood.py
index 8730d07..f825300 100644
--- a/hierarc/Likelihood/LensLikelihood/ddt_dd_gauss_likelihood.py
+++ b/hierarc/Likelihood/LensLikelihood/ddt_dd_gauss_likelihood.py
@@ -25,18 +25,18 @@ def __init__(self, z_lens, z_source, ddt_mean, ddt_sigma, dd_mean, dd_sigma):
         )
         self.num_data = 2
 
-    def log_likelihood(self, ddt, dd, aniso_scaling=None):
+    def log_likelihood(self, ddt, dd, kin_scaling=None):
         """
 
         :param ddt: time-delay distance
         :param dd: angular diameter distance to the deflector
-        :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
+        :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
          dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the
          anisotropy model used to derive the prediction and covariance matrix in the init of this class.
         :return: log likelihood given the single lens analysis
         """
-        if aniso_scaling is not None:
-            dd_ = dd * aniso_scaling[0]
+        if kin_scaling is not None:
+            dd_ = dd * kin_scaling[0]
         else:
             dd_ = dd
         lnlikelihood = (
diff --git a/hierarc/Likelihood/LensLikelihood/ddt_dd_kde_likelihood.py b/hierarc/Likelihood/LensLikelihood/ddt_dd_kde_likelihood.py
index 97ead86..9590cfb 100644
--- a/hierarc/Likelihood/LensLikelihood/ddt_dd_kde_likelihood.py
+++ b/hierarc/Likelihood/LensLikelihood/ddt_dd_kde_likelihood.py
@@ -52,18 +52,18 @@ def __init__(
         self._interpol = interpol
         self.num_data = 2
 
-    def log_likelihood(self, ddt, dd, aniso_scaling=None):
+    def log_likelihood(self, ddt, dd, kin_scaling=None):
         """
 
         :param ddt: time-delay distance
         :param dd: angular diameter distance to the deflector
-        :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
+        :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
          dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the
          anisotropy model used to derive the prediction and covariance matrix in the init of this class.
         :return: log likelihood given the single lens analysis
         """
-        if aniso_scaling is not None:
-            dd_ = dd * aniso_scaling[0]
+        if kin_scaling is not None:
+            dd_ = dd * kin_scaling[0]
         else:
             dd_ = dd
         if self._interpol is True:
diff --git a/hierarc/Likelihood/LensLikelihood/ddt_gauss_kin_likelihood.py b/hierarc/Likelihood/LensLikelihood/ddt_gauss_kin_likelihood.py
index f9cfbdd..d3a5822 100644
--- a/hierarc/Likelihood/LensLikelihood/ddt_gauss_kin_likelihood.py
+++ b/hierarc/Likelihood/LensLikelihood/ddt_gauss_kin_likelihood.py
@@ -55,7 +55,7 @@ def log_likelihood(
         self,
         ddt,
         dd,
-        aniso_scaling=None,
+        kin_scaling=None,
         sigma_v_sys_error=None,
         sigma_v_sys_offset=None,
     ):
@@ -63,7 +63,7 @@ def log_likelihood(
 
         :param ddt: time-delay distance
         :param dd: angular diameter distance to the deflector
-        :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
+        :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
          dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the
          anisotropy model used to derive the prediction and covariance matrix in the init of this class.
         :param sigma_v_sys_error: float (optional) added error on the velocity dispersion measurement in quadrature
@@ -76,19 +76,19 @@ def log_likelihood(
         ) + self._kinlikelihood.log_likelihood(
             ddt,
             dd,
-            aniso_scaling,
+            kin_scaling,
             sigma_v_sys_error=sigma_v_sys_error,
             sigma_v_sys_offset=sigma_v_sys_offset,
         )
         return lnlikelihood
 
-    def sigma_v_prediction(self, ddt, dd, aniso_scaling=1):
+    def sigma_v_prediction(self, ddt, dd, kin_scaling=1):
         """Model prediction mean velocity dispersion vector and model prediction
         covariance matrix.
 
         :param ddt: time-delay distance
         :param dd: angular diameter distance to the deflector
-        :param aniso_scaling: array of size of the velocity dispersion measurement or
+        :param kin_scaling: array of size of the velocity dispersion measurement or
             None, scaling of the predicted dimensionless quantity J (proportional to
             sigma_v^2) of the anisotropy model in the sampling relative to the
             anisotropy model used to derive the prediction and covariance matrix in the
@@ -96,7 +96,7 @@ def sigma_v_prediction(self, ddt, dd, aniso_scaling=1):
         :return: model prediction mean velocity dispersion vector and model prediction
             covariance matrix
         """
-        return self._kinlikelihood.sigma_v_prediction(ddt, dd, aniso_scaling)
+        return self._kinlikelihood.sigma_v_prediction(ddt, dd, kin_scaling)
 
     def sigma_v_measurement(self, sigma_v_sys_error=None, sigma_v_sys_offset=None):
         """
diff --git a/hierarc/Likelihood/LensLikelihood/ddt_hist_kin_likelihood.py b/hierarc/Likelihood/LensLikelihood/ddt_hist_kin_likelihood.py
index ec37141..760dadf 100644
--- a/hierarc/Likelihood/LensLikelihood/ddt_hist_kin_likelihood.py
+++ b/hierarc/Likelihood/LensLikelihood/ddt_hist_kin_likelihood.py
@@ -61,28 +61,28 @@ def __init__(
         )
         self.num_data = self._tdLikelihood.num_data + self._kinlikelihood.num_data
 
-    def log_likelihood(self, ddt, dd, aniso_scaling=None, sigma_v_sys_error=None):
+    def log_likelihood(self, ddt, dd, kin_scaling=None, sigma_v_sys_error=None):
         """
 
         :param ddt: time-delay distance
         :param dd: angular diameter distance to the deflector
-        :param aniso_scaling: numpy array of anisotropy scaling on prediction of Ds/Dds
+        :param kin_scaling: numpy array of anisotropy scaling on prediction of Ds/Dds
         :return: log likelihood given the single lens analysis
         """
         lnlikelihood = self._tdLikelihood.log_likelihood(
             ddt
         ) + self._kinlikelihood.log_likelihood(
-            ddt, dd, aniso_scaling, sigma_v_sys_error=sigma_v_sys_error
+            ddt, dd, kin_scaling, sigma_v_sys_error=sigma_v_sys_error
         )
         return lnlikelihood
 
-    def sigma_v_prediction(self, ddt, dd, aniso_scaling=1):
+    def sigma_v_prediction(self, ddt, dd, kin_scaling=1):
         """Model prediction mean velocity dispersion vector and model prediction
         covariance matrix.
 
         :param ddt: time-delay distance
         :param dd: angular diameter distance to the deflector
-        :param aniso_scaling: array of size of the velocity dispersion measurement or
+        :param kin_scaling: array of size of the velocity dispersion measurement or
             None, scaling of the predicted dimensionless quantity J (proportional to
             sigma_v^2) of the anisotropy model in the sampling relative to the
             anisotropy model used to derive the prediction and covariance matrix in the
@@ -90,7 +90,7 @@ def sigma_v_prediction(self, ddt, dd, aniso_scaling=1):
         :return: model prediction mean velocity dispersion vector and model prediction
             covariance matrix
         """
-        return self._kinlikelihood.sigma_v_prediction(ddt, dd, aniso_scaling)
+        return self._kinlikelihood.sigma_v_prediction(ddt, dd, kin_scaling)
 
     def sigma_v_measurement(self, sigma_v_sys_error=None, sigma_v_sys_offset=None):
         """
diff --git a/hierarc/Likelihood/LensLikelihood/ds_dds_gauss_likelihood.py b/hierarc/Likelihood/LensLikelihood/ds_dds_gauss_likelihood.py
index 03ebbaa..870022a 100644
--- a/hierarc/Likelihood/LensLikelihood/ds_dds_gauss_likelihood.py
+++ b/hierarc/Likelihood/LensLikelihood/ds_dds_gauss_likelihood.py
@@ -19,20 +19,20 @@ def __init__(self, z_lens, z_source, ds_dds_mean, ds_dds_sigma):
         self._ds_dds_sigma2 = ds_dds_sigma**2
         self.num_data = 1
 
-    def log_likelihood(self, ddt, dd, aniso_scaling=None):
+    def log_likelihood(self, ddt, dd, kin_scaling=None):
         """
         Note: kinematics + imaging data can constrain Ds/Dds. The input of Ddt, Dd is transformed here to match Ds/Dds
 
         :param ddt: time-delay distance
         :param dd: angular diameter distance to the deflector
-        :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
+        :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
          dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the
          anisotropy model used to derive the prediction and covariance matrix in the init of this class.
         :return: log likelihood given the single lens analysis
         """
         ds_dds = ddt / dd / (1 + self._z_lens)
-        if aniso_scaling is not None:
-            scaling = aniso_scaling[0]
+        if kin_scaling is not None:
+            scaling = kin_scaling[0]
         else:
             scaling = 1
         ds_dds_ = ds_dds / scaling
diff --git a/hierarc/Likelihood/LensLikelihood/kin_likelihood.py b/hierarc/Likelihood/LensLikelihood/kin_likelihood.py
index e2e13b3..acee56f 100644
--- a/hierarc/Likelihood/LensLikelihood/kin_likelihood.py
+++ b/hierarc/Likelihood/LensLikelihood/kin_likelihood.py
@@ -45,7 +45,7 @@ def log_likelihood(
         self,
         ddt,
         dd,
-        aniso_scaling=None,
+        kin_scaling=None,
         sigma_v_sys_error=None,
         sigma_v_sys_offset=None,
     ):
@@ -54,7 +54,7 @@ def log_likelihood(
 
         :param ddt: time-delay distance
         :param dd: angular diameter distance to the deflector
-        :param aniso_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
+        :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
          dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the
          anisotropy model used to derive the prediction and covariance matrix in the init of this class.
         :param sigma_v_sys_error: float (optional) added error on the velocity dispersion measurement in quadrature
@@ -63,10 +63,10 @@ def log_likelihood(
         :return: log likelihood given the single lens analysis
         """
         ds_dds = np.maximum(ddt / dd / (1 + self._z_lens), 0)
-        if aniso_scaling is None:
+        if kin_scaling is None:
             scaling_ifu = 1
         else:
-            scaling_ifu = aniso_scaling
+            scaling_ifu = kin_scaling
         sigma_v_predict = self.sigma_v_model(ds_dds, scaling_ifu)
         delta = self.sigma_v_measurement_mean(sigma_v_sys_offset) - sigma_v_predict
         cov_error = self.cov_error_measurement(
@@ -98,26 +98,24 @@ def sigma_v_measurement_mean(self, sigma_v_sys_offset=None):
         else:
             return self._sigma_v_measured * (1 + sigma_v_sys_offset)
 
-    def sigma_v_model(self, ds_dds, aniso_scaling=1):
+    def sigma_v_model(self, ds_dds, kin_scaling=1):
         """Model predicted velocity dispersion for the IFU's.
 
         :param ds_dds: Ds/Dds
-        :param aniso_scaling: scaling of the anisotropy affecting sigma_v^2
+        :param kin_scaling: scaling of the anisotropy affecting sigma_v^2
         :return: array of predicted velocity dispersions
         """
-        sigma_v_predict = (
-            np.sqrt(self._j_model * ds_dds * aniso_scaling) * const.c / 1000
-        )
+        sigma_v_predict = np.sqrt(self._j_model * ds_dds * kin_scaling) * const.c / 1000
         return sigma_v_predict
 
-    def cov_error_model(self, ds_dds, aniso_scaling=1):
+    def cov_error_model(self, ds_dds, kin_scaling=1):
         """
 
         :param ds_dds: Ds/Dds
-        :param aniso_scaling: scaling of the anisotropy affecting sigma_v^2
+        :param kin_scaling: scaling of the anisotropy affecting sigma_v^2
         :return: covariance matrix of the error in the predicted model (from mass model uncertainties)
         """
-        scaling_matix = np.outer(np.sqrt(aniso_scaling), np.sqrt(aniso_scaling))
+        scaling_matix = np.outer(np.sqrt(kin_scaling), np.sqrt(kin_scaling))
         return self._error_cov_j_sqrt * scaling_matix * ds_dds * (const.c / 1000) ** 2
 
     def cov_error_measurement(self, sigma_v_sys_error=None):
@@ -134,13 +132,13 @@ def cov_error_measurement(self, sigma_v_sys_error=None):
         else:
             return self._error_cov_measurement
 
-    def sigma_v_prediction(self, ddt, dd, aniso_scaling=1):
+    def sigma_v_prediction(self, ddt, dd, kin_scaling=1):
         """Model prediction mean velocity dispersion vector and model prediction
         covariance matrix.
 
         :param ddt: time-delay distance
         :param dd: angular diameter distance to the deflector
-        :param aniso_scaling: array of size of the velocity dispersion measurement or
+        :param kin_scaling: array of size of the velocity dispersion measurement or
             None, scaling of the predicted dimensionless quantity J (proportional to
             sigma_v^2) of the anisotropy model in the sampling relative to the
             anisotropy model used to derive the prediction and covariance matrix in the
@@ -149,8 +147,8 @@ def sigma_v_prediction(self, ddt, dd, aniso_scaling=1):
             covariance matrix
         """
         ds_dds = np.maximum(ddt / dd / (1 + self._z_lens), 0)
-        sigma_v_predict = self.sigma_v_model(ds_dds, aniso_scaling)
-        cov_error_predict = self.cov_error_model(ds_dds, aniso_scaling)
+        sigma_v_predict = self.sigma_v_model(ds_dds, kin_scaling)
+        cov_error_predict = self.cov_error_model(ds_dds, kin_scaling)
         return sigma_v_predict, cov_error_predict
 
     def sigma_v_measurement(self, sigma_v_sys_error=None, sigma_v_sys_offset=None):
diff --git a/hierarc/Likelihood/anisotropy_scaling.py b/hierarc/Likelihood/anisotropy_scaling.py
deleted file mode 100644
index 0afec08..0000000
--- a/hierarc/Likelihood/anisotropy_scaling.py
+++ /dev/null
@@ -1,156 +0,0 @@
-__author__ = "sibirrer"
-from scipy.interpolate import interp1d, interp2d
-import numpy as np
-
-
-class AnisotropyScalingSingleAperture(object):
-    """Class to manage anisotropy scaling for single slit observation."""
-
-    def __init__(self, ani_param_array, ani_scaling_array):
-        """
-
-        :param ani_param_array: array of anisotropy parameter value
-        :param ani_scaling_array: array with the scaling of J() for single slit
-        """
-        self._evalute_ani = False
-        if ani_param_array is not None and ani_scaling_array is not None:
-            if isinstance(ani_param_array, list):
-                self._dim_scaling = len(ani_param_array)
-            else:
-                self._dim_scaling = 1
-            if self._dim_scaling == 1:
-                self._f_ani = interp1d(
-                    ani_param_array, ani_scaling_array, kind="linear"
-                )
-            elif self._dim_scaling == 2:
-                self._f_ani = interp2d(
-                    ani_param_array[0], ani_param_array[1], ani_scaling_array.T
-                )
-            else:
-                raise ValueError(
-                    "anisotropy scaling with dimension %s not supported."
-                    % self._dim_scaling
-                )
-            self._evalute_ani = True
-
-    def ani_scaling(self, aniso_param_array):
-        """
-
-        :param aniso_param_array: anisotropy parameter array
-        :return: scaling J(a_ani) for single slit
-        """
-        if self._evalute_ani is not True or aniso_param_array is None:
-            return 1
-        if self._dim_scaling == 1:
-            return self._f_ani(aniso_param_array[0])
-        elif self._dim_scaling == 2:
-            return self._f_ani(aniso_param_array[0], aniso_param_array[1])[0]
-
-
-class AnisotropyScalingIFU(object):
-    """Class to manage anisotropy scalings for IFU data."""
-
-    def __init__(
-        self, anisotropy_model="NONE", ani_param_array=None, ani_scaling_array_list=None
-    ):
-        """
-
-        :param anisotropy_model: string, either 'NONE', 'OM' or 'GOM'
-        :param ani_param_array: array of anisotropy parameter value (1d for 'OM' model, 2d for 'GOM' model)
-        :param ani_scaling_array_list: list of array with the scalings of J() for each IFU
-        """
-        self._anisotropy_model = anisotropy_model
-        self._evalute_ani = False
-        if (
-            ani_param_array is not None
-            and ani_scaling_array_list is not None
-            and self._anisotropy_model != "NONE"
-        ):
-            self._evalute_ani = True
-            self._anisotropy_scaling_list = []
-            self._f_ani_list = []
-            for ani_scaling_array in ani_scaling_array_list:
-                self._anisotropy_scaling_list.append(
-                    AnisotropyScalingSingleAperture(
-                        ani_param_array=ani_param_array,
-                        ani_scaling_array=ani_scaling_array,
-                    )
-                )
-
-            if isinstance(ani_param_array, list):
-                self._dim_scaling = len(ani_param_array)
-            else:
-                self._dim_scaling = 1
-            if self._dim_scaling == 1 and anisotropy_model in ["OM", "const"]:
-                self._ani_param_min = np.min(ani_param_array)
-                self._ani_param_max = np.max(ani_param_array)
-            elif self._dim_scaling == 2 and anisotropy_model == "GOM":
-                self._ani_param_min = [min(ani_param_array[0]), min(ani_param_array[1])]
-                self._ani_param_max = [max(ani_param_array[0]), max(ani_param_array[1])]
-            else:
-                raise ValueError(
-                    "anisotropy scaling with dimension %s does not match anisotropy model %s"
-                    % (self._dim_scaling, self._anisotropy_model)
-                )
-
-    def ani_scaling(self, aniso_param_array):
-        """
-
-        :param aniso_param_array: anisotropy parameter array
-        :return: scaling J(a_ani) for the IFU's
-        """
-        if self._evalute_ani is not True or aniso_param_array is None:
-            return [1]
-        scaling_list = []
-        for scaling_class in self._anisotropy_scaling_list:
-            scaling = scaling_class.ani_scaling(aniso_param_array)
-            scaling_list.append(scaling)
-        return np.array(scaling_list)
-
-    def draw_anisotropy(
-        self, a_ani=None, a_ani_sigma=0, beta_inf=None, beta_inf_sigma=0
-    ):
-        """Draw Gaussian distribution and re-sample if outside bounds.
-
-        :param a_ani: mean of the distribution
-        :param a_ani_sigma: std of the distribution
-        :param beta_inf: anisotropy at infinity (relevant for GOM model)
-        :param beta_inf_sigma: std of beta_inf distribution
-        :return: random draw from the distribution
-        """
-        if self._anisotropy_model in ["OM", "const"]:
-            if a_ani < self._ani_param_min or a_ani > self._ani_param_max:
-                raise ValueError(
-                    "anisotropy parameter is out of bounds of the interpolated range!"
-                )
-            # we draw a linear gaussian for 'const' anisotropy and a scaled proportional one for 'OM
-            if self._anisotropy_model == "OM":
-                a_ani_draw = np.random.normal(a_ani, a_ani_sigma * a_ani)
-            else:
-                a_ani_draw = np.random.normal(a_ani, a_ani_sigma)
-            if a_ani_draw < self._ani_param_min or a_ani_draw > self._ani_param_max:
-                return self.draw_anisotropy(a_ani, a_ani_sigma)
-            return np.array([a_ani_draw])
-        elif self._anisotropy_model in ["GOM"]:
-            if (
-                a_ani < self._ani_param_min[0]
-                or a_ani > self._ani_param_max[0]
-                or beta_inf < self._ani_param_min[1]
-                or beta_inf > self._ani_param_max[1]
-            ):
-                raise ValueError(
-                    "anisotropy parameter is out of bounds of the interpolated range!"
-                )
-            a_ani_draw = np.random.normal(a_ani, a_ani_sigma * a_ani)
-            beta_inf_draw = np.random.normal(beta_inf, beta_inf_sigma)
-            if (
-                a_ani_draw < self._ani_param_min[0]
-                or a_ani_draw > self._ani_param_max[0]
-                or beta_inf_draw < self._ani_param_min[1]
-                or beta_inf_draw > self._ani_param_max[1]
-            ):
-                return self.draw_anisotropy(
-                    a_ani, a_ani_sigma, beta_inf, beta_inf_sigma
-                )
-            return np.array([a_ani_draw, beta_inf_draw])
-        return None
diff --git a/hierarc/Likelihood/cosmo_likelihood.py b/hierarc/Likelihood/cosmo_likelihood.py
index f1f9826..5d78870 100644
--- a/hierarc/Likelihood/cosmo_likelihood.py
+++ b/hierarc/Likelihood/cosmo_likelihood.py
@@ -23,10 +23,16 @@ def __init__(
         lambda_mst_sampling=False,
         lambda_mst_distribution="delta",
         anisotropy_sampling=False,
+        gamma_in_sampling=False,
+        gamma_in_distribution="NONE",
+        m2l_sampling=False,
+        m2l_distribution="NONE",
         kappa_ext_sampling=False,
         kappa_ext_distribution="NONE",
         alpha_lambda_sampling=False,
         beta_lambda_sampling=False,
+        alpha_gamma_in_sampling=False,
+        alpha_m2l_sampling=False,
         lambda_ifu_sampling=False,
         lambda_ifu_distribution="NONE",
         sigma_v_systematics=False,
@@ -73,6 +79,10 @@ def __init__(
         kinematic prediction
         :param anisotropy_model: string, specifies the stellar anisotropy model
         :param anisotropy_distribution: string, distribution of the anisotropy parameters
+        :param gamma_in_sampling: bool, if True samples gNFW inner slope parameter
+        :param gamma_in_distribution: string, distribution function of the gamma_in parameter
+        :param m2l_sampling: bool, if True samples a global mass-to-light ratio parameter
+        :param m2l_distribution: string, distribution function of the m2l parameter
         :param sigma_v_systematics: bool, if True samples paramaters relative to systematics in the velocity dispersion
          measurement
         :param sne_apparent_m_sampling: boolean, if True, samples/queries SNe unlensed magnitude distribution
@@ -107,6 +117,12 @@ def __init__(
             lambda_ifu_distribution=lambda_ifu_distribution,
             alpha_lambda_sampling=alpha_lambda_sampling,
             beta_lambda_sampling=beta_lambda_sampling,
+            gamma_in_sampling=gamma_in_sampling,
+            gamma_in_distribution=gamma_in_distribution,
+            m2l_sampling=m2l_sampling,
+            m2l_distribution=m2l_distribution,
+            alpha_gamma_in_sampling=alpha_gamma_in_sampling,
+            alpha_m2l_sampling=alpha_m2l_sampling,
             sne_apparent_m_sampling=sne_apparent_m_sampling,
             sne_distribution=sne_distribution,
             z_apparent_m_anchor=z_apparent_m_anchor,
diff --git a/hierarc/Likelihood/hierarchy_likelihood.py b/hierarc/Likelihood/hierarchy_likelihood.py
index de3a934..533f2c3 100644
--- a/hierarc/Likelihood/hierarchy_likelihood.py
+++ b/hierarc/Likelihood/hierarchy_likelihood.py
@@ -1,12 +1,12 @@
 from hierarc.Likelihood.transformed_cosmography import TransformedCosmography
 from hierarc.Likelihood.LensLikelihood.base_lens_likelihood import LensLikelihoodBase
-from hierarc.Likelihood.anisotropy_scaling import AnisotropyScalingIFU
+from hierarc.Likelihood.parameter_scaling import ParameterScalingIFU
 from hierarc.Util.distribution_util import PDFSampling
 import numpy as np
 import copy
 
 
-class LensLikelihood(TransformedCosmography, LensLikelihoodBase, AnisotropyScalingIFU):
+class LensLikelihood(TransformedCosmography, LensLikelihoodBase, ParameterScalingIFU):
     """Master class containing the likelihood definitions of different analysis for s
     single lens."""
 
@@ -18,8 +18,11 @@ def __init__(
         likelihood_type="TDKin",
         anisotropy_model="NONE",
         ani_param_array=None,
-        ani_scaling_array_list=None,
         ani_scaling_array=None,
+        ani_scaling_array_list=None,
+        gamma_in_array=None,
+        m2l_array=None,
+        param_scaling_grid_list=None,
         num_distribution_draws=50,
         kappa_ext_bias=False,
         kappa_pdf=None,
@@ -42,7 +45,9 @@ def __init__(
          to default prediction. The scaling corresponds to the ani_param_array parameter spacing
          (to generate an interpolation function). A value =1 in ani_scaling_array results in the value stored in the
          provided J() predictions.
-        :param ani_scaling_array_list: list of array with the scalings of J() for each IFU
+        :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.
         :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.
@@ -65,12 +70,32 @@ def __init__(
         TransformedCosmography.__init__(self, z_lens=z_lens, z_source=z_source)
         if ani_scaling_array_list is None and ani_scaling_array is not None:
             ani_scaling_array_list = [ani_scaling_array]
-        AnisotropyScalingIFU.__init__(
-            self,
-            anisotropy_model=anisotropy_model,
-            ani_param_array=ani_param_array,
-            ani_scaling_array_list=ani_scaling_array_list,
-        )
+
+        # AnisotropyScalingIFU.__init__(
+        #     self,
+        #     anisotropy_model=anisotropy_model,
+        #     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 isinstance(ani_param_array, list):
+                param_arrays = ani_param_array + [gamma_in_array, m2l_array]
+            else:
+                param_arrays = [ani_param_array, gamma_in_array, m2l_array]
+            ParameterScalingIFU.__init__(
+                self,
+                anisotropy_model=anisotropy_model,
+                param_arrays=param_arrays,
+                scaling_grid_list=param_scaling_grid_list,
+            )
+        else:
+            ParameterScalingIFU.__init__(
+                self,
+                anisotropy_model=anisotropy_model,
+                param_arrays=ani_param_array,
+                scaling_grid_list=ani_scaling_array_list,
+            )
+
         LensLikelihoodBase.__init__(
             self,
             z_lens=z_lens,
@@ -93,6 +118,8 @@ def __init__(
             self._draw_kappa = False
         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
 
     def lens_log_likelihood(
         self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None
@@ -116,7 +143,7 @@ def lens_log_likelihood(
         delta_lum_dist = self.luminosity_distance_modulus(cosmo, z_apparent_m_anchor)
         # here we effectively change the posteriors of the lens, but rather than changing the instance of the KDE we
         # displace the predicted angular diameter distances in the opposite direction
-        return self.hyper_param_likelihood(
+        a = self.hyper_param_likelihood(
             ddt,
             dd,
             delta_lum_dist,
@@ -125,6 +152,7 @@ def lens_log_likelihood(
             kwargs_source=kwargs_source,
             cosmo=cosmo,
         )
+        return a
 
     def hyper_param_likelihood(
         self,
@@ -221,18 +249,88 @@ def log_likelihood_single(
             kappa_ext=kappa_ext,
             mag_source=mag_source,
         )
-        aniso_param_array = self.draw_anisotropy(**kwargs_kin)
-        aniso_scaling = self.ani_scaling(aniso_param_array)
+        scaling_param_array = self.draw_scaling_params(
+            kwargs_lens=kwargs_lens, **kwargs_kin
+        )
+        kin_scaling = self.param_scaling(scaling_param_array)
 
         lnlikelihood = self.log_likelihood(
             ddt_,
             dd_,
-            aniso_scaling=aniso_scaling,
+            kin_scaling=kin_scaling,
             sigma_v_sys_error=sigma_v_sys_error,
             mu_intrinsic=mag_source_,
         )
         return np.nan_to_num(lnlikelihood)
 
+    def draw_scaling_params(self, kwargs_lens=None, **kwargs_kin):
+        """Draws a realization of the anisotropy parameter scaling from the distribution
+        function.
+
+        :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]])
+        else:
+            return ani_param
+
+    def draw_lens_scaling_params(
+        self,
+        lambda_mst=1,
+        lambda_mst_sigma=0,
+        kappa_ext=0,
+        kappa_ext_sigma=0,
+        gamma_ppn=1,
+        lambda_ifu=1,
+        lambda_ifu_sigma=0,
+        alpha_lambda=0,
+        beta_lambda=0,
+        gamma_in=1,
+        gamma_in_sigma=0,
+        alpha_gamma_in=0,
+        m2l=1,
+        m2l_sigma=0,
+        alpha_m2l=0,
+    ):
+        """Draws a realization of the anisotropy parameter scaling from the
+        distribution.
+
+        :param lambda_mst: MST transform
+        :param lambda_mst_sigma: spread in the distribution
+        :param kappa_ext: external convergence mean in distribution
+        :param kappa_ext_sigma: spread in the distribution
+        :param gamma_ppn: Post-Newtonian parameter
+        :param lambda_ifu: secondary lambda_mst parameter for subset of lenses specified
+            for
+        :param lambda_ifu_sigma: secondary lambda_mst_sigma parameter for subset of
+            lenses specified for
+        :param alpha_lambda: float, linear slope of the lambda_int scaling relation with
+            lens quantity self._lambda_scaling_property
+        :param beta_lambda: float, a second linear slope of the lambda_int scaling
+            relation with lens quantity self._lambda_scaling_property_beta
+        :param gamma_in: inner slope of the NFW profile
+        :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
+        :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(
+                gamma_in + alpha_gamma_in * self._lambda_scaling_property,
+                gamma_in_sigma,
+                m2l + alpha_m2l * self._lambda_scaling_property,
+                m2l_sigma,
+            )
+            return gamma_in_draw, m2l_draw
+        else:
+            return None, None
+
     def angular_diameter_distances(self, cosmo):
         """Time-delay distance Ddt, angular diameter distance to the lens (dd)
 
@@ -312,6 +410,12 @@ def draw_lens(
         lambda_ifu_sigma=0,
         alpha_lambda=0,
         beta_lambda=0,
+        gamma_in=1,
+        gamma_in_sigma=0,
+        alpha_gamma_in=0,
+        m2l=1,
+        m2l_sigma=0,
+        alpha_m2l=0,
     ):
         """Draws a realization of a specific model from the hyper-parameter
         distribution.
@@ -329,6 +433,14 @@ def draw_lens(
             lens quantity self._lambda_scaling_property
         :param beta_lambda: float, a second linear slope of the lambda_int scaling
             relation with lens quantity self._lambda_scaling_property_beta
+        :param gamma_in: inner slope of the NFW profile
+        :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
         :return: draw from the distributions
         """
         if self._mst_ifu is True:
@@ -403,10 +515,12 @@ def sigma_v_measured_vs_predict(self, cosmo, kwargs_lens=None, kwargs_kin=None):
             ddt_, dd_, _ = self.displace_prediction(
                 ddt, dd, gamma_ppn=gamma_ppn, lambda_mst=lambda_mst, kappa_ext=kappa_ext
             )
-            aniso_param_array = self.draw_anisotropy(**kwargs_kin_copy)
-            aniso_scaling = self.ani_scaling(aniso_param_array)
+            scaling_param_array = self.draw_scaling_params(
+                kwargs_lens=kwargs_lens, **kwargs_kin_copy
+            )
+            kin_scaling = self.param_scaling(scaling_param_array)
             sigma_v_predict_i, cov_error_predict_i = self.sigma_v_prediction(
-                ddt_, dd_, aniso_scaling=aniso_scaling
+                ddt_, dd_, kin_scaling=kin_scaling
             )
             sigma_v_predict_mean += sigma_v_predict_i
             cov_error_predict += cov_error_predict_i
diff --git a/hierarc/Likelihood/parameter_scaling.py b/hierarc/Likelihood/parameter_scaling.py
new file mode 100644
index 0000000..3a00f82
--- /dev/null
+++ b/hierarc/Likelihood/parameter_scaling.py
@@ -0,0 +1,208 @@
+__author__ = "sibirrer", "ajshajib"
+
+from scipy.interpolate import interp1d
+from scipy.interpolate import interp2d
+from scipy.interpolate import RegularGridInterpolator
+import numpy as np
+
+
+class ParameterScalingSingleAperture(object):
+    """Class to manage anisotropy scaling for single slit observation."""
+
+    def __init__(self, param_arrays, scaling_grid):
+        """
+
+        :param param_arrays: list of arrays of interpolated parameter values
+        :param scaling_grid: array with the scaling of J() for single slit
+        """
+        self._evalute_scaling = False
+        # check if param arrays is 1d list or 2d list
+        if param_arrays is not None and scaling_grid is not None:
+            if isinstance(param_arrays, list):
+                self._dim_scaling = len(param_arrays)
+            else:
+                self._dim_scaling = 1
+
+            if self._dim_scaling == 1:
+                self._f_ani = interp1d(param_arrays, scaling_grid, kind="linear")
+            elif self._dim_scaling == 2:
+                self._f_ani = interp2d(param_arrays[0], param_arrays[1], scaling_grid.T)
+            else:
+                self._f_ani = RegularGridInterpolator(
+                    tuple(param_arrays),
+                    scaling_grid,
+                )
+            self._evalute_scaling = True
+
+    def param_scaling(self, param_array):
+        """
+
+        :param param_array: anisotropy parameter array
+        :return: scaling J(a_ani) for single slit
+        """
+        if self._evalute_scaling is not True or param_array is None:
+            return 1
+        if self._dim_scaling == 1:
+            return self._f_ani(param_array[0])
+        elif self._dim_scaling == 2:
+            return self._f_ani(param_array[0], param_array[1])[0]
+        else:
+            return self._f_ani(param_array)[0]
+
+
+class ParameterScalingIFU(object):
+    """Class to manage model parameter and anisotropy scalings for IFU data."""
+
+    def __init__(
+        self, anisotropy_model="NONE", param_arrays=None, scaling_grid_list=None
+    ):
+        """
+
+        :param anisotropy_model: string, either 'NONE', 'OM' or 'GOM'
+        :param param_arrays: array of parameter values
+        :param scaling_grid_list: list of array with the scalings of J() for each IFU
+        """
+        self._anisotropy_model = anisotropy_model
+        self._evalute_ani = False
+        if (
+            param_arrays is not None
+            and scaling_grid_list is not None
+            and self._anisotropy_model != "NONE"
+        ):
+            self._evalute_ani = True
+            self._anisotropy_scaling_list = []
+            self._f_ani_list = []
+            for scaling_grid in scaling_grid_list:
+                self._anisotropy_scaling_list.append(
+                    ParameterScalingSingleAperture(param_arrays, scaling_grid)
+                )
+
+            if isinstance(param_arrays, list):
+                self._dim_scaling = len(param_arrays)
+            else:
+                self._dim_scaling = 1
+
+            if anisotropy_model in ["OM", "const"]:
+                if self._dim_scaling == 1:
+                    self._ani_param_min = np.min(param_arrays)
+                    self._ani_param_max = np.max(param_arrays)
+                else:
+                    self._ani_param_min = np.min(param_arrays[0])
+                    self._ani_param_max = np.max(param_arrays[0])
+
+                if self._dim_scaling > 1:
+                    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])
+
+            elif anisotropy_model == "GOM":
+                self._ani_param_min = [min(param_arrays[0]), min(param_arrays[1])]
+                self._ani_param_max = [max(param_arrays[0]), max(param_arrays[1])]
+
+                if self._dim_scaling > 2:
+                    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])
+            else:
+                raise ValueError(
+                    f"Anisotropy model {anisotropy_model} is not recognized!"
+                )
+
+    def param_scaling(self, param_array):
+        """
+
+        :param param_array: parameter array for scaling
+        :return: scaling J(a_ani) for the IFU's
+        """
+        if self._evalute_ani is not True or param_array is None:
+            return [1]
+        scaling_list = []
+        for scaling_class in self._anisotropy_scaling_list:
+            scaling = scaling_class.param_scaling(param_array)
+            scaling_list.append(scaling)
+        return np.array(scaling_list)
+
+    def draw_anisotropy(
+        self, a_ani=None, a_ani_sigma=0, beta_inf=None, beta_inf_sigma=0
+    ):
+        """Draw Gaussian distribution and re-sample if outside bounds.
+
+        :param a_ani: mean of the distribution
+        :param a_ani_sigma: std of the distribution
+        :param beta_inf: anisotropy at infinity (relevant for GOM model)
+        :param beta_inf_sigma: std of beta_inf distribution
+        :return: random draw from the distribution
+        """
+        if self._anisotropy_model in ["OM", "const"]:
+            if a_ani < self._ani_param_min or a_ani > self._ani_param_max:
+                raise ValueError(
+                    "anisotropy parameter is out of bounds of the interpolated range!"
+                )
+            # we draw a linear gaussian for 'const' anisotropy and a scaled proportional one for 'OM
+            if self._anisotropy_model == "OM":
+                a_ani_draw = np.random.normal(a_ani, a_ani_sigma * a_ani)
+            else:
+                a_ani_draw = np.random.normal(a_ani, a_ani_sigma)
+            if a_ani_draw < self._ani_param_min or a_ani_draw > self._ani_param_max:
+                return self.draw_anisotropy(a_ani, a_ani_sigma)
+            return np.array([a_ani_draw])
+        elif self._anisotropy_model in ["GOM"]:
+            if (
+                a_ani < self._ani_param_min[0]
+                or a_ani > self._ani_param_max[0]
+                or beta_inf < self._ani_param_min[1]
+                or beta_inf > self._ani_param_max[1]
+            ):
+                raise ValueError(
+                    "anisotropy parameter is out of bounds of the interpolated range!"
+                )
+            a_ani_draw = np.random.normal(a_ani, a_ani_sigma * a_ani)
+            beta_inf_draw = np.random.normal(beta_inf, beta_inf_sigma)
+            if (
+                a_ani_draw < self._ani_param_min[0]
+                or a_ani_draw > self._ani_param_max[0]
+                or beta_inf_draw < self._ani_param_min[1]
+                or beta_inf_draw > self._ani_param_max[1]
+            ):
+                return self.draw_anisotropy(
+                    a_ani, a_ani_sigma, beta_inf, beta_inf_sigma
+                )
+            return np.array([a_ani_draw, beta_inf_draw])
+        return None
+
+    def draw_lens_parameters(
+        self, gamma_in=None, gamma_in_sigma=0, m2l=None, 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
+        :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 (
+            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)
+
+        return gamma_in_draw, m2l_draw
diff --git a/hierarc/Sampling/ParamManager/lens_param.py b/hierarc/Sampling/ParamManager/lens_param.py
index aefb3c3..60b94de 100644
--- a/hierarc/Sampling/ParamManager/lens_param.py
+++ b/hierarc/Sampling/ParamManager/lens_param.py
@@ -8,12 +8,18 @@ def __init__(
         self,
         lambda_mst_sampling=False,
         lambda_mst_distribution="NONE",
-        kappa_ext_sampling=False,
-        kappa_ext_distribution="NONE",
         lambda_ifu_sampling=False,
         lambda_ifu_distribution="NONE",
+        gamma_in_sampling=False,
+        gamma_in_distribution="NONE",
+        m2l_sampling=False,
+        m2l_distribution="NONE",
+        kappa_ext_sampling=False,
+        kappa_ext_distribution="NONE",
         alpha_lambda_sampling=False,
         beta_lambda_sampling=False,
+        alpha_gamma_in_sampling=False,
+        alpha_m2l_sampling=False,
         kwargs_fixed=None,
         log_scatter=False,
     ):
@@ -21,15 +27,26 @@ def __init__(
 
         :param lambda_mst_sampling: bool, if True adds a global mass-sheet transform parameter in the sampling
         :param lambda_mst_distribution: string, distribution function of the MST transform
-        :param kappa_ext_sampling: bool, if True samples a global external convergence parameter
-        :param kappa_ext_distribution: string, distribution function of the kappa_ext parameter
-        :param lambda_ifu_sampling: bool, if True samples a separate lambda_mst for a second (e.g. IFU) data set
-        independently
+        :param lambda_ifu_sampling: bool, if True samples a separate lambda_mst for a
+            second (e.g. IFU) data set independently
         :param lambda_ifu_distribution: string, distribution function of the lambda_ifu parameter
+        :param gamma_in_sampling: bool, if True samples the inner slope of the GNFW profile
+        :param gamma_in_distribution: string, distribution function of the inner
+            slope of the GNFW profile
+        :param m2l_sampling: bool, if True samples the mass to light ratio of
+            the stars
+        :param m2l_distribution: string, distribution function of the mass to
+            light ratio of the lens
+        :param kappa_ext_sampling: bool, if True samples a global external convergence
+            parameter
+        :param kappa_ext_distribution: string, distribution function of the kappa_ext
+            parameter
         :param alpha_lambda_sampling: bool, if True samples a parameter alpha_lambda, which scales lambda_mst linearly
-         according to a predefined quantity of the lens
+            according to a predefined quantity of the lens
         :param beta_lambda_sampling: bool, if True samples a parameter beta_lambda, which scales lambda_mst linearly
-         according to a predefined quantity of the lens
+            according to a predefined quantity of the lens
+        :param alpha_gamma_in_sampling: bool, if True samples a parameter alpha_gamma_in, which scales gamma_in linearly
+        :param alpha_m2l_sampling: bool, if True samples a parameter alpha_m2l, which scales m2l linearly
         :param log_scatter: boolean, if True, samples the Gaussian scatter amplitude in log space (and thus flat prior in log)
         :param kwargs_fixed: keyword arguments that are held fixed through the sampling
         """
@@ -37,10 +54,17 @@ def __init__(
         self._lambda_mst_distribution = lambda_mst_distribution
         self._lambda_ifu_sampling = lambda_ifu_sampling
         self._lambda_ifu_distribution = lambda_ifu_distribution
+        self._gamma_in_sampling = gamma_in_sampling
+        self._gamma_in_distribution = gamma_in_distribution
+        self._m2l_sampling = m2l_sampling
+        self._m2l_distribution = m2l_distribution
         self._kappa_ext_sampling = kappa_ext_sampling
         self._kappa_ext_distribution = kappa_ext_distribution
         self._alpha_lambda_sampling = alpha_lambda_sampling
         self._beta_lambda_sampling = beta_lambda_sampling
+        self._alpha_gamma_in_sampling = alpha_gamma_in_sampling
+        self._alpha_m2l_sampling = alpha_m2l_sampling
+
         self._log_scatter = log_scatter
         if kwargs_fixed is None:
             kwargs_fixed = {}
@@ -83,6 +107,36 @@ def param_list(self, latex_style=False):
                             list.append(r"$\sigma(\lambda_{\rm ifu})$")
                     else:
                         list.append("lambda_ifu_sigma")
+        if self._gamma_in_sampling:
+            if "gamma_in" not in self._kwargs_fixed:
+                if latex_style is True:
+                    list.append(r"$\gamma_{\rm in}$")
+                else:
+                    list.append("gamma_in")
+            if self._gamma_in_distribution == "GAUSSIAN":
+                if "gamma_in_sigma" not in self._kwargs_fixed:
+                    if latex_style is True:
+                        if self._log_scatter is True:
+                            list.append(r"$\log_{10}\sigma(\gamma_{\rm in})$")
+                        else:
+                            list.append(r"$\sigma(\gamma_{\rm in})$")
+                    else:
+                        list.append("gamma_in_sigma")
+        if self._m2l_sampling is True:
+            if "m2l" not in self._kwargs_fixed:
+                if latex_style is True:
+                    list.append(r"$\Upsilon_{\rm stars}$")
+                else:
+                    list.append("m2l")
+            if self._m2l_distribution == "GAUSSIAN":
+                if "m2l_sigma" not in self._kwargs_fixed:
+                    if latex_style is True:
+                        if self._log_scatter is True:
+                            list.append(r"$\log_{10}\sigma(\Upsilon_{\rm stars})$")
+                        else:
+                            list.append(r"$\sigma(\Upsilon_{\rm stars})$")
+                    else:
+                        list.append("m2l_sigma")
         if self._kappa_ext_sampling is True:
             if "kappa_ext" not in self._kwargs_fixed:
                 if latex_style is True:
@@ -107,6 +161,18 @@ def param_list(self, latex_style=False):
                     list.append(r"$\beta_{\lambda}$")
                 else:
                     list.append("beta_lambda")
+        if self._alpha_gamma_in_sampling is True:
+            if "alpha_gamma_in" not in self._kwargs_fixed:
+                if latex_style is True:
+                    list.append(r"$\alpha_{\gamma_{\rm in}}$")
+                else:
+                    list.append("alpha_gamma_in")
+        if self._alpha_m2l_sampling is True:
+            if "alpha_m2l" not in self._kwargs_fixed:
+                if latex_style is True:
+                    list.append(r"$\alpha_{\Upsilon_{\rm stars}}$")
+                else:
+                    list.append("alpha_m2l")
         return list
 
     def args2kwargs(self, args, i=0):
@@ -146,6 +212,36 @@ def args2kwargs(self, args, i=0):
                     else:
                         kwargs["lambda_ifu_sigma"] = args[i]
                     i += 1
+        if self._gamma_in_sampling is True:
+            if "gamma_in" in self._kwargs_fixed:
+                kwargs["gamma_in"] = self._kwargs_fixed["gamma_in"]
+            else:
+                kwargs["gamma_in"] = args[i]
+                i += 1
+            if self._gamma_in_distribution == "GAUSSIAN":
+                if "gamma_in_sigma" in self._kwargs_fixed:
+                    kwargs["gamma_in_sigma"] = self._kwargs_fixed["gamma_in_sigma"]
+                else:
+                    if self._log_scatter is True:
+                        kwargs["gamma_in_sigma"] = 10 ** (args[i])
+                    else:
+                        kwargs["gamma_in_sigma"] = args[i]
+                    i += 1
+        if self._m2l_sampling is True:
+            if "m2l" in self._kwargs_fixed:
+                kwargs["m2l"] = self._kwargs_fixed["m2l"]
+            else:
+                kwargs["m2l"] = args[i]
+                i += 1
+            if self._m2l_distribution == "GAUSSIAN":
+                if "m2l_sigma" in self._kwargs_fixed:
+                    kwargs["m2l_sigma"] = self._kwargs_fixed["m2l_sigma"]
+                else:
+                    if self._log_scatter is True:
+                        kwargs["m2l_sigma"] = 10 ** (args[i])
+                    else:
+                        kwargs["m2l_sigma"] = args[i]
+                    i += 1
         if self._kappa_ext_sampling is True:
             if "kappa_ext" in self._kwargs_fixed:
                 kwargs["kappa_ext"] = self._kwargs_fixed["kappa_ext"]
@@ -170,6 +266,19 @@ def args2kwargs(self, args, i=0):
             else:
                 kwargs["beta_lambda"] = args[i]
                 i += 1
+        if self._alpha_gamma_in_sampling is True:
+            if "alpha_gamma_in" in self._kwargs_fixed:
+                kwargs["alpha_gamma_in"] = self._kwargs_fixed["alpha_gamma_in"]
+            else:
+                kwargs["alpha_gamma_in"] = args[i]
+                i += 1
+        if self._alpha_m2l_sampling is True:
+            if "alpha_m2l" in self._kwargs_fixed:
+                kwargs["alpha_m2l"] = self._kwargs_fixed["alpha_m2l"]
+            else:
+                kwargs["alpha_m2l"] = args[i]
+                i += 1
+
         return kwargs, i
 
     def kwargs2args(self, kwargs):
@@ -197,6 +306,24 @@ def kwargs2args(self, kwargs):
                         args.append(np.log10(kwargs["lambda_ifu_sigma"]))
                     else:
                         args.append(kwargs["lambda_ifu_sigma"])
+        if self._gamma_in_sampling is True:
+            if "gamma_in" not in self._kwargs_fixed:
+                args.append(kwargs["gamma_in"])
+            if self._gamma_in_distribution == "GAUSSIAN":
+                if "gamma_in_sigma" not in self._kwargs_fixed:
+                    if self._log_scatter is True:
+                        args.append(np.log10(kwargs["gamma_in_sigma"]))
+                    else:
+                        args.append(kwargs["gamma_in_sigma"])
+        if self._m2l_sampling is True:
+            if "m2l" not in self._kwargs_fixed:
+                args.append(kwargs["m2l"])
+            if self._m2l_distribution == "GAUSSIAN":
+                if "m2l_sigma" not in self._kwargs_fixed:
+                    if self._log_scatter is True:
+                        args.append(np.log10(kwargs["m2l_sigma"]))
+                    else:
+                        args.append(kwargs["m2l_sigma"])
         if self._kappa_ext_sampling is True:
             if "kappa_ext" not in self._kwargs_fixed:
                 args.append(kwargs["kappa_ext"])
@@ -209,4 +336,10 @@ def kwargs2args(self, kwargs):
         if self._beta_lambda_sampling is True:
             if "beta_lambda" not in self._kwargs_fixed:
                 args.append(kwargs["beta_lambda"])
+        if self._alpha_gamma_in_sampling is True:
+            if "alpha_gamma_in" not in self._kwargs_fixed:
+                args.append(kwargs["alpha_gamma_in"])
+        if self._alpha_m2l_sampling is True:
+            if "alpha_m2l" not in self._kwargs_fixed:
+                args.append(kwargs["alpha_m2l"])
         return args
diff --git a/hierarc/Sampling/ParamManager/param_manager.py b/hierarc/Sampling/ParamManager/param_manager.py
index 093cb8a..9a4f3bd 100644
--- a/hierarc/Sampling/ParamManager/param_manager.py
+++ b/hierarc/Sampling/ParamManager/param_manager.py
@@ -16,12 +16,18 @@ def __init__(
         anisotropy_sampling=False,
         anisotropy_model="OM",
         anisotropy_distribution="NONE",
+        gamma_in_sampling=False,
+        gamma_in_distribution="NONE",
+        m2l_sampling=False,
+        m2l_distribution="NONE",
         kappa_ext_sampling=False,
         kappa_ext_distribution="NONE",
         lambda_ifu_sampling=False,
         lambda_ifu_distribution="NONE",
         alpha_lambda_sampling=False,
         beta_lambda_sampling=False,
+        alpha_gamma_in_sampling=False,
+        alpha_m2l_sampling=False,
         sigma_v_systematics=False,
         sne_apparent_m_sampling=False,
         sne_distribution="GAUSSIAN",
@@ -56,6 +62,14 @@ def __init__(
         :param anisotropy_sampling: bool, if True adds a global stellar anisotropy parameter that alters the single lens
         kinematic prediction
         :param anisotropy_distribution: string, indicating the distribution function of the anisotropy model
+        :param gamma_in_sampling: bool, if True samples gNFW inner slope parameter
+        :param gamma_in_distribution: string, distribution function of the gamma_in parameter
+        :param m2l_sampling: bool, if True samples the mass-to-light ratio of the stars
+        :param m2l_distribution: string, distribution function of the m2l parameter
+        :param alpha_gamma_in_sampling: bool, if True samples a parameter alpha_gamma_in, which scales gamma_in linearly
+            according to a predefined quantity of the lens
+        :param alpha_m2l_sampling: bool, if True samples a parameter alpha_m2l, which scales m2l linearly
+            according to a predefined quantity of the lens
         :param sne_apparent_m_sampling: boolean, if True, samples/queries SNe unlensed magnitude distribution
          (not intrinsic magnitudes but apparent!)
         :param sne_distribution: string, apparent non-lensed brightness distribution (in linear space).
@@ -83,10 +97,16 @@ def __init__(
             lambda_mst_distribution=lambda_mst_distribution,
             lambda_ifu_sampling=lambda_ifu_sampling,
             lambda_ifu_distribution=lambda_ifu_distribution,
+            gamma_in_sampling=gamma_in_sampling,
+            gamma_in_distribution=gamma_in_distribution,
+            m2l_sampling=m2l_sampling,
+            m2l_distribution=m2l_distribution,
             kappa_ext_sampling=kappa_ext_sampling,
             kappa_ext_distribution=kappa_ext_distribution,
             alpha_lambda_sampling=alpha_lambda_sampling,
             beta_lambda_sampling=beta_lambda_sampling,
+            alpha_gamma_in_sampling=alpha_gamma_in_sampling,
+            alpha_m2l_sampling=alpha_m2l_sampling,
             log_scatter=log_scatter,
             kwargs_fixed=kwargs_fixed_lens,
         )
diff --git a/requirements.txt b/requirements.txt
index 00e04f0..ab46f10 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -7,7 +7,6 @@ mpmath
 h5py
 matplotlib
 pandas
-lenstronomy>=1.11.3
-# -e git://github.com/sibirrer/lenstronomy.git@main#egg=lenstronomy
+lenstronomy @ git+https://github.com/lenstronomy/lenstronomy@main
 scikit-learn>=0.22
 pandas


Test patch:
diff --git a/test/test_LensPosterior/test_kin_constraints.py b/test/test_LensPosterior/test_kin_constraints.py
index 0e86bc1..c01a900 100644
--- a/test/test_LensPosterior/test_kin_constraints.py
+++ b/test/test_LensPosterior/test_kin_constraints.py
@@ -6,7 +6,7 @@
 import unittest
 
 
-class TestIFUKinPosterior(object):
+class TestKinConstraints(object):
     def setup(self):
         pass
 
diff --git a/test/test_LensPosterior/test_kin_constraints_composite.py b/test/test_LensPosterior/test_kin_constraints_composite.py
new file mode 100644
index 0000000..f855ee2
--- /dev/null
+++ b/test/test_LensPosterior/test_kin_constraints_composite.py
@@ -0,0 +1,361 @@
+from hierarc.LensPosterior.kin_constraints_composite import KinConstraintsComposite
+from lenstronomy.Analysis.kinematics_api import KinematicsAPI
+from hierarc.Likelihood.hierarchy_likelihood import LensLikelihood
+import numpy.testing as npt
+import numpy as np
+import pytest
+import unittest
+from astropy.cosmology import FlatLambdaCDM
+
+
+class TestKinConstraintsComposite(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)
+        m2l_array = np.linspace(0.1, 10, 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,
+            m2l_array=m2l_array,
+            rho0_array=rho0_array,
+            r_s_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
+        )
+
+        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 = {
+            "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)
+        m2l_array = np.linspace(0.1, 10, 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,
+            m2l_array=m2l_array,
+            rho0_array=rho0_array,
+            r_s_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 TestRaise(unittest.TestCase):
+    def test_raise(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)
+            m2l_array = np.linspace(0.1, 10, 5)
+            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,
+                m2l_array=m2l_array,
+                rho0_array=rho0_array,
+                r_s_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
+            )
+            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)
+            m2l_array = np.linspace(0.1, 10, 5)
+
+            rho0_array = 10 ** np.random.normal(8, 0.2, 5) / 1e6
+            r_s_array = np.random.normal(0.1, 0.01, 6)
+
+            KinConstraintsComposite(
+                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,
+                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
+            )
+
+
+if __name__ == "__main__":
+    pytest.main()
diff --git a/test/test_Likelihood/test_LensLikelihood/test_base_lens_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_base_lens_likelihood.py
index 41d057a..436a685 100644
--- a/test/test_Likelihood/test_LensLikelihood/test_base_lens_likelihood.py
+++ b/test/test_Likelihood/test_LensLikelihood/test_base_lens_likelihood.py
@@ -101,7 +101,7 @@ def test_log_likelihood(self):
             )
             print(likelihood_type)
             logl = likelihood.log_likelihood(
-                ddt=1, dd=1, aniso_scaling=None, sigma_v_sys_error=1, mu_intrinsic=1
+                ddt=1, dd=1, kin_scaling=None, sigma_v_sys_error=1, mu_intrinsic=1
             )
             print(logl)
             assert logl > -np.inf
@@ -116,7 +116,7 @@ def test_predictions_measurements(self):
             )
             ddt_measurement = likelihood.ddt_measurement()
             likelihood.sigma_v_measurement(sigma_v_sys_error=0)
-            likelihood.sigma_v_prediction(ddt=1, dd=1, aniso_scaling=1)
+            likelihood.sigma_v_prediction(ddt=1, dd=1, kin_scaling=1)
             assert len(ddt_measurement) == 2
 
 
diff --git a/test/test_Likelihood/test_LensLikelihood/test_ddt_dd_gauss_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_ddt_dd_gauss_likelihood.py
index 5a4f634..2387d28 100644
--- a/test/test_Likelihood/test_LensLikelihood/test_ddt_dd_gauss_likelihood.py
+++ b/test/test_Likelihood/test_LensLikelihood/test_ddt_dd_gauss_likelihood.py
@@ -35,9 +35,9 @@ def test_log_likelihood(self):
         )
         npt.assert_almost_equal(lnlog, -0.5, decimal=5)
 
-        aniso_scaling = [1 + self.dd_sigma / self.dd_mean]
+        kin_scaling = [1 + self.dd_sigma / self.dd_mean]
         lnlog = likelihood.log_likelihood(
-            ddt=self.ddt_mean, dd=self.dd_mean, aniso_scaling=aniso_scaling
+            ddt=self.ddt_mean, dd=self.dd_mean, kin_scaling=kin_scaling
         )
         npt.assert_almost_equal(lnlog, -0.5, decimal=5)
 
diff --git a/test/test_Likelihood/test_LensLikelihood/test_ddt_gauss_kin_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_ddt_gauss_kin_likelihood.py
index 19defa9..088639b 100644
--- a/test/test_Likelihood/test_LensLikelihood/test_ddt_gauss_kin_likelihood.py
+++ b/test/test_Likelihood/test_LensLikelihood/test_ddt_gauss_kin_likelihood.py
@@ -63,11 +63,11 @@ def setup(self):
     def test_log_likelihood(self):
         ddt, dd = 9, 0.9
         lnlog_tot = self.ddt_gauss_kin_likelihood.log_likelihood(
-            ddt, dd, aniso_scaling=None, sigma_v_sys_error=None
+            ddt, dd, kin_scaling=None, sigma_v_sys_error=None
         )
         lnlog_ddt = self.ddt_gauss_likelihood.log_likelihood(ddt, dd)
         lnlog_kin = self.kin_likelihood.log_likelihood(
-            ddt, dd, aniso_scaling=None, sigma_v_sys_error=None
+            ddt, dd, kin_scaling=None, sigma_v_sys_error=None
         )
         npt.assert_almost_equal(lnlog_tot, lnlog_ddt + lnlog_kin, decimal=5)
 
@@ -88,7 +88,7 @@ def test_sigma_v_measurement(self):
             sigma_v_predict,
             error_cov_predict,
         ) = self.ddt_gauss_kin_likelihood.sigma_v_prediction(
-            self.ddt_mean, self.dd_mean, aniso_scaling=1
+            self.ddt_mean, self.dd_mean, kin_scaling=1
         )
         assert sigma_v_predict[0] == self.sigma_v_measurement[0]
         assert error_cov_predict[0, 0] == 0
diff --git a/test/test_Likelihood/test_LensLikelihood/test_ddt_hist_kin_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_ddt_hist_kin_likelihood.py
index 3288f9b..01f9e65 100644
--- a/test/test_Likelihood/test_LensLikelihood/test_ddt_hist_kin_likelihood.py
+++ b/test/test_Likelihood/test_LensLikelihood/test_ddt_hist_kin_likelihood.py
@@ -64,7 +64,7 @@ def setup(self):
     def test_log_likelihood(self):
         logl_max = self._ddt_kin_likelihood.log_likelihood(ddt=self._ddt, dd=self._dd)
         logl = self._ddt_kin_likelihood.log_likelihood(
-            self._ddt, self._dd / (1 + self._sigma) ** 2, aniso_scaling=None
+            self._ddt, self._dd / (1 + self._sigma) ** 2, kin_scaling=None
         )
         npt.assert_almost_equal(logl - logl_max, -self._num_ifu / 2, decimal=5)
 
@@ -85,7 +85,7 @@ def test_sigma_v_measurement(self):
             sigma_v_predict,
             error_cov_predict,
         ) = self._ddt_kin_likelihood.sigma_v_prediction(
-            self._ddt, self._dd, aniso_scaling=1
+            self._ddt, self._dd, kin_scaling=1
         )
         assert sigma_v_predict[0] == self.sigma_v_measurement[0]
         assert error_cov_predict[0, 0] == 0
diff --git a/test/test_Likelihood/test_LensLikelihood/test_ds_dds_gauss_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_ds_dds_gauss_likelihood.py
index 44a51dd..665c9cf 100644
--- a/test/test_Likelihood/test_LensLikelihood/test_ds_dds_gauss_likelihood.py
+++ b/test/test_Likelihood/test_LensLikelihood/test_ds_dds_gauss_likelihood.py
@@ -35,7 +35,7 @@ def test_log_likelihood(self):
         npt.assert_almost_equal(lnlog, -0.5, decimal=5)
 
         aniso_scaling = [1.0 / (1 + self.ds_dds_sigma / self.ds_dds_mean)]
-        lnlog = likelihood.log_likelihood(ddt=ddt, dd=dd, aniso_scaling=aniso_scaling)
+        lnlog = likelihood.log_likelihood(ddt=ddt, dd=dd, kin_scaling=aniso_scaling)
         npt.assert_almost_equal(lnlog, -0.5, decimal=5)
 
 
diff --git a/test/test_Likelihood/test_LensLikelihood/test_kin_likelihood.py b/test/test_Likelihood/test_LensLikelihood/test_kin_likelihood.py
index a52bffe..8fc5108 100644
--- a/test/test_Likelihood/test_LensLikelihood/test_kin_likelihood.py
+++ b/test/test_Likelihood/test_LensLikelihood/test_kin_likelihood.py
@@ -33,12 +33,12 @@ def test_log_likelihood(self):
             error_cov_j_sqrt,
             normalized=False,
         )
-        logl = kin_likelihood.log_likelihood(ddt, dd, aniso_scaling=None)
+        logl = kin_likelihood.log_likelihood(ddt, dd, kin_scaling=None)
         npt.assert_almost_equal(logl, 0, decimal=5)
-        logl = kin_likelihood.log_likelihood(ddt * (0.9**2), dd, aniso_scaling=None)
+        logl = kin_likelihood.log_likelihood(ddt * (0.9**2), dd, kin_scaling=None)
         npt.assert_almost_equal(logl, -num_ifu / 2, decimal=5)
         logl = kin_likelihood.log_likelihood(
-            ddt * (1 - np.sqrt(0.1**2 + 0.1**2)) ** 2, dd, aniso_scaling=None
+            ddt * (1 - np.sqrt(0.1**2 + 0.1**2)) ** 2, dd, kin_scaling=None
         )
         npt.assert_almost_equal(logl, -num_ifu, decimal=5)
 
@@ -50,7 +50,7 @@ def test_log_likelihood(self):
         assert error_cov_measurement_[0, 0] == error_cov_measurement[0, 0]
 
         sigma_v_predict, error_cov_predict = kin_likelihood.sigma_v_prediction(
-            ddt, dd, aniso_scaling=1
+            ddt, dd, kin_scaling=1
         )
         assert sigma_v_predict[0] == sigma_v_measurement[0]
         assert error_cov_predict[0, 0] == 0
@@ -82,7 +82,7 @@ def test_error_systematic(self):
         logl = ifu_likelihood.log_likelihood(
             ddt * (1 - np.sqrt(0.1**2)) ** 2,
             dd,
-            aniso_scaling=1,
+            kin_scaling=1,
             sigma_v_sys_error=0.1,
         )
         npt.assert_almost_equal(logl, -1 / 2.0, decimal=5)
@@ -112,10 +112,10 @@ def test_log_likelihood_marg(self):
             normalized=True,
             sigma_sys_error_include=True,
         )
-        logl = ifu_likelihood.log_likelihood(ddt, dd, aniso_scaling=1)
+        logl = ifu_likelihood.log_likelihood(ddt, dd, kin_scaling=1)
         aniso_scaling = 0.8
         logl_ani = ifu_likelihood.log_likelihood(
-            ddt, dd * aniso_scaling, aniso_scaling=aniso_scaling
+            ddt, dd * aniso_scaling, kin_scaling=aniso_scaling
         )
         npt.assert_almost_equal(logl_ani - logl, 0, decimal=1)
 
@@ -123,14 +123,14 @@ def test_log_likelihood_marg(self):
         # analytical Gaussian calculation in the covariance matrix normalization
         sigma_v_sys_error = 0.1
         logl = ifu_likelihood.log_likelihood(
-            ddt, dd, aniso_scaling=1, sigma_v_sys_error=sigma_v_sys_error
+            ddt, dd, kin_scaling=1, sigma_v_sys_error=sigma_v_sys_error
         )
         l_sum = 0
         num_sample = 10000
         for i in range(num_sample):
             sigma_v_pert = np.random.normal(loc=0, scale=sigma_v_sys_error)
             logl_i = ifu_likelihood.log_likelihood(
-                ddt, dd, aniso_scaling=1, sigma_v_sys_offset=sigma_v_pert
+                ddt, dd, kin_scaling=1, sigma_v_sys_offset=sigma_v_pert
             )
             l_sum += np.exp(logl_i)
         logl_average = np.log(l_sum / num_sample)
@@ -149,7 +149,7 @@ def test_log_likelihood_invert_error(self):
             normalized=True,
             sigma_sys_error_include=True,
         )
-        logl = ifu_likelihood.log_likelihood(ddt=1, dd=1, aniso_scaling=1)
+        logl = ifu_likelihood.log_likelihood(ddt=1, dd=1, kin_scaling=1)
         assert logl == -np.inf
 
 
diff --git a/test/test_Likelihood/test_anisotropy_scaling.py b/test/test_Likelihood/test_anisotropy_scaling.py
deleted file mode 100644
index df1e7bd..0000000
--- a/test/test_Likelihood/test_anisotropy_scaling.py
+++ /dev/null
@@ -1,147 +0,0 @@
-import pytest
-import numpy as np
-import unittest
-from hierarc.Likelihood.anisotropy_scaling import (
-    AnisotropyScalingSingleAperture,
-    AnisotropyScalingIFU,
-)
-
-
-class TestAnisotropyScalingSingleAperture(object):
-    def setup(self):
-        ani_param_array = np.linspace(start=0, stop=1, num=10)
-        ani_scaling_array = ani_param_array * 2
-        self.scaling = AnisotropyScalingSingleAperture(
-            ani_param_array, ani_scaling_array
-        )
-
-        ani_param_array = [
-            np.linspace(start=0, stop=1, num=10),
-            np.linspace(start=1, stop=2, num=5),
-        ]
-        ani_scaling_array = np.outer(ani_param_array[0], ani_param_array[1])
-        self.scaling_2d = AnisotropyScalingSingleAperture(
-            ani_param_array, ani_scaling_array
-        )
-
-    def test_ani_scaling(self):
-        scaling = self.scaling.ani_scaling(aniso_param_array=[1])
-        assert scaling == 2
-
-        scaling = self.scaling.ani_scaling(aniso_param_array=None)
-        assert scaling == 1
-
-        scaling = self.scaling_2d.ani_scaling(aniso_param_array=[1, 2])
-        assert scaling == 2
-
-
-class TestAnisotropyScalingIFU(object):
-    def setup(self):
-        ani_param_array = np.linspace(start=0, stop=1, num=10)
-        ani_scaling_array = ani_param_array * 2
-        self.scaling = AnisotropyScalingIFU(
-            anisotropy_model="OM",
-            ani_param_array=ani_param_array,
-            ani_scaling_array_list=[ani_scaling_array],
-        )
-
-        ani_param_array = [
-            np.linspace(start=0, stop=1, num=10),
-            np.linspace(start=1, stop=2, num=5),
-        ]
-        ani_scaling_array = np.outer(ani_param_array[0], ani_param_array[1])
-        self.scaling_2d = AnisotropyScalingIFU(
-            anisotropy_model="GOM",
-            ani_param_array=ani_param_array,
-            ani_scaling_array_list=[ani_scaling_array],
-        )
-
-    def test_ani_scaling(self):
-        scaling = self.scaling.ani_scaling(aniso_param_array=[1])
-        assert scaling[0] == 2
-
-        scaling = self.scaling.ani_scaling(aniso_param_array=None)
-        assert scaling[0] == 1
-
-        scaling = self.scaling_2d.ani_scaling(aniso_param_array=[1, 2])
-        assert scaling[0] == 2
-
-    def test_draw_anisotropy(self):
-        a_ani = 1
-        beta_inf = 1.5
-        param_draw = self.scaling.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.draw_anisotropy(
-                a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1
-            )
-
-        param_draw = self.scaling_2d.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_2d.draw_anisotropy(
-                a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1
-            )
-
-        scaling = AnisotropyScalingIFU(anisotropy_model="NONE")
-        param_draw = scaling.draw_anisotropy(
-            a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0
-        )
-        assert param_draw is None
-
-
-class TestRaise(unittest.TestCase):
-    def test_raise(self):
-        with self.assertRaises(ValueError):
-            ani_param_array = [
-                np.linspace(start=0, stop=1, num=10),
-                np.linspace(start=1, stop=2, num=5),
-                1,
-            ]
-            ani_scaling_array = np.outer(ani_param_array[0], ani_param_array[1])
-            self.scaling_2d = AnisotropyScalingSingleAperture(
-                ani_param_array, ani_scaling_array
-            )
-
-        with self.assertRaises(ValueError):
-            AnisotropyScalingIFU(
-                anisotropy_model="blabla",
-                ani_param_array=np.array([0, 1]),
-                ani_scaling_array_list=[np.array([0, 1])],
-            )
-
-        with self.assertRaises(ValueError):
-            ani_param_array = np.linspace(start=0, stop=1, num=10)
-            ani_scaling_array = ani_param_array * 2
-            scaling = AnisotropyScalingIFU(
-                anisotropy_model="OM",
-                ani_param_array=ani_param_array,
-                ani_scaling_array_list=[ani_scaling_array],
-            )
-            scaling.draw_anisotropy(
-                a_ani=-1, a_ani_sigma=0, beta_inf=-1, beta_inf_sigma=0
-            )
-
-        with self.assertRaises(ValueError):
-            ani_param_array = [
-                np.linspace(start=0, stop=1, num=10),
-                np.linspace(start=1, stop=2, num=5),
-            ]
-            ani_scaling_array = np.outer(ani_param_array[0], ani_param_array[1])
-            scaling = AnisotropyScalingIFU(
-                anisotropy_model="GOM",
-                ani_param_array=ani_param_array,
-                ani_scaling_array_list=[ani_scaling_array],
-            )
-            scaling.draw_anisotropy(
-                a_ani=-1, a_ani_sigma=0, beta_inf=-1, beta_inf_sigma=0
-            )
-
-
-if __name__ == "__main__":
-    pytest.main()
diff --git a/test/test_Likelihood/test_hierarchy_likelihood.py b/test/test_Likelihood/test_hierarchy_likelihood.py
index 319c3a5..a0ce60e 100644
--- a/test/test_Likelihood/test_hierarchy_likelihood.py
+++ b/test/test_Likelihood/test_hierarchy_likelihood.py
@@ -92,6 +92,50 @@ def setup(self):
             **kwargs_likelihood
         )
 
+        gamma_in_array = np.linspace(start=0.1, stop=2.9, num=5)
+        m2l_array = np.linspace(start=1, stop=10, 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)),
+        )
+        self.likelihood_gamma_in_m2l_list_ani = 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,
+            m2l_array=m2l_array,
+            num_distribution_draws=200,
+            kappa_ext_bias=False,
+            kappa_pdf=None,
+            kappa_bin_edges=None,
+            mst_ifu=False,
+            **kwargs_likelihood
+        )
+
+        self.likelihood_gamma_in_m2l = 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,
+            m2l_array=m2l_array,
+            num_distribution_draws=200,
+            kappa_ext_bias=False,
+            kappa_pdf=None,
+            kappa_bin_edges=None,
+            mst_ifu=False,
+            **kwargs_likelihood
+        )
+
     def test_lens_log_likelihood(self):
         np.random.seed(42)
         kwargs_lens = {
@@ -132,9 +176,36 @@ def test_lens_log_likelihood(self):
         )
         assert ln_inf < -10000000
 
+        ln_inf = self.likelihood_single.lens_log_likelihood(
+            self.cosmo, kwargs_lens=None, kwargs_kin=kwargs_kin
+        )
+        npt.assert_almost_equal(ln_inf, 0.0, decimal=1)
+
+        ln_inf = self.likelihood_single.sigma_v_measured_vs_predict(
+            self.cosmo, kwargs_lens=None, kwargs_kin=None
+        )
+        assert np.all([x is None for x in ln_inf])
+
         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()
+        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,
+        }
+        ln_likelihood = self.likelihood_gamma_in_m2l.lens_log_likelihood(
+            self.cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin
+        )
+        npt.assert_almost_equal(ln_likelihood, -0.0, decimal=1)
+
 
 if __name__ == "__main__":
     pytest.main()
diff --git a/test/test_Likelihood/test_parameter_scaling.py b/test/test_Likelihood/test_parameter_scaling.py
new file mode 100644
index 0000000..dad4f0c
--- /dev/null
+++ b/test/test_Likelihood/test_parameter_scaling.py
@@ -0,0 +1,310 @@
+import pytest
+import numpy as np
+import unittest
+from hierarc.Likelihood.parameter_scaling import (
+    ParameterScalingSingleAperture,
+    ParameterScalingIFU,
+)
+
+
+class TestParameterScalingSingleAperture(object):
+    def setup(self):
+        ani_param_array = np.linspace(start=0, stop=1, num=10)
+        param_scaling_array = ani_param_array * 2
+        self.scaling = ParameterScalingSingleAperture(
+            ani_param_array, param_scaling_array
+        )
+
+        ani_param_array = [
+            np.linspace(start=0, stop=1, num=10),
+            np.linspace(start=1, stop=2, num=5),
+        ]
+        param_scaling_array = np.outer(ani_param_array[0], ani_param_array[1])
+        self.scaling_2d = ParameterScalingSingleAperture(
+            ani_param_array, param_scaling_array
+        )
+
+        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)
+
+        param_arrays = [ani_param_array, gamma_in_array, m2l_array]
+        param_scaling_array = np.multiply.outer(
+            ani_param_array,
+            np.outer(gamma_in_array, m2l_array),
+        )
+        self.scaling_nfw = ParameterScalingSingleAperture(
+            param_arrays, 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,
+            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),
+            ),
+        )
+
+        self.scaling_nfw_2d = ParameterScalingSingleAperture(
+            param_arrays, param_scaling_array
+        )
+
+    def test_param_scaling(self):
+        scaling = self.scaling.param_scaling(param_array=[1])
+        assert scaling == np.array([2])
+
+        scaling = self.scaling.param_scaling(param_array=None)
+        assert scaling == 1
+
+        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_2d.param_scaling(param_array=[1, 2, 2.9, 10])
+        assert scaling == 1 * 2 * 2.9 * 10
+
+
+class TestParameterScalingIFU(object):
+    def setup(self):
+        ani_param_array = np.linspace(start=0, stop=1, num=10)
+        param_scaling_array = ani_param_array * 2
+        self.scaling = ParameterScalingIFU(
+            anisotropy_model="OM",
+            param_arrays=ani_param_array,
+            scaling_grid_list=[param_scaling_array],
+        )
+
+        ani_param_array = [
+            np.linspace(start=0, stop=1, num=10),
+            np.linspace(start=1, stop=2, num=5),
+        ]
+        param_scaling_array = np.outer(ani_param_array[0], ani_param_array[1])
+        self.scaling_2d = ParameterScalingIFU(
+            anisotropy_model="GOM",
+            param_arrays=ani_param_array,
+            scaling_grid_list=[param_scaling_array],
+        )
+
+        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)
+
+        param_arrays = [ani_param_array, gamma_in_array, m2l_array]
+        param_scaling_array = np.multiply.outer(
+            ani_param_array,
+            np.outer(gamma_in_array, m2l_array),
+        )
+        self.scaling_nfw = 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,
+            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),
+            ),
+        )
+        self.scaling_nfw_2d = 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
+
+        scaling = self.scaling.param_scaling(param_array=[1])
+        assert scaling[0] == 2
+
+        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_2d.param_scaling(param_array=[1, 2, 2.9, 10])
+        assert scaling[0] == 1 * 2 * 2.9 * 10
+
+    def test_draw_anisotropy(self):
+        a_ani = 1
+        beta_inf = 1.5
+        param_draw = self.scaling.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.draw_anisotropy(
+                a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1
+            )
+        self.scaling._anisotropy_model = "const"
+        param_draw = self.scaling.draw_anisotropy(
+            a_ani=1, a_ani_sigma=0, beta_inf=beta_inf, beta_inf_sigma=0
+        )
+        assert param_draw[0] == 1
+
+        param_draw = self.scaling_2d.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_2d.draw_anisotropy(
+                a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1
+            )
+
+        param_draw = self.scaling_nfw.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.draw_anisotropy(
+                a_ani=1, a_ani_sigma=1, beta_inf=beta_inf, beta_inf_sigma=1
+            )
+
+        param_draw = self.scaling_nfw_2d.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.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
+        )
+        assert param_draw is None
+
+    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
+        )
+        assert param_draw[0] == 1
+        assert param_draw[1] == 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
+            )
+
+        param_draw = self.scaling_nfw_2d.draw_lens_parameters(
+            gamma_in=1, gamma_in_sigma=0, m2l=5, m2l_sigma=0
+        )
+        assert param_draw[0] == 1
+        assert param_draw[1] == 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
+            )
+
+
+class TestRaise(unittest.TestCase):
+    def test_raise(self):
+        with self.assertRaises(ValueError):
+            ParameterScalingIFU(
+                anisotropy_model="blabla",
+                param_arrays=np.array([0, 1]),
+                scaling_grid_list=[np.array([0, 1])],
+            )
+
+        ani_param_array = np.linspace(start=0, stop=1, num=10)
+        ani_scaling_array = ani_param_array * 2
+        scaling = ParameterScalingIFU(
+            anisotropy_model="OM",
+            param_arrays=ani_param_array,
+            scaling_grid_list=[ani_scaling_array],
+        )
+        with self.assertRaises(ValueError):
+            scaling.draw_anisotropy(
+                a_ani=-1, a_ani_sigma=0, beta_inf=-1, beta_inf_sigma=0
+            )
+
+        ani_param_array = [
+            np.linspace(start=0, stop=1, num=10),
+            np.linspace(start=1, stop=2, num=5),
+        ]
+        ani_scaling_array = np.outer(ani_param_array[0], ani_param_array[1])
+        scaling = ParameterScalingIFU(
+            anisotropy_model="GOM",
+            param_arrays=ani_param_array,
+            scaling_grid_list=[ani_scaling_array],
+        )
+        with self.assertRaises(ValueError):
+            scaling.draw_anisotropy(
+                a_ani=0.5, a_ani_sigma=0, beta_inf=-1, beta_inf_sigma=0
+            )
+
+        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)
+
+        param_arrays = [ani_param_array, gamma_in_array, m2l_array]
+        param_scaling_array = np.multiply.outer(
+            ani_param_array,
+            np.multiply.outer(gamma_in_array, m2l_array),
+        )
+        print(param_scaling_array.shape)
+        scaling = ParameterScalingIFU(
+            anisotropy_model="OM",
+            param_arrays=param_arrays,
+            scaling_grid_list=[param_scaling_array],
+        )
+        with self.assertRaises(ValueError):
+            scaling.draw_lens_parameters(
+                gamma_in=-1, gamma_in_sigma=0.1, m2l=5, 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
+            )
+
+        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)
+
+        param_arrays = [ani_param_array, gamma_in_array, m2l_array]
+        param_scaling_array = np.multiply.outer(
+            ani_param_array,
+            np.multiply.outer(gamma_in_array, m2l_array),
+        )
+
+        scaling = ParameterScalingIFU(
+            anisotropy_model="OM",
+            param_arrays=param_arrays,
+            scaling_grid_list=[param_scaling_array],
+        )
+        with self.assertRaises(ValueError):
+            scaling.draw_lens_parameters(
+                gamma_in=1, gamma_in_sigma=0, m2l=0, m2l_sigma=0
+            )
+
+
+if __name__ == "__main__":
+    pytest.main()
diff --git a/test/test_Sampling/test_ParamManager/test_lens_param.py b/test/test_Sampling/test_ParamManager/test_lens_param.py
index 577438e..5b09406 100644
--- a/test/test_Sampling/test_ParamManager/test_lens_param.py
+++ b/test/test_Sampling/test_ParamManager/test_lens_param.py
@@ -91,3 +91,86 @@ def test_args2kwargs(self):
         kwargs_new, i = self._param_fixed.args2kwargs(args, i=0)
         args_new = self._param_fixed.kwargs2args(kwargs_new)
         npt.assert_almost_equal(args_new, args)
+
+
+class TestLensParamGammaInnerM2l(object):
+    def setup(self):
+        self._param = LensParam(
+            gamma_in_sampling=True,
+            gamma_in_distribution="GAUSSIAN",
+            m2l_sampling=True,
+            m2l_distribution="GAUSSIAN",
+            alpha_gamma_in_sampling=True,
+            alpha_m2l_sampling=True,
+            kwargs_fixed={},
+            log_scatter=False,
+        )
+
+        kwargs_fixed = {
+            "gamma_in": 1,
+            "gamma_in_sigma": 0.1,
+            "m2l": 5,
+            "m2l_sigma": 1,
+            "alpha_gamma_in": 0.1,
+            "alpha_m2l": -0.5,
+        }
+        self._param_fixed = LensParam(
+            gamma_in_sampling=True,
+            gamma_in_distribution="GAUSSIAN",
+            m2l_sampling=True,
+            m2l_distribution="GAUSSIAN",
+            alpha_gamma_in_sampling=True,
+            alpha_m2l_sampling=True,
+            kwargs_fixed=kwargs_fixed,
+            log_scatter=False,
+        )
+        self._param_log_scatter = LensParam(
+            gamma_in_sampling=True,
+            gamma_in_distribution="GAUSSIAN",
+            m2l_sampling=True,
+            m2l_distribution="GAUSSIAN",
+            alpha_gamma_in_sampling=True,
+            alpha_m2l_sampling=True,
+            kwargs_fixed={},
+            log_scatter=True,
+        )
+
+    def test_param_list(self):
+        param_list = self._param.param_list(latex_style=False)
+        assert len(param_list) == 6
+        param_list = self._param.param_list(latex_style=True)
+        assert len(param_list) == 6
+
+        param_list = self._param_log_scatter.param_list(latex_style=False)
+        assert len(param_list) == 6
+        param_list = self._param_log_scatter.param_list(latex_style=True)
+        assert len(param_list) == 6
+
+        param_list = self._param_fixed.param_list(latex_style=False)
+        assert len(param_list) == 0
+        param_list = self._param_fixed.param_list(latex_style=True)
+        assert len(param_list) == 0
+
+    def test_args2kwargs(self):
+        kwargs = {
+            "gamma_in": 1,
+            "gamma_in_sigma": 0.1,
+            "m2l": 5,
+            "m2l_sigma": 1,
+            "alpha_gamma_in": 0.1,
+            "alpha_m2l": -0.5,
+        }
+        args = self._param.kwargs2args(kwargs)
+        kwargs_new, i = self._param.args2kwargs(args, i=0)
+        args_new = self._param.kwargs2args(kwargs_new)
+        npt.assert_almost_equal(args_new, args)
+
+        args = self._param_log_scatter.kwargs2args(kwargs)
+        kwargs_new, i = self._param_log_scatter.args2kwargs(args, i=0)
+        args_new = self._param_log_scatter.kwargs2args(kwargs_new)
+        npt.assert_almost_equal(args_new, args)
+
+        args = self._param_fixed.kwargs2args(kwargs)
+        kwargs_new, i = self._param_fixed.args2kwargs(args, i=0)
+        args_new = self._param_fixed.kwargs2args(kwargs_new)
+        npt.assert_almost_equal(args_new, args)