Simulation

We can bring together the detector and source modelling to calculate the expected number of neutrino events and run simulations.

[1]:
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.patches import Circle
from matplotlib.collections import PatchCollection
import h5py

Defining a source and detector model

[2]:
from icecube_tools.detector.effective_area import EffectiveArea
from icecube_tools.detector.energy_resolution import EnergyResolution
from icecube_tools.detector.angular_resolution import AngularResolution
from icecube_tools.detector.detector import IceCube
from icecube_tools.source.flux_model import PowerLawFlux
from icecube_tools.source.source_model import DiffuseSource, PointSource
from icecube_tools.detector.r2021 import R2021IRF
[3]:
# Define detector (see detector model notebook for more info)
aeff = EffectiveArea.from_dataset("20210126", "IC86_II")
irf = R2021IRF.from_period("IC86_II")
# IceCube expects an instance of EffectiveAerea, AngularResolution and
# EnergyResolution, optionally the period (here IC86_I)
# R2021IRF inherits from AngularResolution and EnergyResolution
# just to be able to be used as both
detector = IceCube(aeff, irf, irf, "IC86_II")
[4]:
# Define simple sources (see source model notebook for more info)
diff_flux_norm = 3e-21  # Flux normalisation in units of GeV^-1 cm^-2 s^-1 sr^-1
point_flux_norm = 5e-19  # Flux normalisation in units of GeV^-1 cm^-2 s^-1
norm_energy = 1e5  # Energy of normalisation in units of GeV
min_energy = 1e2  # GeV
max_energy = 1e8  # GeV

diff_power_law = PowerLawFlux(diff_flux_norm, norm_energy, 3.7, min_energy, max_energy)
diff_source = DiffuseSource(diff_power_law, z=0.0)

point_power_law = PowerLawFlux(
    point_flux_norm, norm_energy, 2.5, min_energy, max_energy
)
point_source = PointSource(point_power_law, z=0.0, coord=(np.pi, np.deg2rad(30)))
sources = [diff_source, point_source]

Expected number of neutrino events

Sometimes we just want to predict the number of events from sources in a detector without specifying all detector properties or running a simulation. We can do this with the NeutrinoCalculator. For this, we just need a source list and an effective area.

[5]:
from icecube_tools.neutrino_calculator import NeutrinoCalculator, PhiSolver
[6]:
nu_calc = NeutrinoCalculator(sources, aeff)
nu_calc(
    time=1,  # years
    min_energy=min_energy,
    max_energy=max_energy,  # energy range
    min_cosz=-1,
    max_cosz=1,
)  # cos(zenith) range
[6]:
[185.18577999563084, 24.514085562926446]

The calculator returns a list of expected neutrino event numbers, one for each source.

We may also want to do the inverse, and find the PointSource flux normalisation corresponding to an expected number of events. For this there is the PhiSolver class.

[7]:
phi_solver = PhiSolver(
    aeff, norm_energy, min_energy, max_energy, time=1, min_cosz=-1, max_cosz=1
)
phi_norm = phi_solver(Nnu=15, dec=30, index=2.0)  # degrees  # spectral index
phi_norm  # GeV^-1 cm^-2 s^-1
[7]:
9.368172433902621e-19

Set up and run simulation

[8]:
from icecube_tools.simulator import Simulator, TimeDependentSimulator
[9]:
# Set up simulation
simulator = Simulator(sources, detector)
simulator.time = 1  # year

# Run simulation
simulator.run(show_progress=True, seed=42)
INFO:icecube_tools.simulator:Instantiating simulation.
INFO:icecube_tools.simulator:Random N.
INFO:icecube_tools.simulator:Simulating source 0
INFO:icecube_tools.simulator:Done sampling the spectrum
INFO:icecube_tools.simulator:Sampling angular uncertainty for diffuse source
INFO:icecube_tools.simulator:Simulating source 1
INFO:icecube_tools.simulator:Done sampling the spectrum
INFO:icecube_tools.simulator:Sampling angular uncertainty for point source
INFO:icecube_tools.simulator:Creating array of simulation data

This way, the simulator calculates the expected number of neutrinos from these sources given the observation period, effective area and relevant source properties. We note that we could also run a simulation for a fixed number of neutrinos if we want, simply by passing the optional argument N to simulator.run().

[10]:
# simulator.run(N=10, show_progress=True, seed=42)
[11]:
# Save to file
simulator.save("data/sim_output_86_II.h5")

# Load
with h5py.File("data/sim_output_86_II.h5", "r") as f:
    true_energy = f["true_energy"][()]
    reco_energy = f["reco_energy"][()]
    ra = f["ra"][()]
    dec = f["dec"][()]
    ang_err = f["ang_err"][()]
    source_label = f["source_label"][()]
