Skip to content

Commit

Permalink
Addressed several comments by reviewers
Browse files Browse the repository at this point in the history
  • Loading branch information
moralejo committed Nov 15, 2024
1 parent 94e764a commit eee78fd
Showing 1 changed file with 47 additions and 19 deletions.
66 changes: 47 additions & 19 deletions notebooks/LST1_observation_simulator.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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"
]
},
{
Expand Down Expand Up @@ -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"
]
},
Expand All @@ -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)"
]
Expand Down Expand Up @@ -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)"
]
Expand Down Expand Up @@ -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)"
]
},
{
Expand Down Expand Up @@ -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",
Expand All @@ -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",
Expand All @@ -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": []
Expand Down

0 comments on commit eee78fd

Please sign in to comment.