<jupyter_start><jupyter_text>Notebook for testing the reconstrcution of LST1 dataThis Notebook shows the steps followed within cta-lstchain to reconstruct event data from a corsika-simtelarray file to a dl2 data level, with reconstructed energy, direction and gamma/hadron separation.To run this notebook you will need the last version of cta-lstchain: *git clone https://github.com/cta-observatory/cta-lstchain*        Also, you will need the datafiles from **cta-lstchain-extra**:git clone https://github.com/cta-observatory/cta-lstchain-extra Some imports...<jupyter_code>import joblib
from ctapipe.utils import get_dataset_path
from lstchain.io.config import get_standard_config
from copy import copy
import matplotlib.pyplot as plt
import os
import pandas as pd
import lstchain.reco.r0_to_dl1 as calib
import lstchain.reco.dl1_to_dl2 as reco
import lstchain.visualization.plot_dl2 as lstplot

%matplotlib inline
plt.rcParams['figure.figsize'] = (20, 15)
plt.rcParams['font.size'] = 14<jupyter_output>/home/luca.foffano/workspace/software/miniconda3/envs/lst/lib/python3.7/site-packages/corsikaio/subblocks/dtypes.py:20: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'.
  return np.dtype(dict(**dt))<jupyter_text>Import a simtelarray file of MC simulated data:<jupyter_code>#datafile = get_dataset_path('gamma_test_large.simtel.gz')

# Local files can be also imported.
datafile = "/fefs/aswg/data/mc/DL0/20190415/gamma/south_pointing/gamma_20deg_180deg_run9___cta-prod3-demo-2147m-LaPalma-baseline-mono_off0.4.simtel.gz"<jupyter_output><empty_output><jupyter_text>Calibration and parameters extraction        We have a function for extracting all the image information:        - Charge extraction    - Cleaning    - Hillas Parameters    - Timing parametersWe store this information in a custom container which we write into a **hdf5** file.<jupyter_code>#Limit the number of events to analyze, to speed the examples, if desired.

config = get_standard_config()
config['max_events'] = 1000

calib.r0_to_dl1(datafile, custom_config=config)
output_filename = 'dl1_' + os.path.basename(datafile).split('.simtel')[0] + '.h5'

# If the file exists, it complains

# To speed up and get better results, you can use directly the DL1 file available in 
#output_filename = '/fefs/aswg/data/mc/DL1/20190415/gamma/south_pointing/20200514_v0.5.1_v01/dl1_gamma_south_pointing_20200514_v0.5.1_v01_DL1_testing.h5'<jupyter_output>USING FILTERS:  Filters(complevel=5, complib='blosc:zstd', shuffle=True, bitshuffle=False, fletcher32=True, least_significant_digit=None)
0
100
200
300
400
500<jupyter_text>We can read this file as a pandas dataframe:<jupyter_code># The results are better if we use the complete DL1 file:

output_filename = '/fefs/aswg/data/mc/DL1/20190415/gamma/south_pointing/20200514_v0.5.1_v01/dl1_gamma_south_pointing_20200514_v0.5.1_v01_DL1_testing.h5'


#Take a look at the parameters extracted from the simtelarray file:

from lstchain.tests.test_lstchain import dl1_file, dl1_params_lstcam_key

data = pd.read_hdf(output_filename, key=dl1_params_lstcam_key)
print(data.keys())
data.head()<jupyter_output><empty_output><jupyter_text>We can take a look at some of these features...  Energy<jupyter_code>plt.rcParams['figure.figsize'] = (20, 15)
plt.rcParams['font.size'] = 14
#We can plot some distributions:

plt.hist(data['log_mc_energy'],bins=100) #Log10(Energy) in MeV
plt.xlabel('$log_{10} mcEnergy$ (MeV)')<jupyter_output><empty_output><jupyter_text>Intensity<jupyter_code>plt.hist(data['log_intensity'],bins=100) #Log10(intensity(size))
plt.xlabel('$log_{10} intensity$')<jupyter_output><empty_output><jupyter_text>Width<jupyter_code>plt.hist(data['width'],bins=100) #Width in deg.
plt.xlabel('width (º)')<jupyter_output>/home/luca.foffano/workspace/software/miniconda3/envs/lst/lib/python3.7/site-packages/numpy/lib/histograms.py:839: RuntimeWarning: invalid value encountered in greater_equal
  keep = (tmp_a >= first_edge)
