diff --git a/deerlab/__init__.py b/deerlab/__init__.py
index 43cd8a84d..5e715d4c9 100644
--- a/deerlab/__init__.py
+++ b/deerlab/__init__.py
@@ -15,7 +15,7 @@ def __getattr__(name):
 from .deerload import deerload
 from .selregparam import selregparam
 from .dipolarkernel import dipolarkernel
-from .dipolarbackground import dipolarbackground
+from .dipolarbackground import dipolarbackground,dipolarbackgroundmodel
 from .dipolarmodel import dipolarmodel,ExperimentInfo, dipolarpenalty, ex_4pdeer,ex_3pdeer,ex_rev5pdeer,ex_fwd5pdeer,ex_ridme,ex_sifter,ex_dqc
 from .solvers import snlls, fnnls, cvxnnls
 from .regoperator import regoperator
diff --git a/deerlab/bg_models.py b/deerlab/bg_models.py
index 0345b0a89..15b33d4a4 100644
--- a/deerlab/bg_models.py
+++ b/deerlab/bg_models.py
@@ -7,12 +7,12 @@
 import math as m
 from numpy import pi
 import inspect
+from scipy.special import gamma, hyp2f1, sici
+from deerlab.constants import *
 from copy import deepcopy as _deepcopy
 from deerlab.dipolarkernel import dipolarkernel
 from deerlab.utils import formatted_table
 from deerlab.model import Model
-from scipy.special import gamma, hyp2f1, sici
-from deerlab.constants import *
 
 #---------------------------------------------------------------------------------------
 def hyp2f1_repro(a,b,c,z): 
diff --git a/deerlab/dipolarbackground.py b/deerlab/dipolarbackground.py
index 908927cd3..070ca8014 100644
--- a/deerlab/dipolarbackground.py
+++ b/deerlab/dipolarbackground.py
@@ -1,7 +1,7 @@
 # dipolarbackground.py - Multipathway background generator
 # --------------------------------------------------------
 # This file is a part of DeerLab. License is MIT (see LICENSE.md).
-# Copyright(c) 2019-2023: Luis Fabregas, Stefan Stoll and other contributors.
+# Copyright(c) 2019-2025: Luis Fabregas, Stefan Stoll, Hugo Karas and other contributors.
 
 import numpy as np
 import inspect
@@ -206,3 +206,243 @@ def basis(t,lam):
                 B *= basis(tdip)
 
     return B
