Merge pull request #132 from vuillaut/magic_sens

Add MAGIC sensitivity data and plot
cta-observatory · Nov 17, 2020 · 68bf82d · 68bf82d
2 parents fc76451 + fffb1c2
commit 68bf82d
Show file tree

Hide file tree

Showing 8 changed files with 203 additions and 20 deletions.
diff --git a/ctaplot/ana/ana.py b/ctaplot/ana/ana.py
@@ -7,6 +7,7 @@
 import numpy as np
 from ..io import dataset as ds
 from scipy.stats import binned_statistic, norm
+from astropy.io.ascii import read
 
 _relative_scaling_method = 's1'
 
@@ -42,6 +43,7 @@
            'bias_per_bin',
            'percentile_confidence_interval',
            'logbin_mean',
+           'get_magic_sensitivity',
            ]
 
 
@@ -1135,3 +1137,16 @@ def bias_per_energy(simu, reco, energy, relative_scaling_method=None):
 
     return bias_per_bin(simu, reco, energy, relative_scaling_method=relative_scaling_method, bins=energy_bin)
 
+
+def get_magic_sensitivity():
+    """
+    Load MAGIC differential sensitivity data from file `magic_sensitivity_2014.ecsv`.
+    Extracted from table A.7 in Aleksić, Jelena, et al. "The major upgrade of the MAGIC telescopes,
+    Part II: A performance study using observations of the Crab Nebula." Astroparticle Physics 72 (2016): 76-94.,
+    DOI: 10.1016/j.astropartphys.2015.02.005'
+
+    Returns
+    -------
+    `astropy.table.table.QTable`
+    """
+    return read(ds.get('magic_sensitivity_2014.ecsv'))
diff --git a/ctaplot/ana/tests/test_ana.py b/ctaplot/ana/tests/test_ana.py
@@ -3,6 +3,7 @@
 
 np.random.seed(42)
 
+
 def test_logspace_decades_nbin():
     ca = ana.logspace_decades_nbin(0.1, 10, n=9)
     assert len(ca) == 19
@@ -56,7 +57,7 @@ def test_resolution():
 
     # Test bias
     bias = np.random.rand() * 100
-    y_reco_bias = np.random.normal(loc=loc+bias, scale=scale, size=size)
+    y_reco_bias = np.random.normal(loc=loc + bias, scale=scale, size=size)
 
     assert np.isclose(ana.resolution(y_true, y_reco_bias,
                                      bias_correction=True,
@@ -91,7 +92,7 @@ def test_resolution_per_bin():
         np.testing.assert_allclose(res[:, 0], scale / ana.relative_scaling(y_true, y_reco, scaling).mean(), rtol=1e-1)
 
     bias = 50
-    y_reco = np.random.normal(loc=loc+bias, scale=scale, size=size)
+    y_reco = np.random.normal(loc=loc + bias, scale=scale, size=size)
     bins, res = ana.resolution_per_bin(x, y_true, y_reco, bias_correction=True)
     np.testing.assert_allclose(res[:, 0], scale / ana.relative_scaling(y_true, y_reco).mean(), rtol=1e-1)
 
@@ -147,16 +148,16 @@ def test_distance2d_resolution():
     simu_y = np.ones(size)
     # reconstructed positions on a circle around true position
     t = np.random.rand(size) * np.pi * 2
-    reco_x = 1 + 3*np.cos(t)
-    reco_y = 1 + 3*np.sin(t)
+    reco_x = 1 + 3 * np.cos(t)
+    reco_y = 1 + 3 * np.sin(t)
     res, err_min, err_max = ana.distance2d_resolution(reco_x, reco_y, simu_x, simu_y,
                                                       percentile=68.27, confidence_level=0.95, bias_correction=False)
 
     np.testing.assert_equal(res, 3)
 
     # with different bias on X and Y:
-    reco_x = 5 + 2*np.cos(t)
-    reco_y = 7 + 2*np.sin(t)
+    reco_x = 5 + 2 * np.cos(t)
+    reco_y = 7 + 2 * np.sin(t)
     res, err_min, err_max = ana.distance2d_resolution(reco_x, reco_y, simu_x, simu_y,
                                                       percentile=68.27, confidence_level=0.95, bias_correction=True)
 
@@ -169,16 +170,16 @@ def test_distance2d_resolution_per_bin():
     simu_x = np.ones(size)
     simu_y = np.ones(size)
     t = np.random.rand(size) * np.pi * 2
-    reco_x = 3*np.cos(t)
-    reco_y = 3*np.sin(t)
+    reco_x = 3 * np.cos(t)
+    reco_y = 3 * np.sin(t)
 
     bin, res = ana.distance2d_resolution_per_bin(x, reco_x, reco_y, simu_x, simu_y, bins=10, bias_correction=True)
 
     np.testing.assert_allclose(res, 3, rtol=1e-1)
 
 
 def test_angular_resolution():
-    size=10000
+    size = 10000
     simu_alt = np.random.rand(size)
     simu_az = np.random.rand(size)
     scale = 0.01
@@ -215,7 +216,6 @@ def test_angular_resolution_small_angles():
                                rtol=1e-1,
                                )
 
-
     simu_az = np.random.rand(size)
     simu_alt = np.zeros(size)
     reco_alt = simu_alt
@@ -225,10 +225,12 @@ def test_angular_resolution_small_angles():
                                rtol=1e-1,
                                )
 
