diff --git a/notebooks/LST1_observation_simulator.ipynb b/notebooks/LST1_observation_simulator.ipynb index 824d0bde1..919ab3748 100644 --- a/notebooks/LST1_observation_simulator.ipynb +++ b/notebooks/LST1_observation_simulator.ipynb @@ -19,7 +19,8 @@ "from pyirf.statistics import li_ma_significance\n", "from scipy.stats import moyal, norm, skewnorm\n", "from pathlib import Path\n", - "from lstchain.io.io import get_resource_path" + "from lstchain.io.io import get_resource_path\n", + "from lstchain import version" ] }, { @@ -156,6 +157,8 @@ "# For standard wobble offset (0.4 deg) reasonable values are alpha=0.333 (3 off regions) above 0.2 TeV, \n", "# and alpha=1 below 0.2 TeV. For testing sensitivity with the standard definition, set alpha=0.2\n", "\n", + "# Note: this setting will be overriden in pulsar mode! (see below)\n", + "\n", "alpha = 0.333 # 1" ] }, @@ -176,12 +179,12 @@ "outputs": [], "source": [ "pulsar_mode = False # Set to True to activate it\n", + "on_phase_interval = 0.043 # Crab P1: [-0.017, 0.026] Phase range for integrating the signal\n", + "off_phase_interval = 0.35 # [0.52 - 0.87] Phase range for estimating the background\n", "\n", - "on_phase_interval = 0.043 # Crab P1: [-0.017, 0.026] \n", - "off_phase_interval = 0.35 # [0.52 - 0.87]\n", - "\n", - "alpha = on_phase_interval / off_phase_interval\n", - "print(f'alpha = {alpha:.4f}')\n", + "if pulsar_mode:\n", + " alpha = on_phase_interval / off_phase_interval\n", + " print(f'alpha = {alpha:.4f}')\n", "\n", "# The spectrum is interpreted as average flux in full period (i.e. not just in the on-phase)" ] @@ -239,7 +242,6 @@ "else:\n", " print(f'Source angular radius within which {cut_efficiency:.1%} of the emission is contained: {source_radius}')\n", "\n", - "\n", " background_increase_factor = (theta_cut**2 + \n", " (source_radius.to_value(u.deg)*np.ones_like(theta_cut))**2) / (theta_cut**2)" ] @@ -585,7 +587,9 @@ "num_realizations_b = 100\n", "# Number of realizations (random numbers taken from skewnorm or moyal) for E-migration simulation\n", "\n", - "num_realizations = num_realizations_a * num_realizations_b" + "num_realizations = num_realizations_a * num_realizations_b\n", + "\n", + "np.random.seed(0)" ] }, { @@ -897,7 +901,7 @@ " SED_stat_error = SED*relative_stat_excess_error[displayed_points]\n", " SED_total_error = SED*total_relative_excess_error[displayed_points]\n", "\n", - " plt.figure(figsize=(10,6))\n", + " fig = plt.figure(figsize=(10,6))\n", " intrinsic_SED = (finebincenters*u.TeV)**2 * intrinsic_dFdE(finebincenters*u.TeV)\n", "\n", " if redshift > 0:\n", @@ -916,6 +920,26 @@ " fmt='o', markersize=4, \n", " label='Observed, stat-only')\n", "\n", + " ax = fig.axes[0]\n", + " \n", + " str = 'Source-independent analysis'\n", + " str += f'\\nZenith angle = {zenith[zd_bin]:.1f} degrees'\n", + " str += f'\\nCut efficiencies: {cut_efficiency:.0%}, {cut_efficiency:.0%} ($\\\\theta$, g/h)'\n", + " str += f'\\nEffective observation time = {effective_obs_time.to(u.h):.1f}'\n", + " if pulsar_mode:\n", + " str += '\\nPulsar mode'\n", + " str += f'\\nOn/Off exposure ratio $\\\\alpha$={alpha:.3f}'\n", + " if source_radius == 0:\n", + " str += '\\nPoint-like source'\n", + " else:\n", + " str += f'\\nSource radius: {source_radius:.2f}$^\\\\circ$ ({cut_efficiency:.0%} containment)'\n", + " str += f'\\n\\ncta-lstchain {version.version}'\n", + " \n", + " ax.text(0.05, 0.38, str, \n", + " transform=ax.transAxes, fontsize=12, verticalalignment='top',\n", + " bbox=dict(facecolor='white', alpha=0.8))\n", + "\n", + " \n", " plt.yscale('log')\n", " plt.xscale('log')\n", " plt.xlabel('Energy (TeV)', fontsize=14)\n", @@ -924,38 +948,42 @@ " np.nanmax((SED+SED_total_error)).value*5)\n", "\n", " plt.xlim(mean_etrue_vs_ereco[displayed_points][0]/3, mean_etrue_vs_ereco[displayed_points][-1]*3)\n", + "\n", " plt.legend()\n", " plt.grid()\n", + " plt.show()\n", + " \n", " \n", - " plt.figure(figsize=(10,2))\n", + " plt.figure(figsize=(10,3))\n", " if not pulsar_mode:\n", " plt.scatter(mean_etrue_vs_ereco[displayed_points], SED_total_error/SED, \n", - " s=8, label=f'Observed, stat + backg syst ({backg_systematics_uncertainty:.1%})')\n", + " s=12, label=f'Observed, stat + backg syst ({backg_systematics_uncertainty:.1%})')\n", " plt.scatter(mean_etrue_vs_ereco[displayed_points], SED_stat_error/SED, \n", - " s=8, label='Observed, stat-only')\n", + " s=12, label='Observed, stat-only')\n", " plt.xscale('log')\n", " plt.yscale('log')\n", " plt.xlabel('Energy (TeV)', fontsize=14)\n", " plt.ylabel('Relative SED \\n uncertainty', fontsize=14)\n", " plt.grid()\n", + " plt.legend()\n", " plt.show()\n", - " \n", + " \n", "else:\n", " print(\"Nothing to display. No valid flux points!\")" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "52", + "cell_type": "markdown", + "id": "c927ed19", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "### NOTE: in reality, other systematics from data-MC discrepancies (in Aeff, energy-dependent) may make the point-to point fluctuations larger than shown above. Systematic errors of around a few percent (independent in each bin) are to be expected.\n" + ] }, { "cell_type": "code", "execution_count": null, - "id": "53", + "id": "f3089dde", "metadata": {}, "outputs": [], "source": []