[12]:
"""
for i in [1.5, 2.0, 2.5, 3.0, 3.5, 3.7]:
    norm_energy = 1e5 # Energy of normalisation in units of GeV
    min_energy = 1e2 # GeV
    max_energy = 1e8 # GeV
    phi_solver = PhiSolver(aeff, norm_energy, min_energy, max_energy,
                           time=1, min_cosz=-1, max_cosz=1)
    phi_norm = phi_solver(Nnu=2000,
                          dec=30, # degrees
                          index=i) # spectral index
    phi_norm # GeV^-1 cm^-2 s^-1
    point_power_law = PowerLawFlux(phi_norm, norm_energy, i,
                                   min_energy, max_energy)
    point_source = PointSource(point_power_law, z=0., coord=(np.pi, np.deg2rad(30)))
    sources = [point_source]
    simulator = Simulator(sources, detector)
    simulator.time = 1 # year

    # Run simulation
    simulator.run(show_progress=True, seed=42)
    simulator.save(f"data/sim_output_{i:.1f}.h5")
"""
[12]:
'\nfor i in [1.5, 2.0, 2.5, 3.0, 3.5, 3.7]:\n    norm_energy = 1e5 # Energy of normalisation in units of GeV\n    min_energy = 1e2 # GeV\n    max_energy = 1e8 # GeV\n    phi_solver = PhiSolver(aeff, norm_energy, min_energy, max_energy, \n                           time=1, min_cosz=-1, max_cosz=1)\n    phi_norm = phi_solver(Nnu=2000, \n                          dec=30, # degrees\n                          index=i) # spectral index\n    phi_norm # GeV^-1 cm^-2 s^-1\n    point_power_law = PowerLawFlux(phi_norm, norm_energy, i, \n                                   min_energy, max_energy)\n    point_source = PointSource(point_power_law, z=0., coord=(np.pi, np.deg2rad(30)))\n    sources = [point_source]\n    simulator = Simulator(sources, detector)\n    simulator.time = 1 # year\n\n    # Run simulation\n    simulator.run(show_progress=True, seed=42)\n    simulator.save(f"data/sim_output_{i:.1f}.h5")\n'
[13]:
# Plot energies
bins = np.geomspace(1e2, max_energy)
fig, ax = plt.subplots()
ax.hist(true_energy, bins=bins, alpha=0.7, label="E_true")
ax.hist(reco_energy, bins=bins, alpha=0.7, label="E_reco")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("E [GeV]")
ax.legend()
[13]:
<matplotlib.legend.Legend at 0x7fdc04e49be0>
../_images/notebooks_simulation_18_1.png
[14]:
# Plot directions
ps_sel = source_label == 1

fig, ax = plt.subplots(subplot_kw={"projection": "aitoff"})
fig.set_size_inches((12, 7))

circles = []
for r, d, a in zip(ra[~ps_sel], dec[~ps_sel], ang_err[~ps_sel]):
    circle = Circle((r - np.pi, d), radius=np.deg2rad(a))
    circles.append(circle)
df_nu = PatchCollection(circles)

circles = []
for r, d, a in zip(ra[ps_sel], dec[ps_sel], ang_err[ps_sel]):
    circle = Circle((r - np.pi, d), radius=np.deg2rad(a))
    circles.append(circle)
ps_nu = PatchCollection(circles, color="r")

ax.add_collection(df_nu)
ax.add_collection(ps_nu)

ax.grid()
../_images/notebooks_simulation_19_0.png

Time dependent simulation

We can simulate an observation campaign spanning multiple data periods of IceCube through a “meta class” TimeDependentSimulator:

[15]:
tsim = TimeDependentSimulator(["IC86_I", "IC86_II"], sources)
/opt/hostedtoolcache/Python/3.9.13/x64/lib/python3.9/site-packages/icecube_tools/detector/r2021.py:48: RuntimeWarning: divide by zero encountered in log10
  self.dataset[:, 6:-1] = np.log10(self.dataset[:, 6:-1])
INFO:icecube_tools.simulator:Instantiating simulation.
INFO:icecube_tools.simulator:Instantiating simulation.
WARNING:icecube_tools.simulator:Need to set simulation times, defaults to 1 year each.

We need to set simulation times for all periods. Since for past periods the simulation time shouldn’t be larger than the actual observation time (that is time span - down time of the detector) we need to take care, or rather, we let the class Uptime take care:

[16]:
from icecube_tools.utils.data import Uptime

It lets us calculate the actual observation time through, e.g. IC86_II, vs the time covered:

[17]:
uptime = Uptime()
uptime.time_obs("IC86_II"), uptime.time_span("IC86_II")
[17]:
(<Quantity 0.90890393 yr>, <Quantity 1.01572353 yr>)

We can further define a start and end time of an observation and let it calculate the observation time in each period. Viable possible options are - start and end time in MJD - start time in MJD and duration in years - end time in MJD and duration in years

If the start time is before the earliest period (IC_40), the start time will be set to the earliest possible date.

If the end time is past the last period (IC86_II), then we get an according message and simulate into the future.

