<jupyter_start><jupyter_text>__Notebook for testing the reconstrcution of LST1 data.__This 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**If you have ctapipe already installed in a conda environment:**source activate cta-devpython setup.py install**If you don't have ctapipe installed:**conda env create -f environment.ymlsource activate cta-devpython setup.py installAlso, you will need the datafiles from **cta-lstchain-extra:**git clone https://github.com/misabelber/cta-lstchain-extra     **Some imports...**<jupyter_code>from sklearn.externals import joblib
from ctapipe.utils import get_dataset_path
import matplotlib.pyplot as plt
import os
import pandas as pd
import lstchain.reco.dl0_to_dl1 as calib
import lstchain.reco.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><empty_output><jupyter_text>**Import a simtelarray file of MC simulated data:**<jupyter_code>try:
    datafile = "/home/queenmab/DATA/LST1/Gamma/gamma_20deg_0deg_run8___cta-prod3-lapalma-2147m-LaPalma-FlashCam.simtel.gz"
except:
    datafile = get_dataset_path('gamma_test_large.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>calib.max_events = 1000 #Limit the number of events to analyze, to speed the examples
calib.r0_to_dl1(datafile)
output_filename = 'dl1_' + os.path.basename(datafile).split('.')[0] + '.h5'<jupyter_output>0
100
200
300<jupyter_text>We can read this file as a **pandas dataframe**:<jupyter_code>#Take a look at the parameters extracted from the simtelarray file:
data = pd.read_hdf(output_filename,key='events/LSTCam')
print(data.keys())
data<jupyter_output><empty_output><jupyter_text>**We can take a look at some of these features...**        **Energy**<jupyter_code>#We can plot some distributions:
plt.hist(data['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['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><empty_output><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>**Cengter 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>from sklearn.externals import joblib

try:
    path_models = "/home/queenmab/DATA/LST1/Models"

    fileE = path_models + "/RFreg_Energy.sav"                                                            
    fileD = path_models + "/RFreg_Disp.sav"                                                              
    fileH = path_models + "/RFcls_GH.sav"                                                                
                                                                                                              
    RFreg_Energy = joblib.load(fileE)                                                                         
    RFreg_Disp = joblib.load(fileD)                                                                           
    RFcls_GH = joblib.load(fileH)
except:
    path_models = "../../cta-lstchain-extra/reco/models"
    
    fileE = path_models + "/RFreg_Energy.sav"                                                            
    fileD = path_models + "/RFreg_Disp.sav"                                                              
    fileH = path_models + "/RFcls_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>#Define the features to train the RF:
features = ['intensity',                                                                                  
                'time_gradient',                                                                              
                'width',                                                                                      
                'length',                                                                                     
                'wl',                                                                                        
                'phi',                                                                                        
                'psi']<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:
data_dl2 = reco.ApplyModels(data,features,RFcls_GH,RFreg_Energy,RFreg_Disp)
print(data_dl2.keys())<jupyter_output>Index(['disp', 'event_id', 'gps_time', 'hadroness', 'intensity', 'intercept',
       'kurtosis', 'length', 'mc_alt', 'mc_alt_tel', 'mc_az', 'mc_az_tel',
       'mc_core_distance', 'mc_core_x', 'mc_core_y', 'mc_energy',
       'mc_h_first_int', 'mc_type', 'mc_x_max', 'obs_id', 'phi', 'psi', 'r',
       'skewness', 'src_x', 'src_y', 'time_gradient', 'width', 'wl', 'x', 'y',
       'e_rec', 'disp_rec', 'src_x_rec', 'src_y_rec', 'hadro_rec'],
      dtype='object')<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>#Lets take a look to the results
lstplot.plot_features(data_dl2)<jupyter_output><empty_output><jupyter_text>**Energy reconstruction**<jupyter_code>lstplot.plot_e(data_dl2)<jupyter_output>0.09856435356902216 0.6137690435379151<jupyter_text>**Disp reconstruction**<jupyter_code>lstplot.plot_disp(data_dl2)<jupyter_output>0.08472789618800787 0.24503858392238181<jupyter_text>**Reconstructed source position in camera coordinates**<jupyter_code>lstplot.plot_pos(data_dl2)
<font size="4">
**RHadroness**
plt.hist(data_dl2['hadro_rec'])
plt.xlabel('Hadroness')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1063858512218172, 'fraction_numerical': 0.020521128662076414, 'mean_word_length': 2.2778761061946904, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 66, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 52, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '410067', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2328, 'n_tokens_neox': 2139, 'n_words': 441, 'path': 'notebooks/lstchain.ipynb', 'repo_name': 'pawel21/cta-lstchain', 'size': '7407'}