diff --git a/.gitignore b/.gitignore
index b68aa37..c569ea6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -13,3 +13,9 @@ Pipfile
 **/.venv
 **/.vim
 pyvenv.cfg
+
+build/
+dist/
+
+build/
+dist/
diff --git a/mla/analysis.py b/mla/analysis.py
index 423bad8..1ae94e1 100644
--- a/mla/analysis.py
+++ b/mla/analysis.py
@@ -29,9 +29,12 @@ class SingleSourceLLHAnalysis:
     sob_term_factories: List[SoBTermFactory]
     data_handler_source: Tuple[DataHandler, PointSource]
     _sob_term_factories: List[SoBTermFactory] = field(init=False, repr=False)
-    _data_handler_source: Tuple[DataHandler, PointSource] = field(init=False, repr=False)
-    _trial_generator: SingleSourceTrialGenerator = field(init=False, repr=False)
-    _test_statistic_factory: LLHTestStatisticFactory = field(init=False, repr=False)
+    _data_handler_source: Tuple[DataHandler,
+                                PointSource] = field(init=False, repr=False)
+    _trial_generator: SingleSourceTrialGenerator = field(
+        init=False, repr=False)
+    _test_statistic_factory: LLHTestStatisticFactory = field(
+        init=False, repr=False)
 
     def produce_and_minimize(
         self,
@@ -66,10 +69,12 @@ def sob_term_factories(self) -> List[SoBTermFactory]:
         return self._sob_term_factories
 
     @sob_term_factories.setter
-    def sob_term_factories(self, sob_term_factories: List[SoBTermFactory]) -> None:
+    def sob_term_factories(
+            self,
+            sob_term_factories: List[SoBTermFactory]) -> None:
         """Docstring"""
         self._sob_term_factories = sob_term_factories
-        self._test_statistic_factory = LLHTestStatisticFactory(  # pylint: disable=too-many-function-args
+        self._test_statistic_factory = LLHTestStatisticFactory(
             self.config['LLHTestStatisticFactory'],
             self._sob_term_factories,
         )
@@ -80,8 +85,9 @@ def data_handler_source(self) -> Tuple[DataHandler, PointSource]:
         return self._data_handler_source
 
     @data_handler_source.setter
-    def data_handler_source(
-            self, data_handler_source: Tuple[DataHandler, PointSource]) -> None:
+    def data_handler_source(self,
+                            data_handler_source: Tuple[DataHandler,
+                                                       PointSource]) -> None:
         """Docstring"""
         self._data_handler_source = data_handler_source
         self._trial_generator = SingleSourceTrialGenerator(
diff --git a/mla/core.py b/mla/core.py
index d79be03..13508ff 100644
--- a/mla/core.py
+++ b/mla/core.py
@@ -6,4 +6,5 @@
 def generate_default_config(classes: list) -> dict:
     """Docstring"""
     return {
-        c.__name__: c.generate_config() for c in classes if issubclass(c, Configurable)}
+        c.__name__: c.generate_config()
+        for c in classes if issubclass(c, Configurable)}
diff --git a/mla/data_handlers.py b/mla/data_handlers.py
index 8fa2b31..e4a7cb9 100644
--- a/mla/data_handlers.py
+++ b/mla/data_handlers.py
@@ -31,7 +31,10 @@ class DataHandler(configurable.Configurable):
     __metaclass__ = abc.ABCMeta
 
     @abc.abstractmethod
-    def sample_background(self, n: int, rng: np.random.Generator) -> np.ndarray:
+    def sample_background(
+            self,
+            n: int,
+            rng: np.random.Generator) -> np.ndarray:
         """Docstring"""
 
     @abc.abstractmethod
@@ -47,7 +50,8 @@ def evaluate_background_sindec_pdf(self, events: np.ndarray) -> np.ndarray:
         """Docstring"""
 
     @abc.abstractmethod
-    def build_background_sindec_logenergy_histogram(self, bins: np.ndarray) -> np.ndarray:
+    def build_background_sindec_logenergy_histogram(
+            self, bins: np.ndarray) -> np.ndarray:
         """Docstring"""
 
     @abc.abstractmethod
@@ -77,7 +81,10 @@ class NuSourcesDataHandler(DataHandler):
     _livetime: float = field(init=False, repr=False)
     _sin_dec_bins: np.ndarray = field(init=False, repr=False)
 
-    def sample_background(self, n: int, rng: np.random.Generator) -> np.ndarray:
+    def sample_background(
+            self,
+            n: int,
+            rng: np.random.Generator) -> np.ndarray:
         """Docstring"""
         return rng.choice(self._data, n)
 
@@ -105,7 +112,8 @@ def evaluate_background_sindec_pdf(self, events: np.ndarray) -> np.ndarray:
         """
         return (1 / (2 * np.pi)) * self._dec_spline(events['sindec'])
 
-    def build_background_sindec_logenergy_histogram(self, bins: np.ndarray) -> np.ndarray:
+    def build_background_sindec_logenergy_histogram(
+            self, bins: np.ndarray) -> np.ndarray:
         """Docstring"""
         return np.histogram2d(
             self._data['sindec'],
@@ -188,14 +196,18 @@ def _cut_sim_dec(self) -> None:
             self._sim = self._full_sim[close].copy()
 
             self._sim['ow'] /= 2 * np.pi * (np.min([np.sin(
-                self.config['dec_cut_location'] + self.config['dec_bandwidth (rad)']
+                self.config['dec_cut_location'] +
+                self.config['dec_bandwidth (rad)']
             ), 1]) - np.max([np.sin(
-                self.config['dec_cut_location'] - self.config['dec_bandwidth (rad)']
+                self.config['dec_cut_location'] -
+                self.config['dec_bandwidth (rad)']
             ), -1]))
             self._sim['weight'] /= 2 * np.pi * (np.min([np.sin(
-                self.config['dec_cut_location'] + self.config['dec_bandwidth (rad)']
+                self.config['dec_cut_location'] +
+                self.config['dec_bandwidth (rad)']
             ), 1]) - np.max([np.sin(
-                self.config['dec_cut_location'] - self.config['dec_bandwidth (rad)']
+                self.config['dec_cut_location'] -
+                self.config['dec_bandwidth (rad)']
             ), -1]))
         else:
             self._sim = self._full_sim
@@ -208,7 +220,8 @@ def data_grl(self) -> Tuple[np.ndarray, np.ndarray]:
     @data_grl.setter
     def data_grl(self, data_grl: Tuple[np.ndarray, np.ndarray]) -> None:
         """Docstring"""
-        self._sin_dec_bins = np.linspace(-1, 1, 1 + self.config['sin_dec_bins'])
+        self._sin_dec_bins = np.linspace(-1, 1,
+                                         1 + self.config['sin_dec_bins'])
         self._data = data_grl[0].copy()
         self._grl = data_grl[1].copy()
         if 'sindec' not in self._data.dtype.names:
@@ -232,7 +245,8 @@ def data_grl(self, data_grl: Tuple[np.ndarray, np.ndarray]) -> None:
             self._data['sindec'], bins=self._sin_dec_bins, density=True)
         bin_centers = bins[:-1] + np.diff(bins) / 2
 
-        self._dec_spline = Spline(bin_centers, hist, **self.config['dec_spline_kwargs'])
+        self._dec_spline = Spline(
+            bin_centers, hist, **self.config['dec_spline_kwargs'])
 
     @property
     def livetime(self) -> float:
@@ -301,7 +315,8 @@ def background_time_profile(self, profile: GenericProfile) -> None:
         background_grl = self._grl[background_run_mask]
         self._n_background = background_grl['events'].sum()
         self._n_background /= background_grl['livetime'].sum()
-        self._n_background *= self._contained_livetime(*profile.range, background_grl)
+        self._n_background *= self._contained_livetime(
+            *profile.range, background_grl)
         self._background_time_profile = copy.deepcopy(profile)
 
     @property
@@ -362,31 +377,38 @@ def _contained_livetime(
 
         return contained_livetime
 
-    def sample_background(self, n: int, rng: np.random.Generator) -> np.ndarray:
+    def sample_background(
+            self,
+            n: int,
+            rng: np.random.Generator) -> np.ndarray:
         """Docstring"""
         events = super().sample_background(n, rng)
         return self._randomize_times(events, self._background_time_profile)
 
     def sample_signal(self, n: int, rng: np.random.Generator) -> np.ndarray:
         """Docstring"""
-        events = super().sample_signal(n, rng)
-        return self._randomize_times(events, self._signal_time_profile)
+        self.events = super().sample_signal(n, rng)
+        return self._randomize_times(self.events, self._signal_time_profile)
 
     def _randomize_times(
         self,
         events: np.ndarray,
         time_profile: GenericProfile,
     ) -> np.ndarray:
-        grl_start_cdf = time_profile.cdf(self._grl['start'])
-        grl_stop_cdf = time_profile.cdf(self._grl['stop'])
-        valid = np.logical_and(grl_start_cdf < 1, grl_stop_cdf > 0)
-        rates = grl_stop_cdf[valid] - grl_start_cdf[valid]
-
+        self.grl_start_cdf = time_profile.cdf(self._grl['start'])
+        self.grl_stop_cdf = time_profile.cdf(self._grl['stop'])
+        self.valid = np.logical_and(
+            self.grl_start_cdf <= 1,
+            self.grl_stop_cdf >= 0)
+        self.rates = self.grl_stop_cdf[self.valid] - \
+            self.grl_start_cdf[self.valid]
+        # print(self.rates, self.rates.sum())
+        # print(len(events))
         runs = np.random.choice(
-            self._grl[valid],
+            self._grl[self.valid],
             size=len(events),
             replace=True,
-            p=rates / rates.sum(),
+            p=self.rates / self.rates.sum(),
         )
 
         events['time'] = time_profile.inverse_transform_sample(
diff --git a/mla/minimizers.py b/mla/minimizers.py
index f91f973..2ffdb14 100644
--- a/mla/minimizers.py
+++ b/mla/minimizers.py
@@ -33,26 +33,31 @@ class Minimizer(configurable.Configurable):
 
     @abc.abstractmethod
     def __call__(
-            self, fitting_params: Optional[List[str]] = None) -> Tuple[float, np.ndarray]:
+        self, fitting_params: Optional[List[str]] = None
+    ) -> Tuple[float, np.ndarray]:
         """Docstring"""
 
 
 @dataclasses.dataclass
 class GridSearchMinimizer(Minimizer):
     """Docstring"""
+
     def __call__(
-            self, fitting_params: Optional[List[str]] = None) -> Tuple[float, np.ndarray]:
+        self, fitting_params: Optional[List[str]] = None
+    ) -> Tuple[float, np.ndarray]:
         """Docstring"""
         if fitting_params is None:
             fitting_key_idx_map = self.test_statistic.params.key_idx_map
             fitting_bounds = self.test_statistic.params.bounds
         else:
             fitting_key_idx_map = {
-                key: val for key, val in self.test_statistic.params.key_idx_map.items()
+                key: val
+                for key, val in self.test_statistic.params.key_idx_map.items()
                 if key in fitting_params
             }
             fitting_bounds = {
-                key: val for key, val in self.test_statistic.params.bounds.items()
+                key: val
+                for key, val in self.test_statistic.params.bounds.items()
                 if key in fitting_params
             }
 
@@ -72,13 +77,22 @@ def __call__(
         ])
 
         return self._minimize(
-            points[grid_ts_values.argmin()], fitting_key_idx_map, fitting_bounds)
+            points[grid_ts_values.argmin()],
+            fitting_key_idx_map, fitting_bounds)
 
-    def _eval_test_statistic(self, point: np.ndarray, fitting_key_idx_map: dict) -> float:
+    def _eval_test_statistic(
+            self,
+            point: np.ndarray,
+            fitting_key_idx_map: dict) -> float:
         """Docstring"""
-        return self.test_statistic(self._param_values(point, fitting_key_idx_map))
-
-    def _param_values(self, point: np.ndarray, fitting_key_idx_map: dict) -> np.ndarray:
+        return self.test_statistic(
+            self._param_values(
+                point, fitting_key_idx_map))
+
+    def _param_values(
+            self,
+            point: np.ndarray,
+            fitting_key_idx_map: dict) -> np.ndarray:
         """Docstring"""
         param_values = self.test_statistic.params.value_array.copy()
 
diff --git a/mla/params.py b/mla/params.py
index ad7fadf..c349ca6 100644
--- a/mla/params.py
+++ b/mla/params.py
@@ -37,7 +37,10 @@ def names(self) -> List[str]:
         return [*self.key_idx_map]
 
     @classmethod
-    def from_dict(cls, value_dict: dict, bounds: Optional[dict] = None) -> 'Params':
+    def from_dict(
+            cls,
+            value_dict: dict,
+            bounds: Optional[dict] = None) -> 'Params':
         """Docstring"""
         value_array = np.array(list(value_dict.values()))
         key_idx_map = {key: i for i, key in enumerate(value_dict)}
@@ -50,8 +53,11 @@ def from_array(
         bounds: Optional[dict] = None,
     ) -> 'Params':
         """Docstring"""
-        value_array = rf.structured_to_unstructured(named_value_array, copy=True)[0]
-        key_idx_map = {name: i for i, name in enumerate(named_value_array.dtype.names)}
+        value_array = rf.structured_to_unstructured(
+            named_value_array, copy=True)[0]
+        key_idx_map = {
+            name: i for i, name in enumerate(
+                named_value_array.dtype.names)}
         return cls._build_params(value_array, key_idx_map, bounds)
 
     @classmethod
diff --git a/mla/sob_terms.py b/mla/sob_terms.py
index 5a022e2..09a8dd5 100644
--- a/mla/sob_terms.py
+++ b/mla/sob_terms.py
@@ -173,7 +173,10 @@ def calculate_drop_mask(self, events: np.ndarray) -> np.ndarray:
         return 1 / self.background_time_profile.pdf(events['time']) != 0
 
     def generate_params(self) -> tuple:
-        return self.signal_time_profile.params, self.signal_time_profile.param_bounds
+        return (
+               self.signal_time_profile.params,
+               self.signal_time_profile.param_bounds
+               )
 
 
 @dataclasses.dataclass
@@ -233,8 +236,10 @@ def __post_init__(self) -> None:
 
     def __call__(self, params: Params, events: np.ndarray) -> SoBTerm:
         """Docstring"""
-        sin_dec_idx = np.searchsorted(self._sin_dec_bins[:-1], events['sindec'])
-        log_energy_idx = np.searchsorted(self._log_energy_bins[:-1], events['logE'])
+        sin_dec_idx = np.searchsorted(
+            self._sin_dec_bins[:-1], events['sindec'])
+        log_energy_idx = np.searchsorted(
+            self._log_energy_bins[:-1], events['logE'])
 
         spline_idxs, event_spline_idxs = np.unique(
             [sin_dec_idx - 1, log_energy_idx - 1],
@@ -276,7 +281,8 @@ def _init_sob_map(
             An array of signal-over-background values binned in sin(dec) and
             log(energy) for a given gamma.
         """
-        sig_h = self.data_handler.build_signal_sindec_logenergy_histogram(gamma, bins)
+        sig_h = self.data_handler.build_signal_sindec_logenergy_histogram(
+            gamma, bins)
 
         # Normalize histogram by dec band
         sig_h /= np.sum(sig_h, axis=1)[:, None]
@@ -292,7 +298,10 @@ def _init_sob_map(
             good_bins, good_vals = bin_centers[good], ratio[i][good]
 
             # Do a linear interpolation across the energy range
-            spline = Spline(good_bins, good_vals, **self.config['energy_spline_kwargs'])
+            spline = Spline(
+                good_bins,
+                good_vals,
+                **self.config['energy_spline_kwargs'])
 
             # And store the interpolated values
             ratio[i] = spline(bin_centers)
@@ -306,7 +315,8 @@ def _init_spline_map(self) -> List[List[Spline]]:
         """
         bins = np.array([self._sin_dec_bins, self._log_energy_bins])
         bin_centers = bins[1, :-1] + np.diff(bins[1]) / 2
-        bg_h = self.data_handler.build_background_sindec_logenergy_histogram(bins)
+        bg_h = self.data_handler.build_background_sindec_logenergy_histogram(
+            bins)
 
         # Normalize histogram by dec band
         bg_h /= np.sum(bg_h, axis=1)[:, None]
@@ -323,7 +333,8 @@ def _init_spline_map(self) -> List[List[Spline]]:
         transposed_log_sob_maps = np.log(sob_maps.transpose(1, 2, 0))
 
         splines = [[
-            Spline(self._gamma_bins, log_ratios, **self.config['sob_spline_kwargs'])
+            Spline(self._gamma_bins, log_ratios,
+                   **self.config['sob_spline_kwargs'])
             for log_ratios in dec_bin
         ] for dec_bin in transposed_log_sob_maps]
 
diff --git a/mla/sources.py b/mla/sources.py
index ee937d7..686c74f 100644
--- a/mla/sources.py
+++ b/mla/sources.py
@@ -22,13 +22,18 @@
 @dataclasses.dataclass
 class PointSource(configurable.Configurable):
     """Stores a source object name and location"""
+
     def sample(self, size: int = 1) -> tuple:
         """Sample locations.
 
         Args:
             size: number of points to sample
         """
-        return (np.ones(size) * self.config['ra'], np.ones(size) * self.config['dec'])
+        return (
+            np.ones(size) *
+            self.config['ra'],
+            np.ones(size) *
+            self.config['dec'])
 
     def spatial_pdf(self, events: np.ndarray) -> np.ndarray:
         """calculates the signal probability of events.
@@ -72,6 +77,7 @@ def generate_config(cls) -> dict:
 @dataclasses.dataclass
 class GaussianExtendedSource(PointSource):
     """Gaussian Extended Source"""
+
     def sample(self, size: int = 1) -> np.ndarray:
         """Sample locations.
 
diff --git a/mla/threeml/IceCubeLike.py b/mla/threeml/IceCubeLike.py
index a686250..0fc40ac 100644
--- a/mla/threeml/IceCubeLike.py
+++ b/mla/threeml/IceCubeLike.py
@@ -1,3 +1,5 @@
+from __future__ import print_function
+from __future__ import division
 """Docstring"""
 
 __author__ = 'John Evans and Jason Fan'
@@ -9,8 +11,8 @@
 __email__ = 'klfan@terpmail.umd.edu'
 __status__ = 'Development'
 
-from __future__ import print_function
-from __future__ import division
+# from __future__ import print_function
+# from __future__ import division
 from past.utils import old_div
 import collections
 import scipy
@@ -36,15 +38,16 @@
 from mla import sources
 from mla import minimizers
 from mla import trial_generators
-from mla.utility_functions import newton_method, newton_method_multidataset
+from mla.utility_functions import newton_method_multidataset
 
 __all__ = ["NeutrinoPointSource"]
 r"""This IceCube plugin is currently under develop by Kwok Lung Fan"""
 
 
-class NeutrinoPointSource(PointSource,Node):
+class NeutrinoPointSource(PointSource, Node):
     """
-    Class for NeutrinoPointSource. It is inherited from astromodels PointSource class.
+    Class for NeutrinoPointSource.
+    It is inherited from astromodels PointSource class.
     """
 
     def __init__(
@@ -67,10 +70,12 @@ def __init__(
             source_name:Name of the source
             ra: right ascension in degree
             dec: declination in degree
-            spectral_shape: Shape of the spectrum.Check 3ML example for more detail.
+            spectral_shape:
+              Shape of the spectrum.Check 3ML example for more detail.
             l: galactic longitude in degree
             b: galactic   in degree
-            components: Spectral Component.Check 3ML example for more detail.
+            components:
+               Spectral Component.Check 3ML example for more detail.
             sky_position: sky position
             energy_unit: Unit of the energy
         """
@@ -84,7 +89,8 @@ def __init__(
             (ra is not None and dec is not None)
             ^ (l is not None and b is not None)
             ^ (sky_position is not None)
-        ), "You have to provide one and only one specification for the position"
+        ), "You have to provide one and only one \
+           specification for the position"
 
         # Gather the position
 
@@ -102,10 +108,10 @@ def __init__(
                 except (TypeError, ValueError):
 
                     raise AssertionError(
-                        "RA and Dec must be numbers. If you are confused by this message,"
-                        " you are likely using the constructor in the wrong way. Check"
-                        " the documentation."
-                    )
+                        "RA and Dec must be numbers. "
+                        "If you are confused by this message, "
+                        "you are likely using the constructor in the wrong way"
+                        ". Check the documentation.")
 
                 sky_position = SkyDirection(ra=ra, dec=dec)
 
@@ -117,7 +123,8 @@ def __init__(
 
         # Now gather the component(s)
 
-        # We need either a single component, or a list of components, but not both
+        # We need either a single component,
+        # or a list of components, but not both
         # (that's the ^ symbol)
 
         assert (spectral_shape is not None) ^ (components is not None), (
@@ -155,7 +162,8 @@ def __init__(
         # Components in this case have energy as x and differential flux as y
 
         x_unit = energy_unit
-        y_unit = (energy_unit * current_units.area * current_units.time) ** (-1)
+        y_unit = (energy_unit * current_units.area *
+                  current_units.time) ** (-1)
 
         # Now set the units of the components
         for component in list(self._components.values()):
@@ -197,8 +205,8 @@ def call(self, x, tag=None):
                 # Slow version with units
 
                 results = [
-                    component.shape(x) for component in list(self.components.values())
-                ]
+                    component.shape(x) for component in list(
+                        self.components.values())]
 
                 # We need to sum like this (slower) because using
                 # np.sum will not preserve the units (thanks astropy.units)
@@ -207,12 +215,13 @@ def call(self, x, tag=None):
 
             else:
 
-                # Fast version without units, where x is supposed to be in the same
+                # Fast version without units,
+                # where x is supposed to be in the same
                 # units as currently defined in units.get_units()
 
                 results = [
-                    component.shape(x) for component in list(self.components.values())
-                ]
+                    component.shape(x) for component in list(
+                        self.components.values())]
 
                 return np.sum(results, 0)
 
@@ -276,7 +285,8 @@ def __init__(
 
             # Now gather the component(s)
 
-            # We need either a single component, or a list of components, but not both
+            # We need either a single component,
+            # or a list of components, but not both
             # (that's the ^ symbol)
 
             assert (spectral_shape is not None) ^ (components is not None), (
@@ -289,11 +299,11 @@ def __init__(
 
                 components = [SpectralComponent("main", spectral_shape)]
 
-            # Components in this case have energy as x and differential flux as y
+            # Components in this case have energy as x and differential flux as
+            # y
 
-            diff_flux_units = (current_u.energy * current_u.area * current_u.time) ** (
-                -1
-            )
+            diff_flux_units = (current_u.energy *
+                               current_u.area * current_u.time) ** (-1)
 
             # Now set the units of the components
             for component in components:
@@ -343,7 +353,6 @@ def spatial_shape(self):
         return self._spatial_shape
 
     def get_spatially_integrated_flux(self, energies):
-
         """
         Returns total flux of source at the given energy
         :param energies: energies (array or float)
@@ -397,7 +406,8 @@ def __call__(self, lon, lat, energies):
                 # Slow version with units
 
                 # We need to sum like this (slower) because
-                # using np.sum will not preserve the units (thanks astropy.units)
+                # using np.sum will not preserve the units (thanks
+                # astropy.units)
 
                 result = np.zeros((lat.shape[0], energies.shape[0])) * (
                     u.keV ** -1 * u.cm ** -2 * u.second ** -1 * u.degree ** -2
@@ -505,7 +515,8 @@ def _repr__base(self, rich_output=False):
         repr_dict[key]["spectrum"] = collections.OrderedDict()
 
         for component_name, component in list(self.components.items()):
-            repr_dict[key]["spectrum"][component_name] = component.to_dict(minimal=True)
+            repr_dict[key]["spectrum"][component_name] = component.to_dict(
+                minimal=True)
 
         return dict_to_list(repr_dict, rich_output)
 
@@ -527,11 +538,15 @@ def __init__(self, likelihood_model_instance, A=1):
         r"""Constructor of the class"""
         self.model = likelihood_model_instance
         self.norm = A
-        for source_name, source in likelihood_model_instance.point_sources.items():
+        for source_name, source in (
+            likelihood_model_instance.point_sources.items()
+        ):
             if isinstance(source, NeutrinoPointSource):
                 self.neutrinosource = source_name
                 self.point = True
-        for source_name, source in likelihood_model_instance.extended_sources.items():
+        for source_name, source in (
+                likelihood_model_instance.extended_sources.items()
+        ):
             if isinstance(source, NeutrinoExtendedSource):
                 self.neutrinosource = source_name
                 self.point = False
@@ -539,9 +554,8 @@ def __init__(self, likelihood_model_instance, A=1):
     def __call__(self, energy, **kwargs):
         r"""Evaluate spectrum at E"""
         if self.point:
-            return (
-                self.model.point_sources[self.neutrinosource].call(energy) * self.norm
-            )
+            return (self.model.point_sources[self.neutrinosource].call(
+                energy) * self.norm)
         else:
             return (
                 self.model.extended_sources[self.neutrinosource].call(energy)
@@ -553,7 +567,8 @@ def validate(self):
 
     def __str__(self):
         r"""String representation of class"""
-        return "SpectrumConverter class doesn't support string representation now"
+        return "SpectrumConverter class doesn't \
+                support string representation now"
 
     def copy(self):
         r"""Return copy of this class"""
@@ -580,9 +595,11 @@ def __init__(
             name: name for the plugin
             data: data of experiment
             data_handlers: mla.threeml.data_handlers ThreeMLDataHandler object
-            llh: test_statistics.LLHTestStatistic object. Used to evaluate the ts
+            llh: test_statistics.LLHTestStatistic object.
+                 Used to evaluate the ts
             source: injection location(only when need injection)
-            livetime: livetime in days(calculated using livetime within time profile if None)
+            livetime: livetime in days(calculated using
+                      livetime within time profile if None)
             fix_flux_norm: only fit the spectrum shape
             verbose: print the output or not
 
@@ -612,13 +629,14 @@ def __init__(
             config["dec"] = 0
             source = sources.PointSource(config=config)
         self.injected_source = source
-        trial_config = trial_generators.SingleSourceTrialGenerator.generate_config()
+        trial_config = \
+            trial_generators.SingleSourceTrialGenerator.generate_config()
         self.trial_generator = trial_generators.SingleSourceTrialGenerator(
             trial_config, data_handlers, source
         )
-        analysis_config = analysis.SingleSourceLLHAnalysis.generate_default_config(
-            minimizer_class=minimizers.GridSearchMinimizer
-        )
+        analysis_config = \
+            analysis.SingleSourceLLHAnalysis.generate_default_config(
+                minimizer_class=minimizers.GridSearchMinimizer)
         self.analysis = analysis.SingleSourceLLHAnalysis(
             config=analysis_config,
             minimizer_class=minimizers.GridSearchMinimizer,
@@ -642,8 +660,8 @@ def __init__(
         )
         for key in self.test_statistic.sob_terms.keys():
             if isinstance(
-                self.test_statistic.sob_terms[key], sob_terms.ThreeMLPSEnergyTerm
-            ):
+                    self.test_statistic.sob_terms[key],
+                    sob_terms.ThreeMLPSEnergyTerm):
                 self.energyname = key
         return
 
@@ -653,7 +671,9 @@ def set_model(self, likelihood_model_instance):
 
             return
 
-        for source_name, source in likelihood_model_instance.point_sources.items():
+        for source_name, source in (
+            likelihood_model_instance.point_sources.items()
+        ):
             if isinstance(source, NeutrinoPointSource):
                 self.source_name = source_name
                 ra = source.position.get_ra()
@@ -686,13 +706,16 @@ def set_model(self, likelihood_model_instance):
                     self.test_statistic = self.analysis.test_statistic_factory(
                         Params.from_dict({"ns": 0}), self._data
                     )
-        for source_name, source in likelihood_model_instance.extended_sources.items():
+        for source_name, source in (
+            likelihood_model_instance.extended_sources.items()
+        ):
             if isinstance(source, NeutrinoExtendedSource):
                 self.source_name = source_name
                 ra = source.spatial_shape.lon0.value
                 dec = source.spatial_shape.lat0.value
                 sigma = source.spatial_shape.sigma.value
-                if self._ra == ra and self._dec == dec and self._sigma == sigma:
+                if (self._ra == ra and self._dec == dec
+                        and self._sigma == sigma):
                     self.llh_model = likelihood_model_instance
                     self.energy_sob_factory.spectrum = Spectrum(
                         likelihood_model_instance
@@ -811,8 +834,11 @@ class icecube_analysis(PluginPrototype):
     """Docstring"""
 
     def __init__(
-        self, listoficecubelike, newton_flux_norm=False, name="combine", verbose=False
-    ):
+            self,
+            listoficecubelike,
+            newton_flux_norm=False,
+            name="combine",
+            verbose=False):
         """Docstring"""
         nuisance_parameters = {}
         super(icecube_analysis, self).__init__(name, nuisance_parameters)
@@ -849,9 +875,9 @@ def get_log_like(self, verbose=None):
             for i, icecubeobject in enumerate(self.listoficecubelike):
                 sob.append(icecubeobject.test_statistic._calculate_sob())
                 n_drop.append(icecubeobject.test_statistic.n_dropped)
-                fraction.append(
-                    self.totaln * self.dataset_weight[i] / len(icecubeobject.data)
-                )
+                fraction.append(self.totaln *
+                                self.dataset_weight[i] /
+                                len(icecubeobject.data))
             fraction = np.array(fraction)
             # fraction = fraction/fraction.sum()
             ns_ratio = newton_method_multidataset(sob, n_drop, fraction)
@@ -943,7 +969,9 @@ def injection(self, n_signal=0, flux_norm=None, poisson=False):
                 ratio_injection = self.dataset_ratio * n_signal
                 for i, icecubeobject in enumerate(self.listoficecubelike):
                     icecubeobject.trial_generator.config["fixed_ns"] = True
+                    # print('in IceCubelike test injection', n_signal)
                     injection_signal = np.random.poisson(ratio_injection[i])
+                    # print(injection_signal)
                     tempdata = icecubeobject.trial_generator(injection_signal)
                     self.listoficecubelike[i].update_data(tempdata)
             else:
@@ -958,9 +986,10 @@ def cal_injection_ns(self, flux_norm):
         for icecubeobject in self.listoficecubelike:
             time_intergrated = flux_norm * icecubeobject.livetime * 3600 * 24
             tempns = (
-                time_intergrated
-                * icecubeobject.analysis.data_handler_source[0].sim["weight"].sum()
-            )
+                time_intergrated *
+                icecubeobject.analysis
+                .data_handler_source[0]
+                .sim["weight"].sum())
             ns = ns + tempns
         return ns
 
@@ -969,10 +998,11 @@ def cal_injection_fluxnorm(self, ns):
         totalweight = 0
         for icecubeobject in self.listoficecubelike:
             tempweight = (
-                icecubeobject.analysis.data_handler_source[0].sim["weight"].sum()
-                * icecubeobject.livetime
-                * 3600
-                * 24
-            )
+                icecubeobject.analysis
+                .data_handler_source[0]
+                .sim["weight"].sum() *
+                icecubeobject.livetime *
+                3600 *
+                24)
             totalweight = totalweight + tempweight
         return ns / totalweight
diff --git a/mla/threeml/data_handlers.py b/mla/threeml/data_handlers.py
index 8d65bd1..dcb2021 100644
--- a/mla/threeml/data_handlers.py
+++ b/mla/threeml/data_handlers.py
@@ -1,13 +1,13 @@
 """Docstring"""
 
-__author__ = 'John Evans and Jason Fan'
-__copyright__ = 'Copyright 2024'
-__credits__ = ['John Evans', 'Jason Fan', 'Michael Larson']
-__license__ = 'Apache License 2.0'
-__version__ = '1.4.1'
-__maintainer__ = 'Jason Fan'
-__email__ = 'klfan@terpmail.umd.edu'
-__status__ = 'Development'
+__author__ = "John Evans and Jason Fan"
+__copyright__ = "Copyright 2024"
+__credits__ = ["John Evans", "Jason Fan", "Michael Larson"]
+__license__ = "Apache License 2.0"
+__version__ = "1.4.1"
+__maintainer__ = "Jason Fan"
+__email__ = "klfan@terpmail.umd.edu"
+__status__ = "Development"
 
 import dataclasses
 
@@ -49,7 +49,8 @@ def build_signal_energy_histogram(
     ) -> np.ndarray:
         """
         Building the signal energy histogram.
-        Only used when using MC instead of IRF to build signal energy histogram.
+        Only used when using MC instead of IRF
+         to build signal energy histogram.
 
         Args:
             spectrum: signal spectrum
@@ -64,7 +65,8 @@ def build_signal_energy_histogram(
             density=True,
         )[0]
 
-    def cut_reconstructed_sim(self, dec: float, sampling_width: float) -> np.ndarray:
+    def cut_reconstructed_sim(
+            self, dec: float, sampling_width: float) -> np.ndarray:
         """
         Cutting the MC based on reconstructed dec.
         Only use when using MC instead of IRF to build signal energy histogram.
@@ -103,7 +105,8 @@ def injection_spectrum(self) -> spectral.BaseSpectrum:
         return self._injection_spectrum
 
     @injection_spectrum.setter
-    def injection_spectrum(self, inject_spectrum: spectral.BaseSpectrum) -> None:
+    def injection_spectrum(
+            self, inject_spectrum: spectral.BaseSpectrum) -> None:
         """
         Setting the injection spectrum
 
@@ -124,7 +127,8 @@ def injection_spectrum(self, inject_spectrum: spectral.BaseSpectrum) -> None:
 
         self._full_sim["weight"] = (
             self._full_sim["ow"]
-            * (inject_spectrum(self._full_sim["trueE"] * self._flux_unit_conversion))
+            * (inject_spectrum(self._full_sim["trueE"]
+                               * self._flux_unit_conversion))
             * self._flux_unit_conversion
         )
 
diff --git a/mla/threeml/profilellh.py b/mla/threeml/profilellh.py
index 3bfe386..5af5e10 100644
--- a/mla/threeml/profilellh.py
+++ b/mla/threeml/profilellh.py
@@ -1,13 +1,13 @@
 """Docstring"""
 
-__author__ = 'John Evans and Jason Fan'
-__copyright__ = 'Copyright 2024'
-__credits__ = ['John Evans', 'Jason Fan', 'Michael Larson']
-__license__ = 'Apache License 2.0'
-__version__ = '1.4.1'
-__maintainer__ = 'Jason Fan'
-__email__ = 'klfan@terpmail.umd.edu'
-__status__ = 'Development'
+__author__ = "John Evans and Jason Fan"
+__copyright__ = "Copyright 2024"
+__credits__ = ["John Evans", "Jason Fan", "Michael Larson"]
+__license__ = "Apache License 2.0"
+__version__ = "1.4.1"
+__maintainer__ = "Jason Fan"
+__email__ = "klfan@terpmail.umd.edu"
+__status__ = "Development"
 
 from threeML.plugin_prototype import PluginPrototype
 from scipy.interpolate import RegularGridInterpolator
@@ -45,18 +45,16 @@ def __init__(
         super().__init__(name, nuisance_parameters)
         if spline is not None:
             self.spline = spline
-            self.df = None
+            self.df = df  # None
         else:
             self.df = df
             self.par_name = list(df.columns)
             self.par_name.pop()
-            listofpoint = []
-            shape = []
-            for n in self.par_name:
-                points = np.unique(df[n])
-                listofpoint.append(points)
-                shape.append(points.shape[0])
-            llh = np.reshape(df["llh"].values, shape)
+            listofpoint = [np.unique(df[n]) for n in self.par_name]
+            shape = [len(points) for points in listofpoint]
+            sort_idx = np.lexsort(
+                         [df[p].values for p in reversed(self.par_name)])
+            llh = np.reshape(df["llh"].values[sort_idx], shape)
             self.spline = RegularGridInterpolator(
                 listofpoint, llh, bounds_error=False, fill_value=fill_value
             )
@@ -90,10 +88,14 @@ def get_log_like(self) -> float:
 
     def inner_fit(self) -> float:
         """
-        This is used for the profile likelihood. Keeping fixed all parameters in the
-        LikelihoodModel, this method minimize the logLike over the remaining nuisance
-        parameters, i.e., the parameters belonging only to the model for this
-        particular detector. If there are no nuisance parameters, simply return the
+        This is used for the profile likelihood.
+        Keeping fixed all parameters in the
+        LikelihoodModel, this method minimize the logLike
+        over the remaining nuisance
+        parameters, i.e., the parameters belonging
+        only to the model for this
+        particular detector. If there are no nuisance parameters,
+        simply return the
         logLike value.
         """
 
diff --git a/mla/threeml/sob_terms.py b/mla/threeml/sob_terms.py
index ff7ce5a..581c733 100644
--- a/mla/threeml/sob_terms.py
+++ b/mla/threeml/sob_terms.py
@@ -1,13 +1,13 @@
 """Docstring"""
 
-__author__ = 'John Evans and Jason Fan'
-__copyright__ = 'Copyright 2024'
-__credits__ = ['John Evans', 'Jason Fan', 'Michael Larson']
-__license__ = 'Apache License 2.0'
-__version__ = '1.4.1'
-__maintainer__ = 'Jason Fan'
-__email__ = 'klfan@terpmail.umd.edu'
-__status__ = 'Development'
+__author__ = "John Evans and Jason Fan"
+__copyright__ = "Copyright 2024"
+__credits__ = ["John Evans", "Jason Fan", "Michael Larson"]
+__license__ = "Apache License 2.0"
+__version__ = "1.4.1"
+__maintainer__ = "Jason Fan"
+__email__ = "klfan@terpmail.umd.edu"
+__status__ = "Development"
 
 import dataclasses
 
@@ -105,30 +105,38 @@ class ThreeMLPSEnergyTermFactory(ThreeMLBaseEnergyTermFactory):
     def __post_init__(self) -> None:
         """Docstring"""
         if self.config["list_sin_dec_bins"] is None:
-            self._sin_dec_bins = np.linspace(-1, 1, 1 + self.config["sin_dec_bins"])
+            self._sin_dec_bins = np.linspace(
+                -1, 1, 1 + self.config["sin_dec_bins"])
         else:
             self._sin_dec_bins = self.config["list_sin_dec_bins"]
         if self.config["list_log_energy_bins"] is None:
             self._log_energy_bins = np.linspace(
-                *self.config["log_energy_bounds"], 1 + self.config["log_energy_bins"]
-            )
+                *self.config["log_energy_bounds"],
+                1 + self.config["log_energy_bins"])
         else:
             self._log_energy_bins = self.config["list_log_energy_bins"]
-        self.data_handler.reduced_reco_sim = self.data_handler.cut_reconstructed_sim(
-            self.source.location[1],
-            self.data_handler.config["reco_sampling_width"],
-        )
+        self.data_handler.reduced_reco_sim = \
+            self.data_handler.cut_reconstructed_sim(
+                self.source.location[1],
+                self.data_handler.config["reco_sampling_width"],
+            )
         self._unit_scale = self.config["Energy_convesion(ToGeV)"]
-        self._bins = np.array([self._sin_dec_bins, self._log_energy_bins], dtype=object)
+        self._bins = np.array(
+            [self._sin_dec_bins, self._log_energy_bins], dtype=object)
         self._init_bg_sob_map()
         self._build_ow_hist()
 
-    def __call__(self, params: par.Params, events: np.ndarray) -> sob_terms.SoBTerm:
+    def __call__(
+            self,
+            params: par.Params,
+            events: np.ndarray) -> sob_terms.SoBTerm:
         """Docstring"""
         # Get the bin that each event belongs to
-        sin_dec_idx = np.searchsorted(self._sin_dec_bins[:-1], events["sindec"]) - 1
+        sin_dec_idx = np.searchsorted(
+            self._sin_dec_bins[:-1], events["sindec"]) - 1
 
-        log_energy_idx = np.searchsorted(self._log_energy_bins[:-1], events["logE"]) - 1
+        log_energy_idx = np.searchsorted(
+            self._log_energy_bins[:-1], events["logE"]) - 1
 
         return ThreeMLPSEnergyTerm(
             name=self.config["name"],
@@ -150,18 +158,19 @@ def _build_ow_hist(self) -> np.ndarray:
 
     def get_ns(self) -> float:
         """Docstring"""
-        return (self.spectrum(self._ow_ebin) * self._ow_hist).sum() * self._unit_scale
+        return (self.spectrum(self._ow_ebin) *
+                self._ow_hist).sum() * self._unit_scale
 
     def _init_bg_sob_map(self) -> np.ndarray:
         """Docstring"""
         if self.config["mc_bkgweight"] is None:
-            bg_h = self.data_handler.build_background_sindec_logenergy_histogram(
-                self._bins
-            )
+            bg_h = \
+                self.data_handler.build_background_sindec_logenergy_histogram(
+                    self._bins)
         else:
-            bg_h = self.data_handler.build_mcbackground_sindec_logenergy_histogram(
-                self._bins, self.config["mc_bkgweight"]
-            )
+            bg_h = \
+               self.data_handler.build_mcbackground_sindec_logenergy_histogram(
+                    self._bins, self.config["mc_bkgweight"])
             print("using mc background")
         # Normalize histogram by dec band
         bg_h /= np.sum(bg_h, axis=1)[:, None]
@@ -180,10 +189,11 @@ def source(self) -> sources.PointSource:
     def source(self, source: sources.PointSource) -> None:
         """Docstring"""
         self._source = source
-        self.data_handler.reduced_reco_sim = self.data_handler.cut_reconstructed_sim(
-            self.source.location[1],
-            self.data_handler.config["reco_sampling_width"],
-        )
+        self.data_handler.reduced_reco_sim = \
+            self.data_handler.cut_reconstructed_sim(
+                self.source.location[1],
+                self.data_handler.config["reco_sampling_width"],
+            )
 
     def cal_sob_map(self) -> np.ndarray:
         """Creates sob histogram for a given spectrum.
@@ -194,7 +204,8 @@ def cal_sob_map(self) -> np.ndarray:
         sig_h = self.data_handler.build_signal_energy_histogram(
             self.spectrum, self._bins, self._unit_scale
         )
-        bin_centers = self._log_energy_bins[:-1] + np.diff(self._log_energy_bins) / 2
+        bin_centers = self._log_energy_bins[:-1] + \
+            np.diff(self._log_energy_bins) / 2
         # Normalize histogram by dec band
         sig_h /= np.sum(sig_h, axis=1)[:, None]
 
@@ -268,10 +279,11 @@ class ThreeMLPSIRFEnergyTermFactory(ThreeMLPSEnergyTermFactory):
     _spectrum: spectral.BaseSpectrum = dataclasses.field(
         init=False, repr=False, default=spectral.PowerLaw(1e3, 1e-14, -2)
     )
+
     _bg_sob: np.ndarray = dataclasses.field(init=False, repr=False)
     _sin_dec_bins: np.ndarray = dataclasses.field(
-        init=False, repr=False, default=PSTrackv4_sin_dec_bin
-    )
+        init=False, repr=False,
+        default_factory=lambda: PSTrackv4_sin_dec_bin.copy())
     _log_energy_bins: np.ndarray = dataclasses.field(init=False, repr=False)
     _bins: np.ndarray = dataclasses.field(init=False, repr=False)
     _trueebin: np.ndarray = dataclasses.field(init=False, repr=False)
@@ -282,14 +294,18 @@ class ThreeMLPSIRFEnergyTermFactory(ThreeMLPSEnergyTermFactory):
 
     def __post_init__(self) -> None:
         """Docstring"""
+        print("Calling __post_init__")  # or use logging
+        # self._source = self.config.get("source", None)
         if self.config["list_sin_dec_bins"] is None:
-            self._sin_dec_bins = np.linspace(-1, 1, 1 + self.config["sin_dec_bins"])
+            self._sin_dec_bins = np.linspace(
+                                     -1, 1, 1 + self.config["sin_dec_bins"])
         else:
             self._sin_dec_bins = self.config["list_sin_dec_bins"]
+
         if self.config["list_log_energy_bins"] is None:
             self._log_energy_bins = np.linspace(
-                *self.config["log_energy_bounds"], 1 + self.config["log_energy_bins"]
-            )
+                *self.config["log_energy_bounds"],
+                1 + self.config["log_energy_bins"])
         else:
             self._log_energy_bins = self.config["list_log_energy_bins"]
         lower_sindec = np.maximum(
@@ -306,22 +322,31 @@ def __post_init__(self) -> None:
             ),
             1,
         )
-        lower_sindec_index = np.searchsorted(self._sin_dec_bins, lower_sindec) - 1
+        lower_sindec_index = np.searchsorted(
+            self._sin_dec_bins, lower_sindec) - 1
         uppper_sindec_index = np.searchsorted(self._sin_dec_bins, upper_sindec)
-        self._sindec_bounds = np.array([lower_sindec_index, uppper_sindec_index])
-        self._bins = np.array([self._sin_dec_bins, self._log_energy_bins], dtype=object)
+        # print(lower_sindec_index,uppper_sindec_index)
+        self._sindec_bounds = np.array(
+            [lower_sindec_index, uppper_sindec_index])
+        self._bins = np.array(
+            [self._sin_dec_bins, self._log_energy_bins], dtype=object)
         self._truelogebin = self.config["list_truelogebin"]
         self._unit_scale = self.config["Energy_convesion(ToGeV)"]
         self._init_bg_sob_map()
         self._build_ow_hist()
         self._init_irf()
 
-    def __call__(self, params: par.Params, events: np.ndarray) -> sob_terms.SoBTerm:
+    def __call__(
+            self,
+            params: par.Params,
+            events: np.ndarray) -> sob_terms.SoBTerm:
         """Docstring"""
         # Get the bin that each event belongs to
-        sin_dec_idx = np.searchsorted(self._sin_dec_bins[:-1], events["sindec"]) - 1
+        sin_dec_idx = np.searchsorted(
+            self._sin_dec_bins[:-1], events["sindec"]) - 1
 
-        log_energy_idx = np.searchsorted(self._log_energy_bins[:-1], events["logE"]) - 1
+        log_energy_idx = np.searchsorted(
+            self._log_energy_bins[:-1], events["logE"]) - 1
 
         return ThreeMLPSEnergyTerm(
             name=self.config["name"],
@@ -335,13 +360,13 @@ def __call__(self, params: par.Params, events: np.ndarray) -> sob_terms.SoBTerm:
     def _init_bg_sob_map(self) -> None:
         """Docstring"""
         if self.config["mc_bkgweight"] is None:
-            bg_h = self.data_handler.build_background_sindec_logenergy_histogram(
-                self._bins
-            )
+            bg_h = \
+                self.data_handler.build_background_sindec_logenergy_histogram(
+                 self._bins)
         else:
-            bg_h = self.data_handler.build_mcbackground_sindec_logenergy_histogram(
-                self._bins, self.config["mc_bkgweight"]
-            )
+            bg_h = \
+               self.data_handler.build_mcbackground_sindec_logenergy_histogram(
+                self._bins, self.config["mc_bkgweight"])
             print("using mc background")
         # Normalize histogram by dec band
         bg_h /= np.sum(bg_h, axis=1)[:, None]
@@ -362,13 +387,17 @@ def _init_irf(self) -> None:
         )
         self._trueebin = 10 ** (self._truelogebin[:-1])
         sindec_idx = (
-            np.digitize(np.sin(self.data_handler.full_sim["dec"]), self._sin_dec_bins)
-            - 1
-        )
+            np.digitize(
+                np.sin(
+                    self.data_handler.full_sim["dec"]),
+                self._sin_dec_bins) -
+            1)
 
         for i in range(len(self._sin_dec_bins) - 1):
             events_dec = self.data_handler.full_sim[(sindec_idx == i)]
-            loge_idx = np.digitize(events_dec["logE"], self._log_energy_bins) - 1
+            loge_idx = np.digitize(
+                events_dec["logE"],
+                self._log_energy_bins) - 1
 
             for j in range(len(self._log_energy_bins) - 1):
                 events = events_dec[(loge_idx == j)]
@@ -384,19 +413,23 @@ def _init_irf(self) -> None:
                     weights=events["ow"],
                 )
 
-                # Have to pick an "energy" to assign to the bin. That's complicated, since
-                # you'd (in principle) want the flux-weighted average energy, but we don't
-                # have a flux function here. Instead, try just using the minimum energy of
+                # Have to pick an "energy" to assign to the bin.
+                # That's complicated, since
+                # you'd (in principle) want the flux-weighted average energy,
+                # but we don't have a flux function here. Instead, try just
+                # using the minimum energy of
                 # the bin? Should be fine for small enough bins.
                 # self._trueebin[i,j] = np.exp(bins[:-1] + (bins[1] - bins[0]))
-                # emean[i,j] = np.average(events['trueE'], weights=events['ow'])
+                # emean[i,j] = np.average(events['trueE'],
+                # weights=events['ow'])
 
     def build_sig_h(self, spectrum: spectral.BaseSpectrum) -> np.ndarray:
         """Docstring"""
         sig = np.zeros(self._bg_sob.shape)
         flux = spectrum(self._trueebin * self._unit_scale)  # converting unit
         sig[self._sindec_bounds[0]:self._sindec_bounds[1], :] = np.dot(
-            self._irf[self._sindec_bounds[0]:self._sindec_bounds[1], :, :], flux
+            self._irf[self._sindec_bounds[0]:self._sindec_bounds[1], :, :],
+            flux
         )
         sig /= np.sum(sig, axis=1)[:, None]
         return sig
@@ -410,23 +443,22 @@ def source(self) -> sources.PointSource:
     def source(self, source: sources.PointSource) -> None:
         """Docstring"""
         self._source = source
-        lower_sindec = np.maximum(
-            np.sin(
-                self.source.location[1]
-                - self.data_handler.config["reco_sampling_width"]
-            ),
-            -0.99,
-        )
-        upper_sindec = np.minimum(
-            np.sin(
-                self.source.location[1]
-                + self.data_handler.config["reco_sampling_width"]
-            ),
-            1,
-        )
-        lower_sindec_index = np.searchsorted(self._sin_dec_bins, lower_sindec) - 1
-        uppper_sindec_index = np.searchsorted(self._sin_dec_bins, upper_sindec)
-        self._sindec_bounds = np.array([lower_sindec_index, uppper_sindec_index])
+        # flake8 says that this variable is not used. Kept it under the comment
+        # section if we end up having to use it
+        # lower_sindec = np.maximum(
+        #    np.sin(
+        #        self._source.location[1]
+        #        - self.data_handler.config["reco_sampling_width"]
+        #    ),
+        #    -0.99,
+        # )
+        # upper_sindec = np.minimum(
+        #    np.sin(
+        #        self._source.location[1]
+        #        + self.data_handler.config["reco_sampling_width"]
+        #    ),
+        #    1,
+        # )
 
     def cal_sob_map(self) -> np.ndarray:
         """Creates sob histogram for a given spectrum.
diff --git a/mla/time_profiles.py b/mla/time_profiles.py
index 06fb393..f641ada 100644
--- a/mla/time_profiles.py
+++ b/mla/time_profiles.py
@@ -51,7 +51,7 @@ def pdf(self, times: np.ndarray) -> np.ndarray:
         """Get the probability amplitude given a time for this time profile.
 
         Args:
-            times: An array of event times to get the probability amplitude for.
+            times: An array of event times to get the probability amplitude.
 
         Returns:
             A numpy array of probability amplitudes at the given times.
@@ -119,7 +119,10 @@ def cdf(self, times: np.ndarray) -> np.ndarray:
 
     @abc.abstractmethod
     def inverse_transform_sample(
-            self, start_times: np.ndarray, stop_times: np.ndarray) -> np.ndarray:
+            self,
+            start_times: np.ndarray,
+            stop_times: np.ndarray
+    ) -> np.ndarray:
         """Docstring"""
 
     @property
@@ -157,7 +160,8 @@ def default_params(self) -> Dict[str, float]:
     @property
     @abc.abstractmethod
     def param_dtype(self) -> np.dtype:
-        """Returns the parameter names and datatypes formatted for numpy dtypes.
+        """Returns the parameter names and datatypes
+           formatted for numpy dtypes.
         """
 
 
@@ -179,7 +183,8 @@ class GaussProfile(GenericProfile):
         param_dtype (List[Tuple[str, str]]): The numpy dytpe for the fitting
             parameters.
     """
-    scipy_dist: scipy.stats.distributions.rv_frozen = dataclasses.field(init=False)
+    scipy_dist: scipy.stats.distributions.rv_frozen = dataclasses.field(
+        init=False)
     _mean: float = dataclasses.field(init=False, repr=False)
     _sigma: float = dataclasses.field(init=False, repr=False)
     _param_dtype: ClassVar[np.dtype] = np.dtype(
@@ -238,7 +243,8 @@ def x0(self, times: np.ndarray) -> tuple:
         return x0_mean, x0_sigma
 
     def bounds(self, time_profile: GenericProfile) -> List[tuple]:
-        """Returns good bounds for this time profile given another time profile.
+        """Returns good bounds for this time profile
+           given another time profile.
 
         Limits the mean to be within the range of the other profile and limits
         the sigma to be >= 0 and <= the width of the other profile.
@@ -262,7 +268,10 @@ def cdf(self, times: np.ndarray) -> np.ndarray:
         return self.scipy_dist.cdf(times)
 
     def inverse_transform_sample(
-            self, start_times: np.ndarray, stop_times: np.ndarray) -> np.ndarray:
+            self,
+            start_times: np.ndarray,
+            stop_times: np.ndarray
+    ) -> np.ndarray:
         """Docstring"""
         start_cdfs = self.cdf(start_times)
         stop_cdfs = self.cdf(stop_times)
@@ -336,7 +345,9 @@ class UniformProfile(GenericProfile):
 
     def __post_init__(self) -> None:
         """Constructs the time profile."""
-        self._range = (self.config['start'], self.config['start'] + self.config['length'])
+        self._range = (
+            self.config['start'],
+            self.config['start'] + self.config['length'])
 
     def pdf(self, times: np.ndarray) -> np.ndarray:
         """Calculates the probability for each time.
@@ -404,10 +415,13 @@ def bounds(self, time_profile: GenericProfile
 
     def cdf(self, times: np.ndarray) -> np.ndarray:
         """Docstring"""
-        return np.clip((times - self.range[0]) / (self.range[1] - self.range[0]), 0, 1)
+        return np.clip(
+            (times - self.range[0]) / (self.range[1] - self.range[0]), 0, 1)
 
     def inverse_transform_sample(
-            self, start_times: np.ndarray, stop_times: np.ndarray) -> np.ndarray:
+            self,
+            start_times: np.ndarray,
+            stop_times: np.ndarray) -> np.ndarray:
         """Docstring"""
         return np.random.uniform(
             np.maximum(start_times, self.range[0]),
@@ -417,14 +431,17 @@ def inverse_transform_sample(
     @property
     def params(self) -> dict:
         """Docstring"""
-        return {'start': self._range[0], 'length': self._range[1] - self._range[0]}
+        return {
+            'start': self._range[0],
+            'length': self._range[1] - self._range[0]}
 
     @params.setter
     def params(self, params: Params) -> None:
         """Docstring"""
         if 'start' in params:
             self._range = (
-                params['start'], params['start'] + self._range[1] - self._range[0])
+                params['start'],
+                params['start'] + self._range[1] - self._range[0])
         if 'length' in params:
             self._range = (self._range[0], self._range[0] + params['length'])
 
@@ -462,7 +479,8 @@ class CustomProfile(GenericProfile):
 
     Attributes:
         pdf (Callable[[np.array, Tuple[float, float]], np.array]): The
-            distribution function. This function needs to accept an array of bin
+            distribution function.
+            This function needs to accept an array of bin
             centers and a time window as a tuple, and it needs to return an
             array of probability densities at the given bin centers.
         dist (scipy.stats.rv_histogram): The histogrammed version of the
@@ -499,20 +517,19 @@ def dist(
         self._offset = self.config['offset']
 
         if isinstance(self.config['bins'], int):
-            bin_edges = np.linspace(*self.config['range'], self.config['bins'])
+            bin_edges = np.linspace(*self.config['range'],
+                                    self.config['bins'] + 1)
         else:
-            span = self.config['range'][1] - self.config['range'][0]
-            bin_edges = span * np.array(self.config['bins'])
-
+            bin_edges = np.array(self.config['bins'])
         bin_widths = np.diff(bin_edges)
-        bin_centers = bin_edges[:-1] + bin_widths
-        hist = dist(bin_centers, tuple(self.config['range']))
+        bin_centers = bin_edges[:-1] + bin_widths / 2
 
+        hist, _ = np.histogram(bin_centers, bins=bin_edges)
+        hist = hist.astype(float)
         area_under_hist = np.sum(hist * bin_widths)
-        hist *= 1 / area_under_hist
+        hist *= 1. / area_under_hist
         self._exposure = 1 / np.max(hist)
         hist *= bin_widths
-
         self._dist = scipy.stats.rv_histogram((hist, bin_edges))
 
     def pdf(self, times: np.ndarray) -> np.ndarray:
@@ -579,7 +596,9 @@ def cdf(self, times: np.ndarray) -> np.ndarray:
         return self.dist.cdf(times)
 
     def inverse_transform_sample(
-            self, start_times: np.ndarray, stop_times: np.ndarray) -> np.ndarray:
+            self,
+            start_times: np.ndarray,
+            stop_times: np.ndarray) -> np.ndarray:
         """Docstring"""
         start_cdfs = self.cdf(start_times)
         stop_cdfs = self.cdf(stop_times)
@@ -616,7 +635,8 @@ def offset(self, offset: float) -> None:
     @property
     def range(self) -> Tuple[Optional[float], Optional[float]]:
         return (
-            self.config['range'][0] + self.offset, self.config['range'][1] + self.offset)
+            self.config['range'][0] + self.offset,
+            self.config['range'][1] + self.offset)
 
     @property
     def param_dtype(self) -> np.dtype:
diff --git a/mla/trial_generators.py b/mla/trial_generators.py
index 10939a0..629a118 100644
--- a/mla/trial_generators.py
+++ b/mla/trial_generators.py
@@ -41,7 +41,8 @@ def __call__(self, n_signal: float = 0) -> np.ndarray:
         rng = np.random.default_rng(self.config["random_seed"])
         n_background = rng.poisson(self.data_handler.n_background)
         if not self.config["fixed_ns"]:
-            n_signal = rng.poisson(self.data_handler.calculate_n_signal(n_signal))
+            n_signal = rng.poisson(
+                self.data_handler.calculate_n_signal(n_signal))
 
         background = self.data_handler.sample_background(n_background, rng)
         background["ra"] = rng.uniform(0, 2 * np.pi, len(background))
@@ -58,11 +59,13 @@ def __call__(self, n_signal: float = 0) -> np.ndarray:
         # not present in the data events. These include the true direction,
         # energy, and 'oneweight'.
         signal = rf.drop_fields(
-            signal, [n for n in signal.dtype.names if n not in background.dtype.names]
-        )
+            signal, [
+                n for n in signal.dtype.names if (
+                                    n not in background.dtype.names)])
 
         # Combine the signal background events and time-sort them.
-        # Use recfunctions.stack_arrays to prevent numpy from scrambling entry order
+        # Use recfunctions.stack_arrays to prevent numpy from scrambling entry
+        # order
         if background.dtype == signal.dtype:
             return np.concatenate([background, signal])
         else:
diff --git a/mla/utility_functions.py b/mla/utility_functions.py
index a4eb7f4..e2951ff 100644
--- a/mla/utility_functions.py
+++ b/mla/utility_functions.py
@@ -89,12 +89,12 @@ def rotate(
     ra3 = np.atleast_1d(ra3)
     dec3 = np.atleast_1d(dec3)
 
-    if not len(ra1) == len(dec1) == len(ra2) == len(dec2) == len(ra3) == len(dec3):
+    if not len(ra1) == len(dec1) == len(
+            ra2) == len(dec2) == len(ra3) == len(dec3):
         raise IndexError("Arguments must all have the same dimension.")
 
-    cos_alpha = np.cos(ra2 - ra1) * np.cos(dec1) * np.cos(dec2) + np.sin(dec1) * np.sin(
-        dec2
-    )
+    cos_alpha = np.cos(ra2 - ra1) * np.cos(dec1) * \
+        np.cos(dec2) + np.sin(dec1) * np.sin(dec2)
 
     # correct rounding errors
     cos_alpha[cos_alpha > 1] = 1
@@ -134,7 +134,11 @@ def rotate(
     return r_a, dec
 
 
-def angular_distance(src_ra: float, src_dec: float, r_a: float, dec: float) -> float:
+def angular_distance(
+        src_ra: float,
+        src_dec: float,
+        r_a: float,
+        dec: float) -> float:
     """Computes angular distance between source and location.
 
     Args:
@@ -235,7 +239,8 @@ def newton_method_multidataset(
     return x[i + 1]
 
 
-def trimsim(sim: np.ndarray, fraction: float, scaleow: bool = True) -> np.ndarray:
+def trimsim(sim: np.ndarray, fraction: float,
+            scaleow: bool = True) -> np.ndarray:
     """Keep only fraction of the simulation
 
     Args: