diff --git a/cupix/likelihood/CAMB_model.py b/cupix/likelihood/CAMB_model.py
deleted file mode 100644
index 168d9a5..0000000
--- a/cupix/likelihood/CAMB_model.py
+++ /dev/null
@@ -1,132 +0,0 @@
-import numpy as np
-import copy
-import camb
-from lace.cosmo import camb_cosmo
-from lace.cosmo import fit_linP
-from cupix.likelihood import likelihood_parameter
-import types
-
-
-class CAMBModel(object):
-    """Interface between CAMB object and Theory"""
-
-    def __init__(self, zs, cosmo=None):
-        """Setup from CAMB object and list of redshifts"""
-
-        # list of redshifts at which we evaluate linear power
-        self.zs = zs
-
-        # setup CAMB cosmology object
-        if cosmo is None:
-            self.cosmo = camb_cosmo.get_cosmology()
-        else:
-            self.cosmo = cosmo
-        
-        # cache CAMB results when computed
-        self.cached_camb_results = None
-        # cache wavenumbers and linear power (at zs) when computed
-        self.cached_linP_Mpc = None
-        self.cached_linP_params = None
-
-    def get_camb_results(self):
-        """Check if we have called CAMB.get_results yet, to save time.
-        It returns a CAMB.results object."""
-
-        if self.cached_camb_results is None:
-            print("Calling get_camb_results")
-            self.cached_camb_results = camb_cosmo.get_camb_results(
-                self.cosmo, zs=self.zs, camb_kmax_Mpc=200 #, fast_camb=True
-            )
-
-        return self.cached_camb_results
-
-    def get_linP_Mpc(self):
-        """Check if we have already computed linP_Mpc, to save time.
-        It returns (k_Mpc, zs, linP_Mpc)."""
-
-        if self.cached_linP_Mpc is None:
-            camb_results = self.get_camb_results()
-            self.cached_linP_Mpc = camb_cosmo.get_linP_Mpc(
-                pars=self.cosmo, zs=self.zs, camb_results=camb_results
-            )
-
-        return self.cached_linP_Mpc
-
-
-    def get_linP_Mpc_params(self, kp_Mpc):
-        """Get linear power parameters to call emulator, at each z.
-        Amplitude, slope and running around pivot point kp_Mpc."""
-
-        ## Get the P(k) at each z
-        k_Mpc, z, pk_Mpc = self.get_linP_Mpc()
-        # specify wavenumber range to fit
-        kmin_Mpc = 0.5 * kp_Mpc
-        kmax_Mpc = 2.0 * kp_Mpc
-
-        linP_params = []
-        ## Fit the emulator call params
-        for pk_z in pk_Mpc:
-            linP_Mpc = fit_linP.fit_polynomial(
-                kmin_Mpc / kp_Mpc,
-                kmax_Mpc / kp_Mpc,
-                k_Mpc / kp_Mpc,
-                pk_z,
-                deg=2,
-            )
-            # translate the polynomial to our parameters
-            ln_A_p = linP_Mpc[0]
-            Delta2_p = np.exp(ln_A_p) * kp_Mpc**3 / (2 * np.pi**2)
-            n_p = linP_Mpc[1]
-            # note that the curvature is alpha/2
-            alpha_p = 2.0 * linP_Mpc[2]
-            linP_z = {"Delta2_p": Delta2_p, "n_p": n_p, "alpha_p": alpha_p}
-            linP_params.append(linP_z)
-
-        return linP_params
-
-    
-    def get_dkms_dMpc(self):
-        """Return M(z)=H(z)/(1+z) for each z"""
-        dkms_of_zs = []
-        for z in self.zs:    
-            dkms_of_zs.append(camb_cosmo.dkms_dMpc(self.cosmo, z))
-
-        return dkms_of_zs
-    
-    def get_dAA_dMpc(self):
-        """Return conversion factor from inv cMpc to inv Angstroms for each z"""
-        dAA_of_zs = []
-        for z in self.zs:
-            dAA_of_zs.append(camb_cosmo.dAA_dMpc(self.cosmo, z, 1215.67))
-
-        return dAA_of_zs
-    
-    def get_ddeg_dMpc(self):
-        """Return conversion factor from Mpc to deg for each z"""
-        ddeg_of_zs = []
-        for z in self.zs:
-            ddeg_of_zs.append(camb_cosmo.ddeg_dMpc(self.cosmo, z))
-
-        return ddeg_of_zs
-
-    def get_new_model(self, zs, like_params):
-        """For an arbitrary list of like_params, return a new CAMBModel"""
-
-        # store a dictionary with parameters set to input values
-        camb_param_dict = {}
-
-        # loop over list of likelihood parameters own by this object
-        for mypar in self.get_likelihood_parameters():
-            # loop over input parameters
-            for inpar in like_params:
-                if inpar.name == mypar.name:
-                    camb_param_dict[inpar.name] = inpar.value
-                    continue
-
-        # set cosmology object (use fiducial for parameters not provided)
-        new_cosmo = camb_cosmo.get_cosmology_from_dictionary(
-            camb_param_dict, cosmo_fid=self.cosmo
-        )
-
-        return CAMBModel(zs=zs, cosmo=new_cosmo)
-
diff --git a/cupix/likelihood/cosmologies.py b/cupix/likelihood/cosmologies.py
deleted file mode 100644
index 6e6b97b..0000000
--- a/cupix/likelihood/cosmologies.py
+++ /dev/null
@@ -1,164 +0,0 @@
-import os
-import numpy as np
-
-import lace
-from lace.cosmo import camb_cosmo
-
-
-def get_cosmology_from_label(cosmo_label="default"):
-    if cosmo_label == "default":
-        return camb_cosmo.get_cosmology()
-    elif cosmo_label == "low_omch2":
-        return camb_cosmo.get_cosmology(omch2=0.11)
-    elif cosmo_label == "high_omch2":
-        return camb_cosmo.get_cosmology(omch2=0.13)
-    elif cosmo_label == "omch2_0115":
-        return camb_cosmo.get_cosmology(omch2=0.115)
-    elif cosmo_label == "omch2_0125":
-        return camb_cosmo.get_cosmology(omch2=0.125)
-    elif cosmo_label == "mnu_03":
-        return camb_cosmo.get_cosmology(mnu=0.3)
-    elif cosmo_label == "mnu_06":
-        return camb_cosmo.get_cosmology(mnu=0.6)
-    elif cosmo_label == "SHOES":
-        return camb_cosmo.get_cosmology(H0=73.0)
-    else:
-        raise ValueError("implement cosmo_label " + cosmo_label)
-
-
-def set_cosmo(
-    cosmo_label="mpg_central",
-    return_all=False,
-    nyx_version="models_Nyx_Mar2025_with_CGAN_val_3axes",
-):
-    """Set fiducial cosmology
-
-    Parameters
-    ----------
-    cosmo_label : str
-
-    Returns
-    -------
-    cosmo : object
-    """
-    if (cosmo_label[:3] == "mpg") | (cosmo_label[:3] == "nyx"):
-        if cosmo_label[:3] == "mpg":
-            repo = os.path.dirname(lace.__path__[0]) + "/"
-            fname = repo + ("data/sim_suites/Australia20/mpg_emu_cosmo.npy")
-            get_cosmo = camb_cosmo.get_cosmology_from_dictionary
-        elif cosmo_label[:3] == "nyx":
-            fname = (
-                os.environ["NYX_PATH"] + "nyx_emu_cosmo_" + nyx_version + ".npy"
-            )
-            get_cosmo = camb_cosmo.get_Nyx_cosmology
-
-        try:
-            data_cosmo = np.load(fname, allow_pickle=True).item()
-        except:
-            raise ValueError(f"{fname} not found")
-
-        if cosmo_label in data_cosmo.keys():
-            # print(data_cosmo[cosmo_label]["cosmo_params"])
-            cosmo = get_cosmo(data_cosmo[cosmo_label]["cosmo_params"])
-        else:
-            raise ValueError(f"Cosmo not found in {fname} for {cosmo_label}")
-
-    elif cosmo_label == "Planck18":
-        cosmo = camb_cosmo.get_cosmology(
-            H0=67.66,
-            mnu=0.0,
-            omch2=0.119,
-            ombh2=0.0224,
-            omk=0.0,
-            As=2.105e-09,
-            ns=0.9665,
-            nrun=0.0,
-            pivot_scalar=0.05,
-            w=-1,
-        )
-    elif cosmo_label == "Planck18_nyx":
-        cosmo = camb_cosmo.get_cosmology(
-            H0=67.66,
-            mnu=0.0,
-            omch2=0.119,
-            ombh2=0.0224,
-            omk=0.0,
-            As=2.24e-09,
-            ns=0.937,
-            nrun=0.0,
-            pivot_scalar=0.05,
-            w=-1,
-        )
-    elif cosmo_label == "Planck18_mpg":
-        cosmo = camb_cosmo.get_cosmology(
-            H0=67.0,
-            mnu=0.0,
-            omch2=0.119,
-            ombh2=0.022,
-            omk=0.0,
-            As=2.26e-09,
-            ns=0.982,
-            nrun=0.0,
-            pivot_scalar=0.05,
-            w=-1,
-        )
-    elif cosmo_label == "Planck18_h74":
-        cosmo = camb_cosmo.get_cosmology(
-            H0=74,
-            mnu=0.0,
-            omch2=0.119,
-            ombh2=0.0224,
-            omk=0.0,
-            As=2.105e-09,
-            ns=0.9665,
-            nrun=0.0,
-            pivot_scalar=0.05,
-            w=-1,
-        )
-    elif (cosmo_label == "ACCEL2_6144_160") | (cosmo_label == "accel2"):
-        # https://arxiv.org/pdf/2407.04473
-        # Planck15 ΛCDM Planck TT,TE,EE+lowP (approx...)
-        Omegam = 0.31
-        Omegab = 0.0487
-        h = 0.675
-        omch2 = (Omegam - Omegab) * h**2
-        ombh2 = Omegab * h**2
-        cosmo = camb_cosmo.get_cosmology(
-            H0=h * 100,
-            mnu=0.0,
-            omch2=omch2,
-            ombh2=ombh2,
-            omk=0.0,
-            As=np.exp(3.094) / 1e10,  # Planck15 ΛCDM Planck TT,TE,EE+lowP
-            ns=0.96,
-            nrun=0.0,
-            pivot_scalar=0.05,
-            w=-1,
-        )
-    elif (cosmo_label == "Sherwood_2048_40") | (cosmo_label == "sherwood"):
-        # https://academic.oup.com/mnras/article/464/1/897/2236089
-        # Planck13 ΛCDM Planck+WP+highL+BAO
-        Omegam = 0.308
-        Omegab = 0.0482
-        h = 0.678
-        omch2 = (Omegam - Omegab) * h**2
-        ombh2 = Omegab * h**2
-        cosmo = camb_cosmo.get_cosmology(
-            H0=h * 100,
-            mnu=0.0,
-            omch2=omch2,
-            ombh2=ombh2,
-            omk=0.0,
-            As=np.exp(3.0973) / 1e10,  # Planck13 ΛCDM Planck+WP+highL+BAO
-            ns=0.961,
-            nrun=0.0,
-            pivot_scalar=0.05,
-            w=-1,
-        )
-    else:
-        raise ValueError(f"cosmo_label {cosmo_label} not implemented")
-
-    if return_all:
-        return data_cosmo
-    else:
-        return cosmo
diff --git a/cupix/likelihood/fitter.py b/cupix/likelihood/fitter.py
index d5c56ba..989dc16 100644
--- a/cupix/likelihood/fitter.py
+++ b/cupix/likelihood/fitter.py
@@ -10,7 +10,6 @@
 
 # our own modules
 import cupix, lace
-from lace.cosmo import camb_cosmo
 from cupix.utils.utils import create_print_function, purge_chains
 
 
diff --git a/cupix/likelihood/forestflow_emu.py b/cupix/likelihood/forestflow_emu.py
index 8d5d143..c67a292 100644
--- a/cupix/likelihood/forestflow_emu.py
+++ b/cupix/likelihood/forestflow_emu.py
@@ -14,7 +14,6 @@ class FF_emulator():
     def __init__(
             self,
             z,
-            cosmo_param_dict,
             Nrealizations=3000
         ):
 
@@ -27,7 +26,6 @@ def __init__(
             "kF_Mpc"
         ]
         self.z = z
-        self.cosmo_param_dict = cosmo_param_dict 
         self.emulator_label = "forestflow_emu"
         # self.kp_Mpc = kp_Mpc
         # self.kmax_Mpc = 5 # from Forestflow paper plots, could revisit
@@ -103,4 +101,4 @@ def emulate_P3D_params(self, emu_call, zs):
                 arinyo_coeffs[key] = emu_call[key]
         
         return arinyo_coeffs
-        
\ No newline at end of file
+        
diff --git a/cupix/likelihood/likelihood.py b/cupix/likelihood/likelihood.py
index b3b3632..f479ecd 100644
--- a/cupix/likelihood/likelihood.py
+++ b/cupix/likelihood/likelihood.py
@@ -7,8 +7,6 @@
 from scipy.optimize import minimize
 from scipy.linalg import block_diag
 from cupix.likelihood.window_and_rebin import convolve_window, rebin_theta
