diff --git a/sample_plot_config.ini b/sample_plot_config.ini index 736ed57..e9c6923 100644 --- a/sample_plot_config.ini +++ b/sample_plot_config.ini @@ -1,17 +1,56 @@ +[General] +title = nuspacesim_run +rows = 1 +columns = 1 +figsize = (8,7) +pop_up = yes +save_to_file = no +save_as = pdf +default_color = 0 +# as in 'C0' +default_markersize = 3 +default_colormap = viridis +# output_path = ~/Documents/Pictures + +[Spectra] +spectra_histogram = yes + [Geometry] -geom_beta_tr_hist = no -#interaction_height = no (to be implemented) +geom_beta_tr_hist = yes +path_length_to_detector = yes +# interaction_height = no (to be implemented) +# geometry_overview = no (to be implemented) [Taus] -taus_scatter = no -taus_pexit = no -taus_histogram = yes -taus_overview = yes +energy_hists = yes +beta_hist = yes +tau_beta_hist = yes +tau_lorentz_hist = yes +tau_exit_prob_hist = yes +tau_beta_hex = yes +tau_lorentz_hex = yes +tau_exit_prob_hex = yes [EASOpitcal] -eas_optical_scatter = no -eas_optical_histogram = yes +numpes_vs_beta = yes +altdec_vs_numpes = yes +showerEnergy_vs_numPEs = yes +costhetacheff_vs_beta = yes +altDec_vs_costhetacheff = yes +showerEnergy_vs_costhetacheff = yes +numpes_hist = yes +costhetacheff_hist = yes [Detector] -#event_count_at_detector = no (to be implemented) -#detection_prob_vs_beta = no (to be implemented) +# event_count_at_detector = no (to be implemented) +# detection_prob_vs_beta = no (to be implemented) +# detector_overview = no (to be implemented) + +[Overview] +dashboard = yes +taus_density_beta_overview = yes +taus_histograms_overview = yes +taus_pexit_overview = yes +taus_overview = yes +eas_optical_density_overview = yes +eas_optical_histogram_overview = yes diff --git a/src/nuspacesim/apps/cli.py b/src/nuspacesim/apps/cli.py index 3c87076..f78b7dc 100644 --- a/src/nuspacesim/apps/cli.py +++ b/src/nuspacesim/apps/cli.py @@ -48,19 +48,9 @@ """ -import configparser - import click -from .. import NssConfig, SimulationParameters, simulation -from ..compute import compute -from ..config import FileSpectrum, MonoSpectrum, PowerSpectrum -from ..results_table import ResultsTable -from ..utils import plots -from ..utils.plot_function_registry import registry -from ..xml_config import config_from_xml, create_xml - -__all__ = ["create_config", "run", "show_plot"] +from .commands import create_config, run, show_plot @click.group() @@ -72,272 +62,9 @@ def cli(): # # ctx.obj["DEBUG"] = debug -@cli.command() -@click.option( - "-o", "--output", type=click.Path(exists=False), default=None, help="Output file." -) -@click.option( - "-p", - "--plot", - type=click.Choice(list(registry), case_sensitive=False), - multiple=True, - default=[], - help="Available plotting functions. Select multiple plots with multiple uses of -p", -) -@click.option( - "-P", - "--plotconfig", - type=click.Path( - exists=True, - dir_okay=False, - readable=True, - ), - help="Read selected plotting functions and options from the specified INI file", -) -@click.option("--plotall", is_flag=True, help="Show all result plots.") -@click.option( - "-w", - "--write-stages", - is_flag=True, - help="Write intermediate values after each simulation stage.", -) -@click.option( - "-n", - "--no-result-file", - is_flag=True, - help="Do not save the results to an output file.", -) -@click.option( - "--monospectrum", - type=float, - default=None, - help="Mono Energetic Spectrum Log Energy.", -) -@click.option( - "--powerspectrum", - nargs=3, - type=click.Tuple([float, float, float]), - default=None, - help="Power Spectrum index, lower_bound, upper_bound.", -) -@click.argument( - "config_file", - default=None, - type=click.Path(exists=True, dir_okay=False, readable=True), -) -@click.argument("count", type=float, default=0.0) -def run( - config_file: str, - count: float, - monospectrum, - powerspectrum, - no_result_file: bool, - output: str, - plot: list, - plotconfig: str, - plotall: bool, - write_stages: bool, -) -> None: - """Perform the full nuspacesim simulation. - - Main Simulator for nuspacesim. Given a XML configuration file, and - optionally a count of simulated nutrinos, runs nutrino simulation. - - \f - - Parameters - ---------- - config_file: str - XML configuration file for particular simulation particular. - count : int, optional - Number of thrown trajectories. Optionally override value in config_file. - spectrum_type : str, optional - Type of - output: str, optional - Name of the output file holding the simulation results. - plot: list, optional - Plot the simulation results. - plotconfig: str, optional - INI file to select plots for each module, as well as to specifiy global plot settings. - plotall: bool, optional - Plot all the available the simulation results plots. - no_result_file: bool, optional - Disables writing results to output files. - - - Examples - -------- - Command line usage of the run command may work with the following invocation. - - `nuspacesim run sample_input_file.xml 1e5 8 -o my_sim_results.fits` - """ - - # User Inputs - config = config_from_xml(config_file) - - config.simulation.N = int(config.simulation.N if count == 0.0 else count) - - if monospectrum is not None and powerspectrum is not None: - raise RuntimeError("Only one of --monospectrum or --powerspectrum may be used.") - if monospectrum is not None: - config.simulation.spectrum = MonoSpectrum(monospectrum) - if powerspectrum is not None: - config.simulation.spectrum = PowerSpectrum(*powerspectrum) - - plot = ( - list(registry) - if plotall - else read_plot_config(plotconfig) - if plotconfig - else plot - ) - simulation = compute( - config, - verbose=True, - to_plot=plot, - output_file=output, - write_stages=write_stages, - ) - - if not no_result_file: - simulation.write(output, overwrite=True) - - -@cli.command() -@click.option( - "-n", "--numtrajs", type=float, default=100, help="number of trajectories." -) -@click.option( - "--monospectrum", - type=float, - default=None, - help="Mono Energetic Spectrum Log Energy.", -) -@click.option( - "--powerspectrum", - nargs=3, - type=click.Tuple([float, float, float]), - default=None, - help="Power Spectrum index, lower_bound, upper_bound.", -) -@click.argument("filename") -def create_config(filename: str, numtrajs: float, monospectrum, powerspectrum) -> None: - """Generate a configuration file from the given parameters. - - \f - - Parameters - ---------- - filename: str - Name of output xml configuration file. - numtrajs: float, optional - Number of thrown trajectories. Optionally override value in config_file. - - Examples - -------- - Command line usage of the create_config command may work with the following invocation. - - `nuspacesim create_config -n 1e5 sample_input_file.xml` - """ - if monospectrum is not None and powerspectrum is not None: - raise RuntimeError("Only one of --monospectrum or --powerspectrum may be used.") - - spec = MonoSpectrum() - - if monospectrum is not None: - spec = MonoSpectrum(monospectrum) - if powerspectrum is not None: - spec = PowerSpectrum(*powerspectrum) - # spec = FileSpectrum() - - simulation = SimulationParameters(N=int(numtrajs), spectrum=spec) - - create_xml(filename, NssConfig(simulation=simulation)) - - -@cli.command() -@click.option( - "-p", - "--plot", - type=click.Choice(list(registry), case_sensitive=False), - multiple=True, - default=[], - help="Available plotting functions. Select multiple plots with multiple uses of -p", -) -@click.option( - "-P", - "--plotconfig", - type=click.Path( - exists=True, - dir_okay=False, - readable=True, - ), - help="Read selected plotting functions and options from the specified INI file", -) -@click.option("--plotall", is_flag=True, help="Show all result plots.") -@click.argument( - "simulation_file", - default=None, - type=click.Path(exists=True, dir_okay=False, readable=True), -) -def show_plot( - simulation_file: str, - plot: list, - plotconfig: str, - plotall: bool, -) -> None: - """Show predefined plots of results in simulation file. - - \f - - Parameters - ---------- - simulation_file: str - input ResultsTable fits file. - plot: list, optional - Plot the simulation results. - plotconfig: str, optional - INI file to select plots for each module, as well as to specifiy global plot settings. - plotall: bool, optional - Plot all the available the simulation results plots. - - - Examples - -------- - Command line usage of the run command may work with the following invocation. - - `nuspacesim show_plot my_sim_results.fits -p taus_overview` - """ - - sim = ResultsTable.read(simulation_file) - - plot = ( - list(registry) - if plotall - else read_plot_config(plotconfig) - if plotconfig - else plot - ) - - simulation.geometry.region_geometry.show_plot(sim, plot) - simulation.taus.taus.show_plot(sim, plot) - simulation.eas_optical.eas.show_plot(sim, plot) - plots.show_plot(sim, plot) - - -def read_plot_config(filename): - plot_list = [] - cfg = configparser.ConfigParser() - cfg.read(filename) - for sec in cfg.sections()[1:]: - for key in cfg[sec]: - try: - if cfg[sec].getboolean(key): - plot_list.append(key) - except Exception as e: - print(e, "Config file contains non-valid option") - return plot_list - +cli.add_command(run.run) +cli.add_command(create_config.create_config) +cli.add_command(show_plot.show_plot) if __name__ == "__main__": cli() diff --git a/src/nuspacesim/apps/commands/create_config.py b/src/nuspacesim/apps/commands/create_config.py new file mode 100644 index 0000000..89762a9 --- /dev/null +++ b/src/nuspacesim/apps/commands/create_config.py @@ -0,0 +1,57 @@ +import click + +from ... import NssConfig, SimulationParameters +from ...config import MonoSpectrum, PowerSpectrum # FileSpectrum, +from ...xml_config import create_xml + + +@click.command() +@click.option( + "-n", "--numtrajs", type=float, default=100, help="number of trajectories." +) +@click.option( + "--monospectrum", + type=float, + default=None, + help="Mono Energetic Spectrum Log Energy.", +) +@click.option( + "--powerspectrum", + nargs=3, + type=click.Tuple([float, float, float]), + default=None, + help="Power Spectrum index, lower_bound, upper_bound.", +) +@click.argument("filename") +def create_config(filename: str, numtrajs: float, monospectrum, powerspectrum) -> None: + """Generate a configuration file from the given parameters. + + \f + + Parameters + ---------- + filename: str + Name of output xml configuration file. + numtrajs: float, optional + Number of thrown trajectories. Optionally override value in config_file. + + Examples + -------- + Command line usage of the create_config command may work with the following invocation. + + `nuspacesim create_config -n 1e5 sample_input_file.xml` + """ + if monospectrum is not None and powerspectrum is not None: + raise RuntimeError("Only one of --monospectrum or --powerspectrum may be used.") + + spec = MonoSpectrum() + + if monospectrum is not None: + spec = MonoSpectrum(monospectrum) + if powerspectrum is not None: + spec = PowerSpectrum(*powerspectrum) + # spec = FileSpectrum() + + sim = SimulationParameters(N=int(numtrajs), spectrum=spec) + + create_xml(filename, NssConfig(simulation=sim)) diff --git a/src/nuspacesim/apps/commands/run.py b/src/nuspacesim/apps/commands/run.py new file mode 100644 index 0000000..a0ba6a9 --- /dev/null +++ b/src/nuspacesim/apps/commands/run.py @@ -0,0 +1,142 @@ +import click + +from ...compute import compute +from ...config import MonoSpectrum, PowerSpectrum # FileSpectrum, +from ...utils.plot_function_registry import registry +from ...xml_config import config_from_xml + + +@click.command() +@click.option( + "-o", "--output", type=click.Path(exists=False), default=None, help="Output file." +) +@click.option( + "-p", + "--plot", + type=click.Choice(list(registry), case_sensitive=False), + multiple=True, + default=[], + help="Available plotting functions. Select multiple plots with multiple uses of -p", +) +@click.option( + "-ps", + "--plotsettings", + nargs=7, + type=click.Tuple([int, int, bool, str, bool, int, str]), + default=None, + help="Save plot supplied with -p with given file extension, optionally suppress pop_up", +) +@click.option( + "-pc", + "--plotconfig", + type=click.Path( + exists=True, + dir_okay=False, + readable=True, + ), + help="Read selected plotting functions and options from the specified INI file", +) +@click.option("--plotall", is_flag=True, help="Show all result plots.") +@click.option( + "-w", + "--write-stages", + is_flag=True, + help="Write intermediate values after each simulation stage.", +) +@click.option( + "-n", + "--no-result-file", + is_flag=True, + help="Do not save the results to an output file.", +) +@click.option( + "--monospectrum", + type=float, + default=None, + help="Mono Energetic Spectrum Log Energy.", +) +@click.option( + "--powerspectrum", + nargs=3, + type=click.Tuple([float, float, float]), + default=None, + help="Power Spectrum index, lower_bound, upper_bound.", +) +@click.argument( + "config_file", + default=None, + type=click.Path(exists=True, dir_okay=False, readable=True), +) +@click.argument("count", type=float, default=0.0) +def run( + config_file: str, + count: float, + monospectrum, + powerspectrum, + no_result_file: bool, + output: str, + plot: list, + plotsettings, + plotconfig: str, + plotall: bool, + write_stages: bool, +) -> None: + """Perform the full nuspacesim simulation. + + Main Simulator for nuspacesim. Given a XML configuration file, and + optionally a count of simulated nutrinos, runs nutrino simulation. + + \f + + Parameters + ---------- + config_file: str + XML configuration file for particular simulation particular. + count : int, optional + Number of thrown trajectories. Optionally override value in config_file. + spectrum_type : str, optional + Type of + output: str, optional + Name of the output file holding the simulation results. + plot: list, optional + Plot the simulation results. + plotconfig: str, optional + INI file to select plots for each module, as well as to specifiy global plot settings. + plotall: bool, optional + Plot all the available the simulation results plots. + no_result_file: bool, optional + Disables writing results to output files. + + + Examples + -------- + Command line usage of the run command may work with the following invocation. + + `nuspacesim run sample_input_file.xml 1e5 8 -o my_sim_results.fits` + """ + + # User Inputs + config = config_from_xml(config_file) + + config.simulation.N = int(config.simulation.N if count == 0.0 else count) + + if monospectrum is not None and powerspectrum is not None: + raise RuntimeError("Only one of --monospectrum or --powerspectrum may be used.") + if monospectrum is not None: + config.simulation.spectrum = MonoSpectrum(monospectrum) + if powerspectrum is not None: + config.simulation.spectrum = PowerSpectrum(*powerspectrum) + + plot = list(registry) if plotall else plot + print("plotconfig", plotconfig) + + simulation = compute( + config, + verbose=True, + plot=plot, + plot_config=plotconfig, + output_file=output, + write_stages=write_stages, + ) + if not no_result_file: + simulation.write(output, overwrite=True) diff --git a/src/nuspacesim/apps/commands/show_plot.py b/src/nuspacesim/apps/commands/show_plot.py new file mode 100644 index 0000000..e8a38f7 --- /dev/null +++ b/src/nuspacesim/apps/commands/show_plot.py @@ -0,0 +1,105 @@ +# The Clear BSD License +# +# Copyright (c) 2021 Alexander Reustle and the NuSpaceSim Team +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted (subject to the limitations in the disclaimer +# below) provided that the following conditions are met: +# +# * Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# * Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# * Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY +# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND +# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER +# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +import click + +from ... import simulation +from ...results_table import ResultsTable +from ...utils import dashboard_plot +from ...utils.plot_function_registry import registry +from ...utils.plot_wrapper import PlotWrapper + + +@click.command() +@click.option( + "-p", + "--plot", + type=click.Choice(list(registry), case_sensitive=False), + multiple=True, + default=[], + help="Available plotting functions. Select multiple plots with multiple uses of -p", +) +@click.option( + "-pc", + "--plotconfig", + type=click.Path( + exists=True, + dir_okay=False, + readable=True, + ), + help="Read selected plotting functions and options from the specified INI file", +) +@click.option("--plotall", is_flag=True, help="Show all result plots.") +@click.argument( + "simulation_file", + default=None, + type=click.Path(exists=True, dir_okay=False, readable=True), +) +def show_plot( + simulation_file: str, + plot: list, + plotconfig: str, + plotall: bool, +) -> None: + """Show predefined plots of results in simulation file. + + \f + + Parameters + ---------- + simulation_file: str + input ResultsTable fits file. + plot: list, optional + Plot the simulation results. + plotconfig: str, optional + INI file to select plots for each module, as well as to specifiy global plot settings. + plotall: bool, optional + Plot all the available the simulation results plots. + + + Examples + -------- + Command line usage of the run command may work with the following invocation. + + `nuspacesim show_plot my_sim_results.fits -p taus_overview` + """ + + sim = ResultsTable.read(simulation_file) + + plot = list(registry) if plotall else plot + plot_wrapper = PlotWrapper(to_plot=plot, plotconfig=plotconfig) + + simulation.geometry.region_geometry.show_plot(sim, plot_wrapper) + simulation.taus.taus.show_plot(sim, plot_wrapper) + simulation.eas_optical.eas.show_plot(sim, plot_wrapper) + dashboard_plot.show_plot(sim, plot_wrapper) diff --git a/src/nuspacesim/compute.py b/src/nuspacesim/compute.py index 82d4305..af9d906 100644 --- a/src/nuspacesim/compute.py +++ b/src/nuspacesim/compute.py @@ -57,6 +57,7 @@ from .simulation.geometry.region_geometry import RegionGeom from .simulation.spectra.spectra import Spectra from .simulation.taus.taus import Taus +from .utils.plot_wrapper import PlotWrapper __all__ = ["compute"] @@ -64,8 +65,9 @@ def compute( config: NssConfig, verbose: bool = False, - output_file: str = None, - to_plot: list = [], + plot: list = [], + plot_config=None, + output_file=None, write_stages=False, ) -> ResultsTable: r"""Simulate an upward going shower. @@ -105,8 +107,8 @@ def compute( Flag enabling verbose output. output_file: str, optional Name of file to write intermediate stages - to_plot: list, optional - Call the listed plotting functions as appropritate. + plot_wrapper: PlotWrapper, optional + The contained plotting specifications write_stages: bool, optional Enable writing intermediate results to the output_file. @@ -144,6 +146,12 @@ def mc_logv(mcint, mcintgeo, numEvPass, mcunc, method): eas = EAS(config) eas_radio = EASRadio(config) + plot_wrapper = PlotWrapper( + to_plot=plot, + filename=sim.output_file_basename(output_file), + plotconfig=plot_config, + ) + class StagedWriter: """Optionally write intermediate values to file""" @@ -158,23 +166,23 @@ def add_meta(self, *args, **kwargs): sim.write(output_file, overwrite=True) sw = StagedWriter() - logv(f"Running NuSpaceSim with Energy Spectrum ({config.simulation.spectrum})") logv("Computing [green] Geometries.[/]") - beta_tr, thetaArr, pathLenArr = geom(config.simulation.N, store=sw, plot=to_plot) + beta_tr, thetaArr, pathLenArr = geom( + config.simulation.N, store=sw, plot_wrapper=plot_wrapper + ) logv( f"\t[blue]Threw {config.simulation.N} neutrinos. {beta_tr.size} were valid.[/]" ) logv("Computing [green] Energy Spectra.[/]") log_e_nu, mc_spec_norm, spec_weights_sum = spec( - beta_tr.shape[0], store=sw, plot=to_plot + beta_tr.shape[0], store=sw, plot_wrapper=plot_wrapper ) - logv("Computing [green] Taus.[/]") - tauBeta, tauLorentz, tauEnergy, showerEnergy, tauExitProb = tau( - beta_tr, log_e_nu, store=sw, plot=to_plot + tauBeta, tauLorentz, _, showerEnergy, tauExitProb = tau( + beta_tr, log_e_nu, store=sw, plot_wrapper=plot_wrapper ) logv("Computing [green] Decay Altitudes.[/]") @@ -184,11 +192,7 @@ def add_meta(self, *args, **kwargs): logv("Computing [green] EAS Optical Cherenkov light.[/]") numPEs, costhetaChEff = eas( - beta_tr, - altDec, - showerEnergy, - store=sw, - plot=to_plot, + beta_tr, altDec, showerEnergy, store=sw, plot_wrapper=plot_wrapper ) logv("Computing [green] Optical Monte Carlo Integral.[/]") diff --git a/src/nuspacesim/results_table.py b/src/nuspacesim/results_table.py index 449686c..45b489f 100644 --- a/src/nuspacesim/results_table.py +++ b/src/nuspacesim/results_table.py @@ -119,6 +119,13 @@ def add_meta(self, name: str, value: Any, comment: str) -> None: self.meta[name] = (value, comment) + def output_file_basename(self, filename: Union[str, None] = None) -> str: + return ( + filename.rpartition(".")[0] + if filename + else f"nuspacesim_run_{self.meta['simTime'][0]}" + ) + def write(self, filename: Union[str, None] = None, **kwargs) -> None: r"""Write the simulation results to a file. @@ -135,34 +142,20 @@ def write(self, filename: Union[str, None] = None, **kwargs) -> None: ValueError: If the input format value is not one of fits or hdf5. """ + file_basename = self.output_file_basename(filename) if "format" not in kwargs: kwargs["format"] = "fits" if kwargs["format"] == "fits": - - filename = ( - f"nuspacesim_run_{self.meta['simTime'][0]}.fits" - if filename is None - else filename - ) - super().write(filename, **kwargs) - + super().write(file_basename + ".fits", **kwargs) elif kwargs["format"] == "hdf5": - - filename = ( - f"nuspacesim_run_{self.meta['simTime'][0]}.hdf5" - if filename is None - else filename - ) - if "path" not in kwargs: kwargs["path"] = "/" if "overwrite" not in kwargs: kwargs["overwrite"] = True kwargs["serialize_meta"] = True - - super().write(filename, **kwargs) + super().write(file_basename + ".hdf5", **kwargs) else: raise ValueError(f"File output format {format} not in {{ fits, hdf5 }}!") diff --git a/src/nuspacesim/simulation/eas_optical/__init__.py b/src/nuspacesim/simulation/eas_optical/__init__.py index bfff857..b3f9159 100644 --- a/src/nuspacesim/simulation/eas_optical/__init__.py +++ b/src/nuspacesim/simulation/eas_optical/__init__.py @@ -53,7 +53,7 @@ cphotang """ -__all__ = ["EAS", "atmospheric_models"] +__all__ = ["EAS", "atmospheric_models", "local_plots"] -from . import atmospheric_models +from . import atmospheric_models, local_plots from .eas import EAS diff --git a/src/nuspacesim/simulation/eas_optical/eas.py b/src/nuspacesim/simulation/eas_optical/eas.py index e6f5b75..0751a2e 100644 --- a/src/nuspacesim/simulation/eas_optical/eas.py +++ b/src/nuspacesim/simulation/eas_optical/eas.py @@ -36,7 +36,20 @@ from ...config import NssConfig from ...utils import decorators from .cphotang import CphotAng -from .local_plots import eas_optical_density, eas_optical_histogram +from .local_plots import ( + altdec_vs_beta, + altdec_vs_costhetacheff, + altdec_vs_numpes, + costhetacheff_hist, + costhetacheff_vs_beta, + costhetacheff_vs_numpes, + eas_optical_density_overview, + eas_optical_histogram_overview, + numpes_hist, + numpes_vs_beta, + showerenergy_vs_costhetacheff, + showerenergy_vs_numpes, +) __all__ = ["EAS", "show_plot"] @@ -53,7 +66,7 @@ def __init__(self, config: NssConfig): self.CphotAng = CphotAng(self.config.detector.altitude) @decorators.nss_result_store("altDec", "lenDec") - def altDec(self, beta, tauBeta, tauLorentz, u=None, *args, **kwargs): + def altDec(self, beta, tauBeta, tauLorentz, u=None, *_, **kwargs): """ get decay altitude """ @@ -74,9 +87,21 @@ def altDec(self, beta, tauBeta, tauLorentz, u=None, *args, **kwargs): return altDec, lenDec - @decorators.nss_result_plot(eas_optical_density, eas_optical_histogram) + @decorators.nss_result_plot( + numpes_vs_beta, + altdec_vs_numpes, + altdec_vs_beta, + showerenergy_vs_numpes, + costhetacheff_vs_beta, + altdec_vs_costhetacheff, + showerenergy_vs_costhetacheff, + numpes_hist, + costhetacheff_hist, + eas_optical_density_overview, + eas_optical_histogram_overview, + ) @decorators.nss_result_store("numPEs", "costhetaChEff") - def __call__(self, beta, altDec, showerEnergy, *args, **kwargs): + def __call__(self, beta, altDec, showerEnergy, *_, **kwargs): """ Electromagnetic Air Shower operation. """ @@ -115,8 +140,20 @@ def __call__(self, beta, altDec, showerEnergy, *args, **kwargs): return numPEs, costhetaChEff -def show_plot(sim, plot): - plotfs = (eas_optical_density, eas_optical_histogram) +def show_plot(sim, plot_wrapper): + plotfs = ( + numpes_vs_beta, + altdec_vs_numpes, + altdec_vs_beta, + showerenergy_vs_numpes, + costhetacheff_vs_beta, + altdec_vs_costhetacheff, + showerenergy_vs_costhetacheff, + numpes_hist, + costhetacheff_hist, + eas_optical_density_overview, + eas_optical_histogram_overview, + ) inputs = ("beta_rad", "altDec", "showerEnergy") outputs = ("numPEs", "costhetaChEff") - decorators.nss_result_plot_from_file(sim, inputs, outputs, plotfs, plot) + decorators.nss_result_plot_from_file(sim, inputs, outputs, plotfs, plot_wrapper) diff --git a/src/nuspacesim/simulation/eas_optical/local_plots.py b/src/nuspacesim/simulation/eas_optical/local_plots.py index b0905d3..779ba4d 100644 --- a/src/nuspacesim/simulation/eas_optical/local_plots.py +++ b/src/nuspacesim/simulation/eas_optical/local_plots.py @@ -32,84 +32,270 @@ # POSSIBILITY OF SUCH DAMAGE. import numpy as np -from matplotlib import pyplot as plt -from ...utils.plots import hist2d +from ...utils.plots import hexbin, make_labels -def eas_optical_density(inputs, results, *args, **kwargs): - r"""Plot some density plots""" +def numpes_vs_beta(inputs, results, fig, ax, *_, **kwargs): + _, betas, _, _ = inputs + numPEs, _ = results + binning_b = np.arange( + np.min(np.degrees(betas)) - 1, np.max(np.degrees(betas)) + 2, 1 + ) + im = hexbin( + ax, + x=np.degrees(betas), + y=numPEs, + gs=len(binning_b), + logx=True, + logy=True, + cmap=kwargs["cmap"], + ) + make_labels( + fig, + ax, + "$\\beta$ / $^{\\circ}$", + "#PEs", + clabel="Counts", + im=im, + logx=True, + logy=True, + ) - _, betas, altDec, showerEnergy = inputs - numPEs, costhetaChEff = results - fig, ax = plt.subplots(2, 3, figsize=(15, 8), constrained_layout=True) +def altdec_vs_numpes(inputs, results, fig, ax, *_, **kwargs): + _, _, altDec, _ = inputs + numPEs, _ = results - hist2d(fig, ax[0, 0], np.degrees(betas), numPEs, "β", "numPEs", cmap="plasma") - hist2d( + im = hexbin( + ax, + x=altDec, + y=numPEs, + gs=25, + logx=True, + logy=True, + cmap=kwargs["cmap"], + ) + make_labels( fig, - ax[1, 0], - np.degrees(betas), - costhetaChEff, - "β", - "cos(θ_chEff)", - cmap="plasma", + ax, + "$Decay_\\mathrm{Altitude}$ / km", + "#PEs", + clabel="Counts", + im=im, + logx=True, + logy=True, ) - hist2d( - fig, ax[0, 1], altDec, numPEs, "decay altitude (km)", "numPEs", cmap="plasma" + +def altdec_vs_beta(inputs, results, fig, ax, *_, **kwargs): + _, betas, altDec, _ = inputs + binning_b = np.arange( + np.min(np.degrees(betas)) - 1, np.max(np.degrees(betas)) + 2, 1 + ) + im = hexbin( + ax, + x=np.degrees(betas), + y=altDec, + gs=len(binning_b), + logx=True, + logy=True, + cmap=kwargs["cmap"], ) + make_labels( + fig, + ax, + "$\\beta$ / $^{\\circ}$", + "$Decay_\\mathrm{Altitude}$ / km", + clabel="Counts", + im=im, + logx=True, + logy=True, + ) + - hist2d( +def showerenergy_vs_numpes(inputs, results, fig, ax, *args, **kwargs): + _, _, _, showerEnergy = inputs + numPEs, _ = results + + im = hexbin( + ax, + x=showerEnergy, + y=numPEs, + gs=25, + logx=True, + logy=True, + cmap=kwargs["cmap"], + ) + make_labels( fig, - ax[1, 1], - altDec, - costhetaChEff, - "decay altitude (km)", - "cos(θ_chEff)", - cmap="plasma", + ax, + "$E_\\mathrm{shower}$ / 100 PeV", + "#PEs", + clabel="Counts", + im=im, + logx=True, + logy=True, ) - hist2d( + +def costhetacheff_vs_beta(inputs, results, fig, ax, *_, **kwargs): + _, betas, _, _ = inputs + _, costhetaChEff = results + binning_b = np.arange( + np.min(np.degrees(betas)) - 1, np.max(np.degrees(betas)) + 2, 1 + ) + im = hexbin( + ax, + x=np.degrees(betas), + y=costhetaChEff, + gs=len(binning_b), + logx=True, + logy=True, + cmap=kwargs["cmap"], + ) + make_labels( fig, - ax[0, 2], - showerEnergy, - numPEs, - "showerEnergy (100 PeV)", - "numPEs", - cmap="plasma", + ax, + "$\\beta$ / $^{\\circ}$", + "$\\cos(\\theta_\\mathrm{Cherenkov}$ / rad)", + clabel="Counts", + im=im, + logx=True, + logy=True, ) - hist2d( + +def costhetacheff_vs_numpes(inputs, results, fig, ax, *_, **kwargs): + numPEs, costhetaChEff = results + im = hexbin( + ax, + x=numPEs, + y=costhetaChEff, + gs=25, + logx=True, + logy=True, + cmap=kwargs["cmap"], + ) + make_labels( fig, - ax[1, 2], - showerEnergy, - costhetaChEff, - "showerEnergy (100 PeV)", - "cos(θ_chEff)", - cmap="plasma", + ax, + "#PEs", + "$\\cos(\\theta_\\mathrm{Cherenkov}$ / rad)", + clabel="Counts", + im=im, + logx=True, + logy=True, ) - fig.suptitle("EAS Optical Cherenkov properties.") - plt.show() +def altdec_vs_costhetacheff(inputs, results, fig, ax, *_, **kwargs): + _, _, altDec, _ = inputs + _, costhetaChEff = results + im = hexbin( + ax, + x=altDec, + y=costhetaChEff, + gs=25, + logx=True, + logy=True, + cmap=kwargs["cmap"], + ) + make_labels( + fig, + ax, + "$Decay_\\mathrm{Altitude}$ / km", + "$\\cos(\\theta_\\mathrm{Cherenkov}$ / rad)", + clabel="Counts", + im=im, + logx=True, + logy=True, + ) -def eas_optical_histogram(inputs, results, *args, **kwargs): - r"""Plot some histograms""" - # eas_self, betas, altDec, showerEnergy = inputs - numPEs, costhetaChEff = results +def showerenergy_vs_costhetacheff(inputs, results, fig, ax, *_, **kwargs): + _, _, _, showerEnergy = inputs + _, costhetaChEff = results + im = hexbin( + ax, + x=showerEnergy, + y=costhetaChEff, + gs=25, + logx=True, + logy=True, + cmap=kwargs["cmap"], + ) + make_labels( + fig, + ax, + "$E_\\mathrm{shower}$ / 100 PeV", + "$\\cos(\\theta_\\mathrm{Cherenkov}$ / rad)", + clabel="Counts", + im=im, + logx=True, + logy=True, + ) + + +def numpes_hist(inputs, results, fig, ax, *_, **kwargs): + numPEs, _ = results + ax.hist( + x=numPEs, + bins=np.linspace(min(numPEs), max(numPEs), 100), + color=kwargs["color"][0], + ) + make_labels( + fig, + ax, + "#PEs", + "Counts", + ) + - color = "salmon" - alpha = 1 +def costhetacheff_hist(inputs, results, fig, ax, *_, **kwargs): + _, costhetaChEff = results + ax.hist( + x=costhetaChEff, + bins=np.linspace(min(costhetaChEff), max(costhetaChEff), 100), + color=kwargs["color"][0], + ) + make_labels( + fig, + ax, + "$\\cos(\\theta_\\mathrm{Cherenkov}$ / rad)", + "Counts", + ) + + +def eas_optical_density_overview(inputs, results, fig, ax, *args, **kwargs): + r"""Plot some density plots""" - fig, ax = plt.subplots(2, 1, constrained_layout=True) + ax.remove() + fig.set_size_inches(15, 8) + ax = fig.subplots(nrows=2, ncols=3) + fig.suptitle("EAS Optical Cherenkov properties") - ax[0].hist(numPEs, 100, log=True, facecolor=color, alpha=alpha) - ax[0].set_xlabel("log(numPEs)") + numpes_vs_beta(inputs, results, fig, ax[0, 0], *args, **kwargs) + altdec_vs_numpes(inputs, results, fig, ax[0, 1], *args, **kwargs) + showerenergy_vs_numpes(inputs, results, fig, ax[0, 2], *args, **kwargs) + costhetacheff_vs_beta(inputs, results, fig, ax[1, 0], *args, **kwargs) + altdec_vs_costhetacheff(inputs, results, fig, ax[1, 1], *args, **kwargs) + showerenergy_vs_costhetacheff(inputs, results, fig, ax[1, 2], *args, **kwargs) - ax[1].hist(costhetaChEff, 100, log=True, facecolor=color, alpha=alpha) - ax[1].set_xlabel("log(cos(θ_chEff))") + return "eas_optical_density_overview" + +def eas_optical_histogram_overview(inputs, results, fig, ax, *args, **kwargs): + r"""Plot some histograms""" + + # fig = PlotWrapper( + # plot_kwargs, 2, 1, title="EAS Optical Cherenkov property Histograms" + # ) + + ax.remove() + ax = fig.subplots(nrows=2, ncols=1) fig.suptitle("EAS Optical Cherenkov property Histograms") - plt.show() + + numpes_hist(inputs, results, fig, ax[0], *args, **kwargs) + costhetacheff_hist(inputs, results, fig, ax[1], *args, **kwargs) + return "eas_optical_histogram_overview" diff --git a/src/nuspacesim/simulation/geometry/local_plots.py b/src/nuspacesim/simulation/geometry/local_plots.py index 68b014b..28fa2fb 100644 --- a/src/nuspacesim/simulation/geometry/local_plots.py +++ b/src/nuspacesim/simulation/geometry/local_plots.py @@ -32,17 +32,47 @@ # POSSIBILITY OF SUCH DAMAGE. import numpy as np -from matplotlib import pyplot as plt +from ...utils.plots import hexbin, make_labels -def geom_beta_tr_hist(inputs, results, *args, **kwargs): + +def geom_beta_tr_hist(inputs, results, fig, ax, *_, **kwargs): r"""Plot a histgram of beta trajectories.""" - _ = inputs betas, _, _ = results + ax.hist( + x=np.degrees(betas), + bins=np.linspace(min(np.degrees(betas)), max(np.degrees(betas)), 50), + color=kwargs["color"][0], + ) + make_labels( + fig=fig, + ax=ax, + xlabel="Earth emergence angle $\\beta$ / $^{\\circ}$", + ylabel="Counts", + ) + - plt.hist(np.degrees(betas), 50, alpha=0.75) - plt.xlabel("beta_tr (radians)") - plt.ylabel("frequency (counts)") - plt.title(f"Histogram of {betas.size} Beta Angles") - plt.show() +def path_length_to_detector(inputs, results, fig, ax, *_, **kwargs): + r"""Plot a histgram of beta trajectories.""" + + _ = inputs + betas, _, path_lens = results + binning_b = np.arange( + np.min(np.degrees(betas)) - 1, np.max(np.degrees(betas)) + 2, 1 + ) + im = hexbin( + ax, + x=np.degrees(betas), + y=path_lens, + gs=len(binning_b), + cmap=kwargs["cmap"], + ) + make_labels( + fig, + ax, + "Earth emergence angle $\\beta$ / $^{\\circ}$", + "Path length to detector / km", + clabel="Counts", + im=im, + ) diff --git a/src/nuspacesim/simulation/geometry/region_geometry.py b/src/nuspacesim/simulation/geometry/region_geometry.py index 391b096..ab3e773 100644 --- a/src/nuspacesim/simulation/geometry/region_geometry.py +++ b/src/nuspacesim/simulation/geometry/region_geometry.py @@ -34,7 +34,7 @@ import numpy as np from ...utils import decorators -from .local_plots import geom_beta_tr_hist +from .local_plots import geom_beta_tr_hist, path_length_to_detector __all__ = ["RegionGeom"] @@ -235,7 +235,7 @@ def valid_costhetaNSubV(self): def valid_costhetaTrSubV(self): return self.costhetaTrSubV[self.event_mask] - @decorators.nss_result_plot(geom_beta_tr_hist) + @decorators.nss_result_plot(geom_beta_tr_hist, path_length_to_detector) @decorators.nss_result_store("beta_rad", "theta_rad", "path_len") def __call__(self, numtrajs, *args, **kwargs): """Throw numtrajs events and return valid betas.""" @@ -287,8 +287,8 @@ def mcintegral( return mcintegral, mcintegralgeoonly, numEvPass, mcintegraluncert -def show_plot(sim, plot): - plotfs = tuple([geom_beta_tr_hist]) +def show_plot(sim, plot_wrapper): + plotfs = tuple([geom_beta_tr_hist, path_length_to_detector]) inputs = tuple([0]) outputs = ("beta_rad", "theta_rad", "path_len") - decorators.nss_result_plot_from_file(sim, inputs, outputs, plotfs, plot) + decorators.nss_result_plot_from_file(sim, inputs, outputs, plotfs, plot_wrapper) diff --git a/src/nuspacesim/simulation/spectra/local_plots.py b/src/nuspacesim/simulation/spectra/local_plots.py index 39896bd..7539959 100644 --- a/src/nuspacesim/simulation/spectra/local_plots.py +++ b/src/nuspacesim/simulation/spectra/local_plots.py @@ -31,24 +31,25 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -from matplotlib import pyplot as plt +from ...utils.plots import make_labels -def spectra_histogram(inputs, results, *args, **kwargs): +def spectra_histogram(inputs, results, fig, ax, *_, **kwargs): r"""Plot some histograms""" N, spectrum = inputs log_e_nu = results - - color = "g" - fig = plt.figure(figsize=(8, 7), constrained_layout=True) - ax = fig.add_subplot(211) - ax.hist(log_e_nu, 100, log=False, facecolor=color) - ax.set_xlabel(f"log(E_nu) of {N} events") - - ax = fig.add_subplot(212) - ax.hist(log_e_nu, 100, log=True, facecolor=color) - ax.set_xlabel(f"log(E_nu) of {N} events") - - fig.suptitle(f"Energy Spectra Histogram, Log(E_nu)\n {spectrum}") - plt.show() + ax.hist( + x=log_e_nu, + bins=100, + color=kwargs["color"][0], + label=f"Energy spectrum histogram with {N} events\n {spectrum}", + ) + make_labels( + fig, + ax, + "$E_{\\nu_\\tau}$ / $\\log_\\mathrm{10}\\left(\\frac{E}{\\mathrm{eV}}\\right)$", + "Counts", + logy_scale=True, + ) + ax.legend(loc="best") diff --git a/src/nuspacesim/simulation/taus/local_plots.py b/src/nuspacesim/simulation/taus/local_plots.py index 5ac9f7a..60fb4f9 100644 --- a/src/nuspacesim/simulation/taus/local_plots.py +++ b/src/nuspacesim/simulation/taus/local_plots.py @@ -32,208 +32,261 @@ # POSSIBILITY OF SUCH DAMAGE. import numpy as np -from matplotlib import pyplot as plt -from mpl_toolkits.axes_grid1 import make_axes_locatable +from ...utils.plots import get_profile, hexbin, make_labels -def get_profile(x, y, nbins, useStd=True, *args, **kwargs): - if sum(np.isnan(y)) > 0: - # print "Array contains NaN, removing them for profile plot" - x = x[~np.isnan(y)] - y = y[~np.isnan(y)] +r"""Single plots""" - n, _ = np.histogram(x, bins=nbins, **kwargs) - sy, _ = np.histogram(x, bins=nbins, weights=y, **kwargs) - sy2, _ = np.histogram(x, bins=nbins, weights=y * y, **kwargs) - mean = sy / n - std = np.sqrt(sy2 / n - mean * mean) - if not useStd: - std /= np.sqrt(n) - bincenter = (_[1:] + _[:-1]) / 2 - binwidth = bincenter - _[1:] - return bincenter, mean, std, binwidth - - -def taus_density_beta(inputs, results, *args, **kwargs): - r"""Plot some density plots""" - - _, betas, log_e_nu = inputs - tauBeta, tauLorentz, tauEnergy, showerEnergy, tauExitProb = results - - fig, ax = plt.subplots(2, 2, figsize=(10, 8), sharex=True, constrained_layout=True) - - im = None - - def hist2d(ax, x, y, xlab, ylab): - nonlocal im - _, _, _, im = ax.hist2d( - x=np.degrees(x), y=np.log10(y), bins=(50, 50), cmin=1, cmap="jet" - ) - ax.set_xlabel(xlab) - ax.set_ylabel(ylab) - ax.set_title(f"{xlab} vs {ylab}") - - hist2d(ax[0, 0], betas, tauBeta, "β", "log(τ_β)") - hist2d(ax[0, 1], betas, tauLorentz, "β", "log(τ_Lorentz)") - hist2d(ax[1, 0], betas, showerEnergy, "β", "log(showerEnergy)") - hist2d(ax[1, 1], betas, tauExitProb, "β", "log(Exit Probability (τ))") - - fig.colorbar(im, ax=ax, label="Counts", format="%.0e") - - fig.suptitle("Tau interaction properties vs. β_tr Angles") - plt.show() - - -def taus_histogram(inputs, results, *args, **kwargs): - r"""Plot some histograms""" - - _, _, _ = inputs - tauBeta, tauLorentz, tauEnergy, showerEnergy, tauExitProb = results - - color = "c" - alpha = 1 - - fig, ax = plt.subplots(2, 2, constrained_layout=True) - - ax[0, 0].hist(tauBeta, 100, log=True, facecolor=color, alpha=alpha) - ax[0, 0].set_xlabel("log(τ_β)") - ax[0, 1].hist(tauLorentz, 100, log=True, facecolor=color, alpha=alpha) - ax[0, 1].set_xlabel("log(τ_Lorentz)") - ax[1, 0].hist(showerEnergy, 100, log=True, facecolor=color, alpha=alpha) - ax[1, 0].set_xlabel("log(showerEnergy)") - ax[1, 1].hist(tauExitProb, 100, log=True, facecolor=color, alpha=alpha) - ax[1, 1].set_xlabel("log(PExit(τ))") - - fig.suptitle("Tau interaction property Histograms") - plt.show() - - -def taus_pexit(inputs, results, *args, **kwargs): - _, betas, _ = inputs - _, _, _, _, tauExitProb = results - - color = "c" - alpha = 0.1 / np.log10(betas.size) - - fig, ax = plt.subplots(1, 2, constrained_layout=True) - - ax[0].scatter( - x=np.degrees(betas), - y=np.log10(tauExitProb), - s=1, - c=color, - marker=".", - alpha=alpha, - ) - ax[0].set_xlabel("β") - ax[0].set_ylabel("log(PExit(τ))") - ax[0].set_title("β vs Exit Probability.") - - ax[1].hist(tauExitProb, 100, log=True, facecolor=color) - ax[1].set_ylabel("log(frequency)") - ax[1].set_xlabel("PExit(τ)") - - fig.suptitle("Tau Pexit") - plt.show() - - -def taus_overview(inputs, results, *args, **kwargs): - r"""Overview plot for taus""" - - _, betas, log_e_nu = inputs - _, tauLorentz, tauEnergy, showerEnergy, tauExitProb = results +def energy_hists(inputs, results, fig, ax, *_, **kwargs): + _, _, log_e_nu = inputs + _, _, tauEnergy, showerEnergy, _ = results log_e_nu = log_e_nu + 9 tauEnergy = tauEnergy * 1e9 showerEnergy = showerEnergy * 1e9 * 1e8 - - fig, ax = plt.subplots(2, 2, figsize=(10, 8)) - binning_b = np.arange( - np.min(np.degrees(betas)) - 1, np.max(np.degrees(betas)) + 2, 1 - ) binning_e = np.arange( np.round(np.min(np.log10(showerEnergy)), 1) - 0.1, np.round(np.max(log_e_nu), 1) + 0.1, 0.1, ) - ax[0, 0].hist( - x=log_e_nu, bins=binning_e, color="C0", alpha=0.6, label=r"$E_{\nu_\tau}$" + ax.hist( + x=log_e_nu, + bins=binning_e, + color=kwargs["color"][0], + alpha=0.6, + label="$E_{\\nu_\\tau}$", ) - ax[0, 0].hist( + ax.hist( x=np.log10(tauEnergy), bins=binning_e, - color="C1", + color=kwargs["color"][1], alpha=0.6, - label=r"$E_\tau$", + label="$E_\\tau$", ) - ax[0, 0].hist( + ax.hist( x=np.log10(showerEnergy), bins=binning_e, - color="C2", + color=kwargs["color"][2], alpha=0.6, - label=r"$E_\mathrm{shower}$", + label="$E_\\mathrm{shower}$", + ) + make_labels( + fig, + ax, + "Energy / $\\log_\\mathrm{10}\\left(\\frac{E}{\\mathrm{eV}}\\right)$", + "Counts", + logy_scale=True, ) - ax[0, 0].set_yscale("log") - ax[0, 0].legend(loc="upper left") - ax[0, 0].set_xlabel( - r"Energy / $\log_\mathrm{10}\left(\frac{E}{\mathrm{eV}}\right)$" + ax.legend(loc="best") + + +def beta_hist(inputs, _, fig, ax, *__, **kwargs): + _, betas, _ = inputs + binning_b = np.arange( + np.min(np.degrees(betas)) - 1, np.max(np.degrees(betas)) + 2, 1 + ) + ax.hist( + x=np.degrees(betas), + bins=binning_b, + color=kwargs["color"][0], + ) + make_labels( + fig, + ax, + "Earth emergence angle $\\beta$ / $^{\\circ}$", + "Counts", ) - ax[0, 0].set_ylabel(r"Counts") - n, bins, _ = ax[0, 1].hist(x=np.degrees(betas), bins=binning_b) - ax[0, 1].set_xlabel(r"Earth emergence angle $\beta$ / $^{\circ}$") - ax[0, 1].set_ylabel("Counts") - ax[0, 1].set_yscale("log") - ax[0, 1].set_xlim(min(binning_b), max(binning_b)) - binning_y = np.logspace( - np.log10(np.min(tauLorentz)) - 0.2, np.log10(np.max(tauLorentz)) + 0.2, 25 +def tau_beta_hist(_, results, fig, ax, *__, **kwargs): + tauBeta, *_ = results + ax.hist( + x=tauBeta, + bins=np.linspace(min(tauBeta), max(tauBeta), 100), + color=kwargs["color"][0], + ) + make_labels( + fig, + ax, + "$\\tau_\\beta$", + "Counts", ) - _, _, _, im = ax[1, 0].hist2d( - np.degrees(betas), tauLorentz, bins=[binning_b, binning_y], cmin=1 + + +def tau_lorentz_hist(_, results, fig, ax, *__, **kwargs): + _, tauLorentz, *_ = results + ax.hist( + x=tauLorentz, + bins=np.linspace(min(tauLorentz), max(tauLorentz), 100), + color=kwargs["color"][0], + ) + make_labels( + fig, + ax, + "$\\tau_\\mathrm{Lorentz}$", + "Counts", + ) + + +def tau_exit_prob_hist(_, results, fig, ax, *__, **kwargs): + *_, tauExitProb = results + ax.hist( + x=tauExitProb, + bins=np.linspace(min(tauExitProb), max(tauExitProb), 100), + color=kwargs["color"][0], + ) + make_labels( + fig, + ax, + "$P_\\mathrm{Exit}(\\tau)$", + "Counts", + ) + + +def tau_beta_hex(inputs, results, fig, ax, *args, **kwargs): + _, betas, log_e_nu = inputs + tauBeta, tauLorentz, tauEnergy, showerEnergy, tauExitProb = results + binning_b = np.arange( + np.min(np.degrees(betas)) - 1, np.max(np.degrees(betas)) + 2, 1 + ) + im = hexbin( + ax, + x=np.degrees(betas), + y=tauBeta, + gs=len(binning_b), + logy_scale=True, + cmap=kwargs["cmap"], + ) + make_labels( + fig, + ax, + "Earth emergence angle $\\beta$ / $^{\\circ}$", + "$\\tau_\\beta$", + clabel="Counts", + im=im, + logy_scale=True, + ) + + +def tau_lorentz_hex(inputs, results, fig, ax, *args, **kwargs): + _, betas, log_e_nu = inputs + tauBeta, tauLorentz, tauEnergy, showerEnergy, tauExitProb = results + binning_b = np.arange( + np.min(np.degrees(betas)) - 1, np.max(np.degrees(betas)) + 2, 1 + ) + im = hexbin( + ax, + x=np.degrees(betas), + y=tauLorentz, + gs=len(binning_b), + logy_scale=True, + cmap=kwargs["cmap"], ) bincenter, mean, std, binwidth = get_profile( - np.degrees(betas), tauLorentz, 10, useStd=True, **kwargs + np.degrees(betas), tauLorentz, 10, useStd=True ) - ax[1, 0].errorbar( + ax.errorbar( bincenter, mean, - yerr=std, xerr=binwidth, - color="C1", + yerr=std, + color=kwargs["color"][1], fmt=".", - lw=2, - zorder=5, label="Profile", ) - divider = make_axes_locatable(ax[1, 0]) - cax = divider.append_axes("right", size="5%", pad=0.0) - cbar = fig.colorbar(im, cax=cax) - cbar.set_label("Counts") - ax[1, 0].set_yscale("log") - ax[1, 0].set_xlim(min(binning_b), max(binning_b)) - ax[1, 0].legend(loc="upper center") - ax[1, 0].set_xlabel(r"Earth emergence angle $\beta$ / $^{\circ}$") - ax[1, 0].set_ylabel(r"$\tau_\mathrm{Lorentz}$") - - binning_y = np.logspace( - np.log10(np.min(tauExitProb)) - 0.2, np.log10(np.max(tauExitProb)) + 0.2, 25 - ) - _, _, _, im = ax[1, 1].hist2d( - np.degrees(betas), tauExitProb, bins=[binning_b, binning_y], cmin=1 - ) - divider = make_axes_locatable(ax[1, 1]) - cax = divider.append_axes("right", size="5%", pad=0.0) - cbar = fig.colorbar(im, cax=cax) - cbar.set_label("Counts") - ax[1, 1].set_yscale("log") - ax[1, 1].set_xlim(min(binning_b), max(binning_b)) - ax[1, 1].set_xlabel(r"Earth emergence angle $\beta$ / $^{\circ}$") - ax[1, 1].set_ylabel(r"$P_\mathrm{Exit}(\tau)$") + ax.legend(loc="best") + make_labels( + fig, + ax, + "Earth emergence angle $\\beta$ / $^{\\circ}$", + "$\\tau_\\mathrm{Lorentz}$", + clabel="Counts", + im=im, + logy_scale=True, + ) + + +def tau_exit_prob_hex(inputs, results, fig, ax, *args, **kwargs): + _, betas, log_e_nu = inputs + *_, tauExitProb = results + binning_b = np.arange( + np.min(np.degrees(betas)) - 1, np.max(np.degrees(betas)) + 2, 1 + ) + im = hexbin( + ax, + x=np.degrees(betas), + y=tauExitProb, + gs=len(binning_b), + logy_scale=True, + cmap=kwargs["cmap"], + ) + make_labels( + fig, + ax, + "Earth emergence angle $\\beta$ / $^{\\circ}$", + "$P_\\mathrm{Exit}(\\tau)$", + clabel="Counts", + im=im, + logy_scale=True, + ) + + +r"""Multiplot collections""" + + +def taus_density_beta_overview(inputs, results, fig, ax, *args, **kwargs): + r"""Plot some density plots""" + + ax.remove() + fig.set_size_inches(15, 4) + ax = fig.subplots(nrows=1, ncols=3) + fig.suptitle("Tau interaction properties vs. Earth emergence angles $\\beta$") + tau_beta_hex(inputs, results, fig, ax[0], *args, **kwargs) + tau_lorentz_hex(inputs, results, fig, ax[1], *args, **kwargs) + tau_exit_prob_hex(inputs, results, fig, ax[2], *args, **kwargs) + + return "taus_density_beta" + + +def taus_histograms_overview(inputs, results, fig, ax, *args, **kwargs): + r"""Plot some histograms""" + + ax.remove() + fig.set_size_inches(15, 4) + ax = fig.subplots(nrows=1, ncols=3) + fig.suptitle("$\\tau$ interaction property Histograms") + + tau_beta_hist(inputs, results, fig, ax[0], *args, **kwargs) + tau_lorentz_hist(inputs, results, fig, ax[1], *args, **kwargs) + tau_exit_prob_hist(inputs, results, fig, ax[2], *args, **kwargs) + + return "taus_histogram" + + +def taus_pexit_overview(inputs, results, fig, ax, *args, **kwargs): + + ax.remove() + fig.set_size_inches(10, 4) + ax = fig.subplots(nrows=1, ncols=2) + fig.suptitle("$\\tau$ exit probability") + tau_exit_prob_hex(inputs, results, fig, ax[0], *args, **kwargs) + tau_exit_prob_hist(inputs, results, fig, ax[1], *args, **kwargs) + + return "taus_pexit" + + +def taus_overview(inputs, results, fig, ax, *args, **kwargs): + r"""Overview plot for taus""" + + ax.remove() + fig.set_size_inches(10, 8) + ax = fig.subplots(nrows=2, ncols=2) fig.suptitle("Overview of Tau interaction properties") - fig.tight_layout() - plt.show() - # if "pop_up" in kwargs and kwargs.get("pop_up") is True: - # plt.show() - # fig.savefig("taus_overview", format="pdf") + + energy_hists(inputs, results, fig, ax[0, 0], *args, **kwargs) + beta_hist(inputs, results, fig, ax[0, 1], *args, **kwargs) + tau_lorentz_hex(inputs, results, fig, ax[1, 0], *args, **kwargs) + tau_exit_prob_hex(inputs, results, fig, ax[1, 1], *args, **kwargs) + + return "taus_overview" diff --git a/src/nuspacesim/simulation/taus/taus.py b/src/nuspacesim/simulation/taus/taus.py index 7e68bcd..16d7cab 100644 --- a/src/nuspacesim/simulation/taus/taus.py +++ b/src/nuspacesim/simulation/taus/taus.py @@ -40,7 +40,20 @@ from ...utils import decorators from ...utils.cdf import grid_cdf_sampler from ...utils.grid import NssGrid -from .local_plots import taus_density_beta, taus_histogram, taus_overview, taus_pexit +from .local_plots import ( + beta_hist, + energy_hists, + tau_beta_hex, + tau_beta_hist, + tau_exit_prob_hex, + tau_exit_prob_hist, + tau_lorentz_hex, + tau_lorentz_hist, + taus_density_beta_overview, + taus_histograms_overview, + taus_overview, + taus_pexit_overview, +) try: from importlib.resources import as_file, files @@ -143,7 +156,18 @@ def tau_energy(self, betas, log_e_nu, u=None): return E_tau * (10**log_e_nu) @decorators.nss_result_plot( - taus_density_beta, taus_histogram, taus_pexit, taus_overview + energy_hists, + beta_hist, + tau_beta_hist, + tau_lorentz_hist, + tau_exit_prob_hist, + tau_beta_hex, + tau_lorentz_hex, + tau_exit_prob_hex, + taus_density_beta_overview, + taus_histograms_overview, + taus_pexit_overview, + taus_overview, ) @decorators.nss_result_store( "tauBeta", "tauLorentz", "tauEnergy", "showerEnergy", "tauExitProb" @@ -180,8 +204,21 @@ def __call__(self, betas, log_e_nu, *args, **kwargs): return tauBeta, tauLorentz, tauEnergy, showerEnergy, tauExitProb -def show_plot(sim, plot): +def show_plot(sim, plot_wrapper): inputs = ("beta_rad", "log_e_nu") outputs = ("tauBeta", "tauLorentz", "tauEnergy", "showerEnergy", "tauExitProb") - plotfs = (taus_density_beta, taus_histogram, taus_pexit, taus_overview) - decorators.nss_result_plot_from_file(sim, inputs, outputs, plotfs, plot) + plotfs = ( + energy_hists, + beta_hist, + tau_beta_hist, + tau_lorentz_hist, + tau_exit_prob_hist, + tau_beta_hex, + tau_lorentz_hex, + tau_exit_prob_hex, + taus_density_beta_overview, + taus_histograms_overview, + taus_pexit_overview, + taus_overview, + ) + decorators.nss_result_plot_from_file(sim, inputs, outputs, plotfs, plot_wrapper) diff --git a/src/nuspacesim/utils/__init__.py b/src/nuspacesim/utils/__init__.py index d15df7a..48a7b5a 100644 --- a/src/nuspacesim/utils/__init__.py +++ b/src/nuspacesim/utils/__init__.py @@ -43,9 +43,18 @@ interp misc plots - + plot_wrapper """ -__all__ = ["cdf", "decorators", "grid", "interp", "misc", "plots"] +__all__ = [ + "cdf", + "dashboard_plot", + "decorators", + "grid", + "interp", + "misc", + "plots", + "plot_wrapper", +] -from . import cdf, decorators, grid, interp, misc, plots +from . import cdf, dashboard_plot, decorators, grid, interp, misc, plot_wrapper, plots diff --git a/src/nuspacesim/utils/dashboard_plot.py b/src/nuspacesim/utils/dashboard_plot.py new file mode 100644 index 0000000..b5b19ee --- /dev/null +++ b/src/nuspacesim/utils/dashboard_plot.py @@ -0,0 +1,105 @@ +# The Clear BSD License +# Copyright (c) 2021 Alexander Reustle and the NuSpaceSim Team +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted (subject to the limitations in the disclaimer +# below) provided that the following conditions are met: +# +# * Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# * Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# * Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY +# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND +# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER +# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +import numpy as np +from matplotlib import pyplot as plt + +from ..simulation import eas_optical, taus +from . import decorators +from .plots import hist2d + +__all__ = ["dashboard", "show_plot"] + + +def dashboard(sim, _, fig, ax, *args, **kwargs): + """Full dashboard of plots""" + + ax.remove() + fig.set_size_inches(15, 8) + ax = fig.subplots(nrows=3, ncols=4) + fig.suptitle("Nuspacesim Results Dashboard") + + tau_input = None, sim["beta_rad"], sim["log_e_nu"] + tau_results = ( + sim["tauBeta"], + sim["tauLorentz"], + sim["tauEnergy"], + sim["showerEnergy"], + sim["tauExitProb"], + ) + + eas_input = None, sim["beta_rad"], sim["altDec"], sim["showerEnergy"] + eas_results = sim["numPEs"], sim["costhetaChEff"] + + taus.local_plots.energy_hists( + tau_input, tau_results, fig, ax[0, 0], *args, **kwargs + ) + taus.local_plots.tau_exit_prob_hist( + tau_input, tau_results, fig, ax[1, 0], *args, **kwargs + ) + taus.local_plots.tau_lorentz_hex( + tau_input, tau_results, fig, ax[2, 0], *args, **kwargs + ) + + taus.local_plots.beta_hist(tau_input, tau_results, fig, ax[0, 1], *args, **kwargs) + taus.local_plots.tau_exit_prob_hex( + tau_input, tau_results, fig, ax[1, 1], *args, **kwargs + ) + eas_optical.local_plots.altdec_vs_beta( + eas_input, eas_results, fig, ax[2, 1], *args, **kwargs + ) + + eas_optical.local_plots.numpes_hist( + eas_input, eas_results, fig, ax[0, 2], *args, **kwargs + ) + eas_optical.local_plots.numpes_vs_beta( + eas_input, eas_results, fig, ax[1, 2], *args, **kwargs + ) + eas_optical.local_plots.altdec_vs_numpes( + eas_input, eas_results, fig, ax[2, 2], *args, **kwargs + ) + + # eas_input_density(sim, fig, ax[0, 3]) + eas_optical.local_plots.costhetacheff_hist( + eas_input, eas_results, fig, ax[0, 3], *args, **kwargs + ) + eas_optical.local_plots.costhetacheff_vs_numpes( + eas_input, eas_results, fig, ax[1, 3], *args, **kwargs + ) + eas_optical.local_plots.costhetacheff_vs_beta( + eas_input, eas_results, fig, ax[2, 3], *args, **kwargs + ) + + +@decorators.ensure_plot_registry(dashboard) +def show_plot(sim, plot_wrapper): + plot_wrapper(sim, None, (dashboard,)) diff --git a/src/nuspacesim/utils/decorators.py b/src/nuspacesim/utils/decorators.py index 20ca40d..bfd6a53 100644 --- a/src/nuspacesim/utils/decorators.py +++ b/src/nuspacesim/utils/decorators.py @@ -47,7 +47,6 @@ __all__ = ["nss_result_store", "nss_result_store_scalar", "nss_result_plot"] from functools import wraps -from typing import Callable, Iterable, Union def nss_result_store(*names): @@ -193,24 +192,9 @@ def decorator_plot(func): registry.add(plotname) @wraps(func) - def wrapper_f( - *args, plot: Union[None, str, Iterable, Callable] = None, **kwargs - ): + def wrapper_f(*args, plot_wrapper, **kwargs): values = func(*args, **kwargs) - if isinstance(plot, str): - for plotf in plot_fs: - if plotf.__name__ == plot: - plotf(args, values) - elif callable(plot): - plot(args, values) - elif isinstance(plot, Iterable): - if all(isinstance(p, str) for p in plot): - for plotf in plot_fs: - if plotf.__name__ in plot: - plotf(args, values) - elif all(callable(p) for p in plot): - for plotf in plot: - plotf(args, values) + plot_wrapper(args, values, plot_fs, **kwargs) return values return wrapper_f @@ -218,16 +202,15 @@ def wrapper_f( return decorator_plot -def nss_result_plot_from_file(sim, inputs, outputs, plotfs, plot): - +def nss_result_plot_from_file(sim, inputs, outputs, plotfs, plot_wrapper): f_input = tuple() if inputs is None else tuple(sim[i] for i in inputs) results = tuple() if outputs is None else tuple(sim[o] for o in outputs) @nss_result_plot(*plotfs) - def f(*args, **kwargs): + def f(*_, **__): return results - f(None, *f_input, plot=plot) + f(None, *f_input, plot_wrapper=plot_wrapper) def ensure_plot_registry(*plot_fs): diff --git a/src/nuspacesim/utils/plot_wrapper.py b/src/nuspacesim/utils/plot_wrapper.py new file mode 100644 index 0000000..9eed96b --- /dev/null +++ b/src/nuspacesim/utils/plot_wrapper.py @@ -0,0 +1,137 @@ +import configparser +from os.path import exists + +from matplotlib import pyplot as plt + + +def read_plot_config(configfile): + if configfile is None: + return [], {} + if not exists(configfile): + return [], {} + plot_list = [] + cfg = configparser.ConfigParser() + cfg.read(configfile) + plot_kwargs = { + "title": cfg["General"]["title"], + "rows": cfg["General"].getint("rows"), + "columns": cfg["General"].getint("columns"), + "figsize": eval(cfg["General"]["figsize"]), + "save_as": cfg["General"]["save_as"], + "pop_up": cfg["General"].getboolean("pop_up"), + "save_to_file": cfg["General"].getboolean("save_to_file"), + "default_color": cfg["General"].getint("default_color"), + "default_colormap": cfg["General"].get("default_colormap"), + # "output_path": cfg["General"]["output_path"], + } + for sec in cfg.sections()[1:]: + for key in cfg[sec]: + try: + if cfg[sec].getboolean(key): + plot_list.append(key) + except Exception as e: + print(e, "Config file contains non-valid option") + return plot_list, plot_kwargs + + +class PlotWrapper: + """ + The PlotWrapper class produces figures that are uniformly formatted for nuspacesim as set up in sample_plot_config.ini + """ + + def __init__( + self, + to_plot=[], + rows=None, + cols=None, + figsize=None, + title=None, + save_as=None, + pop_up=None, + save_to_file=None, + default_color=None, + default_colormap=None, + filename=None, + output_path=None, + plotconfig=None, + ): + """ + initialize figure + + rows = number of rows of plots + cols = number of cols of plots + default is 1 for single plot, but can be changed to add subplots for making a multiplot + """ + + cfg_list, cfg_args = read_plot_config(plotconfig) + + cfg_args.setdefault("rows", 1) + cfg_args.setdefault("cols", 1) + cfg_args.setdefault("figsize", (8, 7)) + cfg_args.setdefault("title", "nuspacesim_run") + cfg_args.setdefault("save_as", "pdf") + cfg_args.setdefault("pop_up", True) + cfg_args.setdefault("save_to_file", False) + cfg_args.setdefault("default_color", 0) + cfg_args.setdefault("default_colormap", "viridis") + cfg_args.setdefault("filename", "NuSpaceSim") + cfg_args.setdefault("output_path", ".") + + self.to_plot = list(set(list(to_plot) + cfg_list)) + self.rows = rows if rows else cfg_args["rows"] + self.cols = cols if cols else cfg_args["cols"] + self.figsize = figsize if figsize else cfg_args["figsize"] + self.title = title if title else cfg_args["title"] + self.save_as = save_as if save_as else cfg_args["save_as"] + self.pop_up = pop_up if pop_up is not None else cfg_args["pop_up"] + self.save_to_file = ( + save_to_file if save_to_file is not None else cfg_args["save_to_file"] + ) + self.filename = filename if filename else cfg_args["filename"] + self.output_path = output_path if output_path else cfg_args["output_path"] + self.default_colormap = ( + default_colormap if default_colormap else cfg_args["default_colormap"] + ) + self.default_color = [ + f"C{i+(default_color if default_color else cfg_args['default_color'])}" + for i in range(3) + ] + + def artist_params(self): + + return { + "rows": self.rows, + "cols": self.cols, + "figsize": self.figsize, + "title": self.title, + "cmap": self.default_colormap, + "color": self.default_color, + } + + def init_fig(self, append_title=None): + # initialize figure as subplots + fig, ax = plt.subplots(nrows=self.rows, ncols=self.cols, figsize=self.figsize) + fig.suptitle(f"{self.title}\n{append_title}") + fig.tight_layout() + return fig, ax + + def __call__(self, args, values, plot_fs, **kwargs): + def do_plot(p): + fig, ax = self.init_fig(p.__name__) + title = p(args, values, fig, ax, **self.artist_params(), **kwargs) + self.close(fig, title if title else p.__name__) + + # plot_fs_prs = map(lambda p: (p.__name__, p), plot_fs) + to_plot_now = filter(lambda p: p.__name__ in self.to_plot, plot_fs) + for p in to_plot_now: + do_plot(p) + + def close(self, fig, plot_name): + if self.save_to_file: + fig.savefig( + fname=f"{self.filename}_{plot_name}.{self.save_as}", + bbox_inches="tight", + ) + if self.pop_up: + plt.show() + plt.close(fig) diff --git a/src/nuspacesim/utils/plots.py b/src/nuspacesim/utils/plots.py index e94a0a5..d5f4610 100644 --- a/src/nuspacesim/utils/plots.py +++ b/src/nuspacesim/utils/plots.py @@ -31,11 +31,88 @@ # POSSIBILITY OF SUCH DAMAGE. import numpy as np -from matplotlib import pyplot as plt -from . import decorators +__all__ = ["make_labels", "get_profile", "hexbin", "hist2d"] -__all__ = ["dashboard", "energy_histograms", "show_plot"] + +def make_labels( + fig, + ax, + xlabel, + ylabel, + clabel=None, + im=None, + logx=False, + logy=False, + logx_scale=False, + logy_scale=False, +): + + xl = "$\\log_{10}$" + f"({xlabel})" if logx else xlabel + yl = "$\\log_{10}$" + f"({ylabel})" if logy else ylabel + xs = "log" if logx_scale else "linear" + ys = "log" if logy_scale else "linear" + + ax.set_xlabel(xl) + ax.set_ylabel(yl) + ax.set_xscale(xs) + ax.set_yscale(ys) + if clabel is not None: + cbar = fig.colorbar(im, ax=ax, pad=0.0) + cbar.set_label(clabel) + + +def get_profile(x, y, nbins, useStd=True): + + if sum(np.isnan(y)) > 0: + x = x[~np.isnan(y)] + y = y[~np.isnan(y)] + n, _ = np.histogram(x, bins=nbins) + sy, _ = np.histogram(x, bins=nbins, weights=y) + sy2, _ = np.histogram(x, bins=nbins, weights=y * y) + mean = sy / n + std = np.sqrt(sy2 / n - mean * mean) + if not useStd: + std /= np.sqrt(n) + bincenter = (_[1:] + _[:-1]) / 2 + binwidth = _[1:] - bincenter + + return bincenter, mean, std, binwidth + + +def hexbin( + ax, + x, + y, + gs=25, + logx=False, + logy=False, + logx_scale=False, + logy_scale=False, + cmap=None, +): + + xf = np.log10 if logx else lambda q: q + yf = np.log10 if logy else lambda q: q + + xs = "log" if logx_scale else "linear" + ys = "log" if logy_scale else "linear" + + xmask = x > 0 if logx else np.full(x.shape, True) + ymask = y > 0 if logy else np.full(y.shape, True) + m = xmask & ymask + + im = ax.hexbin( + x=xf(x[m]), + y=yf(y[m]), + gridsize=gs, + mincnt=1, + xscale=xs, + yscale=ys, + cmap=cmap, + edgecolors="none", + ) + return im def hist2d(fig, ax, x, y, xlab, ylab, cmap="jet", logx=True, logy=True): @@ -58,240 +135,3 @@ def hist2d(fig, ax, x, y, xlab, ylab, cmap="jet", logx=True, logy=True): ax.set_title(f"{yl} vs {xl}") cbar = fig.colorbar(im, ax=ax, pad=0.0) cbar.set_label("Counts") - - -def dashboard(sim): - """Full dashboard of plots""" - - fig, ax = plt.subplots(3, 4, figsize=(14, 8), constrained_layout=True) - - energy_histograms(sim, fig, ax[0, 0]) - tau_pexit_hist(sim, fig, ax[1, 0]) - tau_lorentz(sim, fig, ax[2, 0]) - - betas_histogram(sim, fig, ax[0, 1]) - tau_pexit_density(sim, fig, ax[1, 1]) - decay_altitude(sim, fig, ax[2, 1]) - - # decay_altitude_hist(sim, fig, ax[0, 2]) - num_photo_electrons_hist(sim, fig, ax[0, 2]) - num_photo_electrons_density(sim, fig, ax[1, 2]) - num_photo_electrons_altitude(sim, fig, ax[2, 2]) - - # eas_input_density(sim, fig, ax[0, 3]) - cherenkov_angle_hist(sim, fig, ax[0, 3]) - eas_results_density(sim, fig, ax[1, 3]) - cherenkov_angle(sim, fig, ax[2, 3]) - - # tau_betas(sim, fig, ax[1, 2]) - - fig.suptitle("Nuspacesim Results Dashboard", size="x-large") - - plt.show() - - -def energy_histograms(sim, fig, ax=None): - - energy_bins = np.arange( - np.round(np.min(np.log10(sim["showerEnergy"]) + 17), 1) - 0.1, - np.round(np.max(sim["log_e_nu"] + 9), 1) + 0.1, - 0.1, - ) - ax.hist( - x=sim["log_e_nu"] + 9, - bins=energy_bins, - color="C0", - alpha=0.6, - label=r"$E_{\nu_\tau}$", - edgecolor="k", - ) - ax.hist( - x=np.log10(sim["tauEnergy"]) + 9, - bins=energy_bins, - color="C1", - alpha=0.6, - label=r"$E_\tau$", - ) - ax.hist( - x=np.log10(sim["showerEnergy"]) + 17, - bins=energy_bins, - color="C2", - alpha=0.6, - label=r"$E_\mathrm{shower}$", - ) - ax.set_yscale("log") - ax.legend(loc="upper left") - ax.set_xlabel(r"Energy / $\log_\mathrm{10}\left(\frac{E}{\mathrm{eV}}\right)$") - ax.set_ylabel(r"Counts") - - -def betas_histogram(sim, fig, ax): - - beta_bins = np.arange( - np.min(np.degrees(sim["beta_rad"])) - 1, - np.max(np.degrees(sim["beta_rad"])) + 2, - 1, - ) - - ax.hist(x=np.degrees(sim["beta_rad"]), bins=beta_bins) - ax.set_xlabel(r"Earth emergence angle $\beta$ / $^{\circ}$") - ax.set_ylabel("Counts") - ax.set_yscale("log") - ax.set_xlim(min(beta_bins), max(beta_bins)) - - -def tau_lorentz(sim, fig, ax): - hist2d( - fig, - ax, - np.degrees(sim["beta_rad"]), - sim["tauLorentz"], - r"$\beta$", - r"$τ_\mathrm{Lorentz}$", - logx=False, - logy=True, - ) - - -def tau_pexit_hist(sim, _, ax): - ax.hist(sim["tauExitProb"], 100, log=True, facecolor="g") - ax.set_ylabel("Counts") - ax.set_xlabel(r"$P_\mathrm{exit}(\tau)$") - - -def tau_pexit_density(sim, fig, ax): - hist2d( - fig, - ax, - np.degrees(sim["beta_rad"]), - sim["tauExitProb"], - r"$\beta$", - r"$P_\mathrm{exit}(\tau)$", - cmap="viridis", - logx=False, - logy=True, - ) - - -def tau_betas(sim, fig, ax): - hist2d( - fig, - ax, - sim["beta_rad"], - sim["tauBeta"], - r"$\beta$", - r"$τ_β$", - logx=False, - logy=True, - ) - - -def decay_altitude_hist(sim, fig, ax): - ax.hist(sim["altDec"], 100, log=True) - ax.set_ylabel("Counts") - ax.set_xlabel("decay_altitude log(km)") - - -def num_photo_electrons_hist(sim, fig, ax): - m = sim["numPEs"] != 0 - ax.hist(np.log(sim["numPEs"][m]), 100, log=False) - ax.set_ylabel("Counts") - ax.set_xlabel("log(numPEs)") - - -def num_photo_electrons_density(sim, fig, ax): - hist2d( - fig, - ax, - sim["numPEs"], - np.degrees(sim["beta_rad"]), - "numPEs", - "β", - "plasma", - logx=True, - logy=False, - ) - - -def num_photo_electrons_altitude(sim, fig, ax): - hist2d( - fig, - ax, - sim["numPEs"], - sim["altDec"], - "numPEs", - "decay_altitude km", - "plasma", - logx=True, - logy=True, - ) - - -def decay_altitude(sim, fig, ax): - hist2d( - fig, - ax, - np.degrees(sim["beta_rad"]), - sim["altDec"], - "β", - "decay_altitude km", - cmap="viridis", - logx=False, - logy=True, - ) - - -def cherenkov_angle_hist(sim, fig, ax): - ax.hist(np.degrees(np.arccos(sim["costhetaChEff"])), 100, log=True) - ax.set_ylabel("Counts") - ax.set_xlabel("θ_chEff") - - -def eas_input_density(sim, fig, ax): - hist2d( - fig, - ax, - np.degrees(sim["beta_rad"]), - sim["altDec"], - "beta_rad", - "decay alt km", - cmap="jet", - logx=False, - logy=True, - ) - - -def cherenkov_angle(sim, fig, ax): - hist2d( - fig, - ax, - np.degrees(np.arccos(sim["costhetaChEff"])), - np.degrees(sim["beta_rad"]), - "θ_chEff", - "β", - cmap="jet", - logx=False, - logy=False, - ) - - -def eas_results_density(sim, fig, ax): - hist2d( - fig, - ax, - np.degrees(np.arccos(sim["costhetaChEff"])), - sim["numPEs"], - "θ_chEff", - "NumPEs", - cmap="jet", - logx=False, - logy=False, - ) - - -@decorators.ensure_plot_registry(dashboard) -def show_plot(sim, plot): - if dashboard.__name__ in plot: - dashboard(sim) - - # dashboard(sim, plot) diff --git a/test/utils/test_decorators.py b/test/utils/test_decorators.py index 2dd46f2..d0450c1 100644 --- a/test/utils/test_decorators.py +++ b/test/utils/test_decorators.py @@ -59,6 +59,14 @@ def test_nss_result_plot(): plot_written = False iA, iB = np.random.randn(2, 128) + plot_kwargs = { + "save_as": "pdf", + "pop_up": False, + "save_to_file": False, + "default_color": 0, + "default_colormap": "viridis", + "filename": "nuspacesim_run", + } def plotter(inputs, results, *args, **kwargs): nonlocal plot_written @@ -66,7 +74,7 @@ def plotter(inputs, results, *args, **kwargs): assert plot_written assert len(inputs) == 2 assert len(results) == 2 - assert len(args) == 0 + assert len(args) == 2 assert len(kwargs) == 0 assert np.array_equal(inputs[0], iA) assert np.array_equal(inputs[1], iB) @@ -84,35 +92,35 @@ def test_base_f(input1, input2): # test plotter is not called without a plot argument assert not plot_written - cA, cB = test_base_f(iA, iB) + cA, cB = test_base_f(iA, iB, plot_kwargs=plot_kwargs) assert not plot_written assert np.all(np.equal(cA, 0.0)) assert np.all(np.equal(cB, 1.0)) # test plotter is called with a callable plot argument plot_written = False - cA, cB = test_base_f(iA, iB, plot=plotter) + cA, cB = test_base_f(iA, iB, plot=plotter, plot_kwargs=plot_kwargs) assert plot_written assert np.all(np.equal(cA, 0.0)) assert np.all(np.equal(cB, 1.0)) # test plotter is called with a string plot argument plot_written = False - cA, cB = test_base_f(iA, iB, plot=plotter.__name__) + cA, cB = test_base_f(iA, iB, plot=plotter.__name__, plot_kwargs=plot_kwargs) assert plot_written assert np.all(np.equal(cA, 0.0)) assert np.all(np.equal(cB, 1.0)) # test plotter is called with a list of callable plot arguments plot_written = False - cA, cB = test_base_f(iA, iB, plot=list([plotter])) + cA, cB = test_base_f(iA, iB, plot=list([plotter]), plot_kwargs=plot_kwargs) assert plot_written assert np.all(np.equal(cA, 0.0)) assert np.all(np.equal(cB, 1.0)) # test plotter is called with a list of string plot arguments plot_written = False - cA, cB = test_base_f(iA, iB, plot=list([plotter.__name__])) + cA, cB = test_base_f(iA, iB, plot=list([plotter.__name__]), plot_kwargs=plot_kwargs) assert plot_written assert np.all(np.equal(cA, 0.0)) assert np.all(np.equal(cB, 1.0))