diff --git a/doc/scripts/SNEWS2.0_rate_table.py b/doc/scripts/SNEWS2.0_rate_table.py
index 9369936d0..9368b01fb 100644
--- a/doc/scripts/SNEWS2.0_rate_table.py
+++ b/doc/scripts/SNEWS2.0_rate_table.py
@@ -1,6 +1,12 @@
 import numpy as np
 import os
-from snewpy import snowglobes
+
+from snewpy.rate_calculator import RateCalculator
+from snewpy.models.ccsn import Bollig_2016, OConnor_2015
+from snewpy.neutrino import MassHierarchy, MixingParameters, ThreeFlavorMixingParameters
+from snewpy.flavor_transformation import AdiabaticMSW
+
+from astropy import units as u
 
 #Select output format, Mathjax or LaTeX
 while True:
@@ -14,59 +20,46 @@
     else:
         print('Please enter 1 or 0.')
 
-home_directory = os.getcwd()
-SNOwGLoBES_path = None  # change to SNOwGLoBES directory if using a custom detector configuration
-SNEWPY_models_base = "Your/SNEWPY/Models/Path/Here"  # local directory containing model input files ("SNEWPY_models")
+models = { 's11.2' : Bollig_2016(progenitor_mass=11.2<<u.Msun), 
+           's27.0' : Bollig_2016(progenitor_mass=27<<u.Msun),  
+           's40' : OConnor_2015(progenitor_mass=40<<u.Msun) }
 
-d = 10  # distance of supernova in kpc
+transformations = { 'AdiabaticMSW_NMO' : AdiabaticMSW(MixingParameters('NORMAL')), 
+                    'AdiabaticMSW_IMO' : AdiabaticMSW(MixingParameters('INVERTED')) }
 
-dets = ["wc100kt30prct", "ar40kt", "halo1", "halo2", "scint20kt", "novaFD", "icecube", "km3net", "ds20", "xent", "lz",
-        "pandax"]
-ref_mass = {"wc100kt30prct": 100, "ar40kt": 40, "halo1": 0.079, "halo2": 1, "scint20kt": 20, "novaFD": 14,
-            "icecube": 51600, "km3net": 69366 * 3, "ds20": 0.0386, "xent": 0.006, "lz": 0.007, "pandax": 0.004}
+detectors = ["wc100kt30prct", "ar40kt", "halo1", "halo2", "scint20kt", "novaFD", 
+             "icecube", "km3net", "ds20", "xent", "lz", "pandax"]
+        
+detector_masses = {"wc100kt30prct": 100, "ar40kt": 40, "halo1": 0.079, "halo2": 1, "scint20kt": 20, "novaFD": 14,
+                   "icecube": 51600, "km3net": 69366 * 3, "ds20": 0.0386, "xent": 0.006, "lz": 0.007, "pandax": 0.004}
 
-models = {}
-models['s11.2'] = {'type': 'Bollig_2016', 'file_name': 's11.2c'}
-models['s27.0'] = {'type': 'Bollig_2016', 'file_name': 's27.0c'}
-models['s40'] = {'type': 'OConnor_2015', 'file_name': 'M1_neutrinos.dat'}
+detector_effects = {'smeared' : True, 'unsmeared' : False}
 
-transformations = ['AdiabaticMSW_NMO', 'AdiabaticMSW_IMO']
+rc=RateCalculator()
 
 total_events = {}
 
-have_data_saved = False
-if (have_data_saved is False):
-    # Running the modules
+energies = np.linspace(0,100,501)<<u.MeV
+distance = 10*u.kpc
+        
+# Running the modules
+for effects in detector_effects:
+    total_events[effects] = {}
     for model in models:
-        total_events[model] = {}
+        times = models[model].get_time()
+        total_events[effects][model] = {}
         for transformation in transformations:
-            total_events[model][transformation] = {}
-            file_name = models[model]['file_name']
-            modeltype = models[model]['type']
-            outfile = modeltype + "_" + model + "_summed_" + transformation
-            model_dir = SNEWPY_models_base + "/" + modeltype + "/"
-
-            tarredfile = snowglobes.generate_fluence(model_dir + file_name, modeltype, transformation, d, outfile)
-            for det in dets:
-                snowglobes.simulate(SNOwGLoBES_path, tarredfile, detector_input = det)
-                tables = snowglobes.collate(SNOwGLoBES_path, tarredfile, skip_plots = True)
-
-                # for our table, interesting number is the smeared total number of events
-                key = "Collated_" + outfile + "_" + det + "_events_smeared_weighted.dat"
-                total_events[model][transformation][det + "smeared"] = 0
-                for j in range(1, len(tables[key]['header'].split())):
-                    total_events[model][transformation][det + "smeared"] += sum(tables[key]['data'][j])
-
-                key = "Collated_" + outfile + "_" + det + "_events_unsmeared_weighted.dat"
-                total_events[model][transformation][det + "unsmeared"] = 0
-                for j in range(1, len(tables[key]['header'].split())):
-                    total_events[model][transformation][det + "unsmeared"] += sum(tables[key]['data'][j])
-
-    os.chdir(home_directory)
-    np.save("SNEWS2.0_whitepaper_table_data.npy", total_events)
-else:
-    total_events = np.load("SNEWS2.0_whitepaper_table_data.npy", allow_pickle = True).tolist()
+            total_events[effects][model][transformation] = {}
+            #get the flux from the model
+            flux = model.get_flux(t=times, E=energies, distance=distance, flavor_xform=transformation)
+            fluence = flux.integrate('time')
+            for detector in detectors:
+                events = rc.run(fluence, detectors[detector], detector_effects=detector_effects[effects])
+                total_events[effects][model][transformation][detector] = sum([chan.integrate_or_sum('energy').array.squeeze().value for chan in events.values()])     
 
+home_directory = os.getcwd()
+os.chdir(home_directory)
+np.savez("SNEWS2.0_whitepaper_table_data.npz", total_events)
 
 # Now lets make the table:
 def round_to_2(x):
@@ -76,11 +69,11 @@ def round_to_2(x):
         return round(x, -int(np.floor(np.log10(np.abs(x)))) + 1)
 
 