+
+def dipolarbackgroundmodel(experiment=None, basis=None, parametrization='reftimes', spins=2, samespins=True, **kwargs):
+    """
+    Construct a dipolar background model for a given dipolar experiment.
+    This model can be used for evaluating and extrapolating the fitted parameters of the dipolar background.
+    Supports both two-spin and multi-spin systems.
+
+    Parameters
+    ----------
+    experiment : :ref:`ExperimentInfo` or None, optional
+        The Experimental information obtained from experiment models (ex_). If the experiment is not specified, a single-pathway background model is constructed with a refocusing time of 0.
+
+    basis : :ref:`Model`, optional
+        The basis model for the intermolecular (background) contribution. If not specified, a background arising from a homogenous 3D distribution of spins is used.
+
+    parametrization : string, optional
+        Parametrization strategy of the dipolar pathway refocusing times. Must be one of the following:
+
+            * ``'reftimes'`` - Each refocusing time is represented individually as a parameter.
+            * ``'delays'`` - The pulse delays are introduced as parameters from which the refocusing times are computed. Requires ``experiment`` to be specified.
+            * ``'shift'`` - A time shift is introduced as a parameter to represent the variability of the refocusing times from their theoretical values. Requires ``experiment`` to be specified.
+
+        The default is ``'reftimes'``.
+
+    spins : integer, optional
+        Number of spins in the system. If not specified, defaults to two-spin systems. For multi-spin
+        systems (``spins>2``), three-spin interaction contributions are accounted for in the background.
+        The default is ``2``.
+
+    samespins : boolean, optional
+        Only relevant when ``spins>2``. Enables the assumption of spectral permutability, i.e., the
+        assumption that all pairwise interactions have the same amplitude. If enabled, one amplitude
+        parameter ``lam{n}`` is used per pathway. If disabled, individual amplitudes ``lam{n}_{q}``
+        are used for each pathway-interaction combination. Enabled by default.
+
+    **kwargs:
+        Additional keyword arguments. If ``experiment`` is not specified, the number of pathways can be specified with the keyword argument ``npathways`` (default is 1).
+
+    Returns
+    -------
+
+    BModel : :ref:`Model`
+
+    Examples
+    --------
+    To construct a dipolar background model from the fit results using the specified :ref:`ExperimentInfo`::
+
+            Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+            Bfit = results.P_scale*results.evaluate(Bfcn,t)
+            Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
+
+    To construct a dipolar background model for a three-spin system::
+
+            Bfcn = dl.dipolarbackgroundmodel(experimentInfo, spins=3)
+            Bfit = results.P_scale*results.evaluate(Bfcn, t)
+    """
+    from deerlab.dipolarmodel import _populate_dipolar_pathways_parameters
+    from deerlab.model import Model
+    from itertools import permutations
+    from math import factorial
+
+    if basis is None:
+        from deerlab.bg_models import bg_hom3d
+        basis = bg_hom3d
+
+    if experiment is None:
+        npathways = kwargs.get('npathways', 1)
+        harmonics = [1]
+        if npathways>1:
+            pulsedelay_names = ['reftime'+str(i+1) for i in range(npathways)]
+        else:
+            pulsedelay_names = ['reftime']
+        exp_labels = list(np.arange(1,npathways+1,dtype=int))
+    else:
+        npathways = experiment.npathways
+        harmonics = np.array(experiment.harmonics)[:npathways]
+        pulsedelay_names = list(inspect.signature(experiment.reftimes).parameters.keys())
+        exp_labels = experiment.labels
+
+    Q = spins*(spins-1)//2
+    
+
+    signature = []
+    basis_signature = basis.signature.copy()
+    lam_in_basis = 'lam' in basis_signature
+    if lam_in_basis:
+        basis_signature.remove('lam')
+    signature.extend(basis_signature)
+
+    if spins == 2:
+        if parametrization=='reftimes':
+            if npathways==1:
+                signature.extend(['mod','reftime'])
+            else:
+                signature.extend(['lam'+str(i) for i in exp_labels])
+                signature.extend(['reftime'+str(i) for i in exp_labels])
+        elif parametrization=='delays':
+            signature.extend(pulsedelay_names)
+            signature.extend(['lam'+str(i) for i in exp_labels])
+        elif parametrization=='shift':
+            signature.extend(['tshift'])
+            signature.extend(['lam'+str(i) for i in exp_labels])
+    else:
+        if parametrization=='reftimes':
+            if samespins:
+                signature.extend(['lam'+str(i) for i in exp_labels])
+            else:
+                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in exp_labels for q in range(Q)])
+            signature.extend(['reftime'+str(i) for i in exp_labels])
+            signature.append('lamu')
+        elif parametrization=='delays':
+            signature.extend(pulsedelay_names)
+            if samespins:
+                signature.extend(['lam'+str(i) for i in exp_labels])
+            else:
+                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in exp_labels for q in range(Q)])
+            signature.append('lamu')
+        elif parametrization=='shift':
+            signature.extend(['tshift'])
+            if samespins:
+                signature.extend(['lam'+str(i) for i in exp_labels])
+            else:
+                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in exp_labels for q in range(Q)])
+            signature.append('lamu')
+
+    # Construct the dipolar background model
+    def bckgrnd_fun(*param:tuple):
+        basis_constants = list(param[0:len(basis_signature)])
+        t = basis_constants[0] # t is always the first basis parameter
+        pathway_params = np.array(param[len(basis_signature):])
+
+        if spins == 2:
+            if parametrization=='reftimes':
+                λs = pathway_params[np.arange(0, len(pathway_params)//2, 1,)]
+                trefs = pathway_params[np.arange(len(pathway_params)//2,len(pathway_params),1)]
+            elif parametrization=='delays':
+                delays = pathway_params[0:len(experiment.delays)]
+                λs = pathway_params[len(delays):]
+                trefs = experiment.reftimes(*delays)
+            elif parametrization=='shift':
+                shift = pathway_params[0]
+                λs = pathway_params[1:]
+                delays = np.array(experiment.delays) + shift
+                trefs = experiment.reftimes(*delays)
+
+            if lam_in_basis:
+                basis_constants.append(λs[0])
+
+            prod = np.ones_like(t)
+            scale = 1
+            for i in range(npathways):
+                scale -= λs[i]
+                if lam_in_basis:
+                    basis_constants[-1] = λs[i]
+                prod *= basis(t-trefs[i],*basis_constants[1:])
+            return scale*prod
+
+        else:
+            # Multi-spin: parse parameters
+            if parametrization=='reftimes':
+                n_lam = npathways if samespins else npathways*Q
+                λs_flat = pathway_params[:n_lam]
+                trefs = pathway_params[n_lam:n_lam+npathways]
+                lamu = pathway_params[-1]
+            elif parametrization=='delays':
+                delays = pathway_params[0:len(experiment.delays)]
+                remaining = pathway_params[len(experiment.delays):]
+                trefs = experiment.reftimes(*delays)
+                λs_flat = remaining[:-1]
+                lamu = remaining[-1]
+            elif parametrization=='shift':
+                shift = pathway_params[0]
+                remaining = pathway_params[1:]
+                delays = np.array(experiment.delays) + shift
+                trefs = experiment.reftimes(*delays)
+                λs_flat = remaining[:-1]
+                lamu = remaining[-1]
+
+            # Build λs as [Q x npathways] array: λs_matrix[q, idx]
+            if samespins:
+                λs_matrix = np.tile(λs_flat, (Q, 1))
+            else:
+                λs_matrix = np.array(λs_flat).reshape(npathways, Q).T
+
+            # Wrap basis into a callable for dipolarbackground
+            if lam_in_basis:
+                def basis_fcn(tdip, lam):
+                    bc = basis_constants.copy()
+                    bc.append(lam)
+                    return basis(tdip, *bc[1:])
+            else:
+                def basis_fcn(tdip):
+                    return basis(tdip, *basis_constants[1:])
+
+            # Construct all multi-spin dipolar pathways
+            paths, threespin_paths = [], []
+            Λ0 = 1.0
+
+            for idx in range(npathways):
+                tref_idx = trefs[idx]
+                δ_idx = harmonics[idx]
+
+                for perm in set(permutations([0]+[None]*(Q-1))):
+                    q = int(np.where(np.array(perm)==0)[0][0])
+                    λp = λs_matrix[q, idx]
+                    λ2k = factorial(Q-1)*λp*lamu**(Q-1)
+                    tref_path = tuple([tref_idx if n==0 else None for n in perm])
+                    δs_path = tuple([δ_idx if n==0 else 0 for n in perm])
+                    paths.append({'amp': λ2k, 'reftime': tref_path, 'harmonic': δs_path})
+                    Λ0 -= λ2k
+
+                    for idx2 in range(idx+1, npathways):
+                        tref_idx2 = trefs[idx2]
+                        δ_idx2 = harmonics[idx2]
+                        for perm2 in set(permutations([0,1]+[None]*(Q-2))):
+                            perm2_arr = np.array(perm2)
+                            q2 = int(np.where(perm2_arr==1)[0][0])
+                            λm = λs_matrix[q2, idx]
+                            λ3k = 2*factorial(Q-2)*λp*λm*lamu**(Q-2)
+                            tref_path2 = tuple([tref_idx if n==0 else tref_idx2 if n==1 else None for n in perm2])
+                            δs_path2 = tuple([δ_idx if n==0 else δ_idx2 if n==1 else 0 for n in perm2])
+                            threespin_paths.append({'amp': λ3k, 'reftime': tref_path2, 'harmonic': δs_path2})
+                            Λ0 -= λ3k
+
+            paths.append({'amp': Λ0})
+            return dipolarbackground(t, paths+threespin_paths, basis_fcn)
+
+    # Define the dipolar background model
+    bckgrnd_model = Model(bckgrnd_fun, signature=signature,constants=['t'])
+    bckgrnd_model = _populate_dipolar_pathways_parameters(bckgrnd_model, npathways, spins=spins, samespins=samespins,
+                                                          experiment=experiment, parametrization=parametrization,
+                                                          Q=Q if spins>2 else None)
+
+    # Copy the basis parameters to the background model
+    for param in basis_signature:
+        if param =='t':
+            continue
+        getattr(bckgrnd_model,param).setas(getattr(basis,param))
+
+    return bckgrnd_model
diff --git a/deerlab/dipolarmodel.py b/deerlab/dipolarmodel.py
index 7076872e7..97becc300 100644
--- a/deerlab/dipolarmodel.py
+++ b/deerlab/dipolarmodel.py
@@ -4,7 +4,8 @@
 # Copyright(c) 2019-2022: Luis Fabregas, Stefan Stoll and other contributors.
 
 import numpy as np
-from deerlab.dipolarkernel import dipolarkernel, dipolarbackground
+from deerlab.dipolarkernel import dipolarkernel
+from deerlab.dipolarbackground import dipolarbackground, dipolarbackgroundmodel
 from deerlab.regoperator import regoperator
 from deerlab.dd_models import freedist
 from deerlab.model import Model,Penalty, link
@@ -16,7 +17,97 @@
 import inspect
 from scipy.stats import multivariate_normal
 
+
 #===============================================================================
+def _populate_dipolar_pathways_parameters(model,npathways,spins=2,samespins=True,experiment=None,parametrization='reftimes',Q=None):
+    """
+    A private function to populate the dipolar pathways parameters in the model.
+
+    Parameters
+    ----------
+    model : :ref:`Model`
+        Model object to be populated with the dipolar pathways parameters.
+    npathways : int
+        Number of dipolar pathways.
+    spins : int
+        Number of spins in the system.
+    samespins : bool    
+        Enables the assumption of spectral permutability.
+    experiment : :ref:`ExperimentInfo`
+        Experimental information obtained from experiment models (``ex_``).
+    parametrization : str
+        Parametrization strategy of the dipolar pathway refocusing times.
+    Q : int
+        Number of interactions in the spins system.
+
+    Returns
+    -------
+    model : :ref:`Model`
+        Model object populated with the dipolar pathways parameters.
+    """
+
+    # Populate the basic information on the dipolar pathways parameters
+    if experiment is not None:
+        labels = experiment.labels
+        pulsedelay_names = inspect.signature(experiment.reftimes).parameters.keys()
+    else:
+        # Otherwise, just label them in order
+        labels = np.arange(1,npathways+1)
+    
+    if spins>2:
+        if Q is None:
+            raise ValueError('If spins>2, the number of interactions Q must be specified.')
+        
+    
+    for n in range(npathways):
+        if spins==2 and npathways==1:
+            # Special case: use modulation depth notation instead of general pathway amplitude
+            getattr(model,f'mod').set(lb=0,ub=1,par0=0.01,description=f'Modulation depth',unit='')
+        else:
+            pairwise = ''
+            if spins>2: 
+                pairwise = ' pairwise'   
+            if spins==2 or samespins:
+                getattr(model,f'lam{labels[n]}').set(lb=0,ub=1,par0=0.01,description=f'Amplitude of{pairwise} pathway #{labels[n]}',unit='')
+            else:
+                for q in range(Q):
+                    getattr(model,f'lam{labels[n]}_{q+1}').set(lb=0,ub=1,par0=0.01,description=f'Amplitude of{pairwise} pathway #{labels[n]} on interaction #{q+1}',unit='')
+        if parametrization=='reftimes':
+            if experiment is None:
+                theoretical_reftime = 0
+                reftime_variability = 20
+            else:
+                theoretical_reftime = experiment.reftimes(*experiment.delays)[n]
+                reftime_variability = 3*experiment.pulselength
+            if spins==2 and npathways==1:
+                # Special case: use reftime notation
+                getattr(model,f'reftime').set(description=f'Refocusing time',unit='μs',
+                    par0 = theoretical_reftime,  
+                    lb   = theoretical_reftime - reftime_variability,
+                    ub   = theoretical_reftime + reftime_variability)
+            else: 
+                # Special case: use reftime notation
+                getattr(model,f'reftime{labels[n]}').set(description=f'Refocusing time of{pairwise} pathway #{labels[n]}',unit='μs',
+                    par0 = theoretical_reftime,
+                    lb   = theoretical_reftime - reftime_variability,
+                    ub   = theoretical_reftime + reftime_variability)
+    if parametrization=='delays': 
+        experimental_delays = experiment.delays
+        delay_variability = 5*experiment.pulselength
+        for n,delay in enumerate(pulsedelay_names):
+            getattr(model,delay).set(description=f'Pulse delay {delay} of the experiment',unit='μs',
+                par0 = experimental_delays[n],
+                lb   = experimental_delays[n] - delay_variability,
+                ub   = experimental_delays[n] + delay_variability)
+    elif parametrization=='shift':
+        variability = 5*experiment.pulselength
+        getattr(model,'tshift').set(par0=0,lb=-variability,ub=variability,description=f'Variability of experimental pulse delays',unit='μs')
+    if spins>2: 
+        getattr(model,f'lamu').set(lb=0,ub=1,par0=0.5,description='Amplitude of unmodulated pairwise pathway',unit='')
+    
+    return model
+
+
 def dipolarmodel(t, r=None, Pmodel=None, Bmodel=bg_hom3d, experiment=None, parametrization='reftimes', npathways=1,  spins=2, harmonics=None,
                     excbandwidth=np.inf, orisel=None, g=[ge,ge], gridsize=1000, minamp=1e-3, samespins=True, triangles=None, interp=True,):
     """
@@ -121,7 +212,13 @@ def dipolarmodel(t, r=None, Pmodel=None, Bmodel=bg_hom3d, experiment=None, param
     Returns
     -------
     Vmodel : :ref:`Model`
-        Dipolar signal model object.
+        Dipolar signal model object. With these additional attributes:
+            * ``t`` - The time axis of the model.
+            * ``r`` - The distance axis of the model, only for two-spin systems.
+            * ``Bmodel`` - The multi-pathway background model object used in the construction of the dipolar signal model.
+            * ``Pmodel`` - The distance distribution model object used in the construction of the dipolar signal model.
+            
+
     """
 
     # If distance axis not specified default to a long range one
@@ -344,53 +441,9 @@ def _Pmultivar(param):
             [getattr(Pmodel,f'chol{j+1}{q+1}').set(lb=-1,ub=1,par0=0.0,unit='nm',description=f'Cholesky factor ℓ{j+1}{q+1}') for j in np.arange(q+1,Q)]
     Nconstants = len(Pmodel._constantsInfo)
 
-    # Populate the basic information on the dipolar pathways parameters
-    for n in range(npathways):
-        if spins==2 and npathways==1:
-            # Special case: use modulation depth notation instead of general pathway amplitude
-            getattr(PathsModel,f'mod').set(lb=0,ub=1,par0=0.01,description=f'Modulation depth',unit='')
-        else:
-            pairwise = ''
-            if spins>2: 
-                pairwise = ' pairwise'   
-            if spins==2 or samespins:
-                getattr(PathsModel,f'lam{labels[n]}').set(lb=0,ub=1,par0=0.01,description=f'Amplitude of{pairwise} pathway #{labels[n]}',unit='')
-            else:
-                for q in range(Q):
-                    getattr(PathsModel,f'lam{labels[n]}_{q+1}').set(lb=0,ub=1,par0=0.01,description=f'Amplitude of{pairwise} pathway #{labels[n]} on interaction #{q+1}',unit='')
-        if parametrization=='reftimes':
-            if experiment is None:
-                theoretical_reftime = 0
-                reftime_variability = 20
-            else:
-                theoretical_reftime = experiment.reftimes(*experiment.delays)[n]
-                reftime_variability = 3*experiment.pulselength
-            if spins==2 and npathways==1:
-                # Special case: use reftime notation
-                getattr(PathsModel,f'reftime').set(description=f'Refocusing time',unit='μs',
-                    par0 = theoretical_reftime,  
-                    lb   = theoretical_reftime - reftime_variability,
-                    ub   = theoretical_reftime + reftime_variability)
-            else: 
-                # Special case: use reftime notation
-                getattr(PathsModel,f'reftime{labels[n]}').set(description=f'Refocusing time of{pairwise} pathway #{labels[n]}',unit='μs',
-                    par0 = theoretical_reftime,
-                    lb   = theoretical_reftime - reftime_variability,
-                    ub   = theoretical_reftime + reftime_variability)
-    if parametrization=='delays': 
-        experimental_delays = experiment.delays
-        delay_variability = 5*experiment.pulselength
-        for n,delay in enumerate(pulsedelay_names):
-            getattr(PathsModel,delay).set(description=f'Pulse delay {delay} of the experiment',unit='μs',
-                par0 = experimental_delays[n],
-                lb   = experimental_delays[n] - delay_variability,
-                ub   = experimental_delays[n] + delay_variability)
-    elif parametrization=='shift':
-        variability = 5*experiment.pulselength
-        getattr(PathsModel,'tshift').set(par0=0,lb=-variability,ub=variability,description=f'Variability of experimental pulse delays',unit='μs')
-    if spins>2: 
-        getattr(PathsModel,f'lamu').set(lb=0,ub=1,par0=0.5,description='Amplitude of unmodulated pairwise pathway',unit='')
-        
+    PathsModel = _populate_dipolar_pathways_parameters(
+        PathsModel,npathways,spins=spins,samespins=samespins,experiment=experiment,parametrization=parametrization,Q=Q)
+            
     # Construct the signature of the dipolar signal model function
     signature = []
     parameters,linearparam = [],[]
@@ -584,6 +637,16 @@ def Vmultispin_nonlinear_fcn(*nonlin):
                     max(t)+max(reftimes) + 3*experiment.pulselength + dt,
                     dt)
 
+    # 
+    if Bmodel is not None:
+        DipolarSignal.Bmodel = dipolarbackgroundmodel(experiment=experiment, basis=Bmodel, parametrization=parametrization, spins=spins, samespins=samespins, npathways=npathways)
+    
+    if Pmodel is not None:
+        DipolarSignal.Pmodel = Pmodel
+    DipolarSignal.t = t
+    if spins==2:
+        DipolarSignal.r = r
+
     return DipolarSignal
 #===============================================================================
 
diff --git a/deerlab/fit.py b/deerlab/fit.py
index 71d262310..c1c27da30 100644
--- a/deerlab/fit.py
+++ b/deerlab/fit.py
@@ -548,9 +548,15 @@ def _scale(x):
     if len(noiselvl)==1: 
         noiselvl = noiselvl[0]
     
+
     # Generate FitResult object from all the dictionaries
     fitresult = FitResult({**FitResult_param_,**FitResult_paramuq_, **FitResult_dict,'penweights':penweights,'noiselvl':noiselvl,'paramlist':_paramlist, '_param_idx':param_idx}) 
 
+    # Add background if it is part of the model
+    if hasattr(model,'Bmodel'):
+        fitresult['bg'] = fitresult.evaluate(model.Bmodel,model.t)
+        fitresult['bgUncert'] = fitresult.propagate(model.Bmodel,model.t, samples=bootstrap)
+
     fitresult._summary = _print_fitresults(fitresult,model)
 
     return fitresult
diff --git a/docsrc/source/reference.rst b/docsrc/source/reference.rst
index 652289e6a..ebb9df56c 100644
--- a/docsrc/source/reference.rst
+++ b/docsrc/source/reference.rst
@@ -60,6 +60,7 @@ Reference Index
     deerload
     dipolarkernel 
     dipolarbackground
+    dipolarbackgroundmodel
     fftspec
     distancerange
 
diff --git a/examples/basic/ex_bootstrapping.py b/examples/basic/ex_bootstrapping.py
index 1e227ad1c..41565f1be 100644
--- a/examples/basic/ex_bootstrapping.py
+++ b/examples/basic/ex_bootstrapping.py
@@ -47,7 +47,8 @@
 r = np.linspace(2,5,100) # nm
 
 # Construct the model
-Vmodel = dl.dipolarmodel(t,r, experiment = dl.ex_4pdeer(tau1,tau2, pathways=[1]))
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t,r, experiment = experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp,bootstrap=20,bootcores=4)
@@ -69,9 +70,9 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,conc,reftime: results.P_scale*(1-mod)*dl.bg_hom3d(t-reftime,conc,mod)
-Bfit = results.evaluate(Bfcn)
-Bci = results.propagate(Bfcn).ci(95)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfit = results.P_scale*results.evaluate(Bfcn,t)
+Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'
diff --git a/examples/basic/ex_fitting_4pdeer.py b/examples/basic/ex_fitting_4pdeer.py
index 38368f5eb..ea29e35bd 100644
--- a/examples/basic/ex_fitting_4pdeer.py
+++ b/examples/basic/ex_fitting_4pdeer.py
@@ -35,7 +35,8 @@
 r = np.arange(2.5,5,0.01) # nm
 
 # Construct the model
-Vmodel = dl.dipolarmodel(t,r, experiment = dl.ex_4pdeer(tau1,tau2, pathways=[1]))
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t,r, experiment = experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp)
@@ -54,9 +55,9 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,conc,reftime: results.P_scale*(1-mod)*dl.bg_hom3d(t-reftime,conc,mod)
-Bfit = results.evaluate(Bfcn)
-Bci = results.propagate(Bfcn).ci(95)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfit = results.P_scale*results.evaluate(Bfcn,t)
+Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'
diff --git a/examples/basic/ex_fitting_4pdeer_compactness.py b/examples/basic/ex_fitting_4pdeer_compactness.py
index b94f782bb..863eed1f4 100644
--- a/examples/basic/ex_fitting_4pdeer_compactness.py
+++ b/examples/basic/ex_fitting_4pdeer_compactness.py
@@ -38,7 +38,8 @@
 r = np.arange(2,5,0.025) # nm
 
 # Construct the model
-Vmodel = dl.dipolarmodel(t, r, experiment = dl.ex_4pdeer(tau1,tau2, pathways=[1]))
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t, r, experiment = experimentInfo)
 compactness = dl.dipolarpenalty(Pmodel=None, r=r, type='compactness')
 
 # Fit the model to the data