-import lace
-from lace.cosmo import camb_cosmo
 from cupix.utils.utils import is_number_string
 from cupix.likelihood.likelihood_parameter import LikelihoodParameter, likeparam_from_dict
 import warnings
diff --git a/cupix/likelihood/likelihood_parameter.py b/cupix/likelihood/likelihood_parameter.py
index 630b4b1..6fd28b9 100644
--- a/cupix/likelihood/likelihood_parameter.py
+++ b/cupix/likelihood/likelihood_parameter.py
@@ -83,6 +83,8 @@ def get_new_parameter(self, value_in_cube):
 
 def likeparam_from_dict(params_dict):
     # takes a dictionary of parameter names and values, and returns a list of LikelihoodParameter objects
+    if params_dict is None:
+        return []
     like_params = []
     for name, value in params_dict.items():
         like_params.append(LikelihoodParameter(name=name, min_value=-1000, max_value=1000, value=value))
@@ -91,6 +93,8 @@ def likeparam_from_dict(params_dict):
 
 def dict_from_likeparam(like_params):
     # takes a list of LikelihoodParameter objects, and returns a dictionary of parameter names and values
+    if like_params is None:
+        return {}
     params_dict = {}
     for param in like_params:
         params_dict[param.name] = param.value
@@ -120,4 +124,4 @@ def par_index(like_param_list, par_name):
     for i, param in enumerate(like_param_list):
         if param.name == par_name:
             return i
