Towards using Pyirf

benchmarks/DL3/benchmarks_DL3_IRFs_and_sensitivity.ipynb

@@ -253,7 +253,7 @@
     "\n",
     "def plot_energy_response(ax, rmf, label, energy_range, **kwargs):\n",
     "    \"\"\"Simple function to plot the energy resolution (Gammapy)\"\"\"\n",
-    "    energy = rmf.e_true.nodes\n",
+    "    energy = rmf.e_true.center\n",
     "    bias = rmf.get_bias(energy)\n",
     "    resol = rmf.get_resolution(energy)\n",
     "    \n",
@@ -512,7 +512,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.8"
+   "version": "3.7.4"
   }
  },
  "nbformat": 4,

environment.yml

@@ -4,8 +4,8 @@ channels:
   - conda-forge
   - cta-observatory
 dependencies:
-  - ctapipe=0.7.0
-  - gammapy=0.8
+  - ctapipe
+  - gammapy
   - ipython
   - nbsphinx
   - numpydoc
@@ -14,3 +14,6 @@ dependencies:
   - sphinx-automodapi
   - sphinx_rtd_theme
   - vitables
+  - pip:
+    - Cython
+    - pyirf

protopipe/aux/example_config_files/protopipe/performance.yaml

@@ -1,11 +1,13 @@
+
 general:
  # Directory with input data file
- # [...] = your analysis local full path OUTSIDE the Vagrant box
+ # [...] = your analysis local path
+ # Please refer to directory structure shown at Lugano
  indir: '[...]/data/DL2'
  # Template name for input file
- template_input_file: 'DL2_{}_{}_merged.h5' # filled with mode and particle type
+ template_input_file: 'dl2_{}_{}_merged.h5' # filled with mode and particle type
  # Directory for output files
- outdir: '[...]/data/DL3'
+ outdir: '[...]/performance'
  # Output table name
  output_table_name: 'table_best_cutoff'
 
@@ -18,7 +20,10 @@ analysis:
 
  # Normalisation between ON and OFF regions
  alpha: 0.2
-
+
+ # Radius to use for calculating bg rate
+ max_bg_radius: 1.
+
  # Minimimal significance
  min_sigma: 5
 
@@ -42,30 +47,31 @@ analysis:
 
 particle_information:
  gamma:
-  n_events_per_file: 1000000  #  10**5 * 10
+  num_use: 10
+  num_showers: 10**5
   e_min: 0.003
   e_max: 330
   gen_radius: 1400
-  diff_cone: 0
-  gen_gamma: 2
+  gen_gamma: -2
+  diff_cone: 0 
 
  proton:
-  n_events_per_file: 4000000  #  2 * 10**5 * 20
+  num_use: 20
+  num_showers: 2 * 10**5
   e_min: 0.004
   e_max: 600
   gen_radius: 1900
+  gen_gamma: -2
   diff_cone: 10
-  gen_gamma: 2
-  offset_cut: 1.
 
  electron:
-  n_events_per_file: 2000000  #  10**5 * 20
+  num_use: 20
+  num_showers: 10**5
   e_min: 0.003
   e_max: 330
   gen_radius: 1900
+  gen_gamma: -2
   diff_cone: 10
-  gen_gamma: 2
-  offset_cut: 1.
 
 column_definition:
  # Column name for true energy
@@ -76,4 +82,4 @@ column_definition:
  classification_output:
   name: 'gammaness'
   range: [0, 1]
- angular_distance_to_the_src: 'xi'
+ angular_distance_to_the_src: 'xi'

protopipe/perf/cut_optimisation.py

@@ -14,10 +14,8 @@
 class CutsApplicator(object):
     """
     Apply best cut and angular cut to events.
-
     Apply cuts to gamma, proton and electrons that will be further used for
     performance estimation (irf, sensitivity, etc.).
-
     Parameters
     ----------
     config: `dict`
@@ -54,7 +52,6 @@ def apply_cuts(self):
     def apply_cuts_on_data(self, data):
         """
         Flag particle passing angular cut and the best cutoff