@@ -58,9 +59,9 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,conc,reftime: results.P_scale*(1-mod)*dl.bg_hom3d(t-reftime,conc,mod)
-Bfit = results.evaluate(Bfcn)
-Bci = results.propagate(Bfcn).ci(95)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfit = results.P_scale*results.evaluate(Bfcn,t)
+Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'
diff --git a/examples/basic/ex_fitting_4pdeer_gauss.py b/examples/basic/ex_fitting_4pdeer_gauss.py
index e110613e9..cc4536254 100644
--- a/examples/basic/ex_fitting_4pdeer_gauss.py
+++ b/examples/basic/ex_fitting_4pdeer_gauss.py
@@ -37,7 +37,8 @@
 
 # Construct the model
 Pmodel= dl.dd_gauss2
-Vmodel = dl.dipolarmodel(t,r,Pmodel, experiment=dl.ex_4pdeer(tau1,tau2, pathways=[1]))
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t,r,Pmodel, experiment=experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp,reg=False)
@@ -59,9 +60,9 @@
 Pci50 = Puncert.ci(50)/scale
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,conc,reftime: scale*(1-mod)*dl.bg_hom3d(t-reftime,conc,mod)
-Bfit = results.evaluate(Bfcn)
-Bci = results.propagate(Bfcn).ci(95)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfit = scale*results.evaluate(Bfcn,t)
+Bci = scale*results.propagate(Bfcn,t).ci(95)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'
diff --git a/examples/basic/ex_fitting_5pdeer.py b/examples/basic/ex_fitting_5pdeer.py
index 432d31292..39a056e06 100644
--- a/examples/basic/ex_fitting_5pdeer.py
+++ b/examples/basic/ex_fitting_5pdeer.py
@@ -67,9 +67,9 @@
 Pfit =  Pfit
 
 # Extract the unmodulated contribution