We can of course bypass this time setting and operate directly on the instances of Simulator, for example if we’d want to build up large statistics for subsequent likelihood analysis.

[18]:
times = uptime.find_obs_time(start=55569, duration=3)
times
INFO:icecube_tools.utils.data:End time outside of provided data set, sending an owl to Professor Trelawney
[18]:
{'IC86_I': 1.2241905309188956, 'IC86_II': 1.5222023283463137}

The returned dictionary can be used to set the simulation times for an instance of TimeDependentSimulator:

[19]:
tsim = TimeDependentSimulator(["IC86_I", "IC86_II"], sources)
tsim.time = times
/opt/hostedtoolcache/Python/3.9.13/x64/lib/python3.9/site-packages/icecube_tools/detector/r2021.py:48: RuntimeWarning: divide by zero encountered in log10
  self.dataset[:, 6:-1] = np.log10(self.dataset[:, 6:-1])
INFO:icecube_tools.simulator:Instantiating simulation.
INFO:icecube_tools.simulator:Instantiating simulation.
WARNING:icecube_tools.simulator:Need to set simulation times, defaults to 1 year each.
INFO:icecube_tools.simulator:Set simulation times

The simulation is started by calling run(), results can be saved by save(file_prefix), with the filename being {file_prefix}_{p}.h5 with period p.

[20]:
tsim.run()
tsim.save("data", "test_sim")
INFO:icecube_tools.simulator:Simulating period IC86_I.
INFO:icecube_tools.simulator:Random N.
INFO:icecube_tools.simulator:Simulating source 0
INFO:icecube_tools.simulator:Done sampling the spectrum
INFO:icecube_tools.simulator:Sampling angular uncertainty for diffuse source
INFO:icecube_tools.simulator:Simulating source 1
INFO:icecube_tools.simulator:Done sampling the spectrum
INFO:icecube_tools.simulator:Sampling angular uncertainty for point source
INFO:icecube_tools.simulator:Creating array of simulation data
INFO:icecube_tools.simulator:Simulating period IC86_II.
INFO:icecube_tools.simulator:Random N.
INFO:icecube_tools.simulator:Simulating source 0
INFO:icecube_tools.simulator:Done sampling the spectrum
INFO:icecube_tools.simulator:Sampling angular uncertainty for diffuse source
INFO:icecube_tools.simulator:Simulating source 1
INFO:icecube_tools.simulator:Done sampling the spectrum
INFO:icecube_tools.simulator:Sampling angular uncertainty for point source
INFO:icecube_tools.simulator:Creating array of simulation data
[20]:
{'IC86_I': 'data/p_IC86_I_test_sim.h5',
 'IC86_II': 'data/p_IC86_II_test_sim.h5'}
[21]:
"""
for i in [1.5, 2.0, 2.5, 3.0, 3.5, 3.7]:
    norm_energy = 1e5 # Energy of normalisation in units of GeV
    min_energy = 1e2 # GeV
    max_energy = 1e8 # GeV
    phi_solver = PhiSolver(aeff, norm_energy, min_energy, max_energy,
                           time=1, min_cosz=-1, max_cosz=1)
    phi_norm = phi_solver(Nnu=2000,
                          dec=30, # degrees
                          index=i) # spectral index
    phi_norm # GeV^-1 cm^-2 s^-1
    point_power_law = PowerLawFlux(phi_norm, norm_energy, i,
                                   min_energy, max_energy)
    point_source = PointSource(point_power_law, z=0., coord=(np.pi, np.deg2rad(30)))
    sources = [point_source]
    tsim = TimeDependentSimulator(["IC86_I", "IC86_II"], sources)
    for sim in tsim.simulators.values():
        sim.time = 1
    tsim.run()
    tsim.save(f"index_{i}")
"""
[21]:
'\nfor i in [1.5, 2.0, 2.5, 3.0, 3.5, 3.7]:\n    norm_energy = 1e5 # Energy of normalisation in units of GeV\n    min_energy = 1e2 # GeV\n    max_energy = 1e8 # GeV\n    phi_solver = PhiSolver(aeff, norm_energy, min_energy, max_energy, \n                           time=1, min_cosz=-1, max_cosz=1)\n    phi_norm = phi_solver(Nnu=2000, \n                          dec=30, # degrees\n                          index=i) # spectral index\n    phi_norm # GeV^-1 cm^-2 s^-1\n    point_power_law = PowerLawFlux(phi_norm, norm_energy, i, \n                                   min_energy, max_energy)\n    point_source = PointSource(point_power_law, z=0., coord=(np.pi, np.deg2rad(30)))\n    sources = [point_source]\n    tsim = TimeDependentSimulator(["IC86_I", "IC86_II"], sources)\n    for sim in tsim.simulators.values():\n        sim.time = 1\n    tsim.run()\n    tsim.save(f"index_{i}")\n'
[ ]: