Deprecate generate, simulate and collate

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!
+#Just copy & paste the output into LaTeX!

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.

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"]

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)
+        if detector == 'all':
+             detector = list(self.detectors)        
+
+        return self.read_detector(detector,material).run(flux, detector_effects=detector_effects)

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)