-Bfcn = lambda lam1,lam5,reftime1,reftime5,conc: results.P_scale*(1-lam1-lam5)*dl.bg_hom3d(t-reftime1,conc,lam1)*dl.bg_hom3d(t-reftime5,conc,lam5)
-Bfit = results.evaluate(Bfcn)
-Bci = results.propagate(Bfcn).ci(95)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfit = results.P_scale*results.evaluate(Bfcn,t)
+Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'
diff --git a/examples/basic/ex_fitting_ridme.py b/examples/basic/ex_fitting_ridme.py
index f58ccfa6b..357f1d402 100644
--- a/examples/basic/ex_fitting_ridme.py
+++ b/examples/basic/ex_fitting_ridme.py
@@ -36,8 +36,8 @@
 r = np.linspace(1.5,6,50) # nm
 
 # Construct the model
-experimentmodel = dl.ex_ridme(tau1,tau2, pathways=[1])
-Vmodel = dl.dipolarmodel(t,r,Bmodel=dl.bg_strexp, experiment =experimentmodel)
+experimentInfo = dl.ex_ridme(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t,r,Bmodel=dl.bg_strexp, experiment =experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp)
@@ -56,9 +56,9 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,decay,stretch,reftime: results.P_scale*(1-mod)*dl.bg_strexp(t-reftime,decay,stretch)
-Bfit = results.evaluate(Bfcn)
-Bci = results.propagate(Bfcn).ci(95)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo,dl.bg_strexp)
+Bfit = results.P_scale*results.evaluate(Bfcn,t)
+Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'
diff --git a/examples/basic/ex_restraints_4pdeer.py b/examples/basic/ex_restraints_4pdeer.py
index d3c71a103..f24764b85 100644
--- a/examples/basic/ex_restraints_4pdeer.py
+++ b/examples/basic/ex_restraints_4pdeer.py
@@ -42,7 +42,8 @@
 r = np.arange(2,6,0.05) # nm
 
 # Construct dipolar model