-
         Parameters
         ----------
         data: `pandas.DataFrame`
@@ -103,7 +100,6 @@ def apply_cuts_on_data(self, data):
 class CutsDiagnostic(object):
     """
     Class used to get some diagnostic related to the optimal working point.
-
     Parameters
     ----------
     config: `dict`
@@ -274,7 +270,6 @@ class CutsOptimisation(object):
     """
     Class used to find best cutoff to obtain minimal
     sensitivity in a given amount of time.
-
     Parameters
     ----------
     config: `dict`
@@ -290,7 +285,6 @@ def __init__(self, config, evt_dict, verbose_level=0):
     def weight_events(self, model_dict, colname_mc_energy):
         """
         Add a weight column to the files, in order to scale simulated data to reality.
-
         Parameters
         ----------
         model_dict: dict
@@ -320,7 +314,7 @@ def compute_weight(self, energy, particle, model):
             omega_simu = 1.
 
         nsimu = conf_part['n_simulated']
-        index_simu = conf_part['gen_gamma']
+        index_simu = np.absolute(conf_part['gen_gamma'])
         emin = conf_part['e_min'] * u.TeV
         emax = conf_part['e_max'] * u.TeV
         amplitude = 1. * u.Unit('1 / (cm2 s TeV)')
@@ -349,7 +343,6 @@ def find_best_cutoff(self, energy_values, angular_values):
         of energy and theta square cut. Correct the number of events
         according to the ON region which correspond to the angular cut applied to
         the gamma-ray events.
-
         Parameters
         ----------
         energy_values: `astropy.Quantity`
@@ -411,10 +404,10 @@ def find_best_cutoff(self, energy_values, angular_values):
                 # Correct number of background due to acceptance
                 acceptance_g = 2 * np.pi * (1 - np.cos(th_cut.to('rad').value))
                 acceptance_p = 2 * np.pi * (
-                        1 - np.cos(self.config['particle_information']['proton']['offset_cut'] * u.deg.to('rad'))
+                        1 - np.cos(self.config["analysis"]["max_bg_radius"] * u.deg.to('rad'))
                 )
                 acceptance_e = 2 * np.pi * (
-                        1 - np.cos(self.config['particle_information']['electron']['offset_cut'] * u.deg.to('rad'))
+                        1 - np.cos(self.config["analysis"]["max_bg_radius"] * u.deg.to('rad'))
                 )
 
                 # Add corrected weight taking into angular cuts applied to gamma-rays
@@ -691,6 +684,4 @@ def write_results(self, outdir, outfile, format, overwrite=True):
         t['angular_cut'] = Column(angular_cut, unit='TeV')
         for feature in feature_to_save:
             t[feature[0]] = Column(res_to_save[feature[0]])
-        t.write(os.path.join(outdir, outfile), format=format, overwrite=overwrite)
-
-
+        t.write(os.path.join(outdir, outfile), format=format, overwrite=overwrite)

protopipe/perf/irf_maker.py

@@ -16,7 +16,6 @@
 class BkgData(object):
     """
     Class storing background data in a NDDataArray object.
-
     It's a bit hacky, but gammapy sensitivity estimator does not take individual IRF,
     it takes a CTAPerf object. So i'm emulating that... We need a bkg format!!!
     """
@@ -44,7 +43,6 @@ def __init__(self, bkg, aeff, rmf):
 class SensitivityMaker(object):
     """
     Class which estimate sensitivity with IRF
-
     Parameters
     ----------
     config: `dict`