-    raise ValueError("Parameter with name " + par_name + " not found in list of LikelihoodParameter objects.")
\ No newline at end of file
+    raise ValueError("Parameter with name " + par_name + " not found in list of LikelihoodParameter objects.")
diff --git a/cupix/likelihood/lya_theory.py b/cupix/likelihood/lya_theory.py
index b485cee..c002345 100644
--- a/cupix/likelihood/lya_theory.py
+++ b/cupix/likelihood/lya_theory.py
@@ -1,8 +1,7 @@
 import numpy as np
 import copy
-from lace.cosmo import camb_cosmo
+from lace.cosmo import cosmology
 from cupix.likelihood.forestflow_emu import FF_emulator
-from cupix.likelihood import CAMB_model
 from cupix.likelihood.lyaP3D import LyaP3D
 from cupix.likelihood.likelihood_parameter import likeparam_from_dict, LikelihoodParameter, dict_from_likeparam, format_like_params_dict
 import sys
@@ -10,25 +9,6 @@
 from forestflow import priors
 from astropy.io import fits
 
-def set_theory(
-    zs, cosmo_dict=None, default_theory='best_fit_arinyo_from_p1d', p3d_label='arinyo', emulator_label=None, k_unit='iAA', verbose=False
-):
-    """Set theory"""
-
-
-    # set theory
-    theory = Theory(
-        zs = zs,
-        cosmo_dict = cosmo_dict,
-        default_lya_theory=default_theory,
-        p3d_label=p3d_label,
-        emulator_label = emulator_label,
-        k_unit = k_unit,
-        verbose=verbose
-    )
-
-    return theory
-
 
 class Theory(object):
     """Translator between the likelihood object and the emulator. This object
@@ -38,7 +18,7 @@ class Theory(object):
     def __init__(
         self,
         zs,
-        cosmo_dict=None,
+        fid_cosmo=None,
         default_lya_theory='best_fit_arinyo_from_p1d',
         p3d_label='arinyo',
         emulator_label=None,
@@ -48,7 +28,7 @@ def __init__(
         """Setup object to compute predictions for the 1D power spectrum.
         Inputs:
             - zs: redshifts that will be evaluated
-            - cosmo: dictionary of cosmology parameters used for unit conversions, setting the linear power spectrum, etc
+            - fid_cosmo: fiducial Cosmology object, used for unit conversions, setting the linear power spectrum, etc
             - default_lya_theory: string tag to specify which default set of likelihood parameters to use. Options are 'best_fit_arinyo_from_p1d' to use the best-fit Arinyo coefficients from the p1d paper, 'best_fit_igm_from_p1d' to use the best-fit IGM parameters from the p1d paper, and 'best_fit_arinyo_from_colore' to use the best-fit Arinyo coefficients from Laura's CF fits.
             - p3d_label: string tag to specify which P3D model to use. Current option is 'arinyo' to use the Arinyo model.
             - emulator_label: string tag to specify which emulator to use for predicting the Arinyo coefficients. Current option is 'forestflow_emu' to use the forestflow emulator.
@@ -58,122 +38,48 @@ def __init__(
         self.zs = zs
         # specify pivot point used in compressed parameters
         self.k_unit = k_unit
-        self.set_cosmo(cosmo_dict)
+        if fid_cosmo is None:
+            self.fid_cosmo = cosmology.Cosmology()
+        else:
+            self.fid_cosmo = fid_cosmo
         self._load_P3D_model(p3d_label=p3d_label)
         self.default_lya_theory = default_lya_theory
         self.set_default_param_values(tag=default_lya_theory)
         if emulator_label == "forestflow_emu":
-            self.emulator = FF_emulator(zs, self.cosmo_dict)
+            self.emulator = FF_emulator(zs)
         elif emulator_label is None:
             assert "igm" not in self.default_lya_theory, "If default_lya_theory includes 'igm', an emulator must be specified to predict the Arinyo parameters. Please specify an emulator_label, for example 'forestflow_emu'."
         else:
             raise ValueError("Emulator label not recognized. Current only option is 'forestflow_emu'.")
 
-    def set_cosmo_dict(self, cosmo_dict_input):
-        # Make a dictionary with some default values, if none are passed
-        # Replace any passed by user-specified inputs)
-        omnuh2 = 0.0006
-        mnu = omnuh2 * 93.14
-        H0 = 67.36
-        omch2 = 0.12
-        ombh2 = 0.02237
-        As = 2.1e-9
-        ns = 0.9649
-        nrun = 0.0
-        w = -1.0
-        omk = 0
-        cosmo_dict = {
-            'H0': H0,
-            'omch2': omch2,
-            'ombh2': ombh2,
-            'mnu': mnu,
-            'omk': omk,
-            'As': As,
-            'ns': ns,
-            'nrun': nrun,
-            'w': w
-        }
-        if cosmo_dict_input is not None:
-            for key in cosmo_dict_input:
-                if key in cosmo_dict:
-                    cosmo_dict[key] = cosmo_dict_input[key] # update default value with user-provided value
-                else:
-                    print(f"Warning: {key} not recognized as a cosmological parameter, ignoring it.")
-        self.cosmo_dict = cosmo_dict
-
-    def set_cosmo(
-        self,
-        cosmo_dict
-    ):
-        """Setup initial cosmology"""
-        # first, set the comology dictionary
-        self.set_cosmo_dict(cosmo_dict)
-        laceCosmo = camb_cosmo.get_cosmology_from_dictionary(self.cosmo_dict)
-        self.cosmo = {}
-        self.cosmo["zs"] = self.zs
-        self.cosmo["cosmo"] = CAMB_model.CAMBModel(
-            zs=self.zs,
-            cosmo=laceCosmo,
-        )
-        self.cosmo["dkms_dMpc_zs"] = self.cosmo["cosmo"].get_dkms_dMpc()
-        self.cosmo["dAA_dMpc_zs"] = self.cosmo["cosmo"].get_dAA_dMpc()
-        self.cosmo["ddeg_dMpc_zs"] = self.cosmo["cosmo"].get_ddeg_dMpc()
-
-    
-    def fixed_background(self, like_params):
-        """Check if any of the input likelihood parameters would change
-        the background expansion of the initial cosmology"""
-        if like_params is None:
-            return True
-        # look for parameters that would change background
-        for par in like_params:
-            if par.name in ["ombh2", "omch2", "H0", "mnu", "cosmomc_theta"]:
-                return False
-
-        return True
-
 
-    def get_unit_conversions(self, zs, like_params, return_blob=False):
+    def get_unit_conversions(self, zs, like_params):
         """return conversion from Mpc to km/s or Mpc to inv Angstroms,
             and from Mpc to deg for transverse separations"""
-        
-        # compute linear power parameters at all redshifts, and H(z) / (1+z)
-        if self.fixed_background(like_params):
+ 
+        # convert list of LikelihoodParameters to dictionary
+        params_dict = dict_from_likeparam(like_params)
+        print('params_dict', params_dict)
+
+        if self.fid_cosmo.same_background(cosmo_params=params_dict):
             # use background and transfer functions from initial cosmology
+            cosmo = self.fid_cosmo
             if self.verbose:
                 print("recycle transfer function")
-            dkms_dMpc_zs  = []
-            ddeg_dMpc_zs = []
-            dAA_dMpc_zs   = []
-            for z in zs:
-                _ = np.argwhere(self.cosmo["zs"] == z)[0, 0]
-                dkms_dMpc_zs.append(self.cosmo["dkms_dMpc_zs"][_])
-                ddeg_dMpc_zs.append(self.cosmo["ddeg_dMpc_zs"][_])
-                dAA_dMpc_zs.append(self.cosmo["dAA_dMpc_zs"][_])
-            dkms_dMpc_zs  = np.array(dkms_dMpc_zs)
-            ddeg_dMpc_zs = np.array(ddeg_dMpc_zs)
-            dAA_dMpc_zs   = np.array(dAA_dMpc_zs)
-            if return_blob:
-                blob = self.get_blob_fixed_background(like_params)
         else:
-            # setup a new CAMB_model from like_params
-            if self.verbose:
-                print("create new CAMB_model")
-            camb_model = self.cosmo["cosmo"].get_new_model(zs, like_params)
-            dkms_dMpc_zs = camb_model.get_dkms_dMpc()
-            dAA_dMpc_zs = camb_model.get_dAA_dMpc()
-            ddeg_dMpc_zs = camb_model.get_ddeg_dMpc()
-            if return_blob:
-                blob = self.get_blob(camb_model=camb_model)
+            cosmo = cosmology.Cosmology(cosmo_params_dict=params_dict)
+
+        dkms_dMpc_zs  = np.array([cosmo.get_dkms_dMpc(z) for z in zs])
+        ddeg_dMpc_zs = np.array([cosmo.get_ddeg_dMpc(z) for z in zs])
+        dAA_dMpc_zs   = np.array([cosmo.get_dAA_dMpc(z) for z in zs])
+
         if self.k_unit == 'ikms':
             return_conv_k_of_zs = dkms_dMpc_zs
         else:
             return_conv_k_of_zs = dAA_dMpc_zs
-        if return_blob:
-            return return_conv_k_of_zs, ddeg_dMpc_zs, blob
-        else:
-            return return_conv_k_of_zs, ddeg_dMpc_zs
-            
+        return return_conv_k_of_zs, ddeg_dMpc_zs
+
+
     def get_emulator_calls(
         self, iz_choice, theory_inputs
     ):
@@ -185,11 +91,17 @@ def get_emulator_calls(
             if 'T0' in key:
                 sigT_kms = thermal_broadening_kms(theory_inputs[key])
                 z_int = int(key.split('_')[-1])
-                sigT_Mpc = sigT_kms / self.cosmo["dkms_dMpc_zs"][z_int]
+                print('WARNING: check this code (Andreu)')
+                z=self.zs[z_int]
+                # I don't this can assume fiducial cosmo...
+                sigT_Mpc = sigT_kms / self.fid_cosmo.get_dkms_dMpc(z=z)
                 theory_inputs[f'sigT_Mpc_{z_int}'] = sigT_Mpc
             elif 'kF_kms' in key:
                 z_int = int(key.split('_')[-1])
-                kF_Mpc = theory_inputs[key] / self.cosmo["dkms_dMpc_zs"][z_int]
+                print('WARNING: check this code (Andreu)')
+                z=self.zs[z_int]
+                # I don't this can assume fiducial cosmo...
+                kF_Mpc = theory_inputs[key] / self.fid_cosmo.get_dkms_dMpc(z=z)
                 theory_inputs[f'kF_Mpc_{z_int}'] = kF_Mpc
             # later, when varying cosmology, can add here the conversions for Delta* and n* to Deltap and np
 
@@ -269,9 +181,7 @@ def get_px_AA(
     ):
         """Emulate Px in velocity units, for all redshift bins,
         as a function of input likelihood parameters.
-        theta_arcmin is a list of theta bins.
-        It might also return a covariance from the emulator,
-        or a blob with extra information for the fitter."""
+        theta_arcmin is a list of theta bins."""
 
         if verbose is None:
             verbose = self.verbose
@@ -438,7 +348,7 @@ def _load_P3D_model(self, p3d_label='arinyo'):
         """
         if p3d_label == 'arinyo':
             from forestflow.model_p3d_arinyo import ArinyoModel
-            arinyo = ArinyoModel(cosmo=self.cosmo_dict) # set model
+            arinyo = ArinyoModel(fid_cosmo=self.fid_cosmo)
         else:
             sys.exit("Error: no P3D model specified. Please choose a valid P3D model. Current option are 'arinyo'.")
         
@@ -465,4 +375,4 @@ def get_param(self, param_name, iz_choice=None):
         else:
             return self.default_param_dict[param_name]
 
-    
\ No newline at end of file
+    
diff --git a/notebooks/test/example_user_interface_notebook.py b/notebooks/test/example_user_interface_notebook.py
index d81252b..9d68300 100644
--- a/notebooks/test/example_user_interface_notebook.py
+++ b/notebooks/test/example_user_interface_notebook.py
@@ -6,19 +6,19 @@
 #       extension: .py
 #       format_name: percent
 #       format_version: '1.3'
-#       jupytext_version: 1.19.1
+#       jupytext_version: 1.16.4
 #   kernelspec:
-#     display_name: Python 3
+#     display_name: cupix
 #     language: python
-#     name: python3
+#     name: cupix
 # ---
 
 # %%
 import numpy as np
 from cupix.likelihood.likelihood import Likelihood
-from cupix.likelihood.lya_theory import set_theory
+from cupix.likelihood.lya_theory import Theory
 from cupix.likelihood.forestflow_emu import FF_emulator
-from lace.cosmo import camb_cosmo
+from lace.cosmo import cosmology
 import matplotlib.pyplot as plt
 from cupix.likelihood.likelihood_parameter import LikelihoodParameter
 from cupix.px_data.data_DESI_DR2 import DESI_DR2
@@ -35,7 +35,7 @@
 # %%
 data_file = "/global/cfs/cdirs/desi/users/sindhu_s/Lya_Px_measurements/DR2_Px/baseline/binned_out_px-zbins_4-thetabins_9_w_res.hdf5"
 # "../../data/px_measurements/forecast/forecast_ffcentral_cosmo_igm_real_binned_out_px-zbins_4-thetabins_9_w_res_noiseless.hdf5"
-data = DESI_DR2(data_file, kM_max_cut_AA=1, km_max_cut_AA=
+data = DESI_DR2(data_file, kM_max_cut_AA=1, km_max_cut_AA=1.2)
 # kmax_cut determines max of the widely-binned k; finely-binned k will still be fully used for the window matrix application
 
 # %%
@@ -79,18 +79,14 @@ def plot_theta_bins(data, k_M, iz, it_M):
 
 # %%
 z = data.z
-cosmo = {'H0': 67}
+cosmo = cosmology.Cosmology(cosmo_params_dict={'H0': 67})
 # default_theory options are:
 # 'best_fit_arinyo_from_p1d': best fit to the DESI DR1 P1D data from Chaves+2026
 # 'best_fit_igm_from_p1d': same but for IGM parameters
 # 'best_fit_arinyo_from_colore': best fit to xi from colore mocks. Only works for z=2.2, 2.4, 2.6, 2.8
-theory_p1d = set_theory(z, cosmo_dict=cosmo, default_theory='best_fit_igm_from_p1d', p3d_label='arinyo', emulator_label='forestflow_emu', k_unit='iAA', verbose=True)
+theory_p1d = Theory(z, fid_cosmo=cosmo, default_lya_theory='best_fit_igm_from_p1d', p3d_label='arinyo', emulator_label='forestflow_emu', k_unit='iAA', verbose=True)
 # theory_colore = set_theory(z, bkgd_cosmo=cosmo, default_theory='best_fit_arinyo_from_colore', p3d_label='arinyo', emulator_label='forestflow_emu', k_unit='iAA', verbose=True)
 
-# %%
-# check full cosmo dictionary to see what other default parameters were used
-theory_p1d.cosmo_dict
-
 # %%
 theory_p1d.default_param_dict
 
@@ -308,3 +304,5 @@ def plot_theta_bins(data, k_M, iz, it_M):
 likelihood.get_chi2()
 
 # %%
+
+# %%
diff --git a/notebooks/test/likelihood_evals.py b/notebooks/test/likelihood_evals.py
index a175b4f..f9fdc10 100644
--- a/notebooks/test/likelihood_evals.py
+++ b/notebooks/test/likelihood_evals.py
@@ -6,11 +6,11 @@
 #       extension: .py
 #       format_name: percent
 #       format_version: '1.3'
-#       jupytext_version: 1.19.1
+#       jupytext_version: 1.16.4
 #   kernelspec:
-#     display_name: Python 3
+#     display_name: cupix
 #     language: python
-#     name: python3
+#     name: cupix
 # ---
 
 # %% [markdown]
@@ -19,9 +19,9 @@
 # %%
 import numpy as np
 from cupix.likelihood.likelihood import Likelihood
-from cupix.likelihood.lya_theory import set_theory
+from cupix.likelihood.lya_theory import Theory
 from cupix.likelihood.forestflow_emu import FF_emulator
-from lace.cosmo import camb_cosmo
+from lace.cosmo import cosmology
 import matplotlib.pyplot as plt
 from cupix.likelihood.likelihood_parameter import LikelihoodParameter, par_index, dict_from_likeparam
 from cupix.px_data.data_DESI_DR2 import DESI_DR2
@@ -41,14 +41,8 @@
 # mode = 'arinyo' # if the parameters you want to test are Arinyo params
 
 # %%
-forecast_file = f"{cupixpath}/data/px_measurements/forecast//forecast_ffcentral_igm_real_binned_out_px-zbins_4-thetabins_9_w_res_noiseless.hdf5"
-# forecast_file = f"../../data/px_measurements/forecast/forecast_ffcentral_{mode}_real_binned_out_px-zbins_4-thetabins_9_w_res_noiseless.hdf5"
-forecast = DESI_DR2(forecast_file, kmax_cut_AA=1)
-# get default theory from forecast
-with h5.File(forecast_file, 'r') as f:
-    default_theory_label = f['metadata'].attrs['true_lya_theory']
-print(f"Default theory label: {default_theory_label}")
-
+forecast_file = f"../../data/px_measurements/forecast/forecast_ffcentral_igm_real_binned_out_px-zbins_4-thetabins_9_w_res_noiseless.hdf5"
+forecast = DESI_DR2(forecast_file, kM_max_cut_AA=1, km_max_cut_AA=1.2)
 
 # %% [markdown]
 # ### Step 2: Load the emulator with Nrealizations = 1000
@@ -67,7 +61,7 @@
 nrun = 0.0
 w = -1.0
 omk = 0
-cosmo = {
+cosmo_dict = {
     'H0': H0,
     'omch2': omch2,
     'ombh2': ombh2,
@@ -78,8 +72,11 @@
     'nrun': nrun,
     'w': w
 }
+fid_cosmo = cosmology.Cosmology(cosmo_params_dict=cosmo_dict)
+#fid_cosmo = cosmology.Cosmology(label='Planck18')
 
-theory = set_theory(z, bkgd_cosmo=cosmo, default_theory=default_theory_label, p3d_label='arinyo', emulator_label='forestflow_emu', k_unit='iAA', verbose=True)
+theory = Theory(z, fid_cosmo=fid_cosmo, default_lya_theory='best_fit_igm_from_p1d', 
+                p3d_label='arinyo', emulator_label='forestflow_emu', k_unit='iAA', verbose=True)
 
 
 # %% [markdown]
@@ -116,9 +113,8 @@
 
 # %%
 # make sure theory matches forecast data at central value
-theta_A_index = 5
-like_params = None # only use the default cosmo+IGM theory, not like_params as it contains the evaluated Arinyo params
-Px_convolved = like.get_convolved_Px_AA(theta_A=theta_A_index, like_params=like_params)[0] 
+theta_A_index = 8
+Px_convolved = like.get_convolved_Px_AA(theta_A=theta_A_index, like_params=like_params)[0]
 
 # %%
 plt.plot(forecast.k_M_centers_AA, forecast.Px_ZAM[iz_choice, theta_A_index])
@@ -176,7 +172,7 @@
         best = par_vals[i]
         print(f"New best fit found with chi2={chi2_i} at mF={best}")
         best_chi2 = chi2_i
-    
+
 
 # %%
 par_i = par_index(like_params, par_to_test)
@@ -193,3 +189,9 @@
 plt.legend()
 
 # %%
+
+# %%
+
+# %%
+
+# %%