-Vmodel = dl.dipolarmodel(t,r, experiment=dl.ex_4pdeer(tau1,tau2, pathways=[1]))
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t,r, experiment=experimentInfo)
 
 # Fit the model to the data
 fit = dl.fit(Vmodel,Vexp)
diff --git a/examples/intermediate/ex_background_model.py b/examples/intermediate/ex_background_model.py
new file mode 100644
index 000000000..9d47cb60e
--- /dev/null
+++ b/examples/intermediate/ex_background_model.py
@@ -0,0 +1,75 @@
+# %% [markdown]
+""" 
+Generating a background model for a 4-pulse DEER experiment
+-------------------------------------------------------------------------
+
+Since DeerLab v1.2, the background model can be generated for a given experiment using the :func:`deerlab.dipolarbackgroundmodel` function, as is normally automatically calculated in the fit function as `FitResult.bg` and  `FitResult.bgUncert`.
+
+However, it can still be useful to generate the background model separately, and so this example shows how to do that for a 4-pulse DEER experiment.
+
+First we will fit a simple 4-pulse DEER dataset with multiple-pathways. 
+""" 
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import deerlab as dl
+# %%
+# File location
+path = '../data/'
+file = 'example_4pdeer_1.DTA'
+
+# Experimental parameters
+tau1 = 0.3      # First inter-pulse delay, μs
+tau2 = 4.0      # Second inter-pulse delay, μs
+tmin = 0.1      # Start time, μs
+
+# Load the experimental data
+t,Vexp = dl.deerload(path + file)
+
+# Pre-processing
+Vexp = dl.correctphase(Vexp) # Phase correction
+Vexp = Vexp/np.max(Vexp)     # Rescaling (aesthetic)
+t = t - t[0]                  # Account for zerotime
+t = t + tmin    
+# Distance vector
+r = np.arange(2.5,5,0.01) # nm
+
+# Construct the model
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1,2,3])
+Vmodel = dl.dipolarmodel(t,r, experiment = experimentInfo, Bmodel=dl.bg_hom3d)
+
+# Fit the model to the data
+results = dl.fit(Vmodel,Vexp)
+
+# Print results summary
+print(results)
+
+
+# %%
+"""
+Now using DeerLab > v1.2, the background model can be directly plotted from the fit results.
+"""
+plt.figure(figsize=[6,7])
+violet = '#4550e6'
+plt.plot(t,Vexp,'.',color='grey',label='Data')
+plt.plot(t,results.model,linewidth=3,color=violet,label='Fit')
+plt.plot(t,results.bg,'--',linewidth=3,color=violet,label='Unmodulated contribution')
+
+
+# %%
+"""
+Alternatively we could generate the background model separately using the `dl.dipolarbackgroundmodel` function, which takes the experimental information as input.
+"""
+
+bg_model = dl.dipolarbackgroundmodel(experimentInfo, basis = dl.bg_hom3d)
+bg = results.evaluate(bg_model,t)
+bgUncert = results.propagate(bg_model,t)
+
+"""
+In DeerLab <1.2, the multi-pathway background model needs to be generated manually from a function.
+"""
+
+bf_func = lambda lam1,lam2,lam3,reftime1,reftime2,reftime3,conc: results.P_scale*(1-lam1-lam2-lam3)*dl.bg_hom3d(t-reftime1,conc,lam1)*dl.bg_hom3d(t-reftime2,conc,lam2)*dl.bg_hom3d(t-reftime3,conc,lam3)
+bf = results.evaluate(bf_func, t)
+bfUncert = results.propagate(bf_func, t)
diff --git a/examples/intermediate/ex_compactness_with_without.py b/examples/intermediate/ex_compactness_with_without.py
index 0b836934d..39a00d72d 100644
--- a/examples/intermediate/ex_compactness_with_without.py
+++ b/examples/intermediate/ex_compactness_with_without.py
@@ -40,7 +40,8 @@
 r = np.arange(1.5,7,0.05) # nm
 
 # Construct the model
-Vmodel = dl.dipolarmodel(t,r, experiment=dl.ex_4pdeer(tau1,tau2, pathways=[1]))
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t,r, experiment=experimentInfo)
 compactness = dl.dipolarpenalty(Pmodel=None,r=r,type='compactness')
 
 # Fit the model to the data with compactness criterion