/home/luca.foffano/workspace/software/miniconda3/envs/lst/lib/python3.7/site-packages/numpy/lib/histograms.py:840: RuntimeWarning: invalid value encountered in less_equal
  keep &= (tmp_a <= last_edge)<jupyter_text>Length<jupyter_code>plt.hist(data['length'],bins=100) #Length in deg. 
plt.xlabel('length (º)')<jupyter_output><empty_output><jupyter_text>Source position in camera coordinates<jupyter_code>plt.hist2d(data['src_x'],data['src_y'],bins=100)
plt.xlabel('x (m)')
plt.ylabel('y (m)')<jupyter_output><empty_output><jupyter_text>Center of gravity of Hillas ellipse in camera coordinates<jupyter_code>plt.hist2d(data['x'],data['y'],bins=100)
plt.xlabel('x (m)')
plt.ylabel('y (m)')<jupyter_output><empty_output><jupyter_text>Import the Random Forests previously trained<jupyter_code>#path_models = "." # if you are in ~/lstchain/notebooks

# Assuming you are running the notebook in ~/cta-lstchain/notebooks
#path_models = "../../cta-lstchain-extra/reco/models"  

# otherwise, the corresponding RF to load
path_models = "/fefs/aswg/data/models/20190415/south_pointing/20200514_v0.5.1_v01"  

fileE = path_models + "/reg_energy.sav"                                                            
fileD = path_models + "/reg_disp_vector.sav"                                                              
fileH = path_models + "/cls_gh.sav"                                                                

RFreg_Energy = joblib.load(fileE)                                                                         
RFreg_Disp = joblib.load(fileD)                                                                           
RFcls_GH = joblib.load(fileH)<jupyter_output><empty_output><jupyter_text>Features for prediction:<jupyter_code># We are going to take them from the standard configuration

config = get_standard_config()
config['regression_features']<jupyter_output><empty_output><jupyter_text>Now we can apply the Random Forests to the data and reconsctruct the energy, direction and Gamma/Hadron classification<jupyter_code>#Now we have the trained RF that we can use to reconstruct our data into dl2 format:

from lstchain.reco.utils import filter_events

data = filter_events(data)


data_dl2 = reco.apply_models(data, 
                             RFcls_GH, 
                             RFreg_Energy, 
                             RFreg_Disp, 
                             config)
print(data_dl2.keys())<jupyter_output>/home/luca.foffano/workspace/software/miniconda3/envs/lst/lib/python3.7/site-packages/ctapipe/coordinates/__init__.py:50: MissingFrameAttributeWarning: AltAz coordinate has no obstime, assuming same frame
  MissingFrameAttributeWarning,
/home/luca.foffano/workspace/software/miniconda3/envs/lst/lib/python3.7/site-packages/ctapipe/coordinates/__init__.py:58: MissingFrameAttributeWarning: Horizontal coordinate has no location, assuming same frame
  MissingFrameAttributeWarning,<jupyter_text>Plot the results     Features:Since we were using a file only with gamma events, here we can notice the events that has been badly classified as protons.<jupyter_code># Let's take a look to the results


lstplot.plot_features(data_dl2)<jupyter_output><empty_output><jupyter_text>Energy reconstruction<jupyter_code>lstplot.energy_results(data_dl2)<jupyter_output><empty_output><jupyter_text>Disp reconstruction<jupyter_code>lstplot.plot_disp(data_dl2)<jupyter_output>mu = 0.02353091096051025
 sigma = 0.18868411823283618<jupyter_text>Reconstructed source position in camera coordinates<jupyter_code>lstplot.plot_pos(data_dl2)<jupyter_output><empty_output><jupyter_text>RHadroness<jupyter_code># var `hadro_rec` has been deprecated. And we have loaded a file only of gammas...
# In any case, the gammaness is 1 -  hadronnes.

plt.hist(1-data_dl2['gammaness'])
plt.xlabel('Hadroness')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10282743755547102, 'fraction_numerical': 0.035121085330290354, 'mean_word_length': 4.606254442075338, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 69, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '423897', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2891, 'n_tokens_neox': 2638, 'n_words': 553, 'path': 'notebooks/lstchain.ipynb', 'repo_name': 'mireianievas/cta-lstchain', 'size': '7887'}