Merge pull request #61 from johbackm/bugfix

Bugfixes: big soot partilces were lost, gaussian fit in the plot was too wide, split position algorithm now more stable.

pysp2/util/particle_properties.py

@@ -61,7 +61,9 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None):
 
     PkHt_ch1 = input_ds['PkHt_ch1'].values
     PkHt_ch5 = input_ds['PkHt_ch5'].values
-    width = input_ds['PkEnd_ch1'].values - input_ds['PkStart_ch1'].values
+    width_ch1 = input_ds['PkEnd_ch1'].values - input_ds['PkStart_ch1'].values
+    width_ch5 = input_ds['PkEnd_ch5'].values - input_ds['PkStart_ch5'].values
+    width = np.where(np.isnan(width_ch1),width_ch5,width_ch1)
     accepted_incand = width >= Globals.IncanMinWidth
     accepted_incand = np.logical_and(accepted_incand,
                                      input_ds['PkHt_ch2'].values >= Globals.IncanMinPeakHt1)

pysp2/util/peak_fit.py

@@ -509,9 +509,10 @@ def _split_scatter_fit(my_ds, channel):
     num_base_pts_2_avg = 20
     data = my_ds['Data_ch' + str(channel)].values
     V_maxloc = np.argmax(data, axis=1)
-    V_minloc = np.argmin(data, axis=1)
-    invert = V_maxloc < V_minloc
-    data[invert, :] = -data[invert, :]
+    #V_minloc = np.argmin(data, axis=1)
+    #invert = V_maxloc < V_minloc
+    if channel==3:
+        data = -data
     base = np.nanmean(data[:, 0:num_base_pts_2_avg], axis=1)
     V_max = data.max(axis=1)
     conditions = np.logical_and.reduce(((V_max - base) > 1, V_maxloc < len(data), V_maxloc > 0))
@@ -526,7 +527,7 @@ def _split_scatter_fit(my_ds, channel):
     pos_tile = np.tile(pos, (data.shape[1], 1)).T
     counting_up = np.where(np.logical_and(data < 5, counting_up <= pos_tile), counting_up, -1)
     start = counting_up.max(axis=1)
-    fit_coeffs = {'base': base, 'height': height, 'pos': pos, 'start': start, 'inverted': invert}
+    fit_coeffs = {'base': base, 'height': height, 'pos': pos, 'start': start}
 
     return fit_coeffs
 

pysp2/vis/plot_wave.py

@@ -71,7 +71,7 @@ def plot_wave(ds, record_no, chn, plot_fit=True,
         pos = spectra['FtPos_ch' + str(chn)].values
         base = spectra['Base_ch' + str(chn)].values
         width = spectra['PkFWHM_ch' + str(chn)].values
-        Y = _gaus(xspace, amplitude, pos, width/2., base)
+        Y = _gaus(xspace, amplitude, pos, width/2.35482, base)
         ax.plot(xspace, Y)
         ax.text(0.7, 0.5, 'Fit Pos = %3.2f' % pos,
                 transform=ax.transAxes)
@@ -140,6 +140,6 @@ def plot_waves(ds, record_no, plot_fit=True):
         ax.legend(legends[i])
         if i==2:
             lines=ax.get_lines()
-            ax.plot(ds['PkSplitPos_ch3'].isel(event_index=record_no),0,'*',markersize=10,color=lines[0].get_color())
-            ax.plot(ds['PkSplitPos_ch7'].isel(event_index=record_no),0,'*',markersize=10,color=lines[1].get_color())
+            ax.plot(ds['PkSplitPos_ch3'].isel(event_index=record_no),0,'1',markersize=10,color=lines[0].get_color())
+            ax.plot(ds['PkSplitPos_ch7'].isel(event_index=record_no),0,'2',markersize=10,color=lines[1].get_color())
     return display

tests/test_gfit.py

@@ -11,7 +11,7 @@ def test_gaussian_fit():
         np.nanmax(my_binary.PkHt_ch0.values), 98708.92915295, decimal=1)
     np.testing.assert_almost_equal(
         np.nanmax(my_binary.PkHt_ch4.values), 65088.3959945008, decimal=1)
-    # check that there are requal amounts of successful fits for low gain and
+    # check that there are equal amounts of successful fits for low gain and
     # high gain scattering when the peak heighs are large enough
     bl_hg = np.logical_and(my_binary['FtAmp_ch0'] > 30000,
                            my_binary['FtAmp_ch0'] < 50000)
@@ -30,9 +30,9 @@ def test_psds():
     my_psds = pysp2.util.process_psds(my_binary, my_hk, my_ini)
     np.testing.assert_almost_equal(my_psds['NumConcIncan'].max(), 0.95805343)
     np.testing.assert_almost_equal(my_psds['ScatNumEnsemble'].sum(), 254.773995310)
-    np.testing.assert_almost_equal(my_psds['IncanNumEnsemble'].sum(), 32.22939087)
+    np.testing.assert_almost_equal(my_psds['IncanNumEnsemble'].sum(), 32.61006871)
     np.testing.assert_almost_equal(my_psds['ScatMassEnsemble'].sum(), 3.15026266)
-    np.testing.assert_almost_equal(my_psds['IncanMassEnsemble'].sum(), 0.08177226)
+    np.testing.assert_almost_equal(my_psds['IncanMassEnsemble'].sum(), 0.08280955)
     np.testing.assert_almost_equal(my_binary['DeadtimeRelativeBias'].mean(), -0.00023515)
     coeff, beam_profile = pysp2.util.beam_shape(
         my_binary, beam_position_from='peak maximum', Globals=pysp2.util.DMTGlobals())