@@ -69,9 +70,10 @@
     Pci50 = results.PUncert.ci(50)
 
     # Extract the unmodulated contribution
-    Bfcn = lambda mod,conc,reftime: results.P_scale*(1-mod)*dl.bg_hom3d(t-reftime,conc,mod)
-    Bfit = results.evaluate(Bfcn)
-    Bci = results.propagate(Bfcn).ci(95)
+    Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+    Bfit = results.P_scale*results.evaluate(Bfcn,t)
+    Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
+
 
     plt.subplot(2,2,n+1)
     # Plot experimental and fitted data
diff --git a/examples/intermediate/ex_crossing_echoes_masking.py b/examples/intermediate/ex_crossing_echoes_masking.py
index dd60009e2..2defcb860 100644
--- a/examples/intermediate/ex_crossing_echoes_masking.py
+++ b/examples/intermediate/ex_crossing_echoes_masking.py
@@ -110,8 +110,8 @@
 r = np.arange(2,6,0.05) # nm
 
 # Construct the dipolar signal model 
-experiment = dl.ex_4pdeer(tau1,tau2,pathways=[1,2,3])
-Vmodel = dl.dipolarmodel(t,r,experiment=experiment)
+experimentInfo = dl.ex_4pdeer(tau1,tau2,pathways=[1,2,3])
+Vmodel = dl.dipolarmodel(t,r,experiment=experimentInfo)
 
 # Analyze the data while ignoring the crossing echoes
 results = dl.fit(Vmodel,Vexp, mask=mask, noiselvl=noiselevel)
@@ -130,13 +130,8 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-def Vunmodfcn(lam1,lam2,lam3,reftime1,reftime2,reftime3,conc):
-    Lam0 = results.P_scale*(1-lam1-lam2-lam3)
-    Vunmod = Lam0*dl.bg_hom3d(t-reftime1,conc,lam1)
-    Vunmod *= dl.bg_hom3d(t-reftime2,conc,lam2)
-    Vunmod *= dl.bg_hom3d(t-reftime3,conc,lam3)
-    return Vunmod  
-Bfit = results.evaluate(Vunmodfcn)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfit = results.P_scale*results.evaluate(Bfcn,t)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'
diff --git a/examples/intermediate/ex_fitting_4pdeer_pathways.py b/examples/intermediate/ex_fitting_4pdeer_pathways.py
index cef4d058f..b440679f7 100644
--- a/examples/intermediate/ex_fitting_4pdeer_pathways.py
+++ b/examples/intermediate/ex_fitting_4pdeer_pathways.py
@@ -34,8 +34,8 @@
 r = np.arange(2.5,5.5,0.025) # nm
 
 # Construct the model
-experiment = dl.ex_4pdeer(tau1,tau2, pathways=[1,2,3])
-Vmodel = dl.dipolarmodel(t,r,experiment=experiment)
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1,2,3])
+Vmodel = dl.dipolarmodel(t,r,experiment=experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp)
@@ -70,7 +70,7 @@
 lams = [results.lam1, results.lam2, results.lam3]
 reftimes = [results.reftime1, results.reftime2, results.reftime3]
 colors= ['tab:blue',green, red] 
-Vinter = results.P_scale*(1-np.sum(lams))*np.prod([dl.bg_hom3d(t-reftime,results.conc,lam) for lam,reftime in zip(lams,reftimes)],axis=0)
+Vinter = results.P_scale*results.evaluate(dl.dipolarbackgroundmodel(experimentInfo),t)
 for n,(lam,reftime,color) in enumerate(zip(lams,reftimes,colors)):
     Vpath = (1-np.sum(lams) + lam*dl.dipolarkernel(t-reftime,r)@Pfit)*Vinter
     plt.plot(t,Vpath,linewidth=3,label=f'Pathway #{n+1}',color=color)
diff --git a/examples/intermediate/ex_fitting_5pdeer_pathways.py b/examples/intermediate/ex_fitting_5pdeer_pathways.py
index 516e50dce..b1d776ebb 100644
--- a/examples/intermediate/ex_fitting_5pdeer_pathways.py
+++ b/examples/intermediate/ex_fitting_5pdeer_pathways.py
@@ -34,8 +34,8 @@
 r = np.arange(2.5,5.5,0.025) # nm
 
 # Construct the model
-experiment = dl.ex_rev5pdeer(tau1,tau2,tau3, pathways=[1,2,3,5])
-Vmodel = dl.dipolarmodel(t,r,experiment=experiment)
+experimentInfo = dl.ex_rev5pdeer(tau1,tau2,tau3, pathways=[1,2,3,5])
+Vmodel = dl.dipolarmodel(t,r,experiment=experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp)
@@ -72,7 +72,7 @@
 lams = [results.lam1, results.lam2, results.lam3, results.lam5]
 reftimes = [results.reftime1, results.reftime2, results.reftime3, results.reftime5]
 colors= ['tab:blue','tab:orange', red, green] 
-Vinter = results.P_scale*(1-np.sum(lams))*np.prod([dl.bg_hom3d(t-reftime,results.conc,lam) for lam,reftime in zip(lams,reftimes)],axis=0)
+Vinter = results.P_scale*results.evaluate(dl.dipolarbackgroundmodel(experimentInfo),t)
 for (lam,reftime,color,label) in zip(lams,reftimes,colors,labels):
     Vpath = (1-np.sum(lams) + lam*dl.dipolarkernel(t-reftime,r)@Pfit)*Vinter
     plt.plot(t,Vpath,linewidth=3,label=f'Pathway #{label}',color=color)
diff --git a/examples/intermediate/ex_fitting_dqc_pathways.py b/examples/intermediate/ex_fitting_dqc_pathways.py
index 188ba52b5..3e349be0f 100644
--- a/examples/intermediate/ex_fitting_dqc_pathways.py
+++ b/examples/intermediate/ex_fitting_dqc_pathways.py
@@ -37,8 +37,8 @@
 
 # Construct the model
 r = np.arange(2.5,4,0.01) # nm
-experiment = dl.ex_dqc(tau1,tau2,tau3,pathways=[1,2,3])
-Vmodel =  dl.dipolarmodel(t,r,experiment=experiment)
+experimentInfo = dl.ex_dqc(tau1,tau2,tau3,pathways=[1,2,3])
+Vmodel =  dl.dipolarmodel(t,r,experiment=experimentInfo)
 
 # The amplitudes of the second and third pathways must be equal
 Vmodel = dl.link(Vmodel,lam23=['lam2','lam3'])
@@ -56,7 +56,7 @@
 
 # Plot the full detectable range
 tfull = np.arange(-2*tau1,2*tau2-4*tau3,0.008)
-Vmodelext =  dl.dipolarmodel(tfull,r,experiment=experiment)
+Vmodelext =  dl.dipolarmodel(tfull,r,experiment=experimentInfo)
 Vmodelext = dl.link(Vmodelext,lam23=['lam2','lam3'])
 
 # Extract results
@@ -76,7 +76,7 @@
 
 # Plot the individual pathway contributions
 plt.subplot(223)
-Vinter = results.P_scale*(1-np.sum(lams))*np.prod([dl.bg_hom3d(tfull-reftime,results.conc,lam) for lam,reftime in zip(lams,reftimes)],axis=0)
+Vinter = results.P_scale*results.evaluate(dl.dipolarbackgroundmodel(experimentInfo),t)
 for n,(lam,reftime,color) in enumerate(zip(lams,reftimes,colors)):
     Vpath = (1-np.sum(lams) + lam*dl.dipolarkernel(tfull-reftime,r)@Pfit)*Vinter
     plt.plot(tfull,Vpath,label=f'Pathway #{n+1}',color=color)
diff --git a/examples/intermediate/ex_fitting_sparse_4pdeer.py b/examples/intermediate/ex_fitting_sparse_4pdeer.py
index 68adaaa05..687088b38 100644
--- a/examples/intermediate/ex_fitting_sparse_4pdeer.py
+++ b/examples/intermediate/ex_fitting_sparse_4pdeer.py
@@ -41,7 +41,8 @@
 r = np.arange(2,10,0.05) # nm
 
 # Construct the model with the sparse sampled time vector
-Vmodel = dl.dipolarmodel(t,r, experiment = dl.ex_4pdeer(tau1,tau2, pathways=[1]))
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t,r, experiment = experimentInfo)
 compactness = dl.dipolarpenalty(Pmodel=None, r=r, type='compactness')
 
 # Fit the model to the data
@@ -65,9 +66,10 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,conc,reftime: results.P_scale*(1-mod)*dl.bg_hom3d(tuniform-reftime,conc,mod)
-Bfit = results.evaluate(Bfcn)
-Bci = results.propagate(Bfcn).ci(95)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfit = results.P_scale*results.evaluate(Bfcn,t)
+Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
+
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'
diff --git a/test/test_dipolarbackground.py b/test/test_dipolarbackground.py
index f25f7c6e7..73a2420d4 100644
--- a/test/test_dipolarbackground.py
+++ b/test/test_dipolarbackground.py
@@ -1,9 +1,12 @@
 
 import numpy as np
-from deerlab.bg_models import _hom3d, bg_exp, _homfractal
-from deerlab.dipolarbackground import dipolarbackground
+from itertools import permutations
+from math import factorial
+from deerlab.bg_models import _hom3d, bg_exp, _homfractal, bg_hom3d
+from deerlab.dipolarbackground import dipolarbackground, dipolarbackgroundmodel
+from deerlab.dipolarmodel import ex_3pdeer, ex_4pdeer
 from deerlab.utils import assert_docstring
-import pytest 
+import pytest
 
 # Fixtures 
 # -----------------------------------------------------------------------------------
@@ -132,3 +135,171 @@ def test_docstring():
     "Check that the docstring includes all variables and keywords."
     assert_docstring(dipolarbackground)
 # ======================================================================
