diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
index 5b8c7637e2740fdf36d025f3d57050e28a20fc3c..98e51883a91e8cf48beaea002700fdcbab9dbc45 100644
--- a/modules/icing/pirep_goes.py
+++ b/modules/icing/pirep_goes.py
@@ -28,7 +28,6 @@ dir_fmt = '%Y_%m_%d_%j'
 ds_dct = {}
 goes_ds_dct = {}
 
-pirep_file = homedir+'data/pirep/pireps_20180101_20200331.csv'
 
 # --- CLAVRx Radiometric parameters and metadata ------------------------------------------------
 l1b_ds_list = ['temp_10_4um_nom', 'temp_11_0um_nom', 'temp_12_0um_nom', 'temp_13_3um_nom', 'temp_3_75um_nom',
@@ -85,7 +84,7 @@ train_params_night = ['cld_height_acha', 'cld_geo_thick', 'supercooled_cloud_fra
                       'cld_press_acha', 'cld_reff_acha', 'cld_opd_acha', 'cloud_phase', 'cloud_mask']
 
 
-def setup(pirep_file=pirep_file):
+def setup(pirep_file=homedir+'data/pirep/pireps_20180101_20200331.csv'):
     ice_dict, no_ice_dict, neg_ice_dict = pirep_icing(pirep_file)
     return ice_dict, no_ice_dict, neg_ice_dict
 
@@ -1812,6 +1811,9 @@ def tiles_info(filename):
 def analyze(preds_file, truth_file='/Users/tomrink/data/icing_ml/tiles_202109240000_202111212359_l2_test_v3_DAY.h5'):
     h5f = h5py.File(truth_file, 'r')
     nda = h5f['flight_altitude'][:]
+    cld_hgt = h5f['cld_height_acha'][:]
+    cld_dz = h5f['cld_geo_thick'][:]
+    cld_tmp = h5f['cld_temp_acha'][:]
 
     nda[np.logical_and(nda >= 0, nda < 2000)] = 0
     nda[np.logical_and(nda >= 2000, nda < 4000)] = 1
@@ -1823,9 +1825,32 @@ def analyze(preds_file, truth_file='/Users/tomrink/data/icing_ml/tiles_202109240
 
     print(np.sum(nda == 0), np.sum(nda == 1), np.sum(nda == 2), np.sum(nda == 3), np.sum(nda == 4))
     print('---------------------------')
-    print(np.histogram(nda[labels == 0]))
+
+    print('level 0: ')
+    print(np.nanmean(cld_dz[(nda == 0)]), np.nanmean(cld_hgt[(nda == 0)]))
+    print('------------')
+
+    print('level 1: ')
+    print(np.nanmean(cld_dz[(nda == 1)]), np.nanmean(cld_hgt[(nda == 1)]))
+    print('------------')
+
+    print('level 2: ')
+    print(np.nanmean(cld_dz[(nda == 2)]), np.nanmean(cld_hgt[(nda == 2)]))
+    print('------------')
+
+    print('level 3: ')
+    print(np.nanmean(cld_dz[(nda == 3)]), np.nanmean(cld_hgt[(nda == 3)]))
+    print('------------')
+
+    print('level 4: ')
+    print(np.nanmean(cld_dz[(nda == 4)]), np.nanmean(cld_hgt[(nda == 4)]))
+
+    print('----------------------------')
+    print('-----------------------------')
+
+    print('No icing:       ', np.histogram(nda[labels == 0], bins=5)[0])
     print('---------------------------')
-    print(np.histogram(nda[labels == 1]))
+    print('Icing:          ', np.histogram(nda[labels == 1], bins=5)[0])
     print('-----------------------------')
 
     preds = np.where(prob_avg > 0.5, 1, 0)
@@ -1833,17 +1858,68 @@ def analyze(preds_file, truth_file='/Users/tomrink/data/icing_ml/tiles_202109240
     true_ice = (labels == 1) & (preds == 1)
     false_ice = (labels == 0) & (preds == 1)
 
-    print(np.histogram(nda[true_ice]))
+    print('True icing:     ', np.histogram(nda[true_ice], bins=5)[0])
     print('---------------------------')
-    print(np.histogram(nda[false_ice]))
+    print('False icing:    ', np.histogram(nda[false_ice], bins=5)[0])
     print('---------------------------')
 
     true_no_ice = (labels == 0) & (preds == 0)
     false_no_ice = (labels == 1) & (preds == 0)
 
-    print(np.histogram(nda[true_no_ice]))
+    print('True no icing:  ', np.histogram(nda[true_no_ice], bins=5)[0])
     print('---------------------------')
-    print(np.histogram(nda[false_no_ice]))
+    print('False no icing: ', np.histogram(nda[false_no_ice], bins=5)[0])
+    print('---------------------------')
+
+    print('level 0: ')
+    print(np.nanmean(cld_dz[(nda == 0) & false_no_ice]), np.nanmean(cld_hgt[(nda == 0) & false_no_ice]))
+    print(np.nanmean(cld_dz[(nda == 0) & true_no_ice]), np.nanmean(cld_hgt[(nda == 0) & true_no_ice]))
+    print('------------')
+
+    print('level 1: ')
+    print(np.nanmean(cld_dz[(nda == 1) & false_no_ice]), np.nanmean(cld_hgt[(nda == 1) & false_no_ice]))
+    print(np.nanmean(cld_dz[(nda == 1) & true_no_ice]), np.nanmean(cld_hgt[(nda == 1) & true_no_ice]))
+    print('------------')
+
+    print('level 2: ')
+    print(np.nanmean(cld_dz[(nda == 2) & false_no_ice]), np.nanmean(cld_hgt[(nda == 2) & false_no_ice]))
+    print(np.nanmean(cld_dz[(nda == 2) & true_no_ice]), np.nanmean(cld_hgt[(nda == 2) & true_no_ice]))
+    print('------------')
+
+    print('level 3: ')
+    print(np.nanmean(cld_dz[(nda == 3) & false_no_ice]), np.nanmean(cld_hgt[(nda == 3) & false_no_ice]))
+    print(np.nanmean(cld_dz[(nda == 3) & true_no_ice]), np.nanmean(cld_hgt[(nda == 3) & true_no_ice]))
+    print('------------')
+
+    print('level 4: ')
+    print(np.nanmean(cld_dz[(nda == 4) & false_no_ice]), np.nanmean(cld_hgt[(nda == 4) & false_no_ice]))
+    print(np.nanmean(cld_dz[(nda == 4) & true_no_ice]), np.nanmean(cld_hgt[(nda == 4) & true_no_ice]))
+
+    print('----------------------------')
+
+    print('level 0: ')
+    print(np.nanmean(cld_dz[(nda == 0) & false_ice]), np.nanmean(cld_hgt[(nda == 0) & false_ice]))
+    print(np.nanmean(cld_dz[(nda == 0) & true_ice]), np.nanmean(cld_hgt[(nda == 0) & true_ice]))
+    print('------------')
+
+    print('level 1: ')
+    print(np.nanmean(cld_dz[(nda == 1) & false_ice]), np.nanmean(cld_hgt[(nda == 1) & false_ice]))
+    print(np.nanmean(cld_dz[(nda == 1) & true_ice]), np.nanmean(cld_hgt[(nda == 1) & true_ice]))
+    print('------------')
+
+    print('level 2: ')
+    print(np.nanmean(cld_dz[(nda == 2) & false_ice]), np.nanmean(cld_hgt[(nda == 2) & false_ice]))
+    print(np.nanmean(cld_dz[(nda == 2) & true_ice]), np.nanmean(cld_hgt[(nda == 2) & true_ice]))
+    print('------------')
+
+    print('level 3: ')
+    print(np.nanmean(cld_dz[(nda == 3) & false_ice]), np.nanmean(cld_hgt[(nda == 3) & false_ice]))
+    print(np.nanmean(cld_dz[(nda == 3) & true_ice]), np.nanmean(cld_hgt[(nda == 3) & true_ice]))
+    print('------------')
+
+    print('level 4: ')
+    print(np.nanmean(cld_dz[(nda == 4) & false_ice]), np.nanmean(cld_hgt[(nda == 4) & false_ice]))
+    print(np.nanmean(cld_dz[(nda == 4) & true_ice]), np.nanmean(cld_hgt[(nda == 4) & true_ice]))
 
 
 def get_training_parameters(day_night='DAY', l1b_andor_l2='BOTH'):