@@ -158,7 +156,6 @@ def add_sensitivity_to_irf(self):
 class IrfMaker(object):
     """
     Class building IRF for point-like analysis.
-
     Parameters
     ----------
     config: `dict`
@@ -293,7 +290,7 @@ def make_bkg_rate(self):
                 * (
                     1
                     - np.cos(
-                        self.config["particle_information"]["proton"]["offset_cut"]
+                        self.config["analysis"]["max_bg_radius"]
                         * u.deg.to("rad")
                     )
                 )
@@ -304,7 +301,7 @@ def make_bkg_rate(self):
                 * (
                     1
                     - np.cos(
-                        self.config["particle_information"]["electron"]["offset_cut"]
+                        self.config["analysis"]["max_bg_radius"]
                         * u.deg.to("rad")
                     )
                 )
@@ -374,11 +371,11 @@ def make_effective_area(
 
         # Get simulation infos
         cfg_particule = self.config["particle_information"]["gamma"]
-        simu_index = cfg_particule["gen_gamma"]
+        simu_index = np.absolute(cfg_particule["gen_gamma"])
         index = 1.0 - simu_index  # for futur integration
         nsimu_tot = float(cfg_particule["n_files"]) * float(
-            cfg_particule["n_events_per_file"]
-        )
+            cfg_particule["num_showers"]) * float(
+            cfg_particule["num_use"])
         emin_simu = cfg_particule["e_min"]
         emax_simu = cfg_particule["e_max"]
         area_simu = (np.pi * cfg_particule["gen_radius"] ** 2) * u.Unit("m2")
@@ -586,4 +583,4 @@ def _make_edisp_hdu(cls, table_energy, table_migra, table_theta, matrix):
         hdu.header.set(
             "EXTNAME", "ENERGY DISPERSION", "name of this binary table extension "
         )
-        return hdu
+        return hdu

protopipe/scripts/make_performance.py

@@ -64,15 +64,15 @@ def main():
     for particle in ['electron', 'proton']:
         # print('Initial stat: {} {}'.format(len(evt_dict[particle]), particle))
         evt_dict[particle] = evt_dict[particle].query('offset <= {}'.format(
-            cfg['particle_information'][particle]['offset_cut'])
+            cfg['analysis']['max_bg_radius'])
         )
 
     # Add required data in configuration file for future computation
     for particle in particles:
         cfg['particle_information'][particle]['n_files'] = \
             len(np.unique(evt_dict[particle]['obs_id']))
         cfg['particle_information'][particle]['n_simulated'] = \
-            cfg['particle_information'][particle]['n_files'] * cfg['particle_information'][particle]['n_events_per_file']
+            cfg['particle_information'][particle]['n_files'] * cfg['particle_information'][particle]['num_showers'] * cfg['particle_information'][particle]['num_use']
 
     # Define model for the particles
     model_dict = {'gamma': CrabSpectrum('hegra').model,
@@ -182,4 +182,4 @@ def main():
 
 
 if __name__ == '__main__':
-    main()
+    main()

protopipe/scripts/write_dl2.py

@@ -223,8 +223,16 @@ class RecoEvent(tb.IsDescription):
             impact_dict,
         ) in preper.prepare_event(source):
 
+            # True energy
+            true_energy = event.mc.energy.value
+
+            # True direction
+            true_az = event.mc.az
+            true_alt = event.mc.alt
+
             # Angular quantities
             run_array_direction = event.mcheader.run_array_direction
+            pointing_az, pointing_alt = run_array_direction[0], run_array_direction[1]
 
             # Angular separation between true and reco direction
             xi = angular_separation(
@@ -233,8 +241,8 @@ class RecoEvent(tb.IsDescription):
 
             # Angular separation bewteen the center of the camera and the reco direction.
             offset = angular_separation(
-                run_array_direction[0],  # az
-                run_array_direction[1],  # alt
+                pointing_az,
+                pointing_alt,
                 reco_result.az,
                 reco_result.alt,
             )
@@ -348,6 +356,11 @@ class RecoEvent(tb.IsDescription):
                     + n_tels["SST_ASTRI_ASTRICam"]
                     + n_tels["SST_GCT_CHEC"]
                 )
+                reco_event["pointing_az"] = pointing_az.to("deg").value
+                reco_event["pointing_alt"] = pointing_alt.to("deg").value
+                reco_event["true_az"] = true_az.to("deg").value
+                reco_event["true_alt"] = true_alt.to("deg").value
+                reco_event["true_energy"] = true_energy
                 reco_event["reco_energy"] = reco_energy
                 reco_event["reco_alt"] = alt.to("deg").value
                 reco_event["reco_az"] = az.to("deg").value