<jupyter_start><jupyter_text>__cta-lstchain: Notebook for training and storage Random Forests for LST1 data analysis__    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>import numpy as np                                                                                            
import pandas as pd                                                                                                                                                                                      
import matplotlib.pyplot as plt                                                                               
from sklearn.externals import joblib                                                                               
from lstchain.reco import reco_dl1_to_dl2     
from lstchain.visualization import plot_dl2                                                                   
from lstchain.reco import utils
%matplotlib inline
plt.rcParams['figure.figsize'] = (20, 15)
plt.rcParams['font.size'] = 14<jupyter_output><empty_output><jupyter_text>**Get event DL1 files for training.**We need two files, for **gammas** and for **protons**. These gammas are pointlike.<jupyter_code>try:
    PATH_EVENTS = "/home/queenmab/DATA/LST1/Events"
    gammafile = PATH_EVENTS+"/gamma_events_point.h5"
    protonfile = PATH_EVENTS+"/proton_events.h5"

except:
    PATH_EVENTS = "../../cta-lstchain-extra/reco/sample_data/dl1/"
    gammafile = PATH_EVENTS+"gamma_events_point_tiny.h5"
    protonfile = PATH_EVENTS+"proton_events_tiny.h5"<jupyter_output><empty_output><jupyter_text>We read the files as pandas dataframes:<jupyter_code>df_gammas = pd.read_hdf(gammafile)                                                               
df_proton = pd.read_hdf(protonfile)
df_gammas.keys()
df_gammas<jupyter_output><empty_output><jupyter_text>From all the previous information, we choose certain features to train the Random Forests:<jupyter_code>features = ['intensity',                                                                                  
            'time_gradient',                                                                              
            'width',                                                                                      
            'length',                                                                                     
            'wl',                                                                                         
            'phi',                                                                                        
            'psi']<jupyter_output><empty_output><jupyter_text>Now we must split the data into train and test sets. Gamma events will train energy and direction reconstruction, and gamma/hadron separation, but protons are only used for separation.<jupyter_code>train,test = reco_dl1_to_dl2.split_traintest(df_gammas,0.5)
test = test.append(df_proton,ignore_index=True) #Protons are only for testing when trainin Energy/Direction reco.<jupyter_output>/home/queenmab/anaconda3/envs/cta-dev/lib/python3.6/site-packages/pandas/core/frame.py:6211: FutureWarning: Sorting because non-concatenation axis is not aligned. A future version
of pandas will change to not sort by default.

To accept the future behavior, pass 'sort=False'.

To retain the current behavior and silence the warning, pass 'sort=True'.

  sort=sort)<jupyter_text>**Train the Reconstruction:**We train two Random Forest Regressors, from scikit-learn, to reconstruct energy and "disp" of the **test** set.<jupyter_code>RFreg_Energy, RFreg_Disp = reco_dl1_to_dl2.trainRFreco(train,features)<jupyter_output>Given features:  ['intensity', 'time_gradient', 'width', 'length', 'wl', 'phi', 'psi']
Number of events for training:  266413
Training Random Forest Regressor for Energy Reconstruction...
Random Forest trained!
Training Random Forest Regressor for disp Reconstruction...
Random Forest trained!
Done!<jupyter_text>We can now predict the **energy** and **disp** of the test events, and from **disp**, calculate the reconstructed direction.<jupyter_code>#Predict energy and direction for the test set:
test['e_rec'] = RFreg_Energy.predict(test[features])                                                      
test['disp_rec'] = RFreg_Disp.predict(test[features])                                                     
                                                                                                              
test['src_x_rec'],test['src_y_rec'] = utils.disp_to_pos(test['disp_rec'],                                 
                                                       test['x'],                                              
                                                       test['y'],                                              
                                                       test['psi'])