+
+
+# ==================================================================================
+# Tests for dipolarbackgroundmodel
+# ==================================================================================
+
+tau = 0.5
+tau1, tau2 = 0.5, 3.0
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_twospin_single_pathway():
+    "Test dipolarbackgroundmodel for a two-spin system with a single pathway"
+    mod = 0.4
+    tref = tau
+    Bref = (1-mod) * _hom3d(t-tref, conc, mod)
+    experiment = ex_3pdeer(tau, pathways=[1])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d)
+    Bval = Bmodel(t, conc, mod, tref)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_twospin_single_pathway_no_experiment():
+    "Test dipolarbackgroundmodel for a two-spin system with a single pathway with no experiment"
+    mod = 0.4
+    tref = tau
+    Bref = (1-mod) * _hom3d(t-tref, conc, mod)
+    Bmodel = dipolarbackgroundmodel(None, bg_hom3d)
+    Bval = Bmodel(t, conc, mod, tref)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_twospin_multi_pathway():
+    "Test dipolarbackgroundmodel for a two-spin system with multiple pathways"
+    lam1, lam2 = 0.3, 0.1
+    tref1, tref2 = tau1, tau1+tau2
+    Bref = (1-lam1-lam2) * _hom3d(t-tref1, conc, lam1) * _hom3d(t-tref2, conc, lam2)
+    experiment = ex_4pdeer(tau1, tau2, pathways=[1,2])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d)
+    Bval = Bmodel(t, conc, lam1, lam2, tref1, tref2)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_twospin_delays():
+    "Test dipolarbackgroundmodel for a two-spin system with delays parametrization"
+    lam1, lam2 = 0.3, 0.1
+    tref1, tref2 = tau1, tau1+tau2
+    Bref = (1-lam1-lam2) * _hom3d(t-tref1, conc, lam1) * _hom3d(t-tref2, conc, lam2)
+    experiment = ex_4pdeer(tau1, tau2, pathways=[1,2])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d, parametrization='delays')
+    Bval = Bmodel(t, conc, tau1, tau2, lam1, lam2)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_threespin_single_pathway():
+    "Test dipolarbackgroundmodel for a three-spin system with a single pathway"
+    lam, lamu = 0.4, 0.8
+    Q = 3
+    # Construct reference pathways manually
+    Λ0 = 1.0
+    paths_ref = []
+    for perm in set(permutations([0]+[None]*(Q-1))):
+        λ2k = factorial(Q-1)*lam*lamu**(Q-1)
+        tref_path = tuple([tau if n==0 else None for n in perm])
+        δs_path = tuple([1 if n==0 else 0 for n in perm])
+        paths_ref.append({'amp': λ2k, 'reftime': tref_path, 'harmonic': δs_path})
+        Λ0 -= λ2k
+    paths_ref.append({'amp': Λ0})
+    Bref = dipolarbackground(t, paths_ref, lambda td, l: _hom3d(td, conc, l))
+    experiment = ex_3pdeer(tau, pathways=[1])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d, spins=3)
+    Bval = Bmodel(t, conc, lam, tau, lamu)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_threespin_multi_pathway():
+    "Test dipolarbackgroundmodel for a three-spin system with multiple pathways"
+    lam1, lam2, lamu = 0.3, 0.1, 0.7
+    Q = 3
+    harmonics_exp = [1, 1]
+    trefs_exp = [tau1, tau1+tau2]
+
+    # Construct reference pathways manually
+    Λ0 = 1.0
+    paths_ref, threespin_ref = [], []
+    for idx, (tref_idx, δ_idx) in enumerate(zip(trefs_exp, harmonics_exp)):
+        lam_idx = [lam1, lam2][idx]
+        for perm in set(permutations([0]+[None]*(Q-1))):
+            q = int(np.where(np.array(perm)==0)[0][0])
+            λp = lam_idx
+            λ2k = factorial(Q-1)*λp*lamu**(Q-1)
+            paths_ref.append({'amp': λ2k,
+                               'reftime': tuple([tref_idx if n==0 else None for n in perm]),
+                               'harmonic': tuple([δ_idx if n==0 else 0 for n in perm])})
+            Λ0 -= λ2k
+            for idx2 in range(idx+1, 2):
+                tref_idx2 = trefs_exp[idx2]
+                δ_idx2 = harmonics_exp[idx2]
+                for perm2 in set(permutations([0,1]+[None]*(Q-2))):
+                    perm2_arr = np.array(perm2)
+                    q2 = int(np.where(perm2_arr==1)[0][0])
+                    λm = lam_idx
+                    λ3k = 2*factorial(Q-2)*λp*λm*lamu**(Q-2)
+                    threespin_ref.append({'amp': λ3k,
+                                          'reftime': tuple([tref_idx if n==0 else tref_idx2 if n==1 else None for n in perm2]),
+                                          'harmonic': tuple([δ_idx if n==0 else δ_idx2 if n==1 else 0 for n in perm2])})
+                    Λ0 -= λ3k
+    paths_ref.append({'amp': Λ0})
+    Bref = dipolarbackground(t, paths_ref+threespin_ref, lambda td, l: _hom3d(td, conc, l))
+
+    experiment = ex_4pdeer(tau1, tau2, pathways=[1,2])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d, spins=3)
+    Bval = Bmodel(t, conc, lam1, lam2, tau1, tau1+tau2, lamu)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_threespin_nosamespins():
+    "Test dipolarbackgroundmodel for a three-spin system with samespins=False"
+    lam_q = [0.3, 0.4, 0.2]   # different amplitudes per interaction
+    lamu = 0.8
+    Q = 3
+    # Construct reference pathways with per-interaction amplitudes
+    Λ0 = 1.0
+    paths_ref = []
+    for perm in set(permutations([0]+[None]*(Q-1))):
+        q = int(np.where(np.array(perm)==0)[0][0])
+        λ2k = factorial(Q-1)*lam_q[q]*lamu**(Q-1)
+        tref_path = tuple([tau if n==0 else None for n in perm])
+        δs_path = tuple([1 if n==0 else 0 for n in perm])
+        paths_ref.append({'amp': λ2k, 'reftime': tref_path, 'harmonic': δs_path})
+        Λ0 -= λ2k
+    paths_ref.append({'amp': Λ0})
+    Bref = dipolarbackground(t, paths_ref, lambda td, l: _hom3d(td, conc, l))
+
+    experiment = ex_3pdeer(tau, pathways=[1])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d, spins=3, samespins=False)
+    Bval = Bmodel(t, conc, *lam_q, tau, lamu)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_threespin_delays():
+    "Test dipolarbackgroundmodel for a three-spin system with delays parametrization"
+    lam, lamu = 0.4, 0.8
+    Q = 3
+    Λ0 = 1.0
+    paths_ref = []
+    for perm in set(permutations([0]+[None]*(Q-1))):
+        λ2k = factorial(Q-1)*lam*lamu**(Q-1)
+        tref_path = tuple([tau if n==0 else None for n in perm])
+        δs_path = tuple([1 if n==0 else 0 for n in perm])
+        paths_ref.append({'amp': λ2k, 'reftime': tref_path, 'harmonic': δs_path})
+        Λ0 -= λ2k
+    paths_ref.append({'amp': Λ0})
+    Bref = dipolarbackground(t, paths_ref, lambda td, l: _hom3d(td, conc, l))
+
+    experiment = ex_3pdeer(tau, pathways=[1])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d, spins=3, parametrization='delays')
+    # signature: t(const), conc, delay1, lam1, lamu
+    Bval = Bmodel(t, conc, tau, lam, lamu)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
diff --git a/test/test_dipolarmodel.py b/test/test_dipolarmodel.py
index 7d655dd80..84fbb2d9c 100644
--- a/test/test_dipolarmodel.py
+++ b/test/test_dipolarmodel.py
@@ -4,7 +4,7 @@
 from deerlab.model import Model
 from deerlab.fit import fit
 from deerlab.dipolarmodel import ExperimentInfo,dipolarpenalty, dipolarmodel, ex_4pdeer, ex_3pdeer,ex_fwd5pdeer, ex_rev5pdeer, ex_sifter, ex_ridme, ex_dqc
-from deerlab import dd_gauss,dd_gauss2,bg_hom3d,bg_exp
+from deerlab import dd_gauss,dd_gauss2,bg_hom3d,bg_exp,UQResult
 from deerlab.utils import assert_docstring
 from pytest import fixture 
 
@@ -153,6 +153,7 @@ def test_phenomenological_Bmodel(V1path_phenoB):
     Vsim = Vmodel(mod=0.3,reftime=0.0,mean=3,std=0.2,decay=0.1,scale=1e5)
 
     assert np.allclose(Vsim,V1path_phenoB)
+   
 # ======================================================================
 
 # ======================================================================
@@ -164,6 +165,7 @@ def test_no_Bmodel(V1path_noB):
     Vsim = Vmodel(mod=0.3,reftime=0.0,mean=3,std=0.2,scale=1e5)
 
     assert np.allclose(Vsim,V1path_noB)
+    assert not hasattr(Vmodel,'Bmodel')
 # ======================================================================
 
 # ======================================================================
@@ -177,6 +179,19 @@ def test_model_1pathways(V1path):
     assert np.allclose(Vsim,V1path)
 # ======================================================================
 
+# ======================================================================
+def test_model_1pathways_background(V1path): 
+    "Check that the model with one dipolar pathway is correct"
+
+    Vmodel = dipolarmodel(t,r,dd_gauss,bg_hom3d,npathways=1)
+    
+    Vsim = Vmodel(mod=0.3,reftime=0.0,conc=50,mean=3,std=0.2,scale=1e5)
+
+    assert np.allclose(Vsim,V1path)
+    assert hasattr(Vmodel,'Bmodel') and Vmodel.Bmodel is not None
+    assert isinstance(Vmodel.Bmodel,Model)
+# ======================================================================
+
 # ======================================================================
 def test_model_2pathways(V2path): 
     "Check that the model with two dipolar pathways is correct"
@@ -212,6 +227,8 @@ def test_fit_1pathways(V1path):
     result = fit(Vmodel,V1path,ftol=1e-4)
 
     assert np.allclose(result.model,V1path)
+    assert hasattr(result,'bg') and result.bg is not None and isinstance(result.bg,np.ndarray)
+    assert hasattr(result,'bgUncert') and result.bgUncert is not None and isinstance(result.bgUncert,UQResult)
 # ======================================================================
 
 # ======================================================================
@@ -464,7 +481,8 @@ def test_ex_fwd5pdeer_fit(Vfwd5pdeer,Pr):
     result = fit(Vmodel,Vfwd5pdeer,ftol=1e-4)
 
     assert np.allclose(Vfwd5pdeer,result.model,rtol=1e-3) and ovl(result.P/1e5,Pr)>0.975
-# ======================================================================
+# ========
+# ==============================================================
 
 # ======================================================================
 def test_ex_fwd5pdeer_pulselength():