+
 def test_bias_empty():
     x = np.empty(0)
     assert ana.bias(x, x) == 0
 
+
 def test_bias_per_bin():
     size = 100000
     simu = np.ones(size)
@@ -244,4 +246,12 @@ def test_bias_per_energy():
     reco = np.random.normal(loc=2, scale=0.5, size=size)
     energy = np.logspace(-2, 2, size)
     bins, bias = ana.bias_per_energy(simu, reco, energy)
-    np.testing.assert_allclose(bias, 1, rtol=1e-1)
+    np.testing.assert_allclose(bias, 1, rtol=1e-1)
+
+
+def test_get_magic_sensitivity():
+    from astropy.table.table import QTable
+    import astropy.units as u
+    table = ana.get_magic_sensitivity()
+    assert type(table) is QTable
+    assert table['e_min'][0] == 63 * u.GeV
diff --git a/ctaplot/io/dataset.py b/ctaplot/io/dataset.py
@@ -89,6 +89,7 @@
                   'HESS_Impact_Energy_Bias_Std_Mono.txt',
                   'HESS_Impact_Effective_Area_Safe_Mono.txt',
                   'HESS_Impact_Energy_Bias_Stereo.txt',
+                  'magic_sensitivity_2014.ecsv'
                   ]
 
 

diff --git a/ctaplot/plots/plots.py b/ctaplot/plots/plots.py
@@ -14,6 +14,7 @@
 from sklearn import metrics
 from sklearn.multiclass import LabelBinarizer
 from ..io.dataset import load_any_resource
+import astropy.units as u
 
 
 __all__ = ['plot_resolution',
@@ -62,6 +63,7 @@
            'plot_roc_curve_multiclass',
            'plot_roc_curve_gammaness_per_energy',
            'plot_gammaness_distribution',
+           'plot_sensitivity_magic_performance'
            ]
 
 
@@ -2177,3 +2179,49 @@ def plot_gammaness_distribution(mc_type, gammaness, ax=None, **kwargs):
     ax.set_xlabel('gammaness')
     ax.legend()
     return ax
+
+
+def plot_sensitivity_magic_performance(key='lima_5off', ax=None, **kwargs):
+    """
+
+    Parameters
+    ----------
+    key: string
+        'lima_1off': LiMa 1 off position
+        'lima_3off': LiMa 3 off positions
+        'lima_5off': LiMa 5 off positions
+        'snr': Nex/sqrt(Nbkg)
+    ax: `matplotlib.pyplot.axis`
+    kwargs: kwargs for `matplotlib.pyplot.errorbar`
+
+    Returns
+    -------
+    `matplotlib.pyplot.axis`
+    """
+
+    ax = plt.gca() if ax is None else ax
+
+    magic_table = ana.get_magic_sensitivity()
+
+    magic_table['e_err_lo'] = magic_table['e_center'] - magic_table['e_min']
+    magic_table['e_err_hi'] = magic_table['e_max'] - magic_table['e_center']
+
+    kwargs.setdefault('ls', '')
+    kwargs.setdefault('label', f'MAGIC {key}')
+
+    k = 'sensitivity_' + key
+    ax.errorbar(
+        magic_table['e_center'].to_value(u.TeV),
+        y=(magic_table['e_center'] ** 2 * magic_table[k]).to_value(u.Unit('erg cm-2 s-1')),
+        xerr=[magic_table['e_err_lo'].to_value(u.TeV), magic_table['e_err_hi'].to_value(u.TeV)],
+        yerr=(magic_table['e_center'] ** 2 * magic_table[f'{k}_err']).to_value(u.Unit('erg cm-2 s-1')),
+        **kwargs
+    )
+
+    ax.set_xscale('log')
+    ax.set_yscale('log')
+    ax.grid(True, which='both')
+    ax.legend()
+
+    return ax
+
diff --git a/ctaplot/plots/tests/test_plots.py b/ctaplot/plots/tests/test_plots.py
@@ -4,6 +4,7 @@
 
 np.random.seed(42)
 
+
 def plot_energy_distribution():
     SimuE = np.random.rand(100)
     RecoE = np.random.rand(10)
@@ -269,4 +270,9 @@ def test_plot_gammaness_distribution():
     nb_events = 1000
     mc_type = np.random.choice([0, 1, 2, 3], size=nb_events)
     gammaness = np.random.rand(nb_events)
-    plots.plot_gammaness_distribution(mc_type, gammaness)
+    plots.plot_gammaness_distribution(mc_type, gammaness)
+
+
+def test_plot_sensitivity_magic_performance():
+    ax = plots.plot_sensitivity_magic_performance(key='lima_5off')
+    plots.plot_sensitivity_magic_performance(key='lima_3off', ax=ax, color='black', ls='--')
diff --git a/environment.yml b/environment.yml
@@ -11,8 +11,8 @@ dependencies:
   - numpy
   - python >= 3.6
   - pandas
-  - scikit-learn
-  - scipy
-  - nbsphinx
+  - pyaml
   - pytables
   - pytest
+  - scikit-learn
+  - scipy
diff --git a/setup.py b/setup.py
@@ -37,11 +37,7 @@ def get_property(prop, project):
         'scikit-learn',
         'jupyter',
         'ipywidgets',
-        'recommonmark',
-        'sphinx>=1.4',
-        'nbsphinx',
-        'sphinx_rtd_theme',
-        'nbconvert',
+        'pyaml',
     ],
     tests_require=['pytest'],
     classifiers=[

diff --git a/share/iact_sensitivities/magic_sensitivity_2014.ecsv b/share/iact_sensitivities/magic_sensitivity_2014.ecsv
@@ -0,0 +1,107 @@
+# %ECSV 0.9
+# ---
+# datatype:
+# - {name: e_max, unit: GeV, datatype: float64}
+# - {name: e_min, unit: GeV, datatype: float64}
+# - {name: background_rate, unit: 1 / min, datatype: float64}
+# - {name: background_rate_err, unit: 1 / min, datatype: float64}
+# - {name: gamma_rate, unit: 1 / min, datatype: float64}
+# - {name: gamma_rate_err, unit: 1 / min, datatype: float64}
+# - {name: sensitivity_lima_1off, unit: 1 / (cm2 s TeV), datatype: float64}
+# - {name: sensitivity_lima_1off_err, unit: 1 / (cm2 s TeV), datatype: float64}
+# - {name: sensitivity_lima_3off, unit: 1 / (cm2 s TeV), datatype: float64}
+# - {name: sensitivity_lima_3off_err, unit: 1 / (cm2 s TeV), datatype: float64}
+# - {name: sensitivity_lima_5off, unit: 1 / (cm2 s TeV), datatype: float64}
+# - {name: sensitivity_lima_5off_err, unit: 1 / (cm2 s TeV), datatype: float64}
+# - {name: sensitivity_snr, unit: 1 / (cm2 s TeV), datatype: float64}
+# - {name: sensitivity_snr_err, unit: 1 / (cm2 s TeV), datatype: float64}
+# - {name: e_center, unit: GeV, datatype: float64}
+# - {name: e_err_lo, unit: GeV, datatype: float64}
+# - {name: e_err_hi, unit: GeV, datatype: float64}
+# meta: !!omap
+# - {COMMENT: "Table extracted from table A.7 in Aleksi\u0107, Jelena, et al. \"The major upgrade of the MAGIC telescopes, Part II: A\
+#     \ performance study using observations of the Crab Nebula.\" Astroparticle Physics 72 (2016): 76-94., DOI: 10.1016/j.astropartphys.2015.02.005\
+#     \ and multiplied by the CRAB spectrum from  \"Measurement of the Crab Nebula spectrum over three decades in energy with the MAGIC\
+#     \ telescopes\", Aleksi\u0107 2015, DOI:10.1016/j.jheap.2015.01.002"}
+# - __serialized_columns__:
+#     background_rate:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: !astropy.units.Unit {unit: 1 / min}
+#       value: !astropy.table.SerializedColumn {name: background_rate}
+#     background_rate_err:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: !astropy.units.Unit {unit: 1 / min}
+#       value: !astropy.table.SerializedColumn {name: background_rate_err}
+#     e_center:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: &id001 !astropy.units.Unit {unit: GeV}
+#       value: !astropy.table.SerializedColumn {name: e_center}
+#     e_err_hi:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: &id002 !astropy.units.Unit {unit: GeV}
+#       value: !astropy.table.SerializedColumn {name: e_err_hi}
+#     e_err_lo:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: *id001
+#       value: !astropy.table.SerializedColumn {name: e_err_lo}
+#     e_max:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: *id002
+#       value: !astropy.table.SerializedColumn {name: e_max}
+#     e_min:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: *id002
+#       value: !astropy.table.SerializedColumn {name: e_min}
+#     gamma_rate:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: !astropy.units.Unit {unit: 1 / min}
+#       value: !astropy.table.SerializedColumn {name: gamma_rate}
+#     gamma_rate_err:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: !astropy.units.Unit {unit: 1 / min}
+#       value: !astropy.table.SerializedColumn {name: gamma_rate_err}
+#     sensitivity_lima_1off:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: &id003 !astropy.units.Unit {unit: 1 / (cm2 s TeV)}
+#       value: !astropy.table.SerializedColumn {name: sensitivity_lima_1off}
+#     sensitivity_lima_1off_err:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: *id003
+#       value: !astropy.table.SerializedColumn {name: sensitivity_lima_1off_err}
+#     sensitivity_lima_3off:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: *id003
+#       value: !astropy.table.SerializedColumn {name: sensitivity_lima_3off}
+#     sensitivity_lima_3off_err:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: *id003
+#       value: !astropy.table.SerializedColumn {name: sensitivity_lima_3off_err}
+#     sensitivity_lima_5off:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: *id003
+#       value: !astropy.table.SerializedColumn {name: sensitivity_lima_5off}
+#     sensitivity_lima_5off_err:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: *id003
+#       value: !astropy.table.SerializedColumn {name: sensitivity_lima_5off_err}
+#     sensitivity_snr:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: !astropy.units.Unit {unit: 1 / (cm2 s TeV)}
+#       value: !astropy.table.SerializedColumn {name: sensitivity_snr}
+#     sensitivity_snr_err:
+#       __class__: astropy.units.quantity.Quantity
+#       unit: !astropy.units.Unit {unit: 1 / (cm2 s TeV)}
+#       value: !astropy.table.SerializedColumn {name: sensitivity_snr_err}
+# schema: astropy-2.0
+e_max e_min background_rate background_rate_err gamma_rate gamma_rate_err sensitivity_lima_1off sensitivity_lima_1off_err sensitivity_lima_3off sensitivity_lima_3off_err sensitivity_lima_5off sensitivity_lima_5off_err sensitivity_snr sensitivity_snr_err e_center e_err_lo e_err_hi
+100.0 63.0 4.06 0.08 3.01 0.13 7.601958508074096e-10 3.455435685488226e-11 6.1333983417416e-10 2.5915767641161692e-11 5.874240665329984e-10 2.5915767641161692e-11 7.3e-10 3e-11 79.37253933193772 16.372539331937716 20.627460668062284
+158.0 100.0 2.41 0.06 4.29 0.12 1.650780087464995e-10 4.845056860484262e-12 1.339312146433864e-10 3.806830390380492e-12 1.270097048426946e-10 3.460754900345902e-12 1.37e-10 5e-12 125.69805089976535 25.69805089976535 32.30194910023465
+251.0 158.0 0.54 0.03 3.37 0.08 3.910257867777991e-11 1.325511141619658e-12 3.154716517054787e-11 1.0604089132957266e-12 2.982400068644231e-11 1.0604089132957266e-12 3.05e-11 1.3000000000000001e-12 199.14316458266902 41.143164582669016 51.856835417330984
+398.0 251.0 0.066 0.01 1.36 0.05 1.3547361416709163e-11 9.67668672622083e-13 1.0450821664318497e-11 7.741349380976664e-13 9.821837027114142e-12 7.257515044665623e-13 9.300000000000001e-12 8e-13 316.066448709761 65.066448709761 81.933551290239
+631.0 398.0 0.027 0.006 1.22 0.04 3.5556881058906414e-12 3.047732662191978e-13 2.7260275478494916e-12 3.047732662191978e-13 2.556709066616604e-12 2.539777218493315e-13 2.3e-12 3e-13 501.13670789516107 103.13670789516107 129.86329210483893
+1000.0 631.0 0.0133 0.0018 0.88 0.04 1.2365835413508895e-12 6.806881878995721e-14 9.30273856796082e-13 5.105161409246791e-14 8.678774395719544e-13 6.239641722412745e-14 7.2e-13 6e-14 794.35508432942 163.35508432942004 205.64491567057996
+1585.0 1000.0 0.0059 0.0007 0.58 0.03 4.5104238659926285e-13 2.0005912308838273e-14 3.2736947414462623e-13 1.6368473707231312e-14 3.0190740393337754e-13 1.6368473707231312e-14 2.3000000000000003e-13 1.8e-14 1258.967831201417 258.9678312014171 326.0321687985829
+2512.0 1585.0 0.0027 0.0005 0.3 0.02 2.1199333369145533e-13 1.115754387849765e-14 1.4504807042046943e-13 1.0599666684572767e-14 1.3165901776627226e-13 1.0041789490647884e-14 9e-14 1e-14 1995.374651537901 410.37465153790095 516.625348462099
+3981.0 2512.0 0.002 0.0005 0.166 0.016 1.015262173688021e-13 8.187598174903395e-15 7.041334430416919e-14 6.550078539922716e-15 6.22257461292658e-14 6.550078539922716e-15 4.1e-14 7e-15 3162.3206668521143 650.3206668521143 818.6793331478857
+6310.0 3981.0 0.0014 0.0003 0.093 0.012 4.689996060266254e-14 4.59803535320221e-15 3.126664040177503e-14 3.2186247472415472e-15 2.804801565453348e-14 3.2186247472415472e-15 1.7e-14 3e-15 5011.996608139315 1030.996608139315 1298.003391860685
+10000.0 6310.0 0.0046 0.0015 0.06 0.01 2.717581659166208e-14 3.705793171590284e-15 1.9764230248481512e-14 3.705793171590284e-15 1.8528965857951418e-14 2.470528781060189e-15 1.3e-14 3e-15 7943.550843294201 1633.5508432942006 2056.4491567057994