diff --git a/dingo/core/dataset.py b/dingo/core/dataset.py
index 49f9718b8..62893e765 100644
--- a/dingo/core/dataset.py
+++ b/dingo/core/dataset.py
@@ -4,6 +4,7 @@
 import numpy as np
 import pandas as pd
 
+from dingo.core.utils.logging_utils import logger
 from dingo.core.utils.misc import get_version
 
 
@@ -136,7 +137,7 @@ def __init__(
             self.from_dictionary(dictionary)
 
     def to_file(self, file_name: str, mode: str = "w"):
-        print("Saving dataset to " + str(file_name))
+        logger.info("Saving dataset to " + str(file_name))
         save_dict = {
             k: v
             for k, v in vars(self).items()
@@ -150,9 +151,9 @@ def to_file(self, file_name: str, mode: str = "w"):
                 f.attrs["dataset_type"] = self.dataset_type
 
     def from_file(self, file_name: str):
-        print(f"Loading dataset from {str(file_name)}.")
+        logger.info(f"Loading dataset from {str(file_name)}.")
         if self._leave_on_disk_keys:
-            print(f"Omitting data keys {self._leave_on_disk_keys}.")
+            logger.info(f"Omitting data keys {self._leave_on_disk_keys}.")
 
         with h5py.File(file_name, "r") as f:
             loaded_dict = recursive_hdf5_load(
diff --git a/dingo/core/posterior_models/base_model.py b/dingo/core/posterior_models/base_model.py
index 68508d199..5df191c3d 100755
--- a/dingo/core/posterior_models/base_model.py
+++ b/dingo/core/posterior_models/base_model.py
@@ -10,6 +10,7 @@
 
 import torch
 import dingo.core.utils as utils
+from dingo.core.utils.logging_utils import logger
 from torch.utils.data import Dataset
 import time
 import numpy as np
@@ -198,7 +199,7 @@ def network_to_device(self, device):
         #     raise NotImplementedError('This needs testing!')
         #     # dim = 0 [512, ...] -> [256, ...], [256, ...] on 2 GPUs
         #     self.network = torch.nn.DataParallel(self.network)
-        print(f"Putting posterior model to device {self.device}.")
+        logger.info(f"Putting posterior model to device {self.device}.")
         self.network.to(self.device)
 
     def initialize_optimizer_and_scheduler(self):
@@ -394,7 +395,7 @@ def train(
 
         if test_only:
             test_loss = test_epoch(self, test_loader)
-            print(f"test loss: {test_loss:.3f}")
+            logger.info(f"test loss: {test_loss:.3f}")
 
         else:
             while not runtime_limits.limits_exceeded(self.epoch):
@@ -403,24 +404,24 @@ def train(
                 # Training
                 lr = utils.get_lr(self.optimizer)
                 with threadpool_limits(limits=1, user_api="blas"):
-                    print(f"\nStart training epoch {self.epoch} with lr {lr}")
+                    logger.info(f"Start training epoch {self.epoch} with lr {lr}")
                     time_start = time.time()
                     train_loss = train_epoch(self, train_loader)
                     train_time = time.time() - time_start
 
-                    print(
+                    logger.info(
                         "Done. This took {:2.0f}:{:2.0f} min.".format(
                             *divmod(train_time, 60)
                         )
                     )
 
                     # Testing
-                    print(f"Start testing epoch {self.epoch}")
+                    logger.info(f"Start testing epoch {self.epoch}")
                     time_start = time.time()
                     test_loss = test_epoch(self, test_loader)
                     test_time = time.time() - time_start
 
-                    print(
+                    logger.info(
                         "Done. This took {:2.0f}:{:2.0f} min.".format(
                             *divmod(time.time() - time_start, 60)
                         )
@@ -449,7 +450,7 @@ def train(
                             }
                         )
                     except ImportError:
-                        print("wandb not installed. Skipping logging to wandb.")
+                        logger.warning("wandb not installed. Skipping logging to wandb.")
 
                 if early_stopping is not None:
                     # Whether to use train or test loss
@@ -464,9 +465,9 @@ def train(
                             join(train_dir, "best_model.pt"), save_training_info=False
                         )
                     if early_stopping.early_stop:
-                        print("Early stopping")
+                        logger.info("Early stopping")
                         break
-                print(f"Finished training epoch {self.epoch}.\n")
+                logger.info(f"Finished training epoch {self.epoch}.")
 
 
 def train_epoch(pm, dataloader):
diff --git a/dingo/core/result.py b/dingo/core/result.py
index 67017271a..026441de4 100644
--- a/dingo/core/result.py
+++ b/dingo/core/result.py
@@ -16,6 +16,7 @@
 from bilby.core.prior import Constraint, DeltaFunction, PriorDict
 
 from dingo.core.dataset import DingoDataset
+from dingo.core.utils.logging_utils import logger
 from dingo.core.density import train_unconditional_density_estimator
 from dingo.core.utils.misc import recursive_check_dicts_are_equal
 from dingo.core.utils.plotting import get_latex_labels, plot_corner_multi
@@ -167,23 +168,24 @@ def reset_event(self, event_dataset):
         ):
             # This is really just for notification. Actions are only taken if the
             # event metadata differ.
-            print("\nNew event data differ from existing.")
+            logger.warning("New event data differ from existing.")
         self.context = context
 
         if self.event_metadata is not None and self.event_metadata != event_metadata:
-            print("Changes")
-            print("=======")
+            logger.warning("Changes")
+            logger.warning("=======")
             old_minus_new = dict(freeze(self.event_metadata) - freeze(event_metadata))
-            print("Old event metadata:")
+            logger.warning("Old event metadata:")
             for k in sorted(old_minus_new):
-                print(f"  {k}:  {self.event_metadata[k]}")
+                logger.warning(f"  {k}:  {self.event_metadata[k]}")
 
             new_minus_old = dict(freeze(event_metadata) - freeze(self.event_metadata))
-            print("New event metadata:")
+            logger.warning("New event metadata:")
             if self.importance_sampling_metadata.get("updates") is None:
                 self.importance_sampling_metadata["updates"] = {}
             for k in sorted(new_minus_old):
-                print(f"  {k}:  {event_metadata[k]}")
+                logger.warning(f"  {k}:  {event_metadata[k]}")
+                self.importance_sampling_metadata["updates"][k] = event_metadata[k]
                 self.importance_sampling_metadata["updates"][k] = event_metadata[k]
 
             self._rebuild_domain(verbose=True)
@@ -312,12 +314,12 @@ def importance_sample(self, num_processes: int = 1, **likelihood_kwargs):
         valid_samples = np.isfinite(log_prior + delta_log_prob_target)
         theta = theta.iloc[valid_samples]
 
-        print(f"Calculating {len(theta)} likelihoods.")
+        logger.info(f"Calculating {len(theta)} likelihoods.")
         t0 = time.time()
         log_likelihood = self.likelihood.log_likelihood_multi(
             theta, num_processes=num_processes
         )
-        print(f"Done. This took {time.time() - t0:.2f} seconds.")
+        logger.info(f"Done. This took {time.time() - t0:.2f} seconds.")
 
         self.log_noise_evidence = self.likelihood.log_Zn
         self.samples["log_prior"] = log_prior
@@ -494,7 +496,7 @@ def rejection_sample(
         rng = np.random.default_rng(random_state)
         weights = self.samples["weights"].to_numpy(dtype=float)
 
-        print(
+        logger.info(
             f"Rejection sampling: {self.num_samples} samples, "
             f"ESS = {self.effective_sample_size:.0f} "
             f"(efficiency = {100 * self.sample_efficiency:.1f}%)"
@@ -526,7 +528,7 @@ def rejection_sample(
         unweighted = unweighted.drop(
             columns=["weights", "log_prob", "delta_log_prob_target"], errors="ignore"
         )
-        print(f"Produced {len(unweighted)} unweighted samples.")
+        logger.info(f"Produced {len(unweighted)} unweighted samples.")
         return unweighted
 
     def parameter_subset(self, parameters):
@@ -619,12 +621,12 @@ def print_summary(self):
         Display the number of samples, and (if importance sampling is complete) the log
         evidence and number of effective samples.
         """
-        print("Number of samples:", len(self.samples))
+        logger.info("Number of samples: %d", len(self.samples))
         if self.log_evidence is not None:
-            print(
+            logger.info(
                 f"Log(evidence): {self.log_evidence:.3f} +- {self.log_evidence_std:.3f}"
             )
-            print(
+            logger.info(
                 f"Effective samples {self.n_eff:.1f}: "
                 f"(Sample efficiency = {100 * self.sample_efficiency:.2f}%)"
             )
@@ -815,7 +817,7 @@ def plot_log_probs(self, filename="log_probs.png"):
             plt.tight_layout()
             plt.savefig(filename)
         else:
-            print("Results not importance sampled. Cannot produce log_prob plot.")
+            logger.warning("Results not importance sampled. Cannot produce log_prob plot.")
 
     def plot_weights(self, filename="weights.png"):
         """Make a scatter plot of samples weights vs log proposal."""
@@ -844,7 +846,7 @@ def plot_weights(self, filename="weights.png"):
             plt.tight_layout()
             plt.savefig(filename)
         else:
-            print("Results not importance sampled. Cannot plot weights.")
+            logger.warning("Results not importance sampled. Cannot plot weights.")
 
     def get_all_injection_credible_levels(
         self, keys: list[str] = None, weighted: bool = False
@@ -1020,7 +1022,7 @@ def make_pp_plot(
 
     pvalues = []
     latex_labels = get_latex_labels(results[0].prior)
-    print("Key: KS-test p-value")
+    logger.info("Key: KS-test p-value")
     for ii, key in enumerate(credible_levels):
         pp = np.array(
             [
@@ -1030,7 +1032,7 @@ def make_pp_plot(
         )
         pvalue = scipy.stats.kstest(credible_levels[key], "uniform").pvalue
         pvalues.append(pvalue)
-        print("{}: {}".format(key, pvalue))
+        logger.info("{}: {}".format(key, pvalue))
         label = "{} ({:2.3f})".format(latex_labels[key], pvalue)
         plt.plot(x_values, pp, lines[ii], label=label, **kwargs)
 
@@ -1040,7 +1042,7 @@ def make_pp_plot(
         pvalues=pvalues,
         names=list(credible_levels.keys()),
     )
-    print("Combined p-value: {}".format(pvals.combined_pvalue))
+    logger.info("Combined p-value: {}".format(pvals.combined_pvalue))
 
     if title:
         ax.set_title(
diff --git a/dingo/core/samplers.py b/dingo/core/samplers.py
index efc4fc1fa..ab12ddf17 100644
--- a/dingo/core/samplers.py
+++ b/dingo/core/samplers.py
@@ -11,6 +11,7 @@
 from torchvision.transforms import Compose
 
 from dingo.core.posterior_models import BasePosteriorModel
+from dingo.core.utils.logging_utils import logger
 from dingo.core.result import Result
 from dingo.core.result import DATA_KEYS as RESULT_DATA_KEYS
 from dingo.core.utils import torch_detach_to_cpu, IterationTracker
@@ -157,7 +158,7 @@ def _run_sampler(
             x = [x]
         else:
             if context is not None:
-                print("Unconditional model. Ignoring context.")
+                logger.warning("Unconditional model. Ignoring context.")
             x = []
 
         # For a normalizing flow, we get the log_prob for "free" when sampling,
@@ -205,7 +206,7 @@ def run_sampler(
         """
         self.samples = None
 
-        print(f"Running sampler to generate {num_samples} samples.")
+        logger.info(f"Running sampler to generate {num_samples} samples.")
         t0 = time.time()
         if not self.unconditional_model:
             if self.context is None:
@@ -229,7 +230,7 @@ def run_sampler(
         # correction for t_ref) and represent as DataFrame.
         self._post_process(samples)
         self.samples = pd.DataFrame(samples)
-        print(f"Done. This took {time.time() - t0:.1f} s.")
+        logger.info(f"Done. This took {time.time() - t0:.1f} s.")
         sys.stdout.flush()
 
     def log_prob(self, samples: pd.DataFrame | dict) -> np.ndarray:
@@ -430,8 +431,8 @@ def _run_sampler(
             init_samples = self.init_sampler._run_sampler(num_samples, context)
         else:
             if self.num_iterations == 1:
-                print(
-                    f"Warning: Removing initial outliers, but only carrying out "
+                logger.warning(
+                    f"Removing initial outliers, but only carrying out "
                     f"{self.num_iterations} GNPE iteration. This risks biasing "
                     f"results."
                 )
@@ -515,16 +516,15 @@ def _run_sampler(
                 p: x["extrinsic_parameters"][p] for p in self.gnpe_proxy_parameters
             }
 
-            print(
-                f"it {i}.\tmin pvalue: {self.iteration_tracker.pvalue_min:.3f}"
-                f"\tproxy mean: ",
-                *[f"{torch.mean(v).item():.5f}" for v in proxies.values()],
-                "\tproxy std:",
-                *[f"{torch.std(v).item():.5f}" for v in proxies.values()],
-                "\ttimes:",
-                time_sample_start - start_time,
-                time_sample_end - time_sample_start,
-                time.time() - time_sample_end,
+            proxy_means = " ".join(f"{torch.mean(v).item():.5f}" for v in proxies.values())
+            proxy_stds = " ".join(f"{torch.std(v).item():.5f}" for v in proxies.values())
+            logger.debug(
+                f"it {i}. min pvalue: {self.iteration_tracker.pvalue_min:.3f}"
+                f" proxy mean: {proxy_means}"
+                f" proxy std: {proxy_stds}"
+                f" times: {time_sample_start - start_time:.3f}"
+                f" {time_sample_end - time_sample_start:.3f}"
+                f" {time.time() - time_sample_end:.3f}"
             )
 
         #
diff --git a/dingo/core/utils/condor_utils.py b/dingo/core/utils/condor_utils.py
index 1b0e2bbeb..cc51de919 100644
--- a/dingo/core/utils/condor_utils.py
+++ b/dingo/core/utils/condor_utils.py
@@ -2,6 +2,8 @@
 from os.path import join
 import yaml
 
+from dingo.core.utils.logging_utils import logger
+
 
 def resubmit_condor_job(train_dir, train_settings, epoch):
     """
@@ -12,14 +14,14 @@ def resubmit_condor_job(train_dir, train_settings, epoch):
     :return:
     """
     if 'condor_settings' in train_settings:
-        print('Copying log files')
+        logger.info("Copying log files")
         copy_logfiles(train_dir, epoch=epoch)
 
         if epoch >= train_settings['train_settings']['runtime_limits'][
             'max_epochs_total']:
-            print('Training complete, job will not be resubmitted')
+            logger.info("Training complete, job will not be resubmitted")
         else:
-            print('Training incomplete, resubmitting job.')
+            logger.info("Training incomplete, resubmitting job.")
             create_submission_file_and_submit_job(train_dir)
 
 
@@ -76,7 +78,7 @@ def copy_logfiles(log_dir, epoch, name='info', suffixes=('.err','.log','.out')):
         try:
             copyfile(src, dest)
         except:
-            print('Could not copy ' + src)
+            logger.warning("Could not copy " + src)
 
 
 if __name__ == '__main__':
diff --git a/dingo/core/utils/logging_utils.py b/dingo/core/utils/logging_utils.py
index e626fce0b..dd21cd426 100644
--- a/dingo/core/utils/logging_utils.py
+++ b/dingo/core/utils/logging_utils.py
@@ -2,12 +2,6 @@
 import os
 import sys
 
-"""
-Utility functions
-"""
-
-# This is not currently used. We should improve our logging throughout.
-
 
 def check_directory_exists_and_if_not_mkdir(directory, logger):
     """Checks if the given directory exists and creates it if it does not exist
@@ -16,8 +10,6 @@ def check_directory_exists_and_if_not_mkdir(directory, logger):
     ----------
     directory: str
         Name of the directory
-
-    Borrowed from bilby-pipe: https://git.ligo.org/lscsoft/bilby_pipe/
     """
     if not os.path.exists(directory):
         os.makedirs(directory)
@@ -27,7 +19,11 @@ def check_directory_exists_and_if_not_mkdir(directory, logger):
 
 
 def setup_logger(outdir=None, label=None, log_level="INFO"):
-    """Setup logging output: call at the start of the script to use
+    """Setup logging output: call at the start of the script to use.
+
+    Sets up the dingo logger and, if bilby is available, also calls bilby's
+    setup_logger so that both loggers share the same file handler and their
+    entries appear together in the same log file.
 
     Parameters
     ----------
@@ -37,10 +33,7 @@ def setup_logger(outdir=None, label=None, log_level="INFO"):
         ['debug', 'info', 'warning']
         Either a string from the list above, or an integer as specified
         in https://docs.python.org/2/library/logging.html#logging-levels
-
-    Borrowed from bilby-pipe: https://git.ligo.org/lscsoft/bilby_pipe/
     """
-
     if "-v" in sys.argv or "--verbose" in sys.argv:
         log_level = "DEBUG"
 
@@ -52,12 +45,25 @@ def setup_logger(outdir=None, label=None, log_level="INFO"):
     else:
         level = int(log_level)
 
+    # Set up bilby's logger first so we can share its file handler.
+    bilby_logger = None
+    try:
+        from bilby.core.utils.log import setup_logger as bilby_setup_logger
+
+        bilby_setup_logger(outdir=outdir or ".", label=label, log_level=log_level)
+        bilby_logger = logging.getLogger("bilby")
+    except ImportError:
+        pass
+
     logger = logging.getLogger("dingo")
     logger.propagate = False
     logger.setLevel(level)
 
-    streams = [isinstance(h, logging.StreamHandler) for h in logger.handlers]
-    if len(streams) == 0 or not all(streams):
+    # Add a stream handler if not already present.
+    if not any(
+        isinstance(h, logging.StreamHandler) and not isinstance(h, logging.FileHandler)
+        for h in logger.handlers
+    ):
         stream_handler = logging.StreamHandler()
         stream_handler.setFormatter(
             logging.Formatter(
@@ -67,22 +73,34 @@ def setup_logger(outdir=None, label=None, log_level="INFO"):
         stream_handler.setLevel(level)
         logger.addHandler(stream_handler)
 
-    if any([isinstance(h, logging.FileHandler) for h in logger.handlers]) is False:
-        if label:
-            if outdir:
-                check_directory_exists_and_if_not_mkdir(outdir, logger)
-            else:
-                outdir = "."
-            log_file = f"{outdir}/{label}.log"
-            file_handler = logging.FileHandler(log_file)
-            file_handler.setFormatter(
-                logging.Formatter(
-                    "%(asctime)s %(levelname)-8s: %(message)s", datefmt="%H:%M"
+    # Share bilby's file handler so both loggers write to the same log file.
+    if bilby_logger is not None:
+        for handler in bilby_logger.handlers:
+            if isinstance(handler, logging.FileHandler):
+                handler.setFormatter(
+                    logging.Formatter(
+                        "%(asctime)s %(name)s %(levelname)-8s: %(message)s",
+                        datefmt="%H:%M",
+                    )
                 )
+                if not any(isinstance(h, logging.FileHandler) for h in logger.handlers):
+                    logger.addHandler(handler)
+
+    # Fall back to a standalone file handler if bilby is unavailable.
+    if label and not any(isinstance(h, logging.FileHandler) for h in logger.handlers):
+        if outdir:
+            check_directory_exists_and_if_not_mkdir(outdir, logger)
+        else:
+            outdir = "."
+        log_file = f"{outdir}/{label}.log"
+        file_handler = logging.FileHandler(log_file)
+        file_handler.setFormatter(
+            logging.Formatter(
+                "%(asctime)s %(name)s %(levelname)-8s: %(message)s", datefmt="%H:%M"
             )
-
-            file_handler.setLevel(level)
-            logger.addHandler(file_handler)
+        )
+        file_handler.setLevel(level)
+        logger.addHandler(file_handler)
 
     for handler in logger.handlers:
         handler.setLevel(level)
diff --git a/dingo/core/utils/pt_to_hdf5.py b/dingo/core/utils/pt_to_hdf5.py
index 06e6aa96c..c3a0f7e81 100644
--- a/dingo/core/utils/pt_to_hdf5.py
+++ b/dingo/core/utils/pt_to_hdf5.py
@@ -5,6 +5,8 @@
 import json
 import argparse
 
+from dingo.core.utils.logging_utils import logger
+
 
 def parse_args():
     parser = argparse.ArgumentParser(
@@ -32,7 +34,7 @@ def main():
     # This is required for use on CVMFS
     root, ext = os.path.splitext(args.out_file)
     out_file_name = f'{root}_v{args.model_version_number}{ext}'
-    print('Output will be written to', out_file_name)
+    logger.info(f'Output will be written to {out_file_name}')
 
     # Load data into CPU memory since we'll be saving it using CPU libraries
     d = torch.load(args.in_file, map_location=torch.device("cpu"))
diff --git a/dingo/core/utils/trainutils.py b/dingo/core/utils/trainutils.py
index bc7d12154..87de405d0 100644
--- a/dingo/core/utils/trainutils.py
+++ b/dingo/core/utils/trainutils.py
@@ -5,6 +5,8 @@
 import csv
 from typing import Literal
 
+from dingo.core.utils.logging_utils import logger
+
 
 class AvgTracker:
     def __init__(self):
@@ -85,7 +87,7 @@ def __call__(self, val_loss: float):
         elif score < self.best_score + self.delta:
             self.counter += 1
             if self.verbose:
-                print(f"EarlyStopping counter: {self.counter} out of {self.patience}")
+                logger.info(f"EarlyStopping counter: {self.counter} out of {self.patience}")
             if self.counter >= self.patience:
                 self.early_stop = True
             return False
@@ -124,23 +126,26 @@ def get_avg(self):
 
     def print_info(self, batch_idx):
         if batch_idx % self.print_freq == 0:
-            print(
-                "{} Epoch: {} [{}/{} ({:.0f}%)]".format(
+            td, td_avg = self.times["Dataloader"].x, self.times["Dataloader"].get_avg()
+            tn, tn_avg = self.times["Network"].x, self.times["Network"].get_avg()
+            logger.debug(
+                "{} Epoch: {} [{}/{} ({:.0f}%)]  "
+                "Loss: {:.3f} ({:.3f})  "
+                "Time Dataloader: {:.3f} ({:.3f})  "
+                "Time Network: {:.3f} ({:.3f})".format(
                     self.mode,
                     self.epoch,
                     min(batch_idx * self.batch_size, self.len_dataset),
                     self.len_dataset,
                     100.0 * batch_idx * self.batch_size / self.len_dataset,
-                ),
-                end="\t\t",
+                    self.loss,
+                    self.get_avg(),
+                    td,
+                    td_avg,
+                    tn,
+                    tn_avg,
+                )
             )
-            # print loss
-            print(f"Loss: {self.loss:.3f} ({self.get_avg():.3f})", end="\t\t")
-            # print computation times
-            td, td_avg = self.times["Dataloader"].x, self.times["Dataloader"].get_avg()
-            tn, tn_avg = self.times["Network"].x, self.times["Network"].get_avg()
-            print(f"Time Dataloader: {td:.3f} ({td_avg:.3f})", end="\t\t")
-            print(f"Time Network: {tn:.3f} ({tn_avg:.3f})")
 
 
 class RuntimeLimits:
@@ -195,8 +200,8 @@ def limits_exceeded(self, epoch: int = None):
         # check time limit for run
         if self.max_time_per_run is not None:
             if time.time() - self.time_start >= self.max_time_per_run:
-                print(
-                    f"Stop run: Time limit of {self.max_time_per_run} s " f"exceeded."
+                logger.info(
+                    f"Stop run: Time limit of {self.max_time_per_run} s exceeded."
                 )
                 return True
         # check epoch limit for run
@@ -204,14 +209,14 @@ def limits_exceeded(self, epoch: int = None):
             if epoch is None:
                 raise ValueError("epoch required")
             if epoch - self.epoch_start >= self.max_epochs_per_run:
-                print(
+                logger.info(
                     f"Stop run: Epoch limit of {self.max_epochs_per_run} per run reached."
                 )
                 return True
         # check total epoch limit
         if self.max_epochs_total is not None:
             if epoch >= self.max_epochs_total:
-                print(
+                logger.info(
                     f"Stop run: Total epoch limit of {self.max_epochs_total} reached."
                 )
                 return True
@@ -316,13 +321,13 @@ def save_model(pm, log_dir, model_prefix="model", checkpoint_epochs=None):
     """
     # save current model
     model_name = join(log_dir, f"{model_prefix}_latest.pt")
-    print(f"Saving model to {model_name}.", end=" ")
+    logger.info(f"Saving model to {model_name}.")
     pm.save_model(model_name, save_training_info=True)
-    print("Done.")
+    logger.info("Done.")
 
     # potentially copy model to a checkpoint
     if checkpoint_epochs is not None and pm.epoch % checkpoint_epochs == 0:
         model_name_cp = join(log_dir, f"{model_prefix}_{pm.epoch:03d}.pt")
-        print(f"Copy model to checkpoint {model_name_cp}.", end=" ")
+        logger.info(f"Copying model to checkpoint {model_name_cp}.")
         copyfile(model_name, model_name_cp)
-        print("Done.")
+        logger.info("Done.")
diff --git a/dingo/gw/SVD.py b/dingo/gw/SVD.py
index e77b7761d..561d60f9a 100644
--- a/dingo/gw/SVD.py
+++ b/dingo/gw/SVD.py
@@ -3,6 +3,7 @@
 import scipy
 from sklearn.utils.extmath import randomized_svd
 from dingo.core.dataset import DingoDataset
+from dingo.core.utils.logging_utils import logger
 
 class SVDBasis(DingoDataset):
 
@@ -155,15 +156,25 @@ def print_validation_summary(self):
                 if "mismatch" in col:
                     n = int(col.split(sep="=")[-1])
                     mismatches = self.mismatches[col]
-                    print(f"n = {n}")
-                    print("  Mean mismatch = {}".format(np.mean(mismatches)))
-                    print("  Standard deviation = {}".format(np.std(mismatches)))
-                    print("  Max mismatch = {}".format(np.max(mismatches)))
-                    print("  Median mismatch = {}".format(np.median(mismatches)))
-                    print("  Percentiles:")
-                    print("    99    -> {}".format(np.percentile(mismatches, 99)))
-                    print("    99.9  -> {}".format(np.percentile(mismatches, 99.9)))
-                    print("    99.99 -> {}".format(np.percentile(mismatches, 99.99)))
+                    logger.info(
+                        f"n = {n}\n"
+                        "  Mean mismatch = {}\n"
+                        "  Standard deviation = {}\n"
+                        "  Max mismatch = {}\n"
+                        "  Median mismatch = {}\n"
+                        "  Percentiles:\n"
+                        "    99    -> {}\n"
+                        "    99.9  -> {}\n"
+                        "    99.99 -> {}".format(
+                            np.mean(mismatches),
+                            np.std(mismatches),
+                            np.max(mismatches),
+                            np.median(mismatches),
+                            np.percentile(mismatches, 99),
+                            np.percentile(mismatches, 99.9),
+                            np.percentile(mismatches, 99.99),
+                        )
+                    )
 
     def decompress(self, coefficients: np.ndarray):
         """
diff --git a/dingo/gw/conversion/spin_conversion.py b/dingo/gw/conversion/spin_conversion.py
index 9e4b481cf..f2a40b45b 100644
--- a/dingo/gw/conversion/spin_conversion.py
+++ b/dingo/gw/conversion/spin_conversion.py
@@ -8,6 +8,8 @@
 import lalsimulation as LS
 from threadpoolctl import threadpool_limits
 
+from dingo.core.utils.logging_utils import logger
+
 
 DINGO_PE_SPIN_PARAMETERS = (
     "theta_jn",
@@ -225,5 +227,5 @@ def _convert_phase(f_ref, sc_phase_old, sc_phase_new, sample):
 
         return sample_pe_new
     except RuntimeError:
-        print("Failed to convert spins. Saving sample unchanged.")
+        logger.warning("Failed to convert spins. Saving sample unchanged.")
         return sample
diff --git a/dingo/gw/data/data_download.py b/dingo/gw/data/data_download.py
index 31e6b1c92..c0e2f6ed1 100644
--- a/dingo/gw/data/data_download.py
+++ b/dingo/gw/data/data_download.py
@@ -3,6 +3,7 @@
 import pycbc.psd
 import math
 
+from dingo.core.utils.logging_utils import logger
 from dingo.gw.gwutils import get_window
 
 
@@ -42,7 +43,7 @@ def download_psd(det, time_start, time_psd, window, f_s):
     # if strain for PSD data contains nan, shift segment for PSD
     if np.max(np.isnan(psd_strain)):
         dt = math.ceil(np.where(np.isnan(psd_strain))[0][-1] / f_s)
-        print(
+        logger.warning(
             f"Nan encountered in strain data for PSD estimation for detector {det}. "
             f"Shifting strain segment by {dt} seconds."
         )
diff --git a/dingo/gw/data/data_preparation.py b/dingo/gw/data/data_preparation.py
index 0b06b1602..4f355b8b4 100644
--- a/dingo/gw/data/data_preparation.py
+++ b/dingo/gw/data/data_preparation.py
@@ -4,6 +4,7 @@
 from gwpy.timeseries import TimeSeries
 
 from dingo.core.dataset import DingoDataset
+from dingo.core.utils.logging_utils import logger
 from dingo.core.utils.misc import recursive_check_dicts_are_equal
 from dingo.gw.data.data_download import download_raw_data
 from dingo.gw.gwutils import get_window
@@ -45,11 +46,11 @@ def load_raw_data(time_event, settings, event_dataset=None):
                         f"{dataset.settings}"
                     )
             data = vars(dataset)[event]
-            print(f"Data for event at {event} found in {event_dataset}.")
+            logger.info(f"Data for event at {event} found in {event_dataset}.")
             return data
 
     # if this did not work, download the data
-    print(f"Downloading data for event at {event}.")
+    logger.info(f"Downloading data for event at {event}.")
     data = download_raw_data(time_event, **settings)
 
     # optionally save this data to event_dataset
@@ -57,7 +58,7 @@ def load_raw_data(time_event, settings, event_dataset=None):
         dataset = DingoDataset(
             dictionary={event: data, "settings": settings}, data_keys=[event]
         )
-        print(f"Saving data for event at {event} to {event_dataset}.")
+        logger.info(f"Saving data for event at {event} to {event_dataset}.")
         dataset.to_file(event_dataset, mode="a")
 
     return data
diff --git a/dingo/gw/dataset/evaluate_multibanded_domain.py b/dingo/gw/dataset/evaluate_multibanded_domain.py
index b20a2fb43..484603c29 100644
--- a/dingo/gw/dataset/evaluate_multibanded_domain.py
+++ b/dingo/gw/dataset/evaluate_multibanded_domain.py
@@ -4,6 +4,7 @@
 import yaml
 from scipy.interpolate import interp1d
 
+from dingo.core.utils.logging_utils import logger
 from dingo.gw.dataset import generate_parameters_and_polarizations
 from dingo.gw.domains import build_domain, MultibandedFrequencyDomain
 from dingo.gw.gwutils import get_mismatch
@@ -35,13 +36,13 @@ def _evaluate_multibanding_main(
     settings["intrinsic_prior"]["chirp_mass"] = prior["chirp_mass"].minimum
     # Rebuild prior with updated settings.
     prior = build_prior_with_defaults(settings["intrinsic_prior"])
-    print("Prior")
+    logger.info("Prior")
     for k, v in prior.items():
-        print(f"{k}: {v}")
+        logger.info(f"{k}: {v}")
 
     domain = build_domain(settings["domain"])
-    print("\nDomain")
-    print(domain.domain_dict)
+    logger.info("\nDomain")
+    logger.info(domain.domain_dict)
 
     if not isinstance(domain, MultibandedFrequencyDomain):
         raise ValueError("Waveform dataset domain not a MultibandedFrequencyDomain.")
@@ -88,21 +89,21 @@ def _evaluate_multibanding_main(
                 asd_file="aLIGO_ZERO_DET_high_P_asd.txt",
             )
 
-    print("\nMismatches between UFD waveforms and MFD waveforms interpolated to MFD.")
-    print(
+    logger.info("\nMismatches between UFD waveforms and MFD waveforms interpolated to MFD.")
+    logger.info(
         "This is a conservative estimate of the MFD performance when training "
         "networks."
     )
     mismatches = np.concatenate([v for v in mismatches.values()])
-    print(f"num_samples = {num_samples}")
-    print("  Mean mismatch = {}".format(np.mean(mismatches)))
-    print("  Standard deviation = {}".format(np.std(mismatches)))
-    print("  Max mismatch = {}".format(np.max(mismatches)))
-    print("  Median mismatch = {}".format(np.median(mismatches)))
-    print("  Percentiles:")
-    print("    99    -> {}".format(np.percentile(mismatches, 99)))
-    print("    99.9  -> {}".format(np.percentile(mismatches, 99.9)))
-    print("    99.99 -> {}".format(np.percentile(mismatches, 99.99)))
+    logger.info(f"num_samples = {num_samples}")
+    logger.info("  Mean mismatch = {}".format(np.mean(mismatches)))
+    logger.info("  Standard deviation = {}".format(np.std(mismatches)))
+    logger.info("  Max mismatch = {}".format(np.max(mismatches)))
+    logger.info("  Median mismatch = {}".format(np.median(mismatches)))
+    logger.info("  Percentiles:")
+    logger.info("    99    -> {}".format(np.percentile(mismatches, 99)))
+    logger.info("    99.9  -> {}".format(np.percentile(mismatches, 99.9)))
+    logger.info("    99.99 -> {}".format(np.percentile(mismatches, 99.99)))
 
 
 def parse_args():
diff --git a/dingo/gw/dataset/generate_dataset.py b/dingo/gw/dataset/generate_dataset.py
index a52ca9865..4e5886cad 100644
--- a/dingo/gw/dataset/generate_dataset.py
+++ b/dingo/gw/dataset/generate_dataset.py
@@ -24,6 +24,7 @@
     generate_waveforms_parallel,
 )
 from dingo.core.utils.misc import call_func_strict_output_dim
+from dingo.core.utils.logging_utils import logger, setup_logger
 
 
 def generate_parameters_and_polarizations(
@@ -47,7 +48,7 @@ def generate_parameters_and_polarizations(
     pandas DataFrame of parameters
     dictionary of numpy arrays corresponding to waveform polarizations
     """
-    print("Generating dataset of size " + str(num_samples))
+    logger.info("Generating dataset of size " + str(num_samples))
     parameters = pd.DataFrame(prior.sample(num_samples))
 
     if num_processes > 1:
@@ -67,12 +68,12 @@ def generate_parameters_and_polarizations(
         polarizations_ok = {k: v[idx_ok] for k, v in polarizations.items()}
         parameters_ok = parameters.iloc[idx_ok]
         failed_percent = 100 * len(idx_failed) / len(parameters)
-        print(
+        logger.warning(
             f"{len(idx_failed)} out of {len(parameters)} configuration ({failed_percent:.1f}%) failed to generate."
         )
         with pd.option_context("display.max_rows", None, "display.max_columns", None):
-            print(parameters.iloc[idx_failed])
-        print(
+            logger.warning(parameters.iloc[idx_failed].to_string())
+        logger.warning(
             f"Only returning the {len(idx_ok)} successfully generated configurations."
         )
         return parameters_ok, polarizations_ok
@@ -115,7 +116,7 @@ def train_svd_basis(dataset: WaveformDataset, size: int, n_train: int):
     )
     test_parameters.reset_index(drop=True, inplace=True)
 
-    print("Building SVD basis.")
+    logger.info("Building SVD basis.")
     basis = SVDBasis()
     basis.generate_basis(train_data, size)
 
@@ -301,6 +302,8 @@ def _generate_dataset_main(
 
 def main() -> None:
     args = parse_args()
+    out_path = Path(args.out_file)
+    setup_logger(outdir=str(out_path.parent), label=out_path.stem)
     _generate_dataset_main(args.settings_file, args.out_file, args.num_processes)
 
 
diff --git a/dingo/gw/dataset/generate_dataset_dag.py b/dingo/gw/dataset/generate_dataset_dag.py
index 395fa8fa8..718149f72 100644
--- a/dingo/gw/dataset/generate_dataset_dag.py
+++ b/dingo/gw/dataset/generate_dataset_dag.py
@@ -5,6 +5,8 @@
 import yaml
 import copy
 
+from dingo.core.utils.logging_utils import logger, setup_logger
+
 # Fixed file names
 svd_fn = "svd.hdf5"
 settings_part_fn = "settings_part.yaml"
@@ -264,6 +266,7 @@ def create_dag(args, settings):
 
 def main():
     args = parse_args()
+    setup_logger(outdir=args.temp_dir, label="generate_dataset_dag")
 
     # Load settings
     with open(args.settings_file, "r") as f:
@@ -284,7 +287,7 @@ def main():
         pass
 
     dagman.build()
-    print(f"DAG submission file written to {args.submit}.")
+    logger.info(f"DAG submission file written to {args.submit}.")
 
 
 if __name__ == "__main__":
diff --git a/dingo/gw/dataset/utils.py b/dingo/gw/dataset/utils.py
index e8911f87c..5308ef320 100644
--- a/dingo/gw/dataset/utils.py
+++ b/dingo/gw/dataset/utils.py
@@ -1,11 +1,13 @@
 import argparse
 import textwrap
 import copy
+from pathlib import Path
 import pandas as pd
 import numpy as np
 import yaml
 from typing import List
 
+from dingo.core.utils.logging_utils import logger, setup_logger
 from dingo.gw.SVD import SVDBasis
 from dingo.gw.dataset.generate_dataset import train_svd_basis
 from dingo.gw.dataset.waveform_dataset import WaveformDataset
@@ -26,7 +28,7 @@ def merge_datasets(dataset_list: List[WaveformDataset]) -> WaveformDataset:
     WaveformDataset containing the merged data.
     """
 
-    print(f"Merging {len(dataset_list)} datasets into one.")
+    logger.info(f"Merging {len(dataset_list)} datasets into one.")
 
     # This ensures that all the keys are copied into the new dataset. The "extensive"
     # parts of the dataset (parameters, waveforms) will be overwritten by the combined
@@ -87,6 +89,8 @@ def merge_datasets_cli():
         "--settings_file", type=str, help="YAML file containing new dataset settings."
     )
     args = parser.parse_args()
+    out_path = Path(args.out_file)
+    setup_logger(outdir=str(out_path.parent), label=out_path.stem)
 
     dataset_list = []
     for i in range(args.num_parts):
@@ -103,7 +107,7 @@ def merge_datasets_cli():
         merged_dataset.settings = settings
 
     merged_dataset.to_file(args.out_file)
-    print(
+    logger.info(
         f"Complete. New dataset consists of {merged_dataset.settings['num_samples']} "
         f"samples."
     )
@@ -141,6 +145,8 @@ def build_svd_cli():
         "Remainder are used for validation.",
     )
     args = parser.parse_args()
+    out_path = Path(args.out_file)
+    setup_logger(outdir=str(out_path.parent), label=out_path.stem)
 
     dataset = WaveformDataset(file_name=args.dataset_file)
     if args.num_train is None:
@@ -151,7 +157,7 @@ def build_svd_cli():
     basis, n_train, n_test = train_svd_basis(dataset, args.size, n_train)
     # FIXME: This is not an ideal treatment. We should update the waveform generation
     #  to always provide the requested number of waveforms.
-    print(
+    logger.info(
         f"SVD basis trained based on {n_train} waveforms and validated on {n_test} "
         f"waveforms. Note that if this differs from number requested, it will not be "
         f"reflected in the settings file. This is likely due to EOB failure to "
diff --git a/dingo/gw/download_strain_data.py b/dingo/gw/download_strain_data.py
index 336954098..87aa55693 100644
--- a/dingo/gw/download_strain_data.py
+++ b/dingo/gw/download_strain_data.py
@@ -2,6 +2,7 @@
 import pycbc.psd
 from gwpy.timeseries import TimeSeries
 
+from dingo.core.utils.logging_utils import logger
 from dingo.gw.domains import UniformFrequencyDomain
 from dingo.gw.gwutils import (
     get_window,
@@ -107,14 +108,14 @@ def download_strain_data_in_FD(det, time_event, time_segment, time_buffer, windo
         array with the frequency domain strain
     """
     # download strain data
-    print("Downloading strain data for event.", end=" ")
+    logger.info("Downloading strain data for event.")
     event_strain = TimeSeries.fetch_open_data(
         det,
         time_event + time_buffer - time_segment,
         time_event + time_buffer,
         cache=True,
     )
-    print("Done.")
+    logger.info("Done.")
 
     # transform to FD
     if type(window) == dict:
@@ -165,7 +166,7 @@ def download_event_data_in_FD(
     """
     data = {"waveform": {}, "asds": {}}
     for det in detectors:
-        print("Detector {:}:".format(det))
+        logger.info("Detector {:}:".format(det))
 
         data["waveform"][det] = download_strain_data_in_FD(
             det, time_event, time_segment, time_buffer, window
diff --git a/dingo/gw/importance_sampling/diagnostics.py b/dingo/gw/importance_sampling/diagnostics.py
index 39cc95bb5..c178502f8 100644
--- a/dingo/gw/importance_sampling/diagnostics.py
+++ b/dingo/gw/importance_sampling/diagnostics.py
@@ -3,6 +3,7 @@
 import math
 import pandas as pd
 import matplotlib.pyplot as plt
+from dingo.core.utils.logging_utils import logger
 from dingo.core.utils.plotting import plot_corner_multi
 from dingo.gw.result import Result
 
@@ -219,7 +220,7 @@ def plot_diagnostics(
     inds = np.where(weights > threshold)[0]
     theta_new = theta.loc[inds]
     weights_new = weights[inds]
-    print(
+    logger.info(
         f"Generating cornerplot with {len(theta_new)} out of {len(theta)} IS samples."
     )
 
diff --git a/dingo/gw/importance_sampling/importance_weights.py b/dingo/gw/importance_sampling/importance_weights.py
index 7f88fa67e..e6c60a482 100644
--- a/dingo/gw/importance_sampling/importance_weights.py
+++ b/dingo/gw/importance_sampling/importance_weights.py
@@ -9,6 +9,7 @@
 from os.path import dirname, join, isfile, exists
 import argparse
 
+from dingo.core.utils.logging_utils import logger
 from dingo.core.posterior_models import NormalizingFlowPosteriorModel
 from dingo.gw.result import Result
 from dingo.gw.inference.gw_samplers import GWSampler
@@ -100,20 +101,20 @@ def main():
             "path", join(args.outdir, f"nde-{event_name}.pt")
         )
         if isfile(nde_name):
-            print(f"Loading nde at {nde_name} for event {event_name}.")
+            logger.info(f"Loading nde at {nde_name} for event {event_name}.")
             nde = NormalizingFlowPosteriorModel(
                 model_filename=nde_name,
                 device=settings["nde"]["training"]["device"],
                 load_training_info=False,
             )
         else:
-            print(f"Training new nde for event {event_name}.")
+            logger.info(f"Training new nde for event {event_name}.")
             nde = result.train_unconditional_flow(
                 inference_parameters,
                 settings["nde"],
                 train_dir=args.outdir,
             )
-            print(f"Renaming trained nde model to {nde_name}.")
+            logger.info(f"Renaming trained nde model to {nde_name}.")
             rename(join(args.outdir, "model_latest.pt"), nde_name)
 
         # Step 1a: Sample from proposal.
@@ -162,10 +163,10 @@ def main():
     # print(np.std(log_evidences) / np.mean(log_evidences_std))
 
     if synthetic_phase:
-        print(f"Sampling synthetic phase.")
+        logger.info("Sampling synthetic phase.")
         result.sample_synthetic_phase(synthetic_phase_kwargs)
 
-    print(f"Importance sampling.")
+    logger.info("Importance sampling.")
     result.importance_sample(
         num_processes=settings.get("num_processes", 1),
         time_marginalization_kwargs=time_marginalization_kwargs,
@@ -179,7 +180,7 @@ def main():
     diagnostics_dir = join(args.outdir, "IS-diagnostics")
     if not exists(diagnostics_dir):
         makedirs(diagnostics_dir)
-    print("Plotting diagnostics.")
+    logger.info("Plotting diagnostics.")
     plot_diagnostics(
         result,
         diagnostics_dir,
diff --git a/dingo/gw/inference/gw_samplers.py b/dingo/gw/inference/gw_samplers.py
index 2a1b63d89..c7d8f5585 100644
--- a/dingo/gw/inference/gw_samplers.py
+++ b/dingo/gw/inference/gw_samplers.py
@@ -8,6 +8,7 @@
 from torchvision.transforms import Compose
 
 from dingo.core.samplers import Sampler, GNPESampler
+from dingo.core.utils.logging_utils import logger
 from dingo.core.transforms import GetItem, RenameKey
 from dingo.gw.domains import (
     MultibandedFrequencyDomain,
@@ -240,7 +241,7 @@ def _post_process(self, samples: Union[dict, pd.DataFrame], inverse: bool = Fals
             for k, p in prior.items():
                 if isinstance(p, DeltaFunction) and k not in samples:
                     v = p.peak
-                    print(f"Adding fixed parameter {k} = {v} from prior.")
+                    logger.info(f"Adding fixed parameter {k} = {v} from prior.")
                     samples[k] = p.peak * np.ones(num_samples)
         else:
             # Drop non-inference parameters from samples.
@@ -472,8 +473,8 @@ def _initialize_transforms(self):
                 self.gnpe_parameters += transform.input_parameter_names
                 for k, v in transform.kernel.items():
                     self.gnpe_kernel[k] = v
-        print("GNPE parameters: ", self.gnpe_parameters)
-        print("GNPE kernel: ", self.gnpe_kernel)
+        logger.info(f"GNPE parameters: {self.gnpe_parameters}")
+        logger.info(f"GNPE kernel: {self.gnpe_kernel}")
 
         self.transform_pre = Compose(transform_pre)
 
diff --git a/dingo/gw/injection.py b/dingo/gw/injection.py
index 87f7a0660..f1679fb47 100644
--- a/dingo/gw/injection.py
+++ b/dingo/gw/injection.py
@@ -2,6 +2,7 @@
 from bilby.gw.detector import InterferometerList
 from torchvision.transforms import Compose
 
+from dingo.core.utils.logging_utils import logger
 from dingo.gw.noise.asd_dataset import ASDDataset
 from dingo.gw.domains import (
     UniformFrequencyDomain,
@@ -312,7 +313,7 @@ def asd(self, asd):
             if set(asd.asds.keys()) != set(ifo_names):
                 raise KeyError("ASDDataset ifos do not match signal.")
             if asd.domain.domain_dict != self.data_domain.domain_dict:
-                print("Updating ASDDataset domain to match data domain.")
+                logger.info("Updating ASDDataset domain to match data domain.")
                 domain_dict = self.data_domain.domain_dict
                 asd.update_domain(domain_dict)
         elif isinstance(asd, dict):
@@ -402,7 +403,7 @@ def injection(self, theta):
             raise ValueError("self.asd must be set in order to produce injections.")
 
         if self.whiten:
-            print("self.whiten was set to True. Resetting to False.")
+            logger.warning("self.whiten was set to True. Resetting to False.")
             self.whiten = False
 
         data = {}
diff --git a/dingo/gw/likelihood.py b/dingo/gw/likelihood.py
index 49b29941f..f48893ad2 100644
--- a/dingo/gw/likelihood.py
+++ b/dingo/gw/likelihood.py
@@ -22,6 +22,7 @@
 )
 from dingo.gw.domains import build_domain
 from dingo.gw.data.data_preparation import get_event_data_and_domain
+from dingo.core.utils.logging_utils import logger
 
 
 class StationaryGaussianGWLikelihood(GWSignal, Likelihood):
@@ -173,7 +174,7 @@ def __init__(
                 # use endpoint = False for grid, since phase = 0/2pi are equivalent
                 self.phase_grid = np.linspace(0, 2 * np.pi, n_grid, endpoint=False)
             else:
-                print("Using phase marginalization with (2,2) mode approximation.")
+                logger.info("Using phase marginalization with (2,2) mode approximation.")
 
         # Initialize calibration marginalization using the setter from GWSignal.
         self.calibration_marginalization_kwargs = calibration_marginalization_kwargs
@@ -870,15 +871,15 @@ def main():
         try:
             l = likelihood.log_prob(theta)
         except:
-            print(idx)
+            logger.warning(idx)
             l = float("nan")
         log_likelihoods.append(l)
     log_likelihoods = np.array(log_likelihoods)
     log_likelihoods = log_likelihoods[~np.isnan(log_likelihoods)]
-    print(f"mean: {np.mean(log_likelihoods)}")
-    print(f"std: {np.std(log_likelihoods)}")
-    print(f"max: {np.max(log_likelihoods)}")
-    print(f"min: {np.min(log_likelihoods)}")
+    logger.info(f"mean: {np.mean(log_likelihoods)}")
+    logger.info(f"std: {np.std(log_likelihoods)}")
+    logger.info(f"max: {np.max(log_likelihoods)}")
+    logger.info(f"min: {np.min(log_likelihoods)}")
 
 
 if __name__ == "__main__":
diff --git a/dingo/gw/noise/asd_dataset.py b/dingo/gw/noise/asd_dataset.py
index 10e9cf9ef..9e3d1fb02 100644
--- a/dingo/gw/noise/asd_dataset.py
+++ b/dingo/gw/noise/asd_dataset.py
@@ -4,6 +4,7 @@
 
 import numpy as np
 
+from dingo.core.utils.logging_utils import logger
 from dingo.gw.domains import build_domain, UniformFrequencyDomain
 from dingo.gw.domains.base_frequency_domain import BaseFrequencyDomain
 from dingo.gw.gwutils import *
@@ -59,8 +60,8 @@ def __init__(
                     self.gps_times.pop(ifo)
 
         if "window_factor" in self.settings["domain_dict"]:
-            print(
-                "Warning: 'window_factor' is no longer used in ASDDataset. "
+            logger.warning(
+                "'window_factor' is no longer used in ASDDataset. "
                 "Removing from settings."
             )
             self.settings["domain_dict"].pop("window_factor")
@@ -135,14 +136,14 @@ def update_domain(self, domain_update):
             self.domain.domain_dict["type"] == "UniformFrequencyDomain"
             and domain_update["type"] == "MultibandedFrequencyDomain"
         ):
-            print("Updating ASD dataset to MultibandedFrequencyDomain.")
+            logger.info("Updating ASD dataset to MultibandedFrequencyDomain.")
             asd_dataset_decimated = {}
             mfd = build_domain(domain_update)
             ufd = mfd.base_domain
             if not check_domain_compatibility(self.asds, ufd):
                 # If the ASD length is not compatible with the new base UFD,
                 # first truncate it.
-                print(
+                logger.info(
                     f"  Truncating first to new base UniformFrequencyDomain: f_max "
                     f"{self.domain.f_max} Hz -> {ufd.f_max} Hz"
                 )
diff --git a/dingo/gw/noise/generate_dataset.py b/dingo/gw/noise/generate_dataset.py
index e8750e6e2..24f5c34a8 100644
--- a/dingo/gw/noise/generate_dataset.py
+++ b/dingo/gw/noise/generate_dataset.py
@@ -12,6 +12,7 @@
 from dingo.gw.noise.asd_dataset import ASDDataset
 from dingo.gw.noise.generate_dataset_dag import create_dag
 from dingo.gw.noise.utils import merge_datasets, get_time_segments
+from dingo.core.utils.logging_utils import logger, setup_logger
 
 
 def parse_args():
@@ -58,6 +59,7 @@ def generate_dataset():
     Creates and saves an ASD dataset
     """
     args = parse_args()
+    setup_logger(outdir=args.data_dir, label="asd_dataset")
 
     # Load settings
     settings_file = (
@@ -94,11 +96,11 @@ def generate_dataset():
             pass
 
         dagman.build()
-        print(f"DAG submission file written.")
+        logger.info("DAG submission file written.")
 
     else:
 
-        print("Downloading strain data and estimating PSDs...")
+        logger.info("Downloading strain data and estimating PSDs...")
         asd_filename_list = download_and_estimate_psds(
             args.data_dir, settings, time_segments, verbose=args.verbose
         )
@@ -106,7 +108,7 @@ def generate_dataset():
             det: [ASDDataset(asd_file) for asd_file in asd_file_list]
             for det, asd_file_list in asd_filename_list.items()
         }
-        print("Merging single dataset files into one...")
+        logger.info("Merging single dataset files into one...")
         dataset = merge_datasets(asd_dataset_list)
         filename = args.out_name
         if filename is None:
diff --git a/dingo/gw/noise/synthetic/asd_sampling.py b/dingo/gw/noise/synthetic/asd_sampling.py
index 707ca586d..e5bfef1a6 100644
--- a/dingo/gw/noise/synthetic/asd_sampling.py
+++ b/dingo/gw/noise/synthetic/asd_sampling.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 from scipy import stats
+from dingo.core.utils.logging_utils import logger
 from dingo.gw.noise.synthetic.asd_parameterization import fit_broadband_noise
 from dingo.gw.noise.asd_dataset import ASDDataset
 from dingo.gw.noise.synthetic.utils import (
@@ -53,7 +54,7 @@ def fit(self, weights=None):
                         weights=weights,
                     )
                 except np.linalg.LinAlgError:
-                    print(
+                    logger.warning(
                         "Warning: Singular Matrix encountered in spectral KDE. Adding small Gaussian noise..."
                     )
                     perturbed_features = spectral_features[:, i, :] + np.random.normal(
diff --git a/dingo/gw/noise/utils.py b/dingo/gw/noise/utils.py
index 43d1f82ce..2d9fa4578 100644
--- a/dingo/gw/noise/utils.py
+++ b/dingo/gw/noise/utils.py
@@ -12,6 +12,7 @@
 import yaml
 
 from gwpy.table import EventTable
+from dingo.core.utils.logging_utils import logger
 from dingo.gw.noise.asd_dataset import ASDDataset
 
 """
@@ -156,7 +157,7 @@ def merge_datasets(asd_dataset_list):
     merged_dict = {"asds": {}, "gps_times": {}}
 
     for det, asd_list in asd_dataset_list.items():
-        print(f"Merging {len(asd_list)} datasets into one for detector {det}.")
+        logger.info(f"Merging {len(asd_list)} datasets into one for detector {det}.")
         merged_dict["asds"][det] = np.vstack(
             [asd_dataset.asds[det] for asd_dataset in asd_list]
         )
diff --git a/dingo/gw/result.py b/dingo/gw/result.py
index 51616032d..b0499ce28 100644
--- a/dingo/gw/result.py
+++ b/dingo/gw/result.py
@@ -14,6 +14,7 @@
 )
 from dingo.core.multiprocessing import apply_func_with_multiprocessing
 from dingo.core.result import Result as CoreResult
+from dingo.core.utils.logging_utils import logger
 from dingo.gw.conversion import change_spin_conversion_phase
 from dingo.gw.domains import MultibandedFrequencyDomain
 from dingo.gw.domains import build_domain
@@ -199,18 +200,14 @@ def _rebuild_domain(self, verbose=False):
             )
 
             if verbose:
-                print("Rebuilding domain as follows:")
-                print(
-                    yaml.dump(
-                        domain_dict,
-                        default_flow_style=False,
-                        sort_keys=False,
-                    )
+                logger.info(
+                    "Rebuilding domain as follows:\n"
+                    + yaml.dump(domain_dict, default_flow_style=False, sort_keys=False)
                 )
             self.domain = build_domain(domain_dict)
         else:
             if verbose:
-                print("No domain updates found; domain not rebuilt.")
+                logger.info("No domain updates found; domain not rebuilt.")
 
     def _build_prior(self):
         """Build the prior based on model metadata. Called by __init__()."""
@@ -359,7 +356,7 @@ def _build_likelihood(
         if "updates" in self.importance_sampling_metadata:
             if "T" in self.importance_sampling_metadata["updates"]:
                 delta_f_new = 1 / self.importance_sampling_metadata["updates"]["T"]
-                print(
+                logger.info(
                     f'Updating waveform generation delta_f from {wfg_domain_dict["delta_f"]} to {delta_f_new}.'
                 )
                 wfg_domain_dict["delta_f"] = delta_f_new
@@ -452,7 +449,7 @@ def sample_calibration_parameters(self, calibration_sampling_kwargs: dict):
 
         # Sample calibration parameters and calculate log_prob
         num_samples = len(self.samples)
-        print(f"Sampling calibration parameters for {num_samples} samples.")
+        logger.info(f"Sampling calibration parameters for {num_samples} samples.")
 
         delta_log_prob = np.zeros(num_samples)
 
@@ -571,7 +568,7 @@ def sample_synthetic_phase(
             self.synthetic_phase_kwargs.get("num_processes", 1), num_valid_samples // 10
         )
 
-        print(f"Estimating synthetic phase for {num_valid_samples} samples.")
+        logger.info(f"Estimating synthetic phase for {num_valid_samples} samples.")
         t0 = time.time()
 
         if not inverse:
@@ -660,7 +657,7 @@ def sample_synthetic_phase(
             self.samples["log_prob"] = log_prob_array
             del self.samples["phase"]
 
-        print(f"Done. This took {time.time() - t0:.2f} s.")
+        logger.info(f"Done. This took {time.time() - t0:.2f} s.")
 
     def get_samples_bilby_phase(self, num_processes=1):
         """
diff --git a/dingo/gw/training/train_builders.py b/dingo/gw/training/train_builders.py
index ca62856bb..2642eff3d 100755
--- a/dingo/gw/training/train_builders.py
+++ b/dingo/gw/training/train_builders.py
@@ -6,6 +6,7 @@
 from threadpoolctl import threadpool_limits
 from bilby.gw.detector import InterferometerList
 
+from dingo.core.utils.logging_utils import logger
 from dingo.gw.SVD import SVDBasis
 
 from dingo.gw.dataset.waveform_dataset import WaveformDataset
@@ -82,9 +83,9 @@ def set_train_transforms(wfd, data_settings, asd_dataset_path, omit_transforms=N
         List of sub-transforms to omit from the full composition.
     """
 
-    print(f"Setting train transforms.")
+    logger.info("Setting train transforms.")
     if omit_transforms is not None:
-        print("Omitting \n\t" + "\n\t".join([t.__name__ for t in omit_transforms]))
+        logger.info("Omitting: " + ", ".join([t.__name__ for t in omit_transforms]))
 
     # By passing the wfd domain when instantiating the noise dataset, this ensures the
     # domains will match. In particular, it truncates the ASD dataset beyond the new
@@ -142,9 +143,9 @@ def set_train_transforms(wfd, data_settings, asd_dataset_path, omit_transforms=N
     # parameters.
     try:
         standardization_dict = data_settings["standardization"]
-        print("Using previously-calculated parameter standardizations.")
+        logger.info("Using previously-calculated parameter standardizations.")
     except KeyError:
-        print("Calculating new parameter standardizations.")
+        logger.info("Calculating new parameter standardizations.")
         standardization_dict = get_standardization_dict(
             extrinsic_prior_dict,
             wfd,
@@ -258,7 +259,7 @@ def build_svd_for_embedding_network(
         ],
     )
 
-    print("Generating waveforms for embedding network SVD initialization.")
+    logger.info("Generating waveforms for embedding network SVD initialization.")
     time_start = time.time()
     ifos = list(wfd[0]["waveform"].keys())
     waveform_len = len(wfd[0]["waveform"][ifos[0]])
@@ -296,25 +297,25 @@ def build_svd_for_embedding_network(
                 waveforms[ifo][lower : lower + n] = strains[:n]
             if lower + n == num_waveforms:
                 break
-    print(f"...done. This took {time.time() - time_start:.0f} s.")
+    logger.info(f"Done. This took {time.time() - time_start:.0f} s.")
 
     # Reset the standard sharing strategy.
     torch.multiprocessing.set_sharing_strategy(old_sharing_strategy)
 
-    print("Generating SVD basis for ifo:")
+    logger.info("Generating SVD basis for ifo:")
     time_start = time.time()
     basis_dict = {}
     for ifo in ifos:
         basis = SVDBasis()
         basis.generate_basis(waveforms[ifo][:num_training_samples], size)
         basis_dict[ifo] = basis
-        print(f"...{ifo} done.")
-    print(f"...this took {time.time() - time_start:.0f} s.")
+        logger.info(f"...{ifo} done.")
+    logger.info(f"Done. This took {time.time() - time_start:.0f} s.")
 
     if out_dir is not None:
-        print(f"Testing SVD basis matrices.")
+        logger.info("Testing SVD basis matrices.")
         for ifo, basis in basis_dict.items():
-            print(f"...{ifo}:")
+            logger.info(f"...{ifo}:")
             basis.compute_test_mismatches(
                 waveforms[ifo][num_training_samples:],
                 parameters=parameters.iloc[num_training_samples:].reset_index(
@@ -323,19 +324,17 @@ def build_svd_for_embedding_network(
                 verbose=True,
             )
             basis.to_file(os.path.join(out_dir, f"svd_{ifo}.hdf5"))
-    print("Done")
+    logger.info("Done.")
 
     # Return V matrices in standard order. Drop the elements below domain.min_idx,
     # since the neural network expects data truncated below these. The dropped elements
     # should be 0.
-    print(f"Truncating SVD matrices below index {wfd.domain.min_idx}.")
-    print("...V matrix shapes:")
+    logger.info(f"Truncating SVD matrices below index {wfd.domain.min_idx}.")
     V_rb_list = []
     for ifo in data_settings["detectors"]:
         V = basis_dict[ifo].V
         assert np.allclose(V[: wfd.domain.min_idx], 0)
         V = V[wfd.domain.min_idx :]
-        print("      " + str(V.shape))
+        logger.info(f"  V matrix shape for {ifo}: {V.shape}")
         V_rb_list.append(V)
-    print("\n")
     return V_rb_list
diff --git a/dingo/gw/training/train_pipeline.py b/dingo/gw/training/train_pipeline.py
index ac1420b0f..5668c389b 100644
--- a/dingo/gw/training/train_pipeline.py
+++ b/dingo/gw/training/train_pipeline.py
@@ -27,6 +27,7 @@
     build_train_and_test_loaders,
 )
 from dingo.core.utils.trainutils import EarlyStopping
+from dingo.core.utils.logging_utils import logger, setup_logger
 from dingo.gw.dataset import WaveformDataset
 from dingo.core.posterior_models import BasePosteriorModel
 
@@ -61,20 +62,20 @@ def copy_files_to_local(
     if local_dir is not None:
         file_name = file_path.split("/")[-1]
         local_file_path = os.path.join(local_dir, file_name)
-        print(f"Copying file to {local_file_path}")
+        logger.info(f"Copying file to {local_file_path}")
         # Copy file
         start_time = time.time()
         shutil.copy(file_path, local_file_path)
         elapsed_time = time.time() - start_time
-        print("Done. This took {:2.0f}:{:2.0f} min.".format(*divmod(elapsed_time, 60)))
+        logger.info("Done. This took {:2.0f}:{:2.0f} min.".format(*divmod(elapsed_time, 60)))
     elif leave_keys_on_disk and is_condor:
-        print(
-            f"Warning: leave_waveforms_on_disk defaults to True, but local_cache_path is not specified. "
-            f"This means that the waveforms will be loaded during training from {local_file_path} ."
+        logger.warning(
+            f"leave_waveforms_on_disk defaults to True, but local_cache_path is not specified. "
+            f"This means that the waveforms will be loaded during training from {local_file_path}. "
             f"This can lead to unexpected long times for data loading during training due to network traffic. "
             f"To prevent this, specify 'local_cache_path = tmp' in the local settings or set "
             f"leave_waveforms_on_disk = False. However, the latter is not recommended for large datasets since "
-            f"it can lead to memory issues when loading the entire dataset into RAM. "
+            f"it can lead to memory issues when loading the entire dataset into RAM."
         )
 
     return local_file_path
@@ -122,7 +123,7 @@ def prepare_training_new(
 
     if train_settings["model"].get("embedding_kwargs", None):
         # First, build the SVD for seeding the embedding network.
-        print("\nBuilding SVD for initialization of embedding network.")
+        logger.info("Building SVD for initialization of embedding network.")
         initial_weights["V_rb_list"] = build_svd_for_embedding_network(
             wfd,
             train_settings["data"],
@@ -153,9 +154,8 @@ def prepare_training_new(
         "train_settings": train_settings,
     }
 
-    print("\nInitializing new posterior model.")
-    print("Complete settings:")
-    print(yaml.dump(full_settings, default_flow_style=False, sort_keys=False))
+    logger.info("Initializing new posterior model.")
+    logger.info("Complete settings:\n" + yaml.dump(full_settings, default_flow_style=False, sort_keys=False))
 
     pm = build_model_from_kwargs(
         settings=full_settings,
@@ -173,7 +173,7 @@ def prepare_training_new(
                 **local_settings["wandb"],
             )
         except ImportError:
-            print("WandB is enabled but not installed.")
+            logger.warning("WandB is enabled but not installed.")
 
     return pm, wfd
 
@@ -226,7 +226,7 @@ def prepare_training_resume(
                 **local_settings["wandb"],
             )
         except ImportError:
-            print("WandB is enabled but not installed.")
+            logger.warning("WandB is enabled but not installed.")
 
     return pm, wfd
 
@@ -278,7 +278,7 @@ def initialize_stage(
     if not resume:
         # New optimizer and scheduler. If we are resuming, these should have been
         # loaded from the checkpoint.
-        print("Initializing new optimizer and scheduler.")
+        logger.info("Initializing new optimizer and scheduler.")
         pm.optimizer_kwargs = stage["optimizer"]
         pm.scheduler_kwargs = stage["scheduler"]
         pm.initialize_optimizer_and_scheduler()
@@ -295,7 +295,7 @@ def initialize_stage(
             )
     n_grad = get_number_of_model_parameters(pm.network, (True,))
     n_nograd = get_number_of_model_parameters(pm.network, (False,))
-    print(f"Fixed parameters: {n_nograd}\nLearnable parameters: {n_grad}\n")
+    logger.info(f"Fixed parameters: {n_nograd}\nLearnable parameters: {n_grad}")
 
     return train_loader, test_loader
 
@@ -343,14 +343,12 @@ def train_stages(
         stage = stages[n]
 
         if pm.epoch == end_epochs[n] - stage["epochs"]:
-            print(f"\nBeginning training stage {n}. Settings:")
-            print(yaml.dump(stage, default_flow_style=False, sort_keys=False))
+            logger.info(f"Beginning training stage {n}. Settings:\n" + yaml.dump(stage, default_flow_style=False, sort_keys=False))
             train_loader, test_loader = initialize_stage(
                 pm, wfd, stage, local_settings["num_workers"], resume=False
             )
         else:
-            print(f"\nResuming training in stage {n}. Settings:")
-            print(yaml.dump(stage, default_flow_style=False, sort_keys=False))
+            logger.info(f"Resuming training in stage {n}. Settings:\n" + yaml.dump(stage, default_flow_style=False, sort_keys=False))
             train_loader, test_loader = initialize_stage(
                 pm, wfd, stage, local_settings["num_workers"], resume=True
             )
@@ -359,7 +357,7 @@ def train_stages(
             try:
                 early_stopping = EarlyStopping(**stage["early_stopping"])
             except Exception:
-                print(
+                logger.warning(
                     "Early stopping settings invalid. Please pass 'patience', 'delta', 'metric'"
                 )
                 raise
@@ -381,10 +379,10 @@ def train_stages(
 
         if pm.epoch == end_epochs[n]:
             save_file = os.path.join(train_dir, f"model_stage_{n}.pt")
-            print(f"Training stage complete. Saving to {save_file}.")
+            logger.info(f"Training stage complete. Saving to {save_file}.")
             pm.save_model(save_file, save_training_info=True)
         if runtime_limits.local_limits_exceeded(pm.epoch):
-            print("Local runtime limits reached. Ending program.")
+            logger.info("Local runtime limits reached. Ending program.")
             break
 
     if pm.epoch == end_epochs[-1]:
@@ -441,9 +439,10 @@ def train_local():
     args = parse_args()
 
     os.makedirs(args.train_dir, exist_ok=True)
+    setup_logger(outdir=args.train_dir, label="dingo_train")
 
     if args.settings_file is not None:
-        print("Beginning new training run.")
+        logger.info("Beginning new training run.")
         with open(args.settings_file, "r") as fp:
             train_settings = yaml.safe_load(fp)
 
@@ -462,7 +461,7 @@ def train_local():
 
                     local_settings["wandb"]["id"] = wandb.util.generate_id()
                 except ImportError:
-                    print("wandb not installed, cannot generate run id.")
+                    logger.warning("wandb not installed, cannot generate run id.")
             yaml.dump(local_settings, f, default_flow_style=False, sort_keys=False)
 
         pm, wfd = prepare_training_new(train_settings, args.train_dir, local_settings)
@@ -470,7 +469,7 @@ def train_local():
     else:
         if not os.path.isfile(args.checkpoint):
             raise FileNotFoundError(f"Checkpoint not found: {args.checkpoint}")
-        print("Resuming training run.")
+        logger.info("Resuming training run.")
         with open(os.path.join(args.train_dir, "local_settings.yaml"), "r") as f:
             local_settings = yaml.safe_load(f)
         pm, wfd = prepare_training_resume(
@@ -482,11 +481,11 @@ def train_local():
 
     if complete:
         if args.exit_command:
-            print(
+            logger.info(
                 f"All training stages complete. Executing exit command: {args.exit_command}."
             )
             os.system(args.exit_command)
         else:
-            print("All training stages complete.")
+            logger.info("All training stages complete.")
     else:
-        print("Program terminated due to runtime limit.")
+        logger.info("Program terminated due to runtime limit.")
diff --git a/dingo/gw/training/train_pipeline_condor.py b/dingo/gw/training/train_pipeline_condor.py
index 418909d3e..06550f5a2 100644
--- a/dingo/gw/training/train_pipeline_condor.py
+++ b/dingo/gw/training/train_pipeline_condor.py
@@ -4,6 +4,7 @@
 import yaml
 import argparse
 
+from dingo.core.utils.logging_utils import logger, setup_logger
 from dingo.gw.training import (
     prepare_training_new,
     prepare_training_resume,
@@ -61,7 +62,7 @@ def copy_logfiles(log_dir, epoch, name="info", suffixes=(".err", ".log", ".out")
         try:
             copyfile(src, dest)
         except:
-            print("Could not copy " + src)
+            logger.warning("Could not copy " + src)
 
 
 def train_condor():
@@ -78,6 +79,7 @@ def train_condor():
         help="Optional command to execute after completion of training.",
     )
     args = parser.parse_args()
+    setup_logger(outdir=args.train_dir, label="dingo_train")
 
     if not args.start_submission:
 
@@ -90,7 +92,7 @@ def train_condor():
                 raise FileNotFoundError(
                     f"Checkpoint not found: {join(args.train_dir, args.checkpoint)}"
                 )
-            print("Beginning new training run.")
+            logger.info("Beginning new training run.")
             with open(join(args.train_dir, "train_settings.yaml"), "r") as f:
                 train_settings = yaml.safe_load(f)
 
@@ -109,7 +111,7 @@ def train_condor():
 
                         local_settings["wandb"]["id"] = wandb.util.generate_id()
                     except ImportError:
-                        print("wandb not installed, cannot generate run id.")
+                        logger.warning("wandb not installed, cannot generate run id.")
                 yaml.dump(local_settings, f, default_flow_style=False, sort_keys=False)
 
             pm, wfd = prepare_training_new(
@@ -117,7 +119,7 @@ def train_condor():
             )
 
         else:
-            print("Resuming training run.")
+            logger.info("Resuming training run.")
             with open(os.path.join(args.train_dir, "local_settings.yaml"), "r") as f:
                 local_settings = yaml.safe_load(f)
             pm, wfd = prepare_training_resume(
@@ -128,7 +130,7 @@ def train_condor():
 
         complete = train_stages(pm, wfd, args.train_dir, local_settings)
 
-        print("Copying log files")
+        logger.info("Copying log files")
         copy_logfiles(args.train_dir, epoch=pm.epoch)
 
         #
@@ -136,7 +138,7 @@ def train_condor():
         #
 
         if complete:
-            print(
+            logger.info(
                 f"Training complete, job will not be resubmitted. Executing exit command: {args.exit_command}."
             )
             if args.exit_command:
diff --git a/dingo/gw/training/utils.py b/dingo/gw/training/utils.py
index 4cd0732d8..2aeedc964 100644
--- a/dingo/gw/training/utils.py
+++ b/dingo/gw/training/utils.py
@@ -5,6 +5,7 @@
 import yaml
 
 from dingo.core.utils.backward_compatibility import torch_load_with_fallback
+from dingo.core.utils.logging_utils import logger
 
 
 def append_stage():
@@ -25,7 +26,7 @@ def append_stage():
         if k.startswith("stage_")
     ]
     num_stages = len(stages)
-    print(f"Checkpoint training plan consists of {num_stages} stages.")
+    logger.info(f"Checkpoint training plan consists of {num_stages} stages.")
 
     with open(args.stage_settings_file, "r") as f:
         new_stage = yaml.safe_load(f)
@@ -39,21 +40,20 @@ def append_stage():
         current_epoch = d["epoch"]
         stage_epoch = np.sum([s["epochs"] for s in stages[: args.replace]])
         if current_epoch > stage_epoch:
-            print(
-                f"WARNING: Modification to training plan changes a training stage "
+            logger.warning(
+                f"Modification to training plan changes a training stage "
                 f"that has already started. Current model epoch is {current_epoch}. "
                 f"Proceed at your own risk!"
             )
-        print(f"Replacing planned stage {args.replace} with new stage.")
+        logger.info(f"Replacing planned stage {args.replace} with new stage.")
         new_stage_number = args.replace
     else:
-        print(f"Appending new stage to training plan.")
+        logger.info("Appending new stage to training plan.")
         new_stage_number = num_stages
 
     d["metadata"]["train_settings"]["training"][f"stage_{new_stage_number}"] = new_stage
-    print("Summary of new training plan:")
-    print(
-        yaml.dump(
+    logger.info(
+        "Summary of new training plan:\n" + yaml.dump(
             d["metadata"]["train_settings"]["training"],
             default_flow_style=False,
             sort_keys=False,
diff --git a/dingo/gw/transforms/waveform_transforms.py b/dingo/gw/transforms/waveform_transforms.py
index f6d022435..39f884836 100644
--- a/dingo/gw/transforms/waveform_transforms.py
+++ b/dingo/gw/transforms/waveform_transforms.py
@@ -1,6 +1,7 @@
 from typing import Optional
 import numpy as np
 
+from dingo.core.utils.logging_utils import logger
 from dingo.gw.domains import MultibandedFrequencyDomain, UniformFrequencyDomain
 
 
@@ -443,7 +444,7 @@ def __init__(
         self.maximum_frequency = maximum_frequency
 
         if print_output:
-            print(
+            logger.info(
                 f"Transform MaskDataForFrequencyRangeUpdate activated:"
                 f"  Settings: \n"
                 f"    - Minimum_frequency update: {self.minimum_frequency}\n"
diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py
index 43f6a6488..a7a341ee4 100644
--- a/dingo/gw/waveform_generator/waveform_generator.py
+++ b/dingo/gw/waveform_generator/waveform_generator.py
@@ -30,6 +30,7 @@
     TimeDomain,
 )
 from dingo.gw.transforms.waveform_transforms import DecimateAll
+from dingo.core.utils.logging_utils import logger
 
 
 class WaveformGenerator:
@@ -142,12 +143,12 @@ def spin_conversion_phase(self):
     @spin_conversion_phase.setter
     def spin_conversion_phase(self, value):
         if value is None:
-            print(
+            logger.info(
                 "Setting spin_conversion_phase = None. Using phase parameter for "
                 "conversion to cartesian spins."
             )
         else:
-            print(
+            logger.info(
                 f"Setting spin_conversion_phase = {value}. Using this value for the "
                 f"phase parameter for conversion to cartesian spins."
             )
@@ -583,7 +584,7 @@ def generate_FD_waveform(
         # numbers if multibanding is used. If that happens, turn off multibanding to
         # fix this.
         if max(np.max(np.abs(hp.data.data)), np.max(np.abs(hc.data.data))) > 1e-20:
-            print(
+            logger.warning(
                 f"Generation with parameters {parameters_lal} likely numerically "
                 f"unstable due to multibanding, turn off multibanding."
             )
@@ -606,7 +607,7 @@ def generate_FD_waveform(
                 *parameters_lal[lal_dict_idx + 1 :],
             )
             if max(np.max(np.abs(hp.data.data)), np.max(np.abs(hc.data.data))) > 1e-20:
-                print(
+                logger.warning(
                     f"Warning: turning off multibanding for parameters {parameters_lal}"
                     f" likely numerically might not have fixed it, check manually."
                 )
diff --git a/dingo/pipe/dag_creator.py b/dingo/pipe/dag_creator.py
index 077864265..1223a553b 100644
--- a/dingo/pipe/dag_creator.py
+++ b/dingo/pipe/dag_creator.py
@@ -5,7 +5,9 @@
 import copy
 
 from bilby_pipe.job_creation.dag import Dag
-from bilby_pipe.utils import BilbyPipeError, logger
+import logging
+
+from bilby_pipe.utils import BilbyPipeError
 
 from dingo.pipe.utils import _strip_unwanted_submission_keys
 from dingo.pipe.nodes.generation_node import GenerationNode
@@ -15,7 +17,7 @@
 from .nodes.plot_node import PlotNode, PlotPPNode
 from .nodes.sampling_node import SamplingNode
 
-logger.name = "dingo_pipe"
+logger = logging.getLogger("dingo.pipe")
 
 
 def get_trigger_time_list(inputs):
diff --git a/dingo/pipe/data_generation.py b/dingo/pipe/data_generation.py
index c8117027a..e465b3865 100644
--- a/dingo/pipe/data_generation.py
+++ b/dingo/pipe/data_generation.py
@@ -10,12 +10,12 @@
 from bilby.core.prior import PriorDict
 from bilby_pipe.utils import (
     parse_args,
-    logger,
     convert_string_to_dict,
     convert_prior_string_input,
     resolve_filename_with_transfer_fallback,
     BilbyPipeError,
 )
+import logging
 import lalsimulation as LS
 import numpy as np
 
@@ -24,8 +24,9 @@
 from dingo.gw.domains import UniformFrequencyDomain, build_domain_from_model_metadata
 from dingo.gw.injection import Injection
 from dingo.pipe.parser import create_parser
+from dingo.core.utils.logging_utils import setup_logger
 
-logger.name = "dingo_pipe"
+logger = logging.getLogger("dingo.pipe")
 
 
 class DataGenerationInput(BilbyDataGenerationInput):
@@ -534,6 +535,7 @@ def create_generation_parser():
 def main():
     """Data generation main logic"""
     args, unknown_args = parse_args(sys.argv[1:], create_generation_parser())
+    setup_logger(outdir=args.outdir, label=args.label)
     # log_version_information()
     data = DataGenerationInput(args, unknown_args)
     data.save_hdf5()
diff --git a/dingo/pipe/dingo_result.py b/dingo/pipe/dingo_result.py
index 4b0096e13..b06c1dfdd 100644
--- a/dingo/pipe/dingo_result.py
+++ b/dingo/pipe/dingo_result.py
@@ -1,8 +1,11 @@
 import argparse
+import logging
 from pathlib import Path
 
 from dingo.gw.result import Result
 
+logger = logging.getLogger("dingo.pipe")
+
 
 def main():
     parser = argparse.ArgumentParser()
@@ -15,7 +18,7 @@ def main():
     args = parser.parse_args()
 
     if args.merge:
-        print(f"Merging {len(args.result)} parts into complete Result.")
+        logger.info(f"Merging {len(args.result)} parts into complete Result.")
         sub_results = []
         for file_name in args.result:
             sub_results.append(Result(file_name=file_name))
diff --git a/dingo/pipe/importance_sampling.py b/dingo/pipe/importance_sampling.py
index 541e8ebb2..bb195a647 100644
--- a/dingo/pipe/importance_sampling.py
+++ b/dingo/pipe/importance_sampling.py
@@ -8,12 +8,14 @@
 from bilby_pipe.input import Input
 from bilby_pipe.utils import (
     parse_args,
-    logger,
     convert_string_to_dict,
     convert_prior_string_input,
     resolve_filename_with_transfer_fallback,
     BilbyPipeError,
 )
+import logging
+from dingo.core.utils.logging_utils import setup_logger
+logger = logging.getLogger("dingo.pipe")
 
 from dingo.gw.data.event_dataset import EventDataset
 from dingo.gw.domains import MultibandedFrequencyDomain
@@ -21,8 +23,6 @@
 from dingo.pipe.parser import create_parser
 from dingo.gw.result import Result
 
-logger.name = "dingo_pipe"
-
 
 class ImportanceSamplingInput(Input):
     def __init__(self, args, unknown_args):
@@ -285,6 +285,7 @@ def create_sampling_parser():
 def main():
     """Data analysis main logic"""
     args, unknown_args = parse_args(sys.argv[1:], create_sampling_parser())
+    setup_logger(outdir=args.outdir, label=args.label)
     # log_version_information()
     analysis = ImportanceSamplingInput(args, unknown_args)
     analysis.run_sampler()
diff --git a/dingo/pipe/main.py b/dingo/pipe/main.py
index 1c52cedcf..9bdb00bd1 100644
--- a/dingo/pipe/main.py
+++ b/dingo/pipe/main.py
@@ -16,12 +16,15 @@
     convert_string_to_dict,
     get_colored_string,
     get_command_line_arguments,
-    logger,
     parse_args,
     convert_prior_string_input,
     BilbyPipeError,
 )
 
+import logging
+from dingo.core.utils.logging_utils import setup_logger
+logger = logging.getLogger("dingo.pipe")
+
 from .dag_creator import generate_dag
 from .parser import create_parser
 from .utils import dict_to_string
@@ -618,6 +621,8 @@ def main():
     parser = create_parser(top_level=True)
     args, unknown_args = parse_args(get_command_line_arguments(), parser)
 
+    setup_logger(outdir=args.outdir, label=args.label)
+
     importance_sampling_updates, model_args = fill_in_arguments_from_model(args)
     inputs = MainInput(args, unknown_args, importance_sampling_updates)
     write_complete_config_file(parser, args, inputs)
diff --git a/dingo/pipe/nodes/generation_node.py b/dingo/pipe/nodes/generation_node.py
index edbd4d9b2..6f7cfcc57 100644
--- a/dingo/pipe/nodes/generation_node.py
+++ b/dingo/pipe/nodes/generation_node.py
@@ -1,3 +1,4 @@
+import logging
 import os
 
 from bilby_pipe.job_creation.nodes import GenerationNode as BilbyGenerationNode
@@ -5,6 +6,8 @@
 
 from dingo.pipe.utils import _strip_unwanted_submission_keys
 
+logger = logging.getLogger("dingo.pipe")
+
 
 class GenerationNode(BilbyGenerationNode):
 
@@ -101,7 +104,6 @@ def __init__(self, inputs, trigger_time, idx, dag, parent=None, importance_sampl
             for fname in self.extract_paths_from_dict(frame_files)
             if fname.startswith(self.inputs.data_find_urltype)
         ]:
-            from bilby_pipe.utils import logger
             logger.warning(
                 "The following frame files were identified by gwdatafind for this analysis. "
                 "These frames may not be found by the data generation stage as file "
diff --git a/dingo/pipe/nodes/pe_summary_node.py b/dingo/pipe/nodes/pe_summary_node.py
index 6762ece07..876e7a0d1 100644
--- a/dingo/pipe/nodes/pe_summary_node.py
+++ b/dingo/pipe/nodes/pe_summary_node.py
@@ -1,8 +1,10 @@
 import os 
 from bilby_pipe.job_creation.nodes import PESummaryNode as BilbyPESummaryNode
-from bilby_pipe.utils import BilbyPipeError, logger
+import logging
 
-logger.name = "dingo_pipe"
+from bilby_pipe.utils import BilbyPipeError
+
+logger = logging.getLogger("dingo.pipe")
 
 
 class PESummaryNode(BilbyPESummaryNode):
diff --git a/dingo/pipe/plot.py b/dingo/pipe/plot.py
index 7dc8af756..772fa68f9 100644
--- a/dingo/pipe/plot.py
+++ b/dingo/pipe/plot.py
@@ -1,9 +1,12 @@
 from bilby_pipe.bilbyargparser import BilbyArgParser
-from bilby_pipe.utils import get_command_line_arguments, logger, parse_args
+from bilby_pipe.utils import get_command_line_arguments, parse_args
 
+import logging
+
+from dingo.core.utils.logging_utils import setup_logger
 from dingo.gw.result import Result
 
-logger.name = "dingo_pipe"
+logger = logging.getLogger("dingo.pipe")
 
 
 def create_parser():
@@ -46,6 +49,7 @@ def create_parser():
 
 def main():
     args, _ = parse_args(get_command_line_arguments(), create_parser())
+    setup_logger(outdir=args.outdir, label=args.label)
 
     logger.info(f"Generating plots for results file {args.result}")
 
diff --git a/dingo/pipe/pp_test.py b/dingo/pipe/pp_test.py
index 0716b24e6..bbc85502d 100644
--- a/dingo/pipe/pp_test.py
+++ b/dingo/pipe/pp_test.py
@@ -1,12 +1,12 @@
 import argparse
 from pathlib import Path
 
-from bilby_pipe.utils import logger
+import logging
 
 from dingo.core.result import make_pp_plot
 from dingo.gw.result import Result
 
-logger.name = "dingo_pipe"
+logger = logging.getLogger("dingo.pipe")
 
 
 def create_parser():
diff --git a/dingo/pipe/sampling.py b/dingo/pipe/sampling.py
index 9138066bc..e038e7022 100644
--- a/dingo/pipe/sampling.py
+++ b/dingo/pipe/sampling.py
@@ -7,10 +7,12 @@
 from bilby_pipe.input import Input
 from bilby_pipe.utils import (
     parse_args,
-    logger,
     convert_string_to_dict,
     resolve_filename_with_transfer_fallback,
 )
+import logging
+from dingo.core.utils.logging_utils import setup_logger
+logger = logging.getLogger("dingo.pipe")
 
 from dingo.core.posterior_models.build_model import build_model_from_kwargs
 from dingo.gw.data.event_dataset import EventDataset
@@ -216,6 +218,7 @@ def create_sampling_parser():
 def main():
     """Data analysis main logic"""
     args, unknown_args = parse_args(sys.argv[1:], create_sampling_parser())
+    setup_logger(outdir=args.outdir, label=args.label)
     # log_version_information()
     analysis = SamplingInput(args, unknown_args)
     analysis.run_sampler()