SST-1M analysis workflow
Automatic data processing
Raw data from both telescopes are processed automaticaly up to DL3 level
every morning after an observing night on Calculus with Daily analysis
pipeline:
/mnt/nfs22_auger/jurysek/sst1mpipe/sst1mpipe/scripts/daily_analysis.py
Products of automatic processing (all data levels up to DL3) can be found in:
/data/work/analysis/Daily_analysis/daily_data/
Those are sorted by date, source observed, and the telescope used (cs1, cs2 and stereo for coincident events).
All data level directories also contain configuration file and log file for each processed run, so
user can easily check if the processing went OK and which version of ``sst1mpipe``wasused. Random Forest models used in authomatic analysis
can be found in
/data/work/analysis/MC/prod_january_2023/
for both mono and stereo regimes of observation, together with expected performances and sensitivity curves. IRFs used for DL3 production can be found in
/data/work/analysis/IRFs/
DL3 products of the automatic processing can be further analyzed using standard tools such as gammapy, see High level analysis.
Note
As the software development is progressing very fast, automatic processing is usually
done with a slightly obsolete version of sst1mpipe and with a configuration not fully
tunned for atmospheric conditions during particular nights. To obtain paper level
analysis, one should proceed with the steps described in Detailed data processing.
Detailed data processing
All data files are stored on telescope servers cs1/cs2. They are also mirrored on Calculus
(/net/cs1/data/raw/, /net/cs2/data/raw/) and FZU farm (/mnt/sst1m/sst1m_backups/),
where can be processed and analyzed.
To get some feeling of how the event reconstruction and data analysis works, one can take a look at materials from the last SST-1M hands-on analysis session in Prague, December 2023:
Jupyter notebooks on Calculus:
/data/work/analysis/sst1m_analysis_handson_2023/
To understand general logic of the data processing, see scheme of the pipeline for mono and stereo data analysis: Analysis basics.
Configuration files
There are two slightly different configuration files. One is to be used for MC processing (needed only if one wants to get custom RF models / IRFs / performance curves), and one for real data processing. These contain configuration for all reconstruction/analysis steps from R0 to DL3 (or performance evaluation in case of MC).
Default config file for MC:
{
"telescope_coords": {
"tel_001": {
"long_deg": 14.782924,
"lat_deg": 49.911800,
"height_m": 510
},
"tel_002": {
"long_deg": 14.782078,
"lat_deg": 49.913084,
"height_m": 510
}
},
"allowed_tels": null,
"analysis": {
"global_theta_cut": 0.2,
"global_gammaness_cut": 0.6,
"gamma_efficiency": 0.6,
"significance_fit": true,
"log_energy_min_tev":-1,
"log_energy_max_tev": 3,
"n_energy_bins": 21,
"n_bins_migration": 51,
"fov_offset_max_deg": 6,
"nbins_fov_offset": 9,
"fov_offset_bg_max_deg": 8,
"nbins_fov_offset_bg": 17,
"source_offset_max_deg": 1,
"nbins_source_offset": 501,
"gamma_min_simulated_energy_tev": 0.2,
"proton_min_simulated_energy_tev": 0.4,
"n_training_events": 7E5,
"gamma_to_proton_training_ratio": 1.0,
"stereo_reco_weights": "camera_frame_hillas_intensity",
"observation_time_h": 50,
"off_regions": 5,
"bad_pixels": {
"tel_001": [],
"tel_002": []
}
},
"event_selection": {
"camera_frame_hillas_intensity": [0, 1000000],
"camera_frame_hillas_width": [0, 100],
"camera_frame_hillas_length": [0, 100],
"leakage_intensity_width_2": [0, 1]
},
"CameraCalibrator": {
"image_extractor_type" : "LocalPeakWindowSum",
"LocalPeakWindowSum": {
"apply_integration_correction": true,
"window_shift": 3,
"window_width": 7
},
"invalid_pixel_handler_type": "NeighborAverage"
},
"NsbCalibrator": {
"intensity_correction": {
"tel_001": 1.0,
"tel_002": 1.0
}
},
"ImageProcessor": {
"image_cleaner_type": "ImageCleanerSST",
"use_telescope_frame": false,
"TailcutsImageCleaner": {
"picture_threshold_pe": [
["id", 1, 8],
["id", 2, 9]
],
"boundary_threshold_pe": [
["id", 1, 4],
["id", 2, 4]
],
"min_picture_neighbors": [
["id", 1, 2],
["id", 2, 2]
]
},
"ImageCleanerSST": {
"telescope_defaults": {
"tel_001": [
{
"min_nsb_level": 0,
"stdev_scaling": 2.5,
"picture_threshold_pe": 7,
"boundary_threshold_pe": 2.5,
"min_picture_neighbors": 2,
"keep_isolated": false,
"min_time_neighbors" : 1,
"time_limit_ns" : 8,
"only_main_island" : false
}, {
"min_nsb_level": 100000,
"stdev_scaling": 2.5,
"picture_threshold_pe": 8,
"boundary_threshold_pe": 4,
"min_picture_neighbors": 2,
"keep_isolated": false,
"min_time_neighbors" : 1,
"time_limit_ns" : 8,
"only_main_island" : false
}],
"tel_002": [
{
"min_nsb_level": 0,
"stdev_scaling": 2.5,
"picture_threshold_pe": 8,
"boundary_threshold_pe": 2.5,
"min_picture_neighbors": 2,
"keep_isolated": false,
"min_time_neighbors" : 1,
"time_limit_ns" : 8,
"only_main_island" : false
}, {
"min_nsb_level": 100000,
"stdev_scaling": 2.5,
"picture_threshold_pe": 8,
"boundary_threshold_pe": 4,
"min_picture_neighbors": 2,
"keep_isolated": false,
"min_time_neighbors" : 1,
"time_limit_ns" : 8,
"only_main_island" : false
}]
}
}
},
"ShowerProcessor": {
"reconstructor_types": ["HillasReconstructor"]
},
"random_forest_regressor_args": {
"max_depth": 30,
"min_samples_leaf": 10,
"n_jobs": -1,
"n_estimators": 150,
"bootstrap": true,
"criterion": "squared_error",
"max_leaf_nodes": null,
"min_impurity_decrease": 0.0,
"min_samples_split": 10,
"min_weight_fraction_leaf": 0.0,
"oob_score": false,
"random_state": 42,
"verbose": 0,
"warm_start": false
},
"random_forest_classifier_args": {
"max_depth": 30,
"min_samples_leaf": 10,
"n_jobs": -1,
"n_estimators": 100,
"criterion": "gini",
"min_samples_split": 10,
"min_weight_fraction_leaf": 0.0,
"max_leaf_nodes": null,
"min_impurity_decrease": 0.0,
"bootstrap": true,
"oob_score": false,
"random_state": 42,
"verbose": 0,
"warm_start": false,
"class_weight": null
},
"energy_regression_features": [
"log_camera_frame_hillas_intensity",
"camera_frame_hillas_width",
"camera_frame_hillas_length",
"camera_frame_hillas_wl",
"camera_frame_hillas_skewness",
"camera_frame_hillas_kurtosis",
"camera_frame_timing_slope",
"leakage_intensity_width_2",
"camera_frame_hillas_x",
"camera_frame_hillas_y",
"HillasReconstructor_tel_impact_distance",
"HillasReconstructor_h_max"
],
"particle_classification_features": [
"log_camera_frame_hillas_intensity",
"camera_frame_hillas_width",
"camera_frame_hillas_length",
"camera_frame_hillas_wl",
"camera_frame_hillas_skewness",
"camera_frame_hillas_kurtosis",
"camera_frame_timing_slope",
"leakage_intensity_width_2",
"camera_frame_hillas_x",
"camera_frame_hillas_y",
"HillasReconstructor_tel_impact_distance",
"HillasReconstructor_h_max",
"log_reco_energy",
"reco_disp_norm"
],
"disp_method": "disp_norm_sign",
"disp_regression_features": [
"log_camera_frame_hillas_intensity",
"camera_frame_hillas_width",
"camera_frame_hillas_length",
"camera_frame_hillas_wl",
"camera_frame_hillas_skewness",
"camera_frame_hillas_kurtosis",
"camera_frame_timing_slope",
"leakage_intensity_width_2",
"HillasReconstructor_tel_impact_distance",
"HillasReconstructor_h_max"
],
"disp_classification_features": [
"log_camera_frame_hillas_intensity",
"camera_frame_hillas_width",
"camera_frame_hillas_length",
"camera_frame_hillas_wl",
"camera_frame_hillas_skewness",
"camera_frame_hillas_kurtosis",
"camera_frame_timing_slope",
"leakage_intensity_width_2"
]
}
Default config file for data:
{ "array_center_coords": {
"long_deg": 14.7825010,
"lat_deg": 49.9124420,
"height_m": 511.195000
},
"telescope_coords": {
"tel_021": {
"long_deg": 14.782924,
"lat_deg": 49.911800,
"height_m": 510
},
"tel_022": {
"long_deg": 14.782078,
"lat_deg": 49.913084,
"height_m": 510
}
},
"telescope_calibration": {
"tel_021": null,
"tel_022": null,
"bad_calib_px_interpolation": true,
"dynamic_dead_px_interpolation": true
},
"window_transmittance": {
"tel_021": null,
"tel_022": null
},
"telescope_equivalent_focal_length": {
"tel_021": 5.6,
"tel_022": 5.6
},
"stereo": {
"event_matching_method" : "WhiteRabbitClosest",
"SlidingWindow": {
"window_half_width": 1E-3,
"offset_search": {
"time_start": -1E-2,
"time_stop": 1E-2,
"time_step": 2E-5
}
},
"WhiteRabbitClosest": {
"max_time_diff_ns": 1E4
},
"SWATEventIDs": {}
},
"analysis": {
"global_theta_cut": 0.2,
"global_gammaness_cut": 0.75,
"gamma_efficiency": 0.7,
"log_energy_min_tev":-1,
"log_energy_max_tev": 3,
"n_energy_bins": 21,
"gamma_min_simulated_energy_tev": 0.2,
"proton_min_simulated_energy_tev": 0.4,
"n_training_events": 5E5,
"gamma_to_proton_training_ratio": 1.0,
"stereo_reco_weights": "camera_frame_hillas_intensity",
"observation_time_h": 50,
"off_regions": 5,
"bad_pixels": {
"tel_021": [],
"tel_022": []
},
"ped_time_window": 10
},
"event_selection": {
"camera_frame_hillas_intensity": [50, 1000000],
"camera_frame_hillas_width": [0, 100],
"camera_frame_hillas_length": [0, 100],
"leakage_intensity_width_2": [0, 0.7]
},
"CameraCalibrator": {
"image_extractor_type" : "LocalPeakWindowSum",
"LocalPeakWindowSum": {
"apply_integration_correction": true,
"window_shift": 3,
"window_width": 7
},
"invalid_pixel_handler_type": "NeighborAverage"
},
"NsbCalibrator": {
"apply_pixelwise_Vdrop_correction": false,
"apply_global_Vdrop_correction": false,
"intensity_correction": {
"tel_021": 1.0,
"tel_022": 1.0
}
},
"ImageProcessor": {
"image_cleaner_type": "ImageCleanerSST",
"use_telescope_frame": false,
"TailcutsImageCleaner": {
"picture_threshold_pe": [
["id", 21, 8],
["id", 22, 9]
],
"boundary_threshold_pe": [
["id", 21, 4],
["id", 22, 4]
],
"min_picture_neighbors": [
["id", 21, 2],
["id", 22, 2]
]
},
"ImageCleanerSST": {
"min_number_pedestals": 100,
"telescope_defaults": {
"tel_021": [
{
"min_nsb_level": 0,
"stdev_scaling": 2.5,
"picture_threshold_pe": 7,
"boundary_threshold_pe": 2.5,
"min_picture_neighbors": 2,
"keep_isolated": false,
"min_time_neighbors" : 1,
"time_limit_ns" : 8,
"only_main_island" : false
}, {
"min_nsb_level": 2.3526,
"stdev_scaling": 2.5,
"picture_threshold_pe": 8,
"boundary_threshold_pe": 4,
"min_picture_neighbors": 2,
"keep_isolated": false,
"min_time_neighbors" : 1,
"time_limit_ns" : 8,
"only_main_island" : false
}, {
"min_nsb_level": 2.8863,
"stdev_scaling": 2.5,
"picture_threshold_pe": 8,
"boundary_threshold_pe": 4,
"min_picture_neighbors": 2,
"keep_isolated": false,
"min_time_neighbors" : 1,
"time_limit_ns" : 8,
"only_main_island" : false
}],
"tel_022": [
{
"min_nsb_level": 0,
"stdev_scaling": 2.5,
"picture_threshold_pe": 8,
"boundary_threshold_pe": 2.5,
"min_picture_neighbors": 2,
"keep_isolated": false,
"min_time_neighbors" : 1,
"time_limit_ns" : 8,
"only_main_island" : false
}, {
"min_nsb_level": 2.5288,
"stdev_scaling": 2.5,
"picture_threshold_pe": 8,
"boundary_threshold_pe": 4,
"min_picture_neighbors": 2,
"keep_isolated": false,
"min_time_neighbors" : 1,
"time_limit_ns" : 8,
"only_main_island" : false
}, {
"min_nsb_level": 3.2747,
"stdev_scaling": 2.5,
"picture_threshold_pe": 8,
"boundary_threshold_pe": 4,
"min_picture_neighbors": 2,
"keep_isolated": false,
"min_time_neighbors" : 1,
"time_limit_ns" : 8,
"only_main_island" : false
}]
}
}
},
"ShowerProcessor": {
"reconstructor_types": ["HillasReconstructor"]
},
"random_forest_regressor_args": {
"max_depth": 30,
"min_samples_leaf": 10,
"n_jobs": -1,
"n_estimators": 150,
"bootstrap": true,
"criterion": "squared_error",
"max_leaf_nodes": null,
"min_impurity_decrease": 0.0,
"min_samples_split": 10,
"min_weight_fraction_leaf": 0.0,
"oob_score": false,
"random_state": 42,
"verbose": 0,
"warm_start": false
},
"random_forest_classifier_args": {
"max_depth": 30,
"min_samples_leaf": 10,
"n_jobs": -1,
"n_estimators": 100,
"criterion": "gini",
"min_samples_split": 10,
"min_weight_fraction_leaf": 0.0,
"max_leaf_nodes": null,
"min_impurity_decrease": 0.0,
"bootstrap": true,
"oob_score": false,
"random_state": 42,
"verbose": 0,
"warm_start": false,
"class_weight": null
},
"energy_regression_features": [
"log_camera_frame_hillas_intensity",
"camera_frame_hillas_width",
"camera_frame_hillas_length",
"camera_frame_hillas_wl",
"camera_frame_hillas_skewness",
"camera_frame_hillas_kurtosis",
"camera_frame_timing_slope",
"leakage_intensity_width_2",
"camera_frame_hillas_x",
"camera_frame_hillas_y",
"HillasReconstructor_tel_impact_distance",
"HillasReconstructor_h_max"
],
"particle_classification_features": [
"log_camera_frame_hillas_intensity",
"camera_frame_hillas_width",
"camera_frame_hillas_length",
"camera_frame_hillas_wl",
"camera_frame_hillas_skewness",
"camera_frame_hillas_kurtosis",
"camera_frame_timing_slope",
"leakage_intensity_width_2",
"camera_frame_hillas_x",
"camera_frame_hillas_y",
"HillasReconstructor_tel_impact_distance",
"HillasReconstructor_h_max",
"log_reco_energy",
"reco_disp_norm"
],
"disp_method": "disp_norm_sign",
"disp_regression_features": [
"log_camera_frame_hillas_intensity",
"camera_frame_hillas_width",
"camera_frame_hillas_length",
"camera_frame_hillas_wl",
"camera_frame_hillas_skewness",
"camera_frame_hillas_kurtosis",
"camera_frame_timing_slope",
"leakage_intensity_width_2",
"HillasReconstructor_tel_impact_distance",
"HillasReconstructor_h_max"
],
"disp_classification_features": [
"log_camera_frame_hillas_intensity",
"camera_frame_hillas_width",
"camera_frame_hillas_length",
"camera_frame_hillas_wl",
"camera_frame_hillas_skewness",
"camera_frame_hillas_kurtosis",
"camera_frame_timing_slope",
"leakage_intensity_width_2"
],
"mean_charge_to_nsb_rate": {
"tel_021": [
{
"mean_charge_bin_low": 0,
"nsb_rate": 80
}, {
"mean_charge_bin_low": 2.3526,
"nsb_rate": 150
}, {
"mean_charge_bin_low": 2.8863,
"nsb_rate": 270
}
],
"tel_022": [
{
"mean_charge_bin_low": 0,
"nsb_rate": 80
}, {
"mean_charge_bin_low": 2.5288,
"nsb_rate": 150
}, {
"mean_charge_bin_low": 3.2747,
"nsb_rate": 270
}
]
}
}R0 to DL1
To calibrate raw data (R0) or MC (R1) and process them to DL1, one may run sst1mpipe_r0_dl1 script. Inputs are a single raw
.fits.fz data file (containing single telescope data) or .simtel.gz output file of sim_telarray (may contain mono and coincident events
triggering more telescopes).
Output is HDF5 file with a table of DL1 parameters, one output file per run (about 40 seconds of data) and per telescope (DL1 files from individual
telescopes searched for coincidences and merged for stereoscopic reconstruction in the next step DL1 to DL1 stereo).
It applies dc to p.e. calibration on raw R0 waveforms, integrates them, cleans the images (gets rid of the noisy pixels which likely do not
contain any Cherenkov photons), and parametrizes the shower images with Hillas ellipses.
See --help for possible inputs. Some of them, which might not be obvious:
--px-charges- the script stores also distribution of all integrated charges in individual pixels for all events merged. This is useful for further MC/data tunning and to get some impression on the level of NSB in the data.--precise-timestamps- stores also White Rabbit timestamps in the DL1 output with the precision needed for coincident events matching. Keep it on for all data taken after 25th September 2023, when WR was deployed.--pointing-ra/decand--force-pointing- allows to specify the telescope pointing direction. To process data taken after begining of September 2023 it can be ignored (i.e. do not use it for any new data), because the pointing coordinates are being written automaticaly in the fits file header during the datataking and the script understands where to look for it.—-reclean- experimental method of data re-cleaning based on pixel charge variation. For now it needs distribution of pixel charges stored in the first pass of the script (--px-charges). I.e. to apply re-cleaning, one has to run the script for the second time with the—-recleanswitch.
Relevant parts of the config file applied in this analysis step:
telescope_calibration- calibration files based on analysis of dark runs. Should be taken relatively close to the date of observationwindow_transmittance- files with for camera window transmittance correction (measured in the lab and can be kept default)CameraCalibrator- Pulse integration settingsImageProcessor- Settings of image cleaning method and tailcuts. Tailcuts can be set different for different levels of NSBShowerProcessor- Method of shower geometry reconstruction. Only applied if event source contains data from more telescopes, i.e. it’s only relevant for MC in this analysis step.
Random Forest training
Note
In most cases, analyser does not need to train dedicated Random Forest models and this step can be safely skipped using pre-trained RFs
referenced in DL1 to DL2. Training of dedicated RFs is, however, necessary in some performance studies if one wants to use different
configuration for sst1mpipe_r0_dl1 than MC was processed with (e.g. cleaning, peak integration, ..).
Random Forests can be trained on DL1 MC diffuse gammas and diffuse proton files using sst1mpipe_mc_train_rfs script (see
--help for possible inputs). Before running sst1mpipe_mc_train_rfs it is useful to merge many small DL1 files in given MC production (which resulted
from paralelized MC simulations) into a single file per particle with sst1mpipe_merge_hdf5 script to reach satisfactory
statistics for RF training. Outputs are trained models in the scikit.learn format (.sav). There is RF classifier for gamma/hadron
separation, RF regressor for energy reconstruction, and either RF regressor (disp_vector) or RF regressor+classifier (disp_norm_sign)
for arrival direction reconstruction depenting on the method selected (disp_method field in the cfg file).
RF are trained for each telescope, even in case of stereo reconstruction. In stereo, we only use extra stereo features,
which are reconstructed geometricaly, such as HillasReconstructor_h_max and HillasReconstructor_tel_impact_distance.
Then, in DL1 to DL2, reconstruction is performed for each telescope independently, and final reconstructed quantities are
obtained as weighted average of the values for each telescope (except for direction recontruction where MARS-like approach is adopted).
Relevant parts of the config file applied in this analysis step:
Setup of the forests and training procedure
random_forest_regressor_args,random_forest_classifier_argsLists of Random Forest features used for the reconstruction -
energy_regression_features,disp_regression_features,disp_classification_features,particle_classification_features. The very features used for the RF training have to be used later in DL1 to DL2 reconstruction!n_training_events- Total number of events used for individual RF training. I.e. ifn_training_events=200000, 200k diffuse gammas are used for energy regressor and DISP regressor and classifier, and 100k diffuse gammas + 100k diffuse protons is used for particle classifier (ifgamma_to_proton_training_ratio=1).gamma_to_proton_training_ratio- Ratio of gammas and protons in training sample for particle classifier.
DL1 to DL1 stereo
For stereo reconstruction, coincident events have first to be find. In current implementation, tel2 DL1 files are searched
for each tel1 DL1 file to find the closest tel2 event for each tel1 event. Coincident event search results in a new DL1 file containing events from both
telescopes, matched by their event_id. Only coincident events are stored in resulting DL1 files.
This is performed by script sst1mpipe_data_dl1_dl1_stereo (see --help for possible inputs). Input is a single DL1 file from tel1
and a directory with all relevant DL1 files for tel2. Coincidence finder is driven by the config file field stereo. Possible
options are:
SlidingWindow- For analysis of the data without precise White Rabbit timestamps (i.e. taken before 25th September 2023) one needs to use this method. It first searches for the time offset between the two DL1 tables providing maximum number of coindicent events and then selects the closest ones.WhiteRabbitClosest- Works on data with precise WR timestamps in the DL1 table, i.e. all data taken after 25th September 2023. It only finds the closest tel2 event to each tel1 event (precision of WR is high enough to avoid random coincidences for usual trigger rates of the telescopes).SWATEventIDs- After 30th January 2024 the coincident events are tagged by SWAT, providing them with the samearrayEvtNum, resulting in the sameevent_idin the DL1 files. The DL1 events can be then matched just based on theirevent_id.
Note
sst1mpipe_data_dl1_dl1_stereo is not intended to be run on MC, as in MC DL1 the coincident events are already matched by their event_id (mono events are in MC DL1
tables as well, so those can be used for both mono and stereo analysis).
DL1 to DL2
This step uses pre-trained Random Forests to reconstruct parameters of primary gamma-ray photon (gammaness, direction and energy) using Hillas parameters stored in
the DL1 files as features. One can run sst1mpipe_data_dl1_dl2 stript on either mono DL1 files (outputs of sst1mpipe_r0_dl1) for each telescope separately (using RFs for mono reconstruction),
or on stereo DL1 containing coincident events only (outputs of sst1mpipe_data_dl1_dl1_stereo). The script can handle both types of DL1,
but stereo reconstruction has to be requested explicitely using -—stereo switch. RFs trained on MC can be found on Calculus for both mono and stereo
reconstruction and different zenith angles:
/data/work/analysis/MC/prod_january_2023/$SST1MPIPE_VER/models_mono_psf_vaod0.2//data/work/analysis/MC/prod_january_2023/$SST1MPIPE_VER/models_stereo_psf_vaod0.2/
Note
One should always use RF models trained with the same sst1mpipe version that is used for the analysis.
Relevant parts of the config file applied in this analysis step:
Random Forest features used for the reconstruction -
energy_regression_features,disp_regression_features,disp_classification_features,particle_classification_features. These should be the very same features as those used for RF training (check cfg files stored in the directories together with the models)disp_method- Direction reconstruction method used. For now we only usedisp_norm_signwhich requires RF regressor to reconstruct source distance from the image Center of Gravity, and RF classifier to determine on which side along the main axis of the Hillas ellipse the source lies.stereo_reco_weights- Parameter used as a weight for averaging the stereo reconstructed parameters.
DL2 MC to IRFs
Note
In most cases, analyser does not need to produce own Instrument Response Functions and this step can be safely skipped using default IRFs referenced in Automatic data processing. IRF production, however, is necessary in performance studies, or if one uses custom RFs to produce DL2, or applies custom selection cuts in DL2 to DL3 step.
To make IRFs from MC DL2 files, one can run sst1mpipe_mc_make_irfs script, which currently produces only full enclosure IRFs, so it has to be provided with
diffuse protons and diffuse gammas. The script applies event selection cuts defined in the config file (event_selection), including cut on gammaness.
The gammaness cut can be either global (one number independent on energy) or energy dependent (gammaness distribution naturaly depends on energy, so using optimised
energy dependent gammaness cut results in performance improvement). Global gammaness cut can be set in the config file (global_gammaness_cut field), while energy dependent
cuts, must be provided as HDF5 table, where the cuts for individual energy bins are stored, using parameter
--gammaness-cut-dir. These tables can be generated with sst1mpipe_mc_performance (see RF performance and sensitivity). Pre-calculated energy dependent gammaness
cuts are stored on Calculus for mono/stereo and different zenith angles:
/data/work/analysis/MC/prod_january_2023/$SST1MPIPE_VER/performance/*_performance_*
Again, one should make sure that the event selection applied to produce the cuts is the same as for IRFs. The IRF maker creates
some directory structure inside the --output-dir, automaticaly recognizing proper bin in zenith, azimuth, NSB level and gammaness cut applied.
This directory structure should remain untouched for DL2 to DL3 to work properly.
Output IRF files are fully compatible with gammapy and may be read and explored with the use of gammapy funkcionalities:
from gammapy.irf import (
EffectiveAreaTable2D,
PSF3D,
EnergyDispersion2D
Background2D
)
irf_filename = 'SST1M_tel_021_Zen30deg_gcut0.75_irfs.fits'
aeff = EffectiveAreaTable2D.read(irf_filename, hdu="EFFECTIVE AREA")
edisp = EnergyDispersion2D.read(irf_filename, hdu="ENERGY DISPERSION")
psf = PSF3D.read(irf_filename, hdu='POINT SPREAD FUNCTION')
bg_2d = Background2D.read(irf_filename, hdu='BACKGROUND')
DL2 to DL3
sst1mpipe_data_dl2_dl3 is a tool to create DL3 data files from DL2 data files. It is supposed to be provided with a directory with input DL2 files
(typicaly a directory with DL2 for one source observed in one night, but can be run on larger sample as well). It merges the DL2 HDF5 per-run
files into per-wobble DL3 fits files containing only photon lists. It also finds proper IRF based on zenith, azimuth and NSB level for each input DL2 file. It creates
per-night index files needed for further analysis in gammapy. The script applies event selection cuts defined in the config file (event_selection),
including cut on gammaness, so one should make sure that these are the same as used on MC to produce IRFs. Energy dependent gammaness cuts can be used as well, following
the same rules described in DL2 MC to IRFs.
High level analysis
Output DL3 files produced with sst1mpipe_data_dl2_dl3 are fully compatible with gammapy and may be further analyzed using gammapy tools. See e.g.
A typical use case is to run joint gammapy analysis on data from several nights. In such case one has to run create_hdu_indexes script to create
HDU index files indexing all DL3s to be used in the final analysis.
RF performance and sensitivity
TBD