test.keys()<jupyter_output><empty_output><jupyter_text>We use this test events with reconstructed energy and direction to **train the gamma/hadron separation.**We add these two features to the list of features for training:<jupyter_code>features_sep = list(features)                                                                             
features_sep.append('e_rec')                                                                              
features_sep.append('disp_rec')<jupyter_output><empty_output><jupyter_text>We split again the set into **train** and **test**.<jupyter_code>train,test = reco_dl1_to_dl2.split_traintest(test,0.75)<jupyter_output><empty_output><jupyter_text>**Train the gamma/hadron classifier:**Now we train a scikit-learn **RandomForestClassifier** which will separate events in two classes: 0 for **gammas** and 1 for **protons**. We call this parameter **hadroness**.<jupyter_code>#Now we can train the Random Forest Classifier
RFcls_GH = reco_dl1_to_dl2.trainRFsep(train,features_sep)<jupyter_output>Given features:  ['intensity', 'time_gradient', 'width', 'length', 'wl', 'phi', 'psi', 'e_rec', 'disp_rec']
Number of events for training:  365202
Training Random Forest Classifier for Gamma/Hadron separation...
Random Forest trained!
Done!<jupyter_text>Predict the hadroness of the test events:<jupyter_code>test['hadro_rec'] = RFcls_GH.predict(test[features_sep])<jupyter_output><empty_output><jupyter_text>**Save the Random Forests:**We can save these trained RF into files to apply them later on any set of data:<jupyter_code>fileE = "RFreg_Energy.sav"                                                           
fileD = "RFreg_Disp.sav"                                                             
fileH = "RFcls_GH.sav"                                                               
joblib.dump(RFreg_Energy, fileE)                                                                      
joblib.dump(RFreg_Disp, fileD)                                                                        
joblib.dump(RFcls_GH, fileH)<jupyter_output><empty_output><jupyter_text>**Now we can plot some results...****Distribution of features**<jupyter_code>plot_dl2.plot_features(test)<jupyter_output><empty_output><jupyter_text>**Energy reconstruction**<jupyter_code>plot_dl2.plot_e(test)<jupyter_output>0.03707218517239982 0.6655900815918775<jupyter_text>**Disp reconstruction**<jupyter_code>plot_dl2.plot_disp(test)<jupyter_output>0.06304875206486528 0.23996262045022812<jupyter_text>**Source position in camera coordinates**<jupyter_code>plot_dl2.plot_pos(test)<jupyter_output><empty_output><jupyter_text>**Importance of features for Gamma/Hadron separation**<jupyter_code>plot_dl2.plot_importances(RFcls_GH,features_sep)<jupyter_output>Feature importances (gini index)
1. width (0.218487)
2. length (0.149491)
3. disp_rec (0.138147)
4. intensity (0.107225)
5. time_gradient (0.085866)
6. e_rec (0.084655)
7. wl (0.077336)
8. phi (0.076596)
9. psi (0.062197)<jupyter_text>**ROC curve**<jupyter_code>plot_dl2.plot_ROC(RFcls_GH,test,features_sep,-1)<jupyter_output>0.716945659937442<jupyter_text>Mono analysis is very limited, specially at low energies, where it's much more difficult to separate gammas from hadrons without stereo information. If we discard low energy events, we can see that the performance improves.    **For example, we can cut at 500 GeV:**<jupyter_code>cut = 500.

e_cut = np.log10(cut) 
test_cut = test[test['e_rec']>e_cut]
plot_dl2.plot_e(test_cut)
plot_dl2.plot_disp(test_cut)
plot_dl2.plot_pos(test_cut)
plot_dl2.plot_ROC(RFcls_GH,test_cut,features_sep,e_cut)<jupyter_output>0.9364823609226595
{'fraction_non_alphanumeric': 0.08684974324626032, 'fraction_numerical': 0.027461486939048894, 'mean_word_length': 2.0441726129799522, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 70, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 77, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '5596', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2595, 'n_tokens_neox': 2389, 'n_words': 566, 'path': 'notebooks/RF-performance.ipynb', 'repo_name': 'pawel21/cta-lstchain', 'size': '8958'}