-det_maps = {"Super-K": "wc100kt30prct", "Hyper-K": "wc100kt30prct", "IceCube": "icecube", "KM3NeT":"km3net",
-            "LVD": "scint20kt", "KamLAND": "scint20kt", "Borexino": "scint20kt", "JUNO": "scint20kt",
-            "SNO+": "scint20kt", "NO${\\nu}$A": "novaFD", "HALO": "halo1", "HALO-1kT": "halo2", "DUNE": "ar40kt",
-            "MicroBooNe": "ar40kt", "SBND": "ar40kt", "Baksan": "scint20kt", "DarkSide-20k": "ds20", "XENONnT": "xent",
-            "LZ": "lz", "PandaX-4T": "pandax"}
+detector_maps = {"Super-K": "wc100kt30prct", "Hyper-K": "wc100kt30prct", "IceCube": "icecube", "KM3NeT":"km3net",
+                 "LVD": "scint20kt", "KamLAND": "scint20kt", "Borexino": "scint20kt", "JUNO": "scint20kt",
+                 "SNO+": "scint20kt", "NO${\\nu}$A": "novaFD", "HALO": "halo1", "HALO-1kT": "halo2", "DUNE": "ar40kt",
+                 "MicroBooNe": "ar40kt", "SBND": "ar40kt", "Baksan": "scint20kt", "DarkSide-20k": "ds20", "XENONnT": "xent",
+                 "LZ": "lz", "PandaX-4T": "pandax"}
 
 data = {}
 data['Experiment'] = ['Super-K', 'Hyper-K', 'IceCube', 'KM3NeT', 'LVD', 'KamLAND', 'Borexino',
@@ -106,15 +99,15 @@ def round_to_2(x):
 
 for experiment in range(len(data['Experiment'])):
     mass = data['Mass [kt]'][experiment]
-    dettype = det_maps[data['Experiment'][experiment]]
-    base_mass = ref_mass[dettype]
+    detector_type = detector_maps[data['Experiment'][experiment]]
+    base_mass = detector_masses[detector_type]
 
-    counts_LCN = int(total_events['s11.2']['AdiabaticMSW_NMO'][dettype + "smeared"] * mass / base_mass)
-    counts_LCI = int(total_events['s11.2']['AdiabaticMSW_IMO'][dettype + "smeared"] * mass / base_mass)
-    counts_MCN = int(total_events['s27.0']['AdiabaticMSW_NMO'][dettype + "smeared"] * mass / base_mass)
-    counts_MCI = int(total_events['s27.0']['AdiabaticMSW_IMO'][dettype + "smeared"] * mass / base_mass)
-    counts_HCN = int(total_events['s40']['AdiabaticMSW_NMO'][dettype + "smeared"] * mass / base_mass)
-    counts_HCI = int(total_events['s40']['AdiabaticMSW_IMO'][dettype + "smeared"] * mass / base_mass)
+    counts_LCN = int(total_events['smeared']['s11.2']['AdiabaticMSW_NMO'][detector_type] * mass / base_mass)
+    counts_LCI = int(total_events['smeared']['s11.2']['AdiabaticMSW_IMO'][detector_type] * mass / base_mass)
+    counts_MCN = int(total_events['smeared']['s27.0']['AdiabaticMSW_NMO'][detector_type] * mass / base_mass)
+    counts_MCI = int(total_events['smeared']['s27.0']['AdiabaticMSW_IMO'][detector_type] * mass / base_mass)
+    counts_HCN = int(total_events['smeared']['s40']['AdiabaticMSW_NMO'][detector_type] * mass / base_mass)
+    counts_HCI = int(total_events['smeared']['s40']['AdiabaticMSW_IMO'][detector_type] * mass / base_mass)
 
     post = ['', '', '', '', '', '']
     if counts_LCN > 10000:
@@ -149,15 +142,15 @@ def round_to_2(x):
 # the unweighted mass (the entry in SNOwGLoBES), see below for details.  Here we take the
 # effective mass of the s27 normal scenario and discuss the range in the table caption.
 
-dettype = 'icecube'
+detector_type = 'icecube'
 mass = 51600
-data['Mass [kt]'][2] = "~"+str(int(round(mass*total_events['s27.0']['AdiabaticMSW_NMO'][dettype+"smeared"]/
-      total_events['s27.0']['AdiabaticMSW_NMO'][dettype+"unsmeared"], -2)))+"*"
+data['Mass [kt]'][2] = "~"+str(int(round(mass*total_events['smeared']['s27.0']['AdiabaticMSW_NMO'][detector_type]/
+      total_events['unsmeared']['s27.0']['AdiabaticMSW_NMO'][detector_type], -2)))+"*"
 
-dettype = 'km3net'
+detector_type = 'km3net'
 mass = 69366 * 3
-data['Mass [kt]'][3] = "~"+str(int(round(mass*total_events['s27.0']['AdiabaticMSW_NMO'][dettype+"smeared"]/
-      total_events['s27.0']['AdiabaticMSW_NMO'][dettype+"unsmeared"], -1)))+"*"
+data['Mass [kt]'][3] = "~"+str(int(round(mass*total_events['smeared']['s27.0']['AdiabaticMSW_NMO'][detector_type]/
+      total_events['unsmeared']['s27.0']['AdiabaticMSW_NMO'][detector_type], -1)))+"*"
 
 # Formatting the dictionary to be compatible with LaTeX & MathJax (useful for html)
 def dictArray(dictionary):
@@ -208,4 +201,4 @@ def dictArray(dictionary):
     print("\\usepackage{amsmath}\n\\begin{document}")   #loads forever without amsmath package
     print(f"$$\n{dictArray(data)}\n$$")
     print("\\end{document}")
-#Just copy & paste the output into LaTeX!
\ No newline at end of file
+#Just copy & paste the output into LaTeX!
diff --git a/doc/scripts/TimeSeries.py b/doc/scripts/TimeSeries.py
index a9fc21c35..290855bf4 100644
--- a/doc/scripts/TimeSeries.py
+++ b/doc/scripts/TimeSeries.py
@@ -1,18 +1,29 @@
 #!/usr/bin/env python
 
-from snewpy import snowglobes
+from snewpy.models.ccsn import Nakazato_2013
+from snewpy.neutrino import MassHierarchy, MixingParameters
+from snewpy.flavor_transformation import AdiabaticMSW
+from snewpy.rate_calculator import RateCalculator
 
-SNOwGLoBES_path = None  # change to SNOwGLoBES directory if using a custom detector configuration
+import numpy as np
+import astropy.units as u
 
-# arguments for generate_time_series
-model_file = "/path/to/snewpy/models/Nakazato_2013/nakazato-LS220-BH-z0.004-s30.0.fits"
-modeltype = 'Nakazato_2013'
-transformation = 'AdiabaticMSW_NMO'
-d = 10  # Supernova distance in kpc
+model = Nakazato_2013(progenitor_mass=30*u.solMass, revival_time=0*u.ms, metallicity=0.004, eos='LS220')
 
-# Running the modules
-outfile = snowglobes.generate_time_series(model_file, modeltype, transformation, d)
-snowglobes.simulate(SNOwGLoBES_path, outfile, detector_input="icecube")
-snowglobes.collate(SNOwGLoBES_path, outfile)
+transformation = AdiabaticMSW(MixingParameters('NORMAL')) # Desired flavor transformation
+       
+times    = model.get_time()
+energies = np.linspace(0,100,501)<<u.MeV
+distance = 10*u.kpc
+
+flux = model.get_flux(t=times, E=energies, distance=distance, flavor_xform=transformation)
+fluence = flux.integrate('time')
+
+detector = "wc100kt30prct"
+rc = RateCalculator()
+events = rc.run(fluence, detector, detector_effects=True)
+
+with channel in events:
+    filename = f"{model}.{transformation}.{times[0]:.3f}-{times[-1]:.3f},{energies[0]:.3f}-{energies[-1]:.3f},{distance:.3f}.{channel}.npz"
+    events[channel].save(filename)
 
-# An additional, optional argument in simulate() is the detector name, if one wants to only run 1 detector, rather than all of them.
diff --git a/python/snewpy/flux.py b/python/snewpy/flux.py
index 53ce54331..1c329e370 100644
--- a/python/snewpy/flux.py
+++ b/python/snewpy/flux.py
@@ -510,7 +510,7 @@ def plot(flux, projection='energy', styles=None, **kwargs):
         plt.xlabel(f'{projection}, {x.unit._repr_latex_()}')
         plt.ylabel(f'{fP.__class__.__name__}, {x.unit._repr_latex_()}')
         return lines
-
+            
 #some standard container classes that can be used for 
 Flux = Container['1/(MeV*s*m**2)', "d2FdEdT"]
 Fluence = Container[Flux.unit*u.s, "dFdE"]
diff --git a/python/snewpy/rate_calculator.py b/python/snewpy/rate_calculator.py
index 346f7e2eb..037660dff 100644
--- a/python/snewpy/rate_calculator.py
+++ b/python/snewpy/rate_calculator.py
@@ -409,8 +409,9 @@ def run(self, flux:Container, detector:str, material:str=None, detector_effects:
             The incoming neutrino flux (or fluence).
         
         detector: str
-            Name of the detector to calculate the rate. 
+            Name of the detector to calculate the rate, or 'all'.
             Check `RateCalculator.detectors` for the list of options
+            If 'all', all the detector options in SNOwGLoBES will be run            
         
         material: str or None
             Name of the detector material. 
@@ -425,4 +426,7 @@ def run(self, flux:Container, detector:str, material:str=None, detector_effects:
             dict[str, Container]
                 A dictionary with interaction rates (as instances of :class:`snewpy.flux.Container`) for each channel.
         """
-        return self.read_detector(detector,material).run(flux, detector_effects=detector_effects)
\ No newline at end of file
+        if detector == 'all':
+             detector = list(self.detectors)        
+            
+        return self.read_detector(detector,material).run(flux, detector_effects=detector_effects)
diff --git a/python/snewpy/snowglobes.py b/python/snewpy/snowglobes.py
index b15781926..9f258b603 100644
--- a/python/snewpy/snowglobes.py
+++ b/python/snewpy/snowglobes.py
@@ -109,6 +109,53 @@ def _get_model_class(model_type: str):
     except KeyError:
         raise ValueError(f"Model '{model_type}' not found.")
 
+def get_model_class(model_type: str):
+    """Look up model class corresponding to the given model name.
+
+    Parameters
+    ---------
+    model_type : str
+        Model name
+
+    Returns
+    -------
+    Model class corresponding to the given model name
+    """    
+    models_dict = {}
+    modules_list = ["snewpy.models.ccsn", "snewpy.models.presn"]
+    for module_name in modules_list:
+        module = importlib.import_module(module_name)
+        models_dict.update({k:v for k,v in vars(module).items() if isclass(v)})
+
+    try:
+        return models_dict[model_type]
+    except KeyError:
+        raise ValueError(f"Model '{model_type}' not found.")
+
+def get_model_loader(model_type: str):
+    """Look up model class corresponding to the given model name.
+
+    Parameters
+    ---------
+    model_type : str
+        Model name
+
+    Returns
+    -------
+    Model class corresponding to the given model name
+    """    
+    models_dict = {}
+    modules_list = ["snewpy.models.ccsn_loaders", "snewpy.models.presn_loaders"]
+    for module_name in modules_list:
+        module = importlib.import_module(module_name)
+        models_dict.update({k:v for k,v in vars(module).items() if isclass(v)})
+    models_dict['Analytic3Species'] = snewpy.models.ccsn.Analytic3Species
+    
+    try:
+        return models_dict[model_type]
+    except KeyError:
+        raise ValueError(f"Model loader'{model_type}' not found.")
+
 def generate_time_series(model_path, model_type, flavor_transformation, d, output_filename=None, ntbins=30, deltat=None, snmodel_dict={}):
     """Generate time series files in SNOwGLoBES format.
 
@@ -139,14 +186,16 @@ def generate_time_series(model_path, model_type, flavor_transformation, d, outpu
     str
         Path of NumPy archive file with neutrino fluence data.
     """
-    model_class = _get_model_class(model_type)
-
+    
+    warn("generate_time_series is deprecated. Use the SupernovaModel.get_flux() method instead.", DeprecationWarning, stacklevel=2)
+    
+    model_loader = get_model_loader(model_type)
+    model_dir, model_file = os.path.split(os.path.abspath(model_path))
+    snmodel = model_loader(model_path, **snmodel_dict)
+    
     # if flavor_transformation is a string, find the appropriate class
     if isinstance(flavor_transformation, str):
-        flavor_transformation = _get_transformation(flavor_transformation)
-
-    model_dir, model_file = os.path.split(os.path.abspath(model_path))
-    snmodel = model_class(model_path, **snmodel_dict)
+        flavor_transformation = get_transformation(flavor_transformation)
 
     # Subsample the model time. Default to 30 time slices.
     tmin = snmodel.get_time()[0]
@@ -201,18 +250,16 @@ def generate_fluence(model_path, model_type, flavor_transformation, d, output_fi
     str
         Path of NumPy archive file with neutrino fluence data.
     """
-    try:
-        model_class = getattr(snewpy.models.ccsn_loaders, model_type)
-    except AttributeError as e:
-        logging.warn(e)
-        model_class = getattr(snewpy.models.ccsn, model_type)
-
+    
+    warn("generate_fluence is deprecated. Use the SupernovaModel.get_flux() method instead.", DeprecationWarning, stacklevel=2)
+    
+    model_loader = get_model_loader(model_type)
+    model_dir, model_file = os.path.split(os.path.abspath(model_path))
+    snmodel = model_loader(model_path, **snmodel_dict)
+    
     # if flavor_transformation is a string, find the appropriate class
     if isinstance(flavor_transformation, str):
-        flavor_transformation = _get_transformation(flavor_transformation)
-
-    model_dir, model_file = os.path.split(os.path.abspath(model_path))
-    snmodel = model_class(model_path, **snmodel_dict)
+        flavor_transformation = get_transformation(flavor_transformation)
 
     #set the timings up
     #default if inputs are None: full time window of the model
@@ -261,6 +308,9 @@ def simulate(SNOwGLoBESdir, tarball_path, detector_input="all", *, detector_effe
     detector_effects : bool
          Whether to account for detector smearing and efficiency.
     """
+    
+    warn("The simulate function is deprecated. Use the RateCalculator.run() method instead.", DeprecationWarning, stacklevel=2)
+    
     rc = RateCalculator(base_dir=SNOwGLoBESdir)
     if detector_input == 'all':
         detector_input = list(rc.detectors)
@@ -349,6 +399,8 @@ def collate(tarball_path, skip_plots=False, *, smearing=True):
         Dictionary of data tables: One table per time bin; each table contains in the first column the energy bins, in the remaining columns the number of events for each interaction channel in the detector.
     """
 
+    warn("The collate function is deprecated.", DeprecationWarning, stacklevel=2)    
+    
     def aggregate_channels(table, **patterns):
         #rearrange the table to have only channel column
         levels = list(table.columns.names)
diff --git a/python/snewpy/test/_rate_crosscheck_table.py b/python/snewpy/test/_rate_crosscheck_table.py
index 4d75b8297..05cf454d9 100644
--- a/python/snewpy/test/_rate_crosscheck_table.py
+++ b/python/snewpy/test/_rate_crosscheck_table.py
@@ -1,27 +1,27 @@
 #generated by "make_rate_table.py" on 2023-06-02 14:05:49.504799
 rate_table= \
-{('Bollig_2016', 's11.2c', 'AdiabaticMSW_IMO'): {'ar40kt': 2509.8697955822163,
+{('Bollig_2016', 11.2, 'AdiabaticMSW_IMO'): {'ar40kt': 2509.8697955822163,
                                                  'halo1': 3.7370536263087892,
                                                  'halo2': 47.30447628238974,
                                                  'icecube': 332587.17107218195,
                                                  'novaFD': 1955.8019842893325,
                                                  'novaND': 41.91004252048569,
                                                  'wc100kt30prct': 12723.21134007678},
- ('Bollig_2016', 's11.2c', 'AdiabaticMSW_NMO'): {'ar40kt': 2696.396378383758,
+ ('Bollig_2016', 11.2, 'AdiabaticMSW_NMO'): {'ar40kt': 2696.396378383758,
                                                  'halo1': 4.261668358305819,
                                                  'halo2': 53.94516909247873,
                                                  'icecube': 320766.2999953658,
                                                  'novaFD': 1947.9479073661355,
                                                  'novaND': 41.741740872131466,
                                                  'wc100kt30prct': 12645.788228369147},
- ('Bollig_2016', 's27.0c', 'AdiabaticMSW_IMO'): {'ar40kt': 5192.001233928208,
+ ('Bollig_2016', 27.0, 'AdiabaticMSW_IMO'): {'ar40kt': 5192.001233928208,
                                                  'halo1': 8.180298088573792,
                                                  'halo2': 103.54807707055436,
                                                  'icecube': 659449.9649827218,
                                                  'novaFD': 3638.7825760092305,
                                                  'novaND': 77.97391234305493,
                                                  'wc100kt30prct': 23596.070817047243},
- ('Bollig_2016', 's27.0c', 'AdiabaticMSW_NMO'): {'ar40kt': 5520.903869068218,
+ ('Bollig_2016', 27.0, 'AdiabaticMSW_NMO'): {'ar40kt': 5520.903869068218,
                                                  'halo1': 9.105675711861647,
                                                  'halo2': 115.2617178716664,
                                                  'icecube': 661517.4858901916,
diff --git a/python/snewpy/test/simplerate_integrationtest.py b/python/snewpy/test/simplerate_integrationtest.py
index 46f94dd47..e76030a58 100644
--- a/python/snewpy/test/simplerate_integrationtest.py
+++ b/python/snewpy/test/simplerate_integrationtest.py
@@ -2,74 +2,38 @@
 """Integration test based on SNEWS2.0_rate_table_singleexample.py
 """
 import unittest
-from snewpy import snowglobes
-from snewpy import model_path 
-
+from snewpy.rate_calculator import RateCalculator
 from snewpy.models import ccsn
 import astropy.units as u
 
-def preload_model(name:str, **parameters):
-    #initialize the model with given name and parameters
-    model = ccsn.__dict__[name](**parameters)
-    return model
-
-
 class TestSimpleRate(unittest.TestCase):
-
+    
     def test_simplerate(self):
         """Integration test based on SNEWS2.0_rate_table_singleexample.py
         """
-        # Hardcoded paths on GitHub Action runner machines
-        SNOwGLoBES_path = None
-
-        distance = 10  # Supernova distance in kpc
-        detector = "wc100kt30prct" #SNOwGLoBES detector for water Cerenkov
-        modeltype = 'Bollig_2016' # Model type from snewpy.models
-        model = 's11.2c' # Name of model
-        transformation = 'AdiabaticMSW_NMO' # Desired flavor transformation
-
-        # Construct file system path of model file and name of output file
-        model_file_path = f'{model_path}/{modeltype}/{model}'
-        outfile = f'{modeltype}_{model}_{transformation}'
-
-        #make sure the model files are loaded
-        preload_model(modeltype, progenitor_mass=11.2*u.Msun)
-        
+        model = ccsn.Bollig_2016(progenitor_mass=11.2<<u.Msun) # SN model
+        transformation = AdiabaticMSW(MixingParameters('NORMAL')) # Desired flavor transformation
+       
         # Now, do the main work:
         print("Generating fluence files ...")
-        tarredfile = snowglobes.generate_fluence(model_file_path, modeltype, transformation, distance, outfile)
-
-        print("Simulating detector effects with SNOwGLoBES ...")
-        snowglobes.simulate(SNOwGLoBES_path, tarredfile, detector_input=detector, detector_effects=True)
-
-        print("Collating results ...")
-        tables = snowglobes.collate(tarredfile, skip_plots=True, smearing=True)
-
-        # Use results to print the number of events in different interaction channels
-        key = f"Collated_{outfile}_{detector}_events_unsmeared_weighted.dat"
-        total_events = 0
-        for i, channel in enumerate(tables[key]['header'].split()):
-            if i == 0:
-                continue
-            n_events = sum(tables[key]['data'][i])
-            total_events += n_events
-            print(f"{channel:10}: {n_events:.3f} events")
-
-        #Super-K has 32kT inner volume
-        print("Total events in Super-K-like detector:" , 0.32*total_events)
-
-        # Use results to print the number of events in different interaction channels
-        # with efficiency and smearing
-        key = f"Collated_{outfile}_{detector}_events_smeared_weighted.dat"
-        total_events_smeared = 0
-        for i, channel in enumerate(tables[key]['header'].split()):
-            if i == 0:
-                continue
-            n_events = sum(tables[key]['data'][i])
-            total_events_smeared += n_events
-            print(f"{channel:10}: {n_events:.3f} events (smeared)")
+        times    = model.get_time()
+        energies = np.linspace(0,100,501)<<u.MeV
+        distance = 10*u.kpc
+        #get the flux from the model
+        flux = model.get_flux(t=times, E=energies, distance=distance, flavor_xform=transformation)
+        fluence = flux.integrate('time')
+
+        print("Simulating detector effects ...")
+        detector = "wc100kt30prct"
+        rc = RateCalculator()
+
+        events = rc.run(fluence, detector, detector_effects=False)        
+        events_smeared = rc.run(fluence, detector, detector_effects=True)
+        
+        # Compute number of events in all interaction channels
+        total_events  = sum([chan.integrate_or_sum('energy').array.squeeze().value for chan in events.values()])        
+        total_events_smeared  = sum([chan.integrate_or_sum('energy').array.squeeze().value for chan in events_smeared.values()])
 
-        #Super-K has 32kT inner volume
         print("Total events in Super-K-like detector (with smearing):" , 0.32*total_events_smeared)
 
         # We do not use the SNOwGLoBES scaling factors but use other constants so we do not
diff --git a/python/snewpy/test/test_05_snowglobes.py b/python/snewpy/test/test_05_snowglobes.py
index 8dda57b82..a4790fc31 100644
--- a/python/snewpy/test/test_05_snowglobes.py
+++ b/python/snewpy/test/test_05_snowglobes.py
@@ -1,41 +1,48 @@
 import pytest
 from pathlib import Path
 from snewpy.test._rate_crosscheck_table import rate_table
-from snewpy import snowglobes, model_path
+from snewpy import snowglobes
+from snewpy.rate_calculator import RateCalculator
 
-from snewpy.models import ccsn
+from snewpy.models.ccsn import *
 import astropy.units as u
 
 pytestmark=pytest.mark.snowglobes
 
 #get available model parameters from table
 param_values = list(rate_table.keys())
+
+rc = RateCalculator()
 #get available detectors from table
 detectors = list(list(rate_table.values())[0].keys())
 
 #make sure the model files are loaded
-model = ccsn.Bollig_2016
+model = Bollig_2016
 for params in model.get_param_combinations():
     model(**params)
     
-def fluence_calculation(model_name,model_file,transform):
+def fluence_calculation(model_name,model_mass,transform):
     #generating fluence file
-    model_file_path = f'{model_path}/{model_name}/{model_file}'
-    print(model_file_path)
-    return snowglobes.generate_fluence(model_file_path, model_name, transform,d=10)
+    model = snowglobes.get_model_class(model_name)(model_mass<<u.Msun)    
+    times    = model.get_time()
+    energies = np.linspace(0,100,501)<<u.MeV
+    distance = 10*u.kpc
+    flux = model.get_flux(t=times, E=energies, distance=distance, flavor_xform=transformation)
+    fluence = flux.integrate('time')    
+    return fluence
 
-def rates_calculation(fluence_file):
-    tables = snowglobes.simulate(None,fluence_file,detector_input=detectors)
+def rates_calculation(fluence):
+    table = rc.run(fluence, detectors, detector_effects=True)
     result = {}
-    for det,table in tables.items():
-        table = list(table.values())[0] #take the first and the only time bin
-        result[det] = table['weighted']['smeared'].values.sum()
+    for det in table:
+        result[det] += sum([chan.integrate_or_sum('energy').array.squeeze().value for chan in events[det].values()])
+        
     return result
 
 @pytest.mark.parametrize('model_parameters',param_values)
 def test_total_rate_equals_table_value(model_parameters):
-    fluence_file = fluence_calculation(*model_parameters)
-    calculated_rates  = rates_calculation(fluence_file)
+    fluence = fluence_calculation(*model_parameters)
+    calculated_rates  = rates_calculation(fluence)
     for detector in detectors:
         expected = pytest.approx(rate_table[model_parameters][detector], rel=0.01)
         assert calculated_rates[detector] == expected, f"Crosscheck failed for {detector}"
diff --git a/python/snewpy/test/test_presn_rates.py b/python/snewpy/test/test_presn_rates.py
index 0264145c0..83ddabd9b 100644
--- a/python/snewpy/test/test_presn_rates.py
+++ b/python/snewpy/test/test_presn_rates.py
@@ -28,5 +28,5 @@ def test_presn_rate(model_class, model_params, transformation, detector):
     model = model_class(**model_params)
     flux = model.get_flux(T, E, distance=distance, flavor_xform=transformation)
     rate = rc.run(flux, detector='scint20kt', detector_effects=False)['ibd']
-    ibd_events = rate.integrate_or_sum('time').integrate_or_sum('energy').array.squeeze()
-    assert 10<ibd_events<1000
\ No newline at end of file
+    ibd_events = rate.integrate_or_sum('time').integrate_or_sum('energy').array.squeeze().value
+    assert 10<ibd_events<1000
diff --git a/python/snewpy/test/test_rate_calculation.py b/python/snewpy/test/test_rate_calculation.py
index a7249f68b..07a6b6172 100644
--- a/python/snewpy/test/test_rate_calculation.py
+++ b/python/snewpy/test/test_rate_calculation.py
@@ -49,7 +49,7 @@ def test_rate_unsmeared(rc, fluence, detector, expected_total):
     #calculate the event rates
     rates = rc.run(fluence, detector, detector_effects=False)
     #calculate the total IBD event rate 
-    N_total = sum([r.integrate_or_sum('energy').array.value for r in rates.values()])
+    N_total = sum([r.integrate_or_sum('energy').array.squeeze().value for r in rates.values()])
     print(N_total)
     #check the final value
     assert N_total == pytest.approx(expected_total, 0.01)