From ce508606a26d9aa48fd8ccc73e914578a05fc869 Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 3 Apr 2024 19:17:07 +0300 Subject: [PATCH 01/43] added leo_fit() from .leo_fit --- pysp2/util/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysp2/util/__init__.py b/pysp2/util/__init__.py index af8f912..68bba37 100644 --- a/pysp2/util/__init__.py +++ b/pysp2/util/__init__.py @@ -17,4 +17,4 @@ from .DMTGlobals import DMTGlobals from .particle_properties import calc_diams_masses, process_psds from .deadtime import deadtime -from .leo_fit import beam_shape \ No newline at end of file +from .leo_fit import beam_shape,leo_fit From b23968bc0282b45a7a6f0e76042bfa889a25d1a8 Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 4 Apr 2024 11:58:37 +0300 Subject: [PATCH 02/43] unfinnished leo fit code. Not working yet. --- pysp2/util/leo_fit.py | 223 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 181 insertions(+), 42 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index fee4396..6637316 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -1,9 +1,10 @@ import numpy as np -from .peak_fit import _gaus from scipy.optimize import curve_fit +import xarray as xr +from .peak_fit import _gaus, _do_fit_records +#from pysp2.util.peak_fit import _do_fit_records - -def beam_shape(my_binary, beam_position_from='peak maximum', Globals=None): +def beam_shape(my_binary, beam_position_from='split point', Globals=None): """ Calculates the beam shape needed to determine the laser intensity profile @@ -26,7 +27,7 @@ def beam_shape(my_binary, beam_position_from='peak maximum', Globals=None): 'peak maximum' = construct the beam profile arround the maximum peak poistion. The maximum peak position is determied from the peak-height weighted average peak position. - 'split detector' = construct the beam profile around the split position. + 'split point' = construct the beam profile around the split position. The split position is taken from the split detector. Not working yet. Returns @@ -39,26 +40,36 @@ def beam_shape(my_binary, beam_position_from='peak maximum', Globals=None): """ - num_base_pts_2_avg = 5 - bins=my_binary['columns'] + num_base_pts_2_avg = 10 + moving_average_window = 5 + max_amplitude_fraction = 0.03 + bins = my_binary['columns'] + + median_peak_width = my_binary['PkFWHM_ch0'].median().values / 2.35482 + median_peak_pos = my_binary['FtPos_ch0'].median().values + peak_pos_ok = np.logical_and(my_binary['FtPos_ch0'].values >= median_peak_pos - median_peak_width, + my_binary['FtPos_ch0'].values <= median_peak_pos + median_peak_width) - #boolean array for ok particles in the high gain channel + #change this so that it uses scat_reject_key instead. + #scatterin signal ok = True scatter_high_gain_accepted = np.logical_and.reduce(( my_binary['PkFWHM_ch0'].values > Globals.ScatMinWidth, my_binary['PkFWHM_ch0'].values < Globals.ScatMaxWidth, my_binary['PkHt_ch0'].values > Globals.ScatMinPeakHt1, my_binary['PkHt_ch0'].values < Globals.ScatMaxPeakHt1, - my_binary['FtPos_ch0'].values < Globals.ScatMaxPeakPos, - my_binary['FtPos_ch0'].values > Globals.ScatMinPeakPos)) - - #boolean array that is True for particles that have not been triggered by - #the high gain incandesence channel + peak_pos_ok)) +# my_binary['FtPos_ch0'].values < Globals.ScatMaxPeakPos, +# my_binary['FtPos_ch0'].values > Globals.ScatMinPeakPos)) + + #no incandesence signal = True no_incand_trigged = my_binary['PkHt_ch1'].values < Globals.IncanMinPeakHt1 - #find good events that only scatter light + #Particles that only scatter light only_scattering_high_gain = np.logical_and.reduce((scatter_high_gain_accepted, no_incand_trigged)) - + #iloc = "event_index" and "index" with the scattering only particle events + iloc = np.argwhere(only_scattering_high_gain).flatten() + print('High gain scattering particles for beam analysis :: ', np.sum(only_scattering_high_gain)) @@ -71,25 +82,39 @@ def beam_shape(my_binary, beam_position_from='peak maximum', Globals=None): my_high_gain_scatterers.dims['columns'])) \ * np.nan - #weighted mean of beam peak position. Weight is scattering amplitude. - high_gain_peak_pos = int( - np.sum(np.multiply(my_high_gain_scatterers['PkPos_ch0'].values, - my_high_gain_scatterers['PkHt_ch0'].values))/ \ - np.sum(my_high_gain_scatterers['PkHt_ch0'].values)) + my_high_gain_profiles_ = np.zeros_like(my_high_gain_profiles)*np.nan + + mean_high_gain_max_peak_pos = int(np.nanmean(my_high_gain_scatterers['PkPos_ch0'].values)) + mean_high_gain_split_pos = int(np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values)) + + #cross to center + c2c = my_high_gain_scatterers['FtPos_ch0'].values - my_high_gain_scatterers['PkSplitPos_ch3'].values + mean_c2c = np.mean(c2c) #mean cross to centre + + high_gain_split_position = my_high_gain_scatterers['PkSplitPos_ch3'].values + + mean_split_pos = np.round(np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values)).astype(np.int32) - #loop through all particle events + #loop through all particle events (THIS CAN BE MADE SMARTER WITH MUCH OF THE CALCULATIONS BEFORE THE LOOP) for i in my_high_gain_scatterers['event_index']: data = my_high_gain_scatterers['Data_ch0'].sel(event_index=i).values #base level base = np.mean(data[0:num_base_pts_2_avg]) #peak height - peak_height = data.max()-base - #peak position + peak_height = data.max() - base + #max peak position peak_pos = data.argmax() - #normalize the profile to range [0,1] + #split position + split_pos = my_high_gain_scatterers['PkSplitPos_ch3'].sel(event_index=i).values + #normalize the profile to range [~0,1] profile = (data - base) / peak_height + #insert the profile as it was recorded (no shifting due to PEAK POSITION or PSD POSITION) + my_high_gain_profiles_[i,:] = profile #distance to the mean beam peak position - peak_difference = high_gain_peak_pos - peak_pos + if beam_position_from == 'peak maximum': + peak_difference = mean_high_gain_max_peak_pos - peak_pos + elif beam_position_from == 'split point': + peak_difference = mean_split_pos - split_pos #insert so that the peak is at the right position (accounts for #particles travelling at different speeds) if peak_difference > 0: @@ -100,21 +125,135 @@ def beam_shape(my_binary, beam_position_from='peak maximum', Globals=None): else: my_high_gain_profiles[i, :] = profile - #get the beam profile - beam_profile = np.nanmean(my_high_gain_profiles, axis=0) - #find values that are lower than 5% of the max value. - low_values = np.argwhere(beam_profile < 0.05) - #fit the gaussian curve to the beginning of the profile only. The tail - #can deviate from zero substantially and is not of interest. - fit_to = low_values[low_values > high_gain_peak_pos].min() - - #initial guess - p0 = np.array([beam_profile.max() - beam_profile.min(), - np.argmax(beam_profile), 20., - np.nanmin(beam_profile)]).astype(float) - #fit gaussian curve - #coeff[amplitude, peakpos, width , baseline] - coeff, var_matrix = curve_fit(_gaus, bins[:fit_to], beam_profile[:fit_to], - p0=p0, method='lm', maxfev=40, ftol=1e-3) - - return coeff, beam_profile \ No newline at end of file + # MOVING AVERAGE OF THE BEAM PROFILE TO FIND THE DISTANCE BETWEEN THE SPLIT POINT AND THE POINT IN THE + #LASER BEAM WHERE PARTICLES CAN START TO EVAPORATE + + #moving average of the beam shape with a window of moving_average_window + moving_high_gain_profile_window = np.lib.stride_tricks.sliding_window_view(my_high_gain_profiles, + moving_average_window, axis=0) + moving_avg_high_gain_profiles_ = np.nanmean(moving_high_gain_profile_window,axis=2) + + moving_avg_high_gain_profiles = np.zeros_like(moving_avg_high_gain_profiles_) * np.nan + moving_avg_high_gain_split_to_leo_pos = np.zeros(moving_avg_high_gain_profiles_.shape[0]) * np.nan + moving_avg_high_gain_max_leo_pos = np.zeros(moving_avg_high_gain_profiles_.shape[0]) * np.nan + + moving_avg_high_gain_max_leo_amplitude_factor = np.zeros(moving_avg_high_gain_profiles_.shape[0]) * np.nan + + for i in range(moving_avg_high_gain_profiles_.shape[0]): + i_profile = moving_avg_high_gain_profiles_[i,:] + i_max = np.nanargmax(i_profile) + i_range = i_profile[i_max] - np.nanmin(i_profile[:i_max]) + moving_avg_high_gain_profiles[i,:] = (i_profile - np.nanmin(i_profile[:i_max])) / i_range + #interpolate here to get the exact position in fraction (not integer) :: which posiiton (float) is the 0.03 cross in + moving_avg_high_gain_max_leo_pos[i] = np.argwhere(moving_avg_high_gain_profiles[i,:] >= max_amplitude_fraction).min()-1 + moving_avg_high_gain_split_to_leo_pos[i] = moving_avg_high_gain_max_leo_pos[i] - high_gain_split_position[i] + moving_avg_high_gain_max_leo_amplitude_factor[i] = 1./ moving_avg_high_gain_profiles[i, np.round(moving_avg_high_gain_max_leo_pos[i]).astype(int)] + + #cleaning up + moving_avg_high_gain_max_leo_amplitude_factor = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, + np.nan, moving_avg_high_gain_max_leo_amplitude_factor) + + moving_avg_high_gain_max_leo_pos = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, + np.nan, moving_avg_high_gain_max_leo_pos) + + + #moving average of beam width + moving_high_gain_beam_width = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['PkFWHM_ch0'].values, + moving_average_window, axis=0) + moving_avg_high_gain_beam_width = np.nanmean(moving_high_gain_beam_width,axis=1) + + + leo_FtMaxPosAmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan #will later be converted to int + leo_FtMaxPosAmpFactor_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor + + leo_PkFWHM_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan + leo_PkFWHM_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_beam_width + + leo_PkPos_ch0 = np.zeros(scatter_high_gain_accepted.shape) * np.nan + leo_PkPos_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_beam_peakpos + + output_ds = my_binary.copy() + output_ds['leo_FtMaxPos_ch0'] = (('event_index'), leo_FtMaxPos_ch0) + output_ds['leo_FtMaxPosAmpFactor_ch0'] = (('event_index'), leo_FtMaxPosAmpFactor_ch0) + output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) + output_ds['leo_PkPos_ch0'] = (('event_index'), leo_PkPos_ch0) + + output_ds['leo_FtMaxPos_ch0'] = output_ds['leo_FtMaxPos_ch0'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") + output_ds['leo_FtMaxPosAmpFactor_ch0'] = output_ds['leo_FtMaxPosAmpFactor_ch0'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") + output_ds['leo_PkFWHM_ch0'] = output_ds['leo_PkFWHM_ch0'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") + output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") + + #c2c time series (data from scattering particles but extrapolted to all paricles) + #leo_PkPos_ch0 (calculated from c2c) --> position of the max amplitude calculated from the split and extrapolated c2c + #leo_FtMaxPos_ch0 (from split) + + return output_ds + +def leo_fit(my_binary,Globals=None): + bins = my_binary['columns'].astype('float').values + #number of points at the beginning to use for base line average + num_base_pts_2_avg = 10 + #max_amplitude_fraction = 0.03 + + # bl_scattering_ok = my_binary['ScatRejectKey'].values == 0 + # bl_only_scattering_particles = np.logical_and(bl_scattering_ok, + # my_binary['PkHt_ch1'].values < Globals.IncanMinPeakHt1) + + # #Particles that only scatter light and ch0 not saturated + # bl_only_scattering_particles_ch0 = np.logical_and(my_binary['PkHt_ch0'].values < Globals.ScatMaxPeakHt1, + # bl_only_scattering_particles) + + #split to peak height difference (in bins) for scattering only particles + split_to_peak_high_gain = my_binary['PkPos_ch0'].values - my_binary['PkSplitPos_ch3'].values + #For particles with inandesence signal, set to NaN since the peak needn't + #be where the laser intensity is the highest, so se to NaN + #split_to_peak_high_gain[~bl_only_scattering_particles_ch0] = np.nan + + #get the information about the gaussian fits + pos = my_binary['leo_PkPos_ch0'].values + #amplitude = my_binary['PkHt_ch0'].values + width = my_binary['leo_PkFWHM_ch0'].values / 2.35482 #2*sqrt(2*log(2)) + data_ch0 = my_binary['Data_ch0'].values + + #mean of the first num_base_pts_2_avg points + #leo_base_ch0 = np.mean(data_ch0[:, 0:num_base_pts_2_avg], axis=1) + #mean of the lowest 3 points + leo_base_ch0 = np.mean(data_ch0[:, 0:num_base_pts_2_avg], axis=1) + + leo_fit_max_pos = my_binary['leo_FtMaxPos_ch0'].astype(int).values + leo_FtMaxPosAmpFactor_ch0 = my_binary['leo_FtMaxPosAmpFactor_ch0'].values + leo_PkHt_ch0 = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan + leo_PkHt_ch0_ = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan + + for i in range(my_binary.dims['event_index']): + #NAGOYA STYLE + fractional_peak_height_ch0 = data_ch0[i, leo_fit_max_pos[i]] - leo_base_ch0[i] + estimated_peak_height_ch0 = fractional_peak_height_ch0 * leo_FtMaxPosAmpFactor_ch0[i] + leo_PkHt_ch0_[i] = estimated_peak_height_ch0 + max_value = data_ch0[i,:].max() - data_ch0[i,:].min() + bins_ = bins[:leo_fit_max_pos[i]] + #signals + data_ch0_ = data_ch0[i, :leo_fit_max_pos[i]] + leo_coeff, var_matrix = curve_fit( + lambda x, a: _gaus(x, a, pos[i], width[i], leo_base_ch0[i]), + bins_[:], data_ch0_[:], p0=[max_value], maxfev=40, + ftol=1e-3, method='lm' ) #, bounds=(0, 1e6)) #, method='lm' + leo_PkHt_ch0[i] = leo_coeff[0] + + output_ds = my_binary.copy() + output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0) + output_ds['leo_FtAmp_ch0_'] = (('index'), leo_PkHt_ch0_) + output_ds['leo_FtMaxPos_ch0'] = (('index'), leo_fit_max_pos) + output_ds['leo_Base_ch0'] = (('index'), leo_base_ch0) + + #my_high_gain_scatterers['FtAmp_ch0'].plot(marker='.',linestyle='') + #my_high_gain_scatterers['leo_FtAmp_ch0_'].plot(marker='+',linestyle='') + #my_high_gain_scatterers['leo_FtAmp_ch0'].plot(marker='x',linestyle='') + + return output_ds + +#my_binary = pysp2.util.beam_shape(my_binary, beam_position_from='peak maximum', Globals=global_settings) From 2b219123e06fcdac0b8a5456974e6da48f510095 Mon Sep 17 00:00:00 2001 From: John Backman Date: Sat, 1 Feb 2025 19:32:55 +0200 Subject: [PATCH 03/43] leo_PkPos_ch0 for all particles using cross-to-centre distance from scattering only particles --- pysp2/util/leo_fit.py | 56 +++++++++++++++++++++++++++---------------- 1 file changed, 36 insertions(+), 20 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 6637316..29c9930 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -45,10 +45,12 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): max_amplitude_fraction = 0.03 bins = my_binary['columns'] - median_peak_width = my_binary['PkFWHM_ch0'].median().values / 2.35482 - median_peak_pos = my_binary['FtPos_ch0'].median().values - peak_pos_ok = np.logical_and(my_binary['FtPos_ch0'].values >= median_peak_pos - median_peak_width, - my_binary['FtPos_ch0'].values <= median_peak_pos + median_peak_width) + # median_peak_width = my_binary['PkFWHM_ch0'].median().values / 2.35482 + # median_peak_pos = my_binary['FtPos_ch0'].median().values + # peak_pos_ok = np.logical_and(my_binary['FtPos_ch0'].values >= median_peak_pos - median_peak_width, + # my_binary['FtPos_ch0'].values <= median_peak_pos + median_peak_width) + + #change this so that it uses scat_reject_key instead. #scatterin signal ok = True @@ -57,12 +59,13 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): my_binary['PkFWHM_ch0'].values < Globals.ScatMaxWidth, my_binary['PkHt_ch0'].values > Globals.ScatMinPeakHt1, my_binary['PkHt_ch0'].values < Globals.ScatMaxPeakHt1, - peak_pos_ok)) -# my_binary['FtPos_ch0'].values < Globals.ScatMaxPeakPos, -# my_binary['FtPos_ch0'].values > Globals.ScatMinPeakPos)) + my_binary['FtPos_ch0'].values < Globals.ScatMaxPeakPos, + my_binary['FtPos_ch0'].values > Globals.ScatMinPeakPos)) #no incandesence signal = True - no_incand_trigged = my_binary['PkHt_ch1'].values < Globals.IncanMinPeakHt1 + no_incand_trigged = np.logical_and( + my_binary['PkHt_ch1'].values < Globals.IncanMinPeakHt1, + my_binary['PkHt_ch5'].values < Globals.IncanMinPeakHt2) #Particles that only scatter light only_scattering_high_gain = np.logical_and.reduce((scatter_high_gain_accepted, @@ -85,16 +88,15 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): my_high_gain_profiles_ = np.zeros_like(my_high_gain_profiles)*np.nan mean_high_gain_max_peak_pos = int(np.nanmean(my_high_gain_scatterers['PkPos_ch0'].values)) - mean_high_gain_split_pos = int(np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values)) + mean_high_gain_split_pos = np.round(np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values)).astype(np.int32) + #mean_high_gain_split_pos = int(np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values)) #cross to center - c2c = my_high_gain_scatterers['FtPos_ch0'].values - my_high_gain_scatterers['PkSplitPos_ch3'].values - mean_c2c = np.mean(c2c) #mean cross to centre + high_gain_c2c = my_high_gain_scatterers['FtPos_ch0'].values - my_high_gain_scatterers['PkSplitPos_ch3'].values + high_gain_mean_c2c = np.mean(high_gain_c2c) #mean cross to centre high_gain_split_position = my_high_gain_scatterers['PkSplitPos_ch3'].values - mean_split_pos = np.round(np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values)).astype(np.int32) - #loop through all particle events (THIS CAN BE MADE SMARTER WITH MUCH OF THE CALCULATIONS BEFORE THE LOOP) for i in my_high_gain_scatterers['event_index']: data = my_high_gain_scatterers['Data_ch0'].sel(event_index=i).values @@ -114,7 +116,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): if beam_position_from == 'peak maximum': peak_difference = mean_high_gain_max_peak_pos - peak_pos elif beam_position_from == 'split point': - peak_difference = mean_split_pos - split_pos + peak_difference = mean_high_gain_split_pos - split_pos #insert so that the peak is at the right position (accounts for #particles travelling at different speeds) if peak_difference > 0: @@ -160,23 +162,37 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #moving average of beam width moving_high_gain_beam_width = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['PkFWHM_ch0'].values, moving_average_window, axis=0) - moving_avg_high_gain_beam_width = np.nanmean(moving_high_gain_beam_width,axis=1) + #moving_avg_high_gain_beam_width = np.nanmedian(moving_high_gain_beam_width,axis=1) + moving_median_high_gain_beam_width = np.nanmedian(moving_high_gain_beam_width,axis=1) + + moving_high_gain_c2c = np.lib.stride_tricks.sliding_window_view(high_gain_c2c, + moving_average_window, axis=0) + moving_median_high_gain_c2c = np.nanmedian(moving_high_gain_c2c,axis=1) leo_FtMaxPosAmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan #will later be converted to int leo_FtMaxPosAmpFactor_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor leo_PkFWHM_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - leo_PkFWHM_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_beam_width + leo_PkFWHM_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_beam_width - leo_PkPos_ch0 = np.zeros(scatter_high_gain_accepted.shape) * np.nan - leo_PkPos_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_beam_peakpos + leo_c2c_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan + leo_c2c_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_c2c + + #leo_PkPos_ch0 = np.zeros(scatter_high_gain_accepted.shape) * np.nan + #leo_PkPos_ch0[iloc[:-moving_average_window+1]] = output_ds = my_binary.copy() - output_ds['leo_FtMaxPos_ch0'] = (('event_index'), leo_FtMaxPos_ch0) + #output_ds['leo_FtMaxPos_ch0'] = (('event_index'), leo_FtMaxPos_ch0) output_ds['leo_FtMaxPosAmpFactor_ch0'] = (('event_index'), leo_FtMaxPosAmpFactor_ch0) output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) - output_ds['leo_PkPos_ch0'] = (('event_index'), leo_PkPos_ch0) + #output_ds['leo_PkPos_ch0'] = (('event_index'), leo_PkPos_ch0) + + # ADD HERE: CALCULATE leo_PkPos_ch0 from split_point + moving_median_high_gain_c2c for all particles + output_ds['leo_PkPos_ch0'] = (('event_index'), leo_c2c_ch0) + output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") + output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'] + my_binary['PkSplitPos_ch3'] output_ds['leo_FtMaxPos_ch0'] = output_ds['leo_FtMaxPos_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") From fd3b4d8e99d18e2992ac2b36b8dfd1a39252deb4 Mon Sep 17 00:00:00 2001 From: John Backman Date: Sat, 1 Feb 2025 20:52:57 +0200 Subject: [PATCH 04/43] small edits and comments --- pysp2/util/leo_fit.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 29c9930..3c6d2d4 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -179,14 +179,11 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): leo_c2c_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_c2c_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_c2c - #leo_PkPos_ch0 = np.zeros(scatter_high_gain_accepted.shape) * np.nan - #leo_PkPos_ch0[iloc[:-moving_average_window+1]] = output_ds = my_binary.copy() - #output_ds['leo_FtMaxPos_ch0'] = (('event_index'), leo_FtMaxPos_ch0) + output_ds['leo_FtMaxPos_ch0'] = (('event_index'), leo_FtMaxPos_ch0) output_ds['leo_FtMaxPosAmpFactor_ch0'] = (('event_index'), leo_FtMaxPosAmpFactor_ch0) output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) - #output_ds['leo_PkPos_ch0'] = (('event_index'), leo_PkPos_ch0) # ADD HERE: CALCULATE leo_PkPos_ch0 from split_point + moving_median_high_gain_c2c for all particles output_ds['leo_PkPos_ch0'] = (('event_index'), leo_c2c_ch0) @@ -202,9 +199,9 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): method="nearest", fill_value="extrapolate") output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") - - #c2c time series (data from scattering particles but extrapolted to all paricles) - #leo_PkPos_ch0 (calculated from c2c) --> position of the max amplitude calculated from the split and extrapolated c2c + #ToDo: + #done: c2c time series (data from scattering particles but extrapolted to all paricles) + #done: leo_PkPos_ch0 (calculated from c2c) --> position of the max amplitude calculated from the split and extrapolated c2c #leo_FtMaxPos_ch0 (from split) return output_ds @@ -224,7 +221,7 @@ def leo_fit(my_binary,Globals=None): # bl_only_scattering_particles) #split to peak height difference (in bins) for scattering only particles - split_to_peak_high_gain = my_binary['PkPos_ch0'].values - my_binary['PkSplitPos_ch3'].values + #split_to_peak_high_gain = my_binary['PkPos_ch0'].values - my_binary['PkSplitPos_ch3'].values #For particles with inandesence signal, set to NaN since the peak needn't #be where the laser intensity is the highest, so se to NaN #split_to_peak_high_gain[~bl_only_scattering_particles_ch0] = np.nan From 7999b2254013fd24e95563f59368d6ca2d78e25f Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 15 May 2025 16:43:40 +0300 Subject: [PATCH 05/43] added leo_PkFWHM, leo_MaxPos_ch0, leo_PkPos_ch0, etc. --- pysp2/util/leo_fit.py | 72 +++++++++++++++++++++++++++++++++---------- 1 file changed, 55 insertions(+), 17 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 3c6d2d4..c44b944 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -51,7 +51,6 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): # my_binary['FtPos_ch0'].values <= median_peak_pos + median_peak_width) - #change this so that it uses scat_reject_key instead. #scatterin signal ok = True scatter_high_gain_accepted = np.logical_and.reduce(( @@ -88,13 +87,13 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): my_high_gain_profiles_ = np.zeros_like(my_high_gain_profiles)*np.nan mean_high_gain_max_peak_pos = int(np.nanmean(my_high_gain_scatterers['PkPos_ch0'].values)) - mean_high_gain_split_pos = np.round(np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values)).astype(np.int32) + mean_high_gain_split_pos_float = np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values) + mean_high_gain_split_pos = np.round(mean_high_gain_split_pos_float).astype(np.int32) #mean_high_gain_split_pos = int(np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values)) #cross to center high_gain_c2c = my_high_gain_scatterers['FtPos_ch0'].values - my_high_gain_scatterers['PkSplitPos_ch3'].values - high_gain_mean_c2c = np.mean(high_gain_c2c) #mean cross to centre - + high_gain_split_position = my_high_gain_scatterers['PkSplitPos_ch3'].values #loop through all particle events (THIS CAN BE MADE SMARTER WITH MUCH OF THE CALCULATIONS BEFORE THE LOOP) @@ -127,13 +126,20 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): else: my_high_gain_profiles[i, :] = profile - # MOVING AVERAGE OF THE BEAM PROFILE TO FIND THE DISTANCE BETWEEN THE SPLIT POINT AND THE POINT IN THE + #MOVING AVERAGE OF THE BEAM PROFILE TO FIND THE DISTANCE BETWEEN THE SPLIT POINT AND THE POINT IN THE #LASER BEAM WHERE PARTICLES CAN START TO EVAPORATE #moving average of the beam shape with a window of moving_average_window moving_high_gain_profile_window = np.lib.stride_tricks.sliding_window_view(my_high_gain_profiles, moving_average_window, axis=0) - moving_avg_high_gain_profiles_ = np.nanmean(moving_high_gain_profile_window,axis=2) +# moving_avg_high_gain_profiles_ = np.nanmean(moving_high_gain_profile_window,axis=2) + moving_avg_high_gain_profiles_ = np.nanmedian(moving_high_gain_profile_window,axis=2) + + # # + # moving_high_gain_FtPos = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['FtPos_ch0'].values, + # moving_average_window, axis=0) + # moving_high_gain_FtPos_ = np.nanmedian(moving_high_gain_FtPos,axis=1) + moving_avg_high_gain_profiles = np.zeros_like(moving_avg_high_gain_profiles_) * np.nan moving_avg_high_gain_split_to_leo_pos = np.zeros(moving_avg_high_gain_profiles_.shape[0]) * np.nan @@ -148,7 +154,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): moving_avg_high_gain_profiles[i,:] = (i_profile - np.nanmin(i_profile[:i_max])) / i_range #interpolate here to get the exact position in fraction (not integer) :: which posiiton (float) is the 0.03 cross in moving_avg_high_gain_max_leo_pos[i] = np.argwhere(moving_avg_high_gain_profiles[i,:] >= max_amplitude_fraction).min()-1 - moving_avg_high_gain_split_to_leo_pos[i] = moving_avg_high_gain_max_leo_pos[i] - high_gain_split_position[i] + moving_avg_high_gain_split_to_leo_pos[i] = moving_avg_high_gain_max_leo_pos[i] - mean_high_gain_split_pos moving_avg_high_gain_max_leo_amplitude_factor[i] = 1./ moving_avg_high_gain_profiles[i, np.round(moving_avg_high_gain_max_leo_pos[i]).astype(int)] #cleaning up @@ -158,6 +164,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): moving_avg_high_gain_max_leo_pos = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, np.nan, moving_avg_high_gain_max_leo_pos) + moving_avg_high_gain_split_to_leo_pos = np.where(moving_avg_high_gain_split_to_leo_pos < -20. , + np.nan, moving_avg_high_gain_split_to_leo_pos) #moving average of beam width moving_high_gain_beam_width = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['PkFWHM_ch0'].values, @@ -167,42 +175,72 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): moving_high_gain_c2c = np.lib.stride_tricks.sliding_window_view(high_gain_c2c, moving_average_window, axis=0) + #JB OK moving_median_high_gain_c2c = np.nanmedian(moving_high_gain_c2c,axis=1) - - leo_FtMaxPosAmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan #will later be converted to int - leo_FtMaxPosAmpFactor_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor + + #JB LATER + #leo_FtMaxPosAmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan #will later be converted to int + #leo_FtMaxPosAmpFactor_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor + #JB OK leo_PkFWHM_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_PkFWHM_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_beam_width + #JB OK leo_c2c_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_c2c_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_c2c + #JB OK + leo_split_to_leo_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan + leo_split_to_leo_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_split_to_leo_pos + + + """ + WHAT IS NEEDED: + * GAUSSIAN WIDTH FOR ALL PARTICLES - FROM PROFILES (DONE) + * PEAK POSITION FOR ALL PARTICLES - FROM C2C (DONE) + * LAST BIN TO USE FOR LEO FIT (from split detector) + """ output_ds = my_binary.copy() - output_ds['leo_FtMaxPos_ch0'] = (('event_index'), leo_FtMaxPos_ch0) - output_ds['leo_FtMaxPosAmpFactor_ch0'] = (('event_index'), leo_FtMaxPosAmpFactor_ch0) +# output_ds['leo_FtMaxPos_ch0'] = (('event_index'), leo_FtMaxPos_ch0) #same as leo_PkPos_ch0? +# output_ds['leo_FtMaxPosAmpFactor_ch0'] = (('event_index'), leo_FtMaxPosAmpFactor_ch0) output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) - # ADD HERE: CALCULATE leo_PkPos_ch0 from split_point + moving_median_high_gain_c2c for all particles + #distance from cross-to-centre (split point to laser maximum intensity). This comes from scattering only partilces output_ds['leo_PkPos_ch0'] = (('event_index'), leo_c2c_ch0) + #interpolate to all particles (including and most importantly rBC particles) output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") + #add split position to cross-to-centre distance to get location of the (estimated) location of the peak maximum + #This is needed for the LEO-fit. output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'] + my_binary['PkSplitPos_ch3'] - output_ds['leo_FtMaxPos_ch0'] = output_ds['leo_FtMaxPos_ch0'].interpolate_na(dim="event_index", - method="nearest", fill_value="extrapolate") - output_ds['leo_FtMaxPosAmpFactor_ch0'] = output_ds['leo_FtMaxPosAmpFactor_ch0'].interpolate_na(dim="event_index", + #First add t_alpha_to_split data + output_ds['leo_MaxPos_ch0'] = (('event_index'), leo_split_to_leo_ch0) + #interpolate to all particles + output_ds['leo_MaxPos_ch0'] = output_ds['leo_MaxPos_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") + #calculate the position at which the leo fits should end based on the split position (of all particles) + output_ds['leo_MaxPos_ch0'] = output_ds['leo_MaxPos_ch0'] + output_ds['PkSplitPos_ch3'] + +#same as leo_PkPos_ch0 +# output_ds['leo_FtMaxPos_ch0'] = output_ds['leo_FtMaxPos_ch0'].interpolate_na(dim="event_index", +# method="nearest", fill_value="extrapolate") +# output_ds['leo_FtMaxPosAmpFactor_ch0'] = output_ds['leo_FtMaxPosAmpFactor_ch0'].interpolate_na(dim="event_index", +# method="nearest", fill_value="extrapolate") output_ds['leo_PkFWHM_ch0'] = output_ds['leo_PkFWHM_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") + + + #ToDo: #done: c2c time series (data from scattering particles but extrapolted to all paricles) #done: leo_PkPos_ch0 (calculated from c2c) --> position of the max amplitude calculated from the split and extrapolated c2c - #leo_FtMaxPos_ch0 (from split) + #done: leo_FtMaxPos_ch0 (from split) return output_ds From 1ef8e85a9a4a42be9328a8c0ef2710913aa97379 Mon Sep 17 00:00:00 2001 From: John Backman Date: Mon, 19 May 2025 11:37:04 +0300 Subject: [PATCH 06/43] changed to leo_EndPos_ch0 for clarity --- pysp2/util/leo_fit.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index c44b944..c9e3b02 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -178,7 +178,6 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #JB OK moving_median_high_gain_c2c = np.nanmedian(moving_high_gain_c2c,axis=1) - #JB LATER #leo_FtMaxPosAmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan #will later be converted to int #leo_FtMaxPosAmpFactor_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor @@ -195,7 +194,6 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): leo_split_to_leo_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_split_to_leo_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_split_to_leo_pos - """ WHAT IS NEEDED: * GAUSSIAN WIDTH FOR ALL PARTICLES - FROM PROFILES (DONE) @@ -218,12 +216,16 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'] + my_binary['PkSplitPos_ch3'] #First add t_alpha_to_split data - output_ds['leo_MaxPos_ch0'] = (('event_index'), leo_split_to_leo_ch0) + output_ds['leo_EndPos_ch0'] = (('event_index'), leo_split_to_leo_ch0) #interpolate to all particles - output_ds['leo_MaxPos_ch0'] = output_ds['leo_MaxPos_ch0'].interpolate_na(dim="event_index", + output_ds['leo_EndPos_ch0'] = output_ds['leo_EndPos_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") #calculate the position at which the leo fits should end based on the split position (of all particles) - output_ds['leo_MaxPos_ch0'] = output_ds['leo_MaxPos_ch0'] + output_ds['PkSplitPos_ch3'] + output_ds['leo_EndPos_ch0'] = output_ds['leo_EndPos_ch0'] + output_ds['PkSplitPos_ch3'] + #cleaning up + output_ds['leo_EndPos_ch0'].values = np.where(output_ds['leo_EndPos_ch0'].values < num_base_pts_2_avg, + np.nan, output_ds['leo_EndPos_ch0'].values) + #same as leo_PkPos_ch0 # output_ds['leo_FtMaxPos_ch0'] = output_ds['leo_FtMaxPos_ch0'].interpolate_na(dim="event_index", @@ -236,11 +238,10 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): method="nearest", fill_value="extrapolate") - #ToDo: #done: c2c time series (data from scattering particles but extrapolted to all paricles) #done: leo_PkPos_ch0 (calculated from c2c) --> position of the max amplitude calculated from the split and extrapolated c2c - #done: leo_FtMaxPos_ch0 (from split) + #done: leo_FtEndPos_ch0 (from split) return output_ds From 0159dfc340eee6aa961cef6df65e6519468041e2 Mon Sep 17 00:00:00 2001 From: John Backman Date: Mon, 19 May 2025 12:29:36 +0300 Subject: [PATCH 07/43] changes to leo_fit() --- pysp2/util/leo_fit.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index c9e3b02..49026be 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -276,23 +276,26 @@ def leo_fit(my_binary,Globals=None): #mean of the lowest 3 points leo_base_ch0 = np.mean(data_ch0[:, 0:num_base_pts_2_avg], axis=1) - leo_fit_max_pos = my_binary['leo_FtMaxPos_ch0'].astype(int).values - leo_FtMaxPosAmpFactor_ch0 = my_binary['leo_FtMaxPosAmpFactor_ch0'].values + leo_fit_max_pos = my_binary['leo_EndPos_ch0'].astype(int).values + #leo_FtMaxPosAmpFactor_ch0 = my_binary['leo_FtMaxPosAmpFactor_ch0'].values leo_PkHt_ch0 = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan leo_PkHt_ch0_ = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan for i in range(my_binary.dims['event_index']): - #NAGOYA STYLE - fractional_peak_height_ch0 = data_ch0[i, leo_fit_max_pos[i]] - leo_base_ch0[i] - estimated_peak_height_ch0 = fractional_peak_height_ch0 * leo_FtMaxPosAmpFactor_ch0[i] - leo_PkHt_ch0_[i] = estimated_peak_height_ch0 + #one factor for the peak height: + #fractional_peak_height_ch0 = data_ch0[i, leo_fit_max_pos[i]] - leo_base_ch0[i] + #estimated_peak_height_ch0 = fractional_peak_height_ch0 * leo_FtMaxPosAmpFactor_ch0[i] + #leo_PkHt_ch0_[i] = estimated_peak_height_ch0 max_value = data_ch0[i,:].max() - data_ch0[i,:].min() bins_ = bins[:leo_fit_max_pos[i]] + if len(bins_) < 2: + leo_PkHt_ch0[i] = np.nan + continue #signals data_ch0_ = data_ch0[i, :leo_fit_max_pos[i]] leo_coeff, var_matrix = curve_fit( lambda x, a: _gaus(x, a, pos[i], width[i], leo_base_ch0[i]), - bins_[:], data_ch0_[:], p0=[max_value], maxfev=40, + bins_[:], data_ch0_[:], p0=[leo_base_ch0[i]+max_value], maxfev=40, ftol=1e-3, method='lm' ) #, bounds=(0, 1e6)) #, method='lm' leo_PkHt_ch0[i] = leo_coeff[0] From 1b461b94100245085020504e255304e196d3ab91 Mon Sep 17 00:00:00 2001 From: John Backman Date: Tue, 20 May 2025 15:41:48 +0300 Subject: [PATCH 08/43] changes to baseline for leo fits... --- pysp2/util/leo_fit.py | 50 +++++++++++++++++++++++++++---------------- 1 file changed, 32 insertions(+), 18 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 49026be..89c54a2 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -42,7 +42,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): num_base_pts_2_avg = 10 moving_average_window = 5 - max_amplitude_fraction = 0.03 + max_amplitude_fraction = 0.05 bins = my_binary['columns'] # median_peak_width = my_binary['PkFWHM_ch0'].median().values / 2.35482 @@ -153,13 +153,13 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): i_range = i_profile[i_max] - np.nanmin(i_profile[:i_max]) moving_avg_high_gain_profiles[i,:] = (i_profile - np.nanmin(i_profile[:i_max])) / i_range #interpolate here to get the exact position in fraction (not integer) :: which posiiton (float) is the 0.03 cross in - moving_avg_high_gain_max_leo_pos[i] = np.argwhere(moving_avg_high_gain_profiles[i,:] >= max_amplitude_fraction).min()-1 + moving_avg_high_gain_max_leo_pos[i] = np.argwhere(moving_avg_high_gain_profiles[i,:] >= max_amplitude_fraction).min()-1 +1 moving_avg_high_gain_split_to_leo_pos[i] = moving_avg_high_gain_max_leo_pos[i] - mean_high_gain_split_pos moving_avg_high_gain_max_leo_amplitude_factor[i] = 1./ moving_avg_high_gain_profiles[i, np.round(moving_avg_high_gain_max_leo_pos[i]).astype(int)] #cleaning up - moving_avg_high_gain_max_leo_amplitude_factor = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, - np.nan, moving_avg_high_gain_max_leo_amplitude_factor) + # moving_avg_high_gain_max_leo_amplitude_factor = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, + # np.nan, moving_avg_high_gain_max_leo_amplitude_factor) moving_avg_high_gain_max_leo_pos = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, np.nan, moving_avg_high_gain_max_leo_pos) @@ -173,6 +173,13 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #moving_avg_high_gain_beam_width = np.nanmedian(moving_high_gain_beam_width,axis=1) moving_median_high_gain_beam_width = np.nanmedian(moving_high_gain_beam_width,axis=1) + #moving average of beam width + + moving_high_gain_base = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['Base_ch0'].values, + moving_average_window, axis=0) + moving_median_high_gain_base = np.nanpercentile(moving_high_gain_base, 33,axis=1) + + moving_high_gain_c2c = np.lib.stride_tricks.sliding_window_view(high_gain_c2c, moving_average_window, axis=0) #JB OK @@ -186,6 +193,9 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): leo_PkFWHM_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_PkFWHM_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_beam_width + leo_Base_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan + leo_Base_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_base + #JB OK leo_c2c_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_c2c_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_c2c @@ -206,6 +216,10 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): # output_ds['leo_FtMaxPosAmpFactor_ch0'] = (('event_index'), leo_FtMaxPosAmpFactor_ch0) output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) + output_ds['leo_Base_ch0'] = (('event_index'), leo_Base_ch0) + output_ds['leo_Base_ch0'] = output_ds['leo_Base_ch0'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") + #distance from cross-to-centre (split point to laser maximum intensity). This comes from scattering only partilces output_ds['leo_PkPos_ch0'] = (('event_index'), leo_c2c_ch0) #interpolate to all particles (including and most importantly rBC particles) @@ -273,13 +287,14 @@ def leo_fit(my_binary,Globals=None): #mean of the first num_base_pts_2_avg points #leo_base_ch0 = np.mean(data_ch0[:, 0:num_base_pts_2_avg], axis=1) - #mean of the lowest 3 points - leo_base_ch0 = np.mean(data_ch0[:, 0:num_base_pts_2_avg], axis=1) - + leo_base_ch0 = my_binary['leo_Base_ch0'].values + #data_ch0_sorted = np.sort(data_ch0[:, 0:num_base_pts_2_avg], axis=1) + #leo_base_ch0 = np.mean(data_ch0_sorted[:, 0:int(num_base_pts_2_avg/3)], axis=1) + leo_fit_max_pos = my_binary['leo_EndPos_ch0'].astype(int).values #leo_FtMaxPosAmpFactor_ch0 = my_binary['leo_FtMaxPosAmpFactor_ch0'].values leo_PkHt_ch0 = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan - leo_PkHt_ch0_ = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan + #leo_PkHt_ch0_ = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan for i in range(my_binary.dims['event_index']): #one factor for the peak height: @@ -287,27 +302,26 @@ def leo_fit(my_binary,Globals=None): #estimated_peak_height_ch0 = fractional_peak_height_ch0 * leo_FtMaxPosAmpFactor_ch0[i] #leo_PkHt_ch0_[i] = estimated_peak_height_ch0 max_value = data_ch0[i,:].max() - data_ch0[i,:].min() - bins_ = bins[:leo_fit_max_pos[i]] + #bins_ = bins[num_base_pts_2_avg:leo_fit_max_pos[i]] + bins_ = bins[leo_fit_max_pos[i]-3:leo_fit_max_pos[i]] if len(bins_) < 2: leo_PkHt_ch0[i] = np.nan continue #signals - data_ch0_ = data_ch0[i, :leo_fit_max_pos[i]] + #data_ch0_ = data_ch0[i, num_base_pts_2_avg:leo_fit_max_pos[i]] + data_ch0_ = data_ch0[i, leo_fit_max_pos[i]-3:leo_fit_max_pos[i]] leo_coeff, var_matrix = curve_fit( lambda x, a: _gaus(x, a, pos[i], width[i], leo_base_ch0[i]), - bins_[:], data_ch0_[:], p0=[leo_base_ch0[i]+max_value], maxfev=40, - ftol=1e-3, method='lm' ) #, bounds=(0, 1e6)) #, method='lm' + bins_[:], data_ch0_[:], p0=[data_ch0[i,:].max()], maxfev=100, + ftol=1e-5, method='lm' ) #, bounds=(0, 1e6)) #, method='lm' leo_PkHt_ch0[i] = leo_coeff[0] output_ds = my_binary.copy() output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0) - output_ds['leo_FtAmp_ch0_'] = (('index'), leo_PkHt_ch0_) - output_ds['leo_FtMaxPos_ch0'] = (('index'), leo_fit_max_pos) + #output_ds['leo_FtAmp_ch0_'] = (('index'), leo_PkHt_ch0_) + #output_ds['leo_FtMaxPos_ch0'] = (('index'), leo_fit_max_pos) output_ds['leo_Base_ch0'] = (('index'), leo_base_ch0) - - #my_high_gain_scatterers['FtAmp_ch0'].plot(marker='.',linestyle='') - #my_high_gain_scatterers['leo_FtAmp_ch0_'].plot(marker='+',linestyle='') - #my_high_gain_scatterers['leo_FtAmp_ch0'].plot(marker='x',linestyle='') + return output_ds From eee8d4ca3257ccef730134620590bc41d97c71bb Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 21 May 2025 15:55:21 +0300 Subject: [PATCH 09/43] improvements --- pysp2/util/leo_fit.py | 37 +++++++++++++++++++++---------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 89c54a2..0b9d1b5 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -42,7 +42,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): num_base_pts_2_avg = 10 moving_average_window = 5 - max_amplitude_fraction = 0.05 + max_amplitude_fraction = 0.033 bins = my_binary['columns'] # median_peak_width = my_binary['PkFWHM_ch0'].median().values / 2.35482 @@ -153,19 +153,25 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): i_range = i_profile[i_max] - np.nanmin(i_profile[:i_max]) moving_avg_high_gain_profiles[i,:] = (i_profile - np.nanmin(i_profile[:i_max])) / i_range #interpolate here to get the exact position in fraction (not integer) :: which posiiton (float) is the 0.03 cross in - moving_avg_high_gain_max_leo_pos[i] = np.argwhere(moving_avg_high_gain_profiles[i,:] >= max_amplitude_fraction).min()-1 +1 + moving_avg_high_gain_max_leo_pos[i] = np.argwhere(moving_avg_high_gain_profiles[i,:] >= max_amplitude_fraction).min()-1 moving_avg_high_gain_split_to_leo_pos[i] = moving_avg_high_gain_max_leo_pos[i] - mean_high_gain_split_pos moving_avg_high_gain_max_leo_amplitude_factor[i] = 1./ moving_avg_high_gain_profiles[i, np.round(moving_avg_high_gain_max_leo_pos[i]).astype(int)] #cleaning up - # moving_avg_high_gain_max_leo_amplitude_factor = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, - # np.nan, moving_avg_high_gain_max_leo_amplitude_factor) - moving_avg_high_gain_max_leo_pos = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, - np.nan, moving_avg_high_gain_max_leo_pos) - - moving_avg_high_gain_split_to_leo_pos = np.where(moving_avg_high_gain_split_to_leo_pos < -20. , + np.nan, moving_avg_high_gain_max_leo_pos) + moving_avg_high_gain_split_to_leo_pos = np.where(moving_avg_high_gain_split_to_leo_pos < -30. , np.nan, moving_avg_high_gain_split_to_leo_pos) + moving_avg_high_gain_max_leo_amplitude_factor = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, + np.nan, moving_avg_high_gain_max_leo_amplitude_factor) + + # leo_pos_unique = np.unique(moving_avg_high_gain_split_to_leo_pos) + # mean_leo_pos_unique = np.zeros_like(leo_pos_unique) + # mean_leo_pos_unique_num = np.zeros_like(leo_pos_unique) + # for i,pos in enumerate(leo_pos_unique): + # bl = moving_avg_high_gain_split_to_leo_pos == pos + # mean_leo_pos_unique[i] = np.nanmedian(moving_avg_high_gain_max_leo_amplitude_factor[bl]) + # mean_leo_pos_unique_num[i] = sum(bl) #moving average of beam width moving_high_gain_beam_width = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['PkFWHM_ch0'].values, @@ -173,20 +179,19 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #moving_avg_high_gain_beam_width = np.nanmedian(moving_high_gain_beam_width,axis=1) moving_median_high_gain_beam_width = np.nanmedian(moving_high_gain_beam_width,axis=1) - #moving average of beam width - + #Moving leo_Base_ch0 moving_high_gain_base = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['Base_ch0'].values, moving_average_window, axis=0) - moving_median_high_gain_base = np.nanpercentile(moving_high_gain_base, 33,axis=1) - + moving_median_high_gain_base = np.nanpercentile(moving_high_gain_base, 10,axis=1) + #JB OK + #Moving c2c high gain moving_high_gain_c2c = np.lib.stride_tricks.sliding_window_view(high_gain_c2c, moving_average_window, axis=0) - #JB OK moving_median_high_gain_c2c = np.nanmedian(moving_high_gain_c2c,axis=1) #JB LATER - #leo_FtMaxPosAmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan #will later be converted to int + #leo_FtMaxPosAmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan #leo_FtMaxPosAmpFactor_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor #JB OK @@ -212,8 +217,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): """ output_ds = my_binary.copy() -# output_ds['leo_FtMaxPos_ch0'] = (('event_index'), leo_FtMaxPos_ch0) #same as leo_PkPos_ch0? -# output_ds['leo_FtMaxPosAmpFactor_ch0'] = (('event_index'), leo_FtMaxPosAmpFactor_ch0) + + output_ds['leo_AmpFactor_ch0'] = (('event_index'), np.zeros_like(leo_PkFWHM_ch0)+np.nanmean(moving_avg_high_gain_max_leo_amplitude_factor)) output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) output_ds['leo_Base_ch0'] = (('event_index'), leo_Base_ch0) From f533374f041c7a1a74d91a8b86e972a164d4378f Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 21 May 2025 16:36:16 +0300 Subject: [PATCH 10/43] leo amplitude factor for peak height added (testing) --- pysp2/util/leo_fit.py | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 0b9d1b5..bb89fb0 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -293,19 +293,15 @@ def leo_fit(my_binary,Globals=None): #mean of the first num_base_pts_2_avg points #leo_base_ch0 = np.mean(data_ch0[:, 0:num_base_pts_2_avg], axis=1) leo_base_ch0 = my_binary['leo_Base_ch0'].values - #data_ch0_sorted = np.sort(data_ch0[:, 0:num_base_pts_2_avg], axis=1) - #leo_base_ch0 = np.mean(data_ch0_sorted[:, 0:int(num_base_pts_2_avg/3)], axis=1) + data_ch0_sorted = np.sort(data_ch0[:, 0:num_base_pts_2_avg], axis=1) + leo_base_ch0_ = np.min(data_ch0_sorted[:, 0:int(num_base_pts_2_avg)], axis=1) leo_fit_max_pos = my_binary['leo_EndPos_ch0'].astype(int).values - #leo_FtMaxPosAmpFactor_ch0 = my_binary['leo_FtMaxPosAmpFactor_ch0'].values + leo_AmpFactor_ch0 = my_binary['leo_AmpFactor_ch0'].values leo_PkHt_ch0 = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan - #leo_PkHt_ch0_ = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan + leo_PkHt_ch0_ = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan for i in range(my_binary.dims['event_index']): - #one factor for the peak height: - #fractional_peak_height_ch0 = data_ch0[i, leo_fit_max_pos[i]] - leo_base_ch0[i] - #estimated_peak_height_ch0 = fractional_peak_height_ch0 * leo_FtMaxPosAmpFactor_ch0[i] - #leo_PkHt_ch0_[i] = estimated_peak_height_ch0 max_value = data_ch0[i,:].max() - data_ch0[i,:].min() #bins_ = bins[num_base_pts_2_avg:leo_fit_max_pos[i]] bins_ = bins[leo_fit_max_pos[i]-3:leo_fit_max_pos[i]] @@ -313,18 +309,18 @@ def leo_fit(my_binary,Globals=None): leo_PkHt_ch0[i] = np.nan continue #signals - #data_ch0_ = data_ch0[i, num_base_pts_2_avg:leo_fit_max_pos[i]] + data_ch0_ = data_ch0[i, num_base_pts_2_avg:leo_fit_max_pos[i]] data_ch0_ = data_ch0[i, leo_fit_max_pos[i]-3:leo_fit_max_pos[i]] leo_coeff, var_matrix = curve_fit( lambda x, a: _gaus(x, a, pos[i], width[i], leo_base_ch0[i]), bins_[:], data_ch0_[:], p0=[data_ch0[i,:].max()], maxfev=100, ftol=1e-5, method='lm' ) #, bounds=(0, 1e6)) #, method='lm' leo_PkHt_ch0[i] = leo_coeff[0] - + leo_PkHt_ch0_[i] = (data_ch0[i, leo_fit_max_pos[i]] - leo_base_ch0_[i]) * leo_AmpFactor_ch0[i] + output_ds = my_binary.copy() - output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0) - #output_ds['leo_FtAmp_ch0_'] = (('index'), leo_PkHt_ch0_) - #output_ds['leo_FtMaxPos_ch0'] = (('index'), leo_fit_max_pos) + output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0_) + output_ds['leo_FtAmp_ch0_'] = (('index'), leo_PkHt_ch0_) output_ds['leo_Base_ch0'] = (('index'), leo_base_ch0) From ba6f7e8ee304199b1b305bae992a9801a469cd04 Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 20 Aug 2025 12:07:15 +0300 Subject: [PATCH 11/43] xr.dims to xr.sizes --- pysp2/util/leo_fit.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index bb89fb0..64575a3 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -80,8 +80,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): event_index = only_scattering_high_gain) #numpy array for the normalized beam profiels - my_high_gain_profiles = np.zeros((my_high_gain_scatterers.dims['index'], - my_high_gain_scatterers.dims['columns'])) \ + my_high_gain_profiles = np.zeros((my_high_gain_scatterers.sizes['index'], + my_high_gain_scatterers.sizes['columns'])) \ * np.nan my_high_gain_profiles_ = np.zeros_like(my_high_gain_profiles)*np.nan @@ -213,7 +213,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): WHAT IS NEEDED: * GAUSSIAN WIDTH FOR ALL PARTICLES - FROM PROFILES (DONE) * PEAK POSITION FOR ALL PARTICLES - FROM C2C (DONE) - * LAST BIN TO USE FOR LEO FIT (from split detector) + * LAST BIN TO USE FOR LEO FIT (from split detector) (DONE) + * do also for LG. """ output_ds = my_binary.copy() @@ -301,7 +302,7 @@ def leo_fit(my_binary,Globals=None): leo_PkHt_ch0 = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan leo_PkHt_ch0_ = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan - for i in range(my_binary.dims['event_index']): + for i in range(my_binary.sizes['event_index']): max_value = data_ch0[i,:].max() - data_ch0[i,:].min() #bins_ = bins[num_base_pts_2_avg:leo_fit_max_pos[i]] bins_ = bins[leo_fit_max_pos[i]-3:leo_fit_max_pos[i]] @@ -322,6 +323,7 @@ def leo_fit(my_binary,Globals=None): output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0_) output_ds['leo_FtAmp_ch0_'] = (('index'), leo_PkHt_ch0_) output_ds['leo_Base_ch0'] = (('index'), leo_base_ch0) + output_ds['leo_Base_ch0_'] = (('index'), leo_base_ch0_) return output_ds From 0f865c141a00c8f50ae460f6bee0dee7fc0baeff Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 20 Aug 2025 12:08:44 +0300 Subject: [PATCH 12/43] force 'fork' method with multiprocessing.Pool --- pysp2/util/peak_fit.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pysp2/util/peak_fit.py b/pysp2/util/peak_fit.py index 9c75d2d..86f38e2 100644 --- a/pysp2/util/peak_fit.py +++ b/pysp2/util/peak_fit.py @@ -1,7 +1,7 @@ import numpy as np import time import dask.bag as db -from multiprocessing import Pool +from multiprocessing import Pool,get_context from itertools import repeat from scipy.optimize import curve_fit @@ -176,17 +176,16 @@ def gaussian_fit(my_ds, config, parallel=False, num_records=None): proc_records = the_bag.map(fit_record).compute() #use multiprocessing.Pool to do the curve fits in parallel elif parallel == 'multiprocessing': - with Pool() as pool: + with get_context("fork").Pool() as pool: proc_records = pool.starmap(_do_fit_records, zip(repeat(my_ds), range(num_records), - repeat(num_trig_pts))) + repeat(num_trig_pts)), chunksize=2000) #else, no parallelism else: proc_records = [] for i in range(num_records): proc_records.append(_do_fit_records(my_ds, i, num_trig_pts)) - FtAmp = np.stack([x[0] for x in proc_records]) FtPos = np.stack([x[1] for x in proc_records]) Base = np.stack([x[2] for x in proc_records]) @@ -604,3 +603,4 @@ def _gaussian_sat_fit(my_ds, record_number): 'height': height, 'error': error} return fit_coeffs + From 67b11a1837715c654baafc8a5168afc7060ebc27 Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 10 Sep 2025 07:54:38 +0300 Subject: [PATCH 13/43] bugfix leo peak amplitude factor --- pysp2/util/leo_fit.py | 37 +++++++++++++++++++++++-------------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 64575a3..56ee7c7 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -1,8 +1,8 @@ import numpy as np from scipy.optimize import curve_fit import xarray as xr -from .peak_fit import _gaus, _do_fit_records -#from pysp2.util.peak_fit import _do_fit_records +#from .peak_fit import _gaus, _do_fit_records +from pysp2.util.peak_fit import _do_fit_records def beam_shape(my_binary, beam_position_from='split point', Globals=None): """ @@ -40,9 +40,9 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): """ - num_base_pts_2_avg = 10 - moving_average_window = 5 - max_amplitude_fraction = 0.033 + num_base_pts_2_avg = 10 #take from globals + moving_average_window = 5 #make an argument out of this + max_amplitude_fraction = 0.033 #make and argument out of this bins = my_binary['columns'] # median_peak_width = my_binary['PkFWHM_ch0'].median().values / 2.35482 @@ -96,7 +96,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): high_gain_split_position = my_high_gain_scatterers['PkSplitPos_ch3'].values - #loop through all particle events (THIS CAN BE MADE SMARTER WITH MUCH OF THE CALCULATIONS BEFORE THE LOOP) + #loop through all particle events (THIS CAN BE MADE SMARTER WITH MUCH OF + #THE CALCULATIONS BEFORE THE LOOP) --> TBD later for i in my_high_gain_scatterers['event_index']: data = my_high_gain_scatterers['Data_ch0'].sel(event_index=i).values #base level @@ -110,7 +111,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #normalize the profile to range [~0,1] profile = (data - base) / peak_height #insert the profile as it was recorded (no shifting due to PEAK POSITION or PSD POSITION) - my_high_gain_profiles_[i,:] = profile + #my_high_gain_profiles_[i,:] = profile #distance to the mean beam peak position if beam_position_from == 'peak maximum': peak_difference = mean_high_gain_max_peak_pos - peak_pos @@ -153,7 +154,10 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): i_range = i_profile[i_max] - np.nanmin(i_profile[:i_max]) moving_avg_high_gain_profiles[i,:] = (i_profile - np.nanmin(i_profile[:i_max])) / i_range #interpolate here to get the exact position in fraction (not integer) :: which posiiton (float) is the 0.03 cross in - moving_avg_high_gain_max_leo_pos[i] = np.argwhere(moving_avg_high_gain_profiles[i,:] >= max_amplitude_fraction).min()-1 + #and skip if it is the where the baseline is calculated + above_max_leo_pos = np.ndarray.flatten(np.argwhere(moving_avg_high_gain_profiles[i,:] >= max_amplitude_fraction)) + moving_avg_high_gain_max_leo_pos[i] = above_max_leo_pos[above_max_leo_pos>num_base_pts_2_avg].min()-1 + moving_avg_high_gain_split_to_leo_pos[i] = moving_avg_high_gain_max_leo_pos[i] - mean_high_gain_split_pos moving_avg_high_gain_max_leo_amplitude_factor[i] = 1./ moving_avg_high_gain_profiles[i, np.round(moving_avg_high_gain_max_leo_pos[i]).astype(int)] @@ -198,6 +202,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): leo_PkFWHM_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_PkFWHM_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_beam_width + #TAKE FROM ACTUAL PARTICLE TRACE WITH INCANDESENCE (NOT NEEDED IN ACTUAL LEO FIT?) leo_Base_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_Base_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_base @@ -215,24 +220,27 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): * PEAK POSITION FOR ALL PARTICLES - FROM C2C (DONE) * LAST BIN TO USE FOR LEO FIT (from split detector) (DONE) * do also for LG. + * TAKE BASELINE SCAT FROM ACTUAL TRACERS IN LEO_FIT """ output_ds = my_binary.copy() - output_ds['leo_AmpFactor_ch0'] = (('event_index'), np.zeros_like(leo_PkFWHM_ch0)+np.nanmean(moving_avg_high_gain_max_leo_amplitude_factor)) + output_ds['leo_AmpFactor_ch0'] = (('event_index'), np.zeros_like(leo_PkFWHM_ch0) + + np.nanmean(moving_avg_high_gain_max_leo_amplitude_factor)) output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) output_ds['leo_Base_ch0'] = (('event_index'), leo_Base_ch0) output_ds['leo_Base_ch0'] = output_ds['leo_Base_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") - #distance from cross-to-centre (split point to laser maximum intensity). This comes from scattering only partilces + #distance from cross-to-centre (split point to laser maximum intensity). + #This comes from scattering only partilces output_ds['leo_PkPos_ch0'] = (('event_index'), leo_c2c_ch0) - #interpolate to all particles (including and most importantly rBC particles) + #interpolate to all particles (including, and most importantly, rBC containing particles) output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") - #add split position to cross-to-centre distance to get location of the (estimated) location of the peak maximum - #This is needed for the LEO-fit. + #add split position to cross-to-centre distance to get location of the (estimated) + #location of the peak maximum. This is needed for the LEO-fit. output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'] + my_binary['PkSplitPos_ch3'] #First add t_alpha_to_split data @@ -240,7 +248,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #interpolate to all particles output_ds['leo_EndPos_ch0'] = output_ds['leo_EndPos_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") - #calculate the position at which the leo fits should end based on the split position (of all particles) + #calculate the position at which the leo fits should end based on the split position + #(of all particles) output_ds['leo_EndPos_ch0'] = output_ds['leo_EndPos_ch0'] + output_ds['PkSplitPos_ch3'] #cleaning up output_ds['leo_EndPos_ch0'].values = np.where(output_ds['leo_EndPos_ch0'].values < num_base_pts_2_avg, From 1196753df5ff974d0b929b44c3838563216d5246 Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 10 Sep 2025 09:41:04 +0300 Subject: [PATCH 14/43] tidying up --- pysp2/util/leo_fit.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 56ee7c7..26162cc 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -2,7 +2,7 @@ from scipy.optimize import curve_fit import xarray as xr #from .peak_fit import _gaus, _do_fit_records -from pysp2.util.peak_fit import _do_fit_records +from pysp2.util.peak_fit import _do_fit_records, _gaus def beam_shape(my_binary, beam_position_from='split point', Globals=None): """ @@ -229,9 +229,9 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): np.nanmean(moving_avg_high_gain_max_leo_amplitude_factor)) output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) - output_ds['leo_Base_ch0'] = (('event_index'), leo_Base_ch0) - output_ds['leo_Base_ch0'] = output_ds['leo_Base_ch0'].interpolate_na(dim="event_index", - method="nearest", fill_value="extrapolate") + #output_ds['leo_Base_ch0'] = (('event_index'), leo_Base_ch0) + #output_ds['leo_Base_ch0'] = output_ds['leo_Base_ch0'].interpolate_na(dim="event_index", + # method="nearest", fill_value="extrapolate") #distance from cross-to-centre (split point to laser maximum intensity). #This comes from scattering only partilces @@ -302,9 +302,9 @@ def leo_fit(my_binary,Globals=None): #mean of the first num_base_pts_2_avg points #leo_base_ch0 = np.mean(data_ch0[:, 0:num_base_pts_2_avg], axis=1) - leo_base_ch0 = my_binary['leo_Base_ch0'].values + #leo_base_ch0 = my_binary['leo_Base_ch0'].values data_ch0_sorted = np.sort(data_ch0[:, 0:num_base_pts_2_avg], axis=1) - leo_base_ch0_ = np.min(data_ch0_sorted[:, 0:int(num_base_pts_2_avg)], axis=1) + leo_base_ch0 = np.min(data_ch0_sorted[:, 0:int(num_base_pts_2_avg)], axis=1) leo_fit_max_pos = my_binary['leo_EndPos_ch0'].astype(int).values leo_AmpFactor_ch0 = my_binary['leo_AmpFactor_ch0'].values @@ -312,7 +312,7 @@ def leo_fit(my_binary,Globals=None): leo_PkHt_ch0_ = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan for i in range(my_binary.sizes['event_index']): - max_value = data_ch0[i,:].max() - data_ch0[i,:].min() + #max_value = data_ch0[i,:].max() - data_ch0[i,:].min() #bins_ = bins[num_base_pts_2_avg:leo_fit_max_pos[i]] bins_ = bins[leo_fit_max_pos[i]-3:leo_fit_max_pos[i]] if len(bins_) < 2: @@ -326,13 +326,13 @@ def leo_fit(my_binary,Globals=None): bins_[:], data_ch0_[:], p0=[data_ch0[i,:].max()], maxfev=100, ftol=1e-5, method='lm' ) #, bounds=(0, 1e6)) #, method='lm' leo_PkHt_ch0[i] = leo_coeff[0] - leo_PkHt_ch0_[i] = (data_ch0[i, leo_fit_max_pos[i]] - leo_base_ch0_[i]) * leo_AmpFactor_ch0[i] + leo_PkHt_ch0_[i] = (data_ch0[i, leo_fit_max_pos[i]] - leo_base_ch0[i]) * leo_AmpFactor_ch0[i] output_ds = my_binary.copy() output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0_) output_ds['leo_FtAmp_ch0_'] = (('index'), leo_PkHt_ch0_) + #output_ds['leo_Base_ch0'] = (('index'), leo_base_ch0) output_ds['leo_Base_ch0'] = (('index'), leo_base_ch0) - output_ds['leo_Base_ch0_'] = (('index'), leo_base_ch0_) return output_ds From eb9175793f9c80823349d0a7f88f9f75c47d6f13 Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 10 Sep 2025 09:41:24 +0300 Subject: [PATCH 15/43] adding incand stats --- pysp2/util/particle_properties.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/pysp2/util/particle_properties.py b/pysp2/util/particle_properties.py index 04ca08f..b4a7b10 100644 --- a/pysp2/util/particle_properties.py +++ b/pysp2/util/particle_properties.py @@ -60,11 +60,13 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None): print("Number of scattering particles rejected for peak pos. = %d" % rejectFtPosTotal) PkHt_ch1 = input_ds['PkHt_ch1'].values + rejectMinIncandTotal = np.sum(PkHt_ch0 < Globals.IncanMinPeakHt1) PkHt_ch5 = input_ds['PkHt_ch5'].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 + rejectIncandWidthTotal = width < Globals.IncanMinWidth accepted_incand = np.logical_and(accepted_incand, input_ds['PkHt_ch2'].values >= Globals.IncanMinPeakHt1) accepted_incand = np.logical_and(accepted_incand, @@ -85,6 +87,15 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None): input_ds['PkPos_ch5'].values > Globals.IncanMinPeakPos, input_ds['PkPos_ch5'].values < Globals.IncanMaxPeakPos)) accepted_incand = np.logical_or(unsat_incand, sat_incand) + numIncandFlag = np.sum(accepted_incand) + + if debug: + print("Number of incandescence particles accepted = %d" % numIncandFlag) + print("Number of incandescence particles rejected for min. peak height = %d" % rejectMinIncandTotal) + print("Number of incandescence particles rejected for peak width = %d" % rejectIncandWidthTotal) + #print("Number of incandescence particles rejected for fat peak = %d" % rejectFatPeakTotal) + #print("Number of incandescence particles rejected for peak pos. = %d" % rejectFtPosTotal) + Scat_not_sat = 1e-18 * (Globals.c0Scat1 + Globals.c1Scat1*PkHt_ch0 + Globals.c2Scat1*PkHt_ch0**2) Scat_sat = 1e-18 * (Globals.c0Scat2 + Globals.c1Scat2*PkHt_ch4 + Globals.c2Scat2*PkHt_ch4**2) From ccc0c11d1492396ef786a5461d4f4eb64d20bcdf Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 10 Sep 2025 16:13:40 +0300 Subject: [PATCH 16/43] started adding low gain stuff --- pysp2/util/leo_fit.py | 169 +++++++++++++++++++++++++++++++----------- 1 file changed, 125 insertions(+), 44 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 26162cc..6179d21 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -60,6 +60,15 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): my_binary['PkHt_ch0'].values < Globals.ScatMaxPeakHt1, my_binary['FtPos_ch0'].values < Globals.ScatMaxPeakPos, my_binary['FtPos_ch0'].values > Globals.ScatMinPeakPos)) + + scatter_low_gain_accepted = np.logical_and.reduce(( + my_binary['PkFWHM_ch0'].values > Globals.ScatMinWidth, + my_binary['PkFWHM_ch0'].values < Globals.ScatMaxWidth, + my_binary['PkHt_ch0'].values > Globals.ScatMinPeakHt2, + my_binary['PkHt_ch0'].values < Globals.ScatMaxPeakHt2, + my_binary['FtPos_ch0'].values < Globals.ScatMaxPeakPos, + my_binary['FtPos_ch0'].values > Globals.ScatMinPeakPos)) + #no incandesence signal = True no_incand_trigged = np.logical_and( @@ -69,63 +78,111 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #Particles that only scatter light only_scattering_high_gain = np.logical_and.reduce((scatter_high_gain_accepted, no_incand_trigged)) + only_scattering_low_gain = np.logical_and.reduce((scatter_low_gain_accepted, + no_incand_trigged)) + #iloc = "event_index" and "index" with the scattering only particle events - iloc = np.argwhere(only_scattering_high_gain).flatten() + iloc_high_gain = np.argwhere(only_scattering_high_gain).flatten() + iloc_low_gain = np.argwhere(only_scattering_high_gain).flatten() print('High gain scattering particles for beam analysis :: ', np.sum(only_scattering_high_gain)) + print('Low gain scattering particles for beam analysis :: ', + np.sum(only_scattering_low_gain)) #make an xarray of the purely scattering particles my_high_gain_scatterers = my_binary.sel(index = only_scattering_high_gain, event_index = only_scattering_high_gain) + my_low_gain_scatterers = my_binary.sel(index = only_scattering_low_gain, + event_index = only_scattering_low_gain) + #numpy array for the normalized beam profiels my_high_gain_profiles = np.zeros((my_high_gain_scatterers.sizes['index'], my_high_gain_scatterers.sizes['columns'])) \ * np.nan + my_low_gain_profiles = np.zeros((my_low_gain_scatterers.sizes['index'], + my_low_gain_scatterers.sizes['columns'])) \ + * np.nan - my_high_gain_profiles_ = np.zeros_like(my_high_gain_profiles)*np.nan + #my_high_gain_profiles_ = np.zeros_like(my_high_gain_profiles)*np.nan mean_high_gain_max_peak_pos = int(np.nanmean(my_high_gain_scatterers['PkPos_ch0'].values)) mean_high_gain_split_pos_float = np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values) mean_high_gain_split_pos = np.round(mean_high_gain_split_pos_float).astype(np.int32) #mean_high_gain_split_pos = int(np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values)) - + mean_low_gain_max_peak_pos = int(np.nanmean(my_low_gain_scatterers['PkPos_ch4'].values)) + mean_low_gain_split_pos_float = np.nanmean(my_low_gain_scatterers['PkSplitPos_ch7'].values) + mean_low_gain_split_pos = np.round(mean_low_gain_split_pos_float).astype(np.int32) + #mean_low_gain_split_pos = int(np.nanmean(my_low_gain_scatterers['PkSplitPos_ch7'].values)) + #cross to center high_gain_c2c = my_high_gain_scatterers['FtPos_ch0'].values - my_high_gain_scatterers['PkSplitPos_ch3'].values - - high_gain_split_position = my_high_gain_scatterers['PkSplitPos_ch3'].values + low_gain_c2c = my_low_gain_scatterers['FtPos_ch4'].values - my_low_gain_scatterers['PkSplitPos_ch7'].values + + #high_gain_split_position = my_high_gain_scatterers['PkSplitPos_ch3'].values + #low_gain_split_position = my_low_gain_scatterers['PkSplitPos_ch7'].values #loop through all particle events (THIS CAN BE MADE SMARTER WITH MUCH OF #THE CALCULATIONS BEFORE THE LOOP) --> TBD later for i in my_high_gain_scatterers['event_index']: - data = my_high_gain_scatterers['Data_ch0'].sel(event_index=i).values + data_ch0 = my_high_gain_scatterers['Data_ch0'].sel(event_index=i).values #base level - base = np.mean(data[0:num_base_pts_2_avg]) + base_ch0 = np.mean(data_ch0[0:num_base_pts_2_avg]) #peak height - peak_height = data.max() - base + peak_height_ch0 = data_ch0.max() - base_ch0 #max peak position - peak_pos = data.argmax() + peak_pos_ch0 = data_ch0.argmax() #split position - split_pos = my_high_gain_scatterers['PkSplitPos_ch3'].sel(event_index=i).values + split_pos_ch3 = my_high_gain_scatterers['PkSplitPos_ch3'].sel(event_index=i).values #normalize the profile to range [~0,1] - profile = (data - base) / peak_height + profile_ch0 = (data_ch0 - base_ch0) / peak_height_ch0 #insert the profile as it was recorded (no shifting due to PEAK POSITION or PSD POSITION) #my_high_gain_profiles_[i,:] = profile #distance to the mean beam peak position if beam_position_from == 'peak maximum': - peak_difference = mean_high_gain_max_peak_pos - peak_pos + peak_difference_ch0 = mean_high_gain_max_peak_pos - peak_pos_ch0 elif beam_position_from == 'split point': - peak_difference = mean_high_gain_split_pos - split_pos + peak_difference_ch0 = mean_high_gain_split_pos - split_pos_ch3 #insert so that the peak is at the right position (accounts for #particles travelling at different speeds) - if peak_difference > 0: - my_high_gain_profiles[i, peak_difference:] = profile[:len(data) - - peak_difference] - elif peak_difference < 0: - my_high_gain_profiles[i, :len(data)+peak_difference] = profile[-peak_difference:] + if peak_difference_ch0 > 0: + my_high_gain_profiles[i, peak_difference_ch0:] = profile_ch0[:len(data_ch0) - + peak_difference_ch0] + elif peak_difference_ch0 < 0: + my_high_gain_profiles[i, :len(data_ch0)+peak_difference_ch0] = profile_ch0[-peak_difference_ch0:] else: - my_high_gain_profiles[i, :] = profile + my_high_gain_profiles[i, :] = profile_ch0 + + for i in my_low_gain_scatterers['event_index']: + data_ch4 = my_high_gain_scatterers['Data_ch4'].sel(event_index=i).values + #base level + base_ch4 = np.mean(data_ch4[0:num_base_pts_2_avg]) + #peak height + peak_height_ch4 = data_ch4.max() - base_ch4 + #max peak position + peak_pos_ch4 = data_ch4.argmax() + #split position + split_pos_ch7 = my_high_gain_scatterers['PkSplitPos_ch7'].sel(event_index=i).values + #normalize the profile to range [~0,1] + profile_ch4 = (data_ch4 - base_ch4) / peak_height_ch4 + #insert the profile as it was recorded (no shifting due to PEAK POSITION or PSD POSITION) + #my_high_gain_profiles_[i,:] = profile + #distance to the mean beam peak position + if beam_position_from == 'peak maximum': + peak_difference_ch4 = mean_low_gain_max_peak_pos - peak_pos_ch4 + elif beam_position_from == 'split point': + peak_difference_ch4 = mean_low_gain_split_pos - split_pos_ch7 + #insert so that the peak is at the right position (accounts for + #particles travelling at different speeds) + if peak_difference_ch4 > 0: + my_low_gain_profiles[i, peak_difference_ch4:] = profile_ch4[:len(data_ch4) - + peak_difference_ch4] + elif peak_difference_ch4 < 0: + my_low_gain_profiles[i, :len(data_ch4)+peak_difference_ch4] = profile_ch4[-peak_difference_ch4:] + else: + my_low_gain_profiles[i, :] = profile_ch4 + #MOVING AVERAGE OF THE BEAM PROFILE TO FIND THE DISTANCE BETWEEN THE SPLIT POINT AND THE POINT IN THE #LASER BEAM WHERE PARTICLES CAN START TO EVAPORATE @@ -133,20 +190,21 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #moving average of the beam shape with a window of moving_average_window moving_high_gain_profile_window = np.lib.stride_tricks.sliding_window_view(my_high_gain_profiles, moving_average_window, axis=0) -# moving_avg_high_gain_profiles_ = np.nanmean(moving_high_gain_profile_window,axis=2) moving_avg_high_gain_profiles_ = np.nanmedian(moving_high_gain_profile_window,axis=2) - - # # - # moving_high_gain_FtPos = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['FtPos_ch0'].values, - # moving_average_window, axis=0) - # moving_high_gain_FtPos_ = np.nanmedian(moving_high_gain_FtPos,axis=1) - + moving_low_gain_profile_window = np.lib.stride_tricks.sliding_window_view(my_low_gain_profiles, + moving_average_window, axis=0) + moving_avg_low_gain_profiles_ = np.nanmedian(moving_low_gain_profile_window,axis=2) moving_avg_high_gain_profiles = np.zeros_like(moving_avg_high_gain_profiles_) * np.nan moving_avg_high_gain_split_to_leo_pos = np.zeros(moving_avg_high_gain_profiles_.shape[0]) * np.nan moving_avg_high_gain_max_leo_pos = np.zeros(moving_avg_high_gain_profiles_.shape[0]) * np.nan + moving_avg_low_gain_profiles = np.zeros_like(moving_avg_low_gain_profiles_) * np.nan + moving_avg_low_gain_split_to_leo_pos = np.zeros(moving_avg_low_gain_profiles_.shape[0]) * np.nan + moving_avg_low_gain_max_leo_pos = np.zeros(moving_avg_low_gain_profiles_.shape[0]) * np.nan + moving_avg_high_gain_max_leo_amplitude_factor = np.zeros(moving_avg_high_gain_profiles_.shape[0]) * np.nan + moving_avg_low_gain_max_leo_amplitude_factor = np.zeros(moving_avg_low_gain_profiles_.shape[0]) * np.nan for i in range(moving_avg_high_gain_profiles_.shape[0]): i_profile = moving_avg_high_gain_profiles_[i,:] @@ -161,6 +219,19 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): moving_avg_high_gain_split_to_leo_pos[i] = moving_avg_high_gain_max_leo_pos[i] - mean_high_gain_split_pos moving_avg_high_gain_max_leo_amplitude_factor[i] = 1./ moving_avg_high_gain_profiles[i, np.round(moving_avg_high_gain_max_leo_pos[i]).astype(int)] + for i in range(moving_avg_low_gain_profiles_.shape[0]): + i_profile = moving_avg_low_gain_profiles_[i,:] + i_max = np.nanargmax(i_profile) + i_range = i_profile[i_max] - np.nanmin(i_profile[:i_max]) + moving_avg_low_gain_profiles[i,:] = (i_profile - np.nanmin(i_profile[:i_max])) / i_range + #interpolate here to get the exact position in fraction (not integer) :: which posiiton (float) is the 0.03 cross in + #and skip if it is the where the baseline is calculated + above_max_leo_pos = np.ndarray.flatten(np.argwhere(moving_avg_low_gain_profiles[i,:] >= max_amplitude_fraction)) + moving_avg_low_gain_max_leo_pos[i] = above_max_leo_pos[above_max_leo_pos>num_base_pts_2_avg].min()-1 + + moving_avg_low_gain_split_to_leo_pos[i] = moving_avg_low_gain_max_leo_pos[i] - mean_low_gain_split_pos + moving_avg_low_gain_max_leo_amplitude_factor[i] = 1./ moving_avg_low_gain_profiles[i, np.round(moving_avg_low_gain_max_leo_pos[i]).astype(int)] + #cleaning up moving_avg_high_gain_max_leo_pos = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, np.nan, moving_avg_high_gain_max_leo_pos) @@ -169,50 +240,60 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): moving_avg_high_gain_max_leo_amplitude_factor = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, np.nan, moving_avg_high_gain_max_leo_amplitude_factor) - # leo_pos_unique = np.unique(moving_avg_high_gain_split_to_leo_pos) - # mean_leo_pos_unique = np.zeros_like(leo_pos_unique) - # mean_leo_pos_unique_num = np.zeros_like(leo_pos_unique) - # for i,pos in enumerate(leo_pos_unique): - # bl = moving_avg_high_gain_split_to_leo_pos == pos - # mean_leo_pos_unique[i] = np.nanmedian(moving_avg_high_gain_max_leo_amplitude_factor[bl]) - # mean_leo_pos_unique_num[i] = sum(bl) + moving_avg_low_gain_max_leo_pos = np.where(moving_avg_low_gain_max_leo_pos < num_base_pts_2_avg, + np.nan, moving_avg_low_gain_max_leo_pos) + moving_avg_low_gain_split_to_leo_pos = np.where(moving_avg_low_gain_split_to_leo_pos < -30. , + np.nan, moving_avg_low_gain_split_to_leo_pos) + moving_avg_low_gain_max_leo_amplitude_factor = np.where(moving_avg_low_gain_max_leo_pos < num_base_pts_2_avg, + np.nan, moving_avg_low_gain_max_leo_amplitude_factor) + #moving average of beam width moving_high_gain_beam_width = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['PkFWHM_ch0'].values, moving_average_window, axis=0) - #moving_avg_high_gain_beam_width = np.nanmedian(moving_high_gain_beam_width,axis=1) + moving_low_gain_beam_width = np.lib.stride_tricks.sliding_window_view(my_low_gain_scatterers['PkFWHM_ch4'].values, + moving_average_window, axis=0) moving_median_high_gain_beam_width = np.nanmedian(moving_high_gain_beam_width,axis=1) + moving_median_low_gain_beam_width = np.nanmedian(moving_low_gain_beam_width,axis=1) - #Moving leo_Base_ch0 + #Moving leo_Base moving_high_gain_base = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['Base_ch0'].values, moving_average_window, axis=0) + moving_low_gain_base = np.lib.stride_tricks.sliding_window_view(my_low_gain_scatterers['Base_ch4'].values, + moving_average_window, axis=0) moving_median_high_gain_base = np.nanpercentile(moving_high_gain_base, 10,axis=1) + moving_median_low_gain_base = np.nanpercentile(moving_low_gain_base, 10,axis=1) - #JB OK #Moving c2c high gain moving_high_gain_c2c = np.lib.stride_tricks.sliding_window_view(high_gain_c2c, moving_average_window, axis=0) + moving_low_gain_c2c = np.lib.stride_tricks.sliding_window_view(low_gain_c2c, + moving_average_window, axis=0) + moving_median_high_gain_c2c = np.nanmedian(moving_high_gain_c2c,axis=1) + moving_median_low_gain_c2c = np.nanmedian(moving_low_gain_c2c,axis=1) #JB LATER #leo_FtMaxPosAmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan #leo_FtMaxPosAmpFactor_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor + + #HIT KOM DU #JB OK leo_PkFWHM_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - leo_PkFWHM_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_beam_width + leo_PkFWHM_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_beam_width #TAKE FROM ACTUAL PARTICLE TRACE WITH INCANDESENCE (NOT NEEDED IN ACTUAL LEO FIT?) leo_Base_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - leo_Base_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_base + leo_Base_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_base #JB OK leo_c2c_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - leo_c2c_ch0[iloc[:-moving_average_window+1]] = moving_median_high_gain_c2c + leo_c2c_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_c2c #JB OK leo_split_to_leo_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - leo_split_to_leo_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_split_to_leo_pos + leo_split_to_leo_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_avg_high_gain_split_to_leo_pos """ WHAT IS NEEDED: @@ -309,7 +390,7 @@ def leo_fit(my_binary,Globals=None): leo_fit_max_pos = my_binary['leo_EndPos_ch0'].astype(int).values leo_AmpFactor_ch0 = my_binary['leo_AmpFactor_ch0'].values leo_PkHt_ch0 = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan - leo_PkHt_ch0_ = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan + leo_PkHt_ch0 = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan for i in range(my_binary.sizes['event_index']): #max_value = data_ch0[i,:].max() - data_ch0[i,:].min() @@ -326,11 +407,11 @@ def leo_fit(my_binary,Globals=None): bins_[:], data_ch0_[:], p0=[data_ch0[i,:].max()], maxfev=100, ftol=1e-5, method='lm' ) #, bounds=(0, 1e6)) #, method='lm' leo_PkHt_ch0[i] = leo_coeff[0] - leo_PkHt_ch0_[i] = (data_ch0[i, leo_fit_max_pos[i]] - leo_base_ch0[i]) * leo_AmpFactor_ch0[i] + leo_PkHt_ch0[i] = (data_ch0[i, leo_fit_max_pos[i]] - leo_base_ch0[i]) * leo_AmpFactor_ch0[i] output_ds = my_binary.copy() - output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0_) - output_ds['leo_FtAmp_ch0_'] = (('index'), leo_PkHt_ch0_) + output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0) + #output_ds['leo_FtAmp_ch0_'] = (('index'), leo_PkHt_ch0_) #output_ds['leo_Base_ch0'] = (('index'), leo_base_ch0) output_ds['leo_Base_ch0'] = (('index'), leo_base_ch0) From 739c71e283ee4ec8423caed7659c4093461a2639 Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 10 Sep 2025 16:51:56 +0300 Subject: [PATCH 17/43] more low gain stuff --- pysp2/util/leo_fit.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 6179d21..2b33d98 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -264,7 +264,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): moving_median_high_gain_base = np.nanpercentile(moving_high_gain_base, 10,axis=1) moving_median_low_gain_base = np.nanpercentile(moving_low_gain_base, 10,axis=1) - #Moving c2c high gain + #Moving cross to centre (c2c) moving_high_gain_c2c = np.lib.stride_tricks.sliding_window_view(high_gain_c2c, moving_average_window, axis=0) moving_low_gain_c2c = np.lib.stride_tricks.sliding_window_view(low_gain_c2c, @@ -277,7 +277,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #leo_FtMaxPosAmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan #leo_FtMaxPosAmpFactor_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor - #HIT KOM DU + #HIT KOM DU_ #JB OK leo_PkFWHM_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan From 0171bc764cd56114cca82854a2cfeee7137e1541 Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 11 Sep 2025 10:26:00 +0300 Subject: [PATCH 18/43] low added to beam_shape --- pysp2/util/leo_fit.py | 53 +++++++++++++++++++++++++------------------ 1 file changed, 31 insertions(+), 22 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 2b33d98..d15eb5e 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -277,23 +277,29 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #leo_FtMaxPosAmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan #leo_FtMaxPosAmpFactor_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor - #HIT KOM DU_ - #JB OK leo_PkFWHM_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_PkFWHM_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_beam_width - + leo_PkFWHM_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan + leo_PkFWHM_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_median_low_gain_beam_width + #TAKE FROM ACTUAL PARTICLE TRACE WITH INCANDESENCE (NOT NEEDED IN ACTUAL LEO FIT?) leo_Base_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_Base_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_base + leo_Base_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan + leo_Base_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_median_low_gain_base #JB OK leo_c2c_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_c2c_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_c2c + leo_c2c_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan + leo_c2c_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_median_low_gain_c2c #JB OK leo_split_to_leo_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_split_to_leo_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_avg_high_gain_split_to_leo_pos + leo_split_to_leo_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan + leo_split_to_leo_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_avg_low_gain_split_to_leo_pos """ WHAT IS NEEDED: @@ -309,50 +315,53 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): output_ds['leo_AmpFactor_ch0'] = (('event_index'), np.zeros_like(leo_PkFWHM_ch0) + np.nanmean(moving_avg_high_gain_max_leo_amplitude_factor)) output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) - - #output_ds['leo_Base_ch0'] = (('event_index'), leo_Base_ch0) - #output_ds['leo_Base_ch0'] = output_ds['leo_Base_ch0'].interpolate_na(dim="event_index", - # method="nearest", fill_value="extrapolate") - + output_ds['leo_AmpFactor_ch4'] = (('event_index'), np.zeros_like(leo_PkFWHM_ch4) + + np.nanmean(moving_avg_low_gain_max_leo_amplitude_factor)) + output_ds['leo_PkFWHM_ch4'] = (('event_index'), leo_PkFWHM_ch4) + + #distance from cross-to-centre (split point to laser maximum intensity). #This comes from scattering only partilces output_ds['leo_PkPos_ch0'] = (('event_index'), leo_c2c_ch0) + output_ds['leo_PkPos_ch4'] = (('event_index'), leo_c2c_ch4) #interpolate to all particles (including, and most importantly, rBC containing particles) output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") + output_ds['leo_PkPos_ch4'] = output_ds['leo_PkPos_ch4'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") #add split position to cross-to-centre distance to get location of the (estimated) #location of the peak maximum. This is needed for the LEO-fit. output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'] + my_binary['PkSplitPos_ch3'] + output_ds['leo_PkPos_ch4'] = output_ds['leo_PkPos_ch4'] + my_binary['PkSplitPos_ch7'] #First add t_alpha_to_split data output_ds['leo_EndPos_ch0'] = (('event_index'), leo_split_to_leo_ch0) + output_ds['leo_EndPos_ch4'] = (('event_index'), leo_split_to_leo_ch4) #interpolate to all particles output_ds['leo_EndPos_ch0'] = output_ds['leo_EndPos_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") + output_ds['leo_EndPos_ch4'] = output_ds['leo_EndPos_ch4'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") + #calculate the position at which the leo fits should end based on the split position #(of all particles) output_ds['leo_EndPos_ch0'] = output_ds['leo_EndPos_ch0'] + output_ds['PkSplitPos_ch3'] + output_ds['leo_EndPos_ch4'] = output_ds['leo_EndPos_ch4'] + output_ds['PkSplitPos_ch7'] #cleaning up output_ds['leo_EndPos_ch0'].values = np.where(output_ds['leo_EndPos_ch0'].values < num_base_pts_2_avg, np.nan, output_ds['leo_EndPos_ch0'].values) - - -#same as leo_PkPos_ch0 -# output_ds['leo_FtMaxPos_ch0'] = output_ds['leo_FtMaxPos_ch0'].interpolate_na(dim="event_index", -# method="nearest", fill_value="extrapolate") -# output_ds['leo_FtMaxPosAmpFactor_ch0'] = output_ds['leo_FtMaxPosAmpFactor_ch0'].interpolate_na(dim="event_index", -# method="nearest", fill_value="extrapolate") + output_ds['leo_EndPos_ch4'].values = np.where(output_ds['leo_EndPos_ch4'].values < num_base_pts_2_avg, + np.nan, output_ds['leo_EndPos_ch4'].values) + output_ds['leo_PkFWHM_ch0'] = output_ds['leo_PkFWHM_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") + output_ds['leo_PkFWHM_ch4'] = output_ds['leo_PkFWHM_ch4'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") output_ds['leo_PkPos_ch0'] = output_ds['leo_PkPos_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") - - - #ToDo: - #done: c2c time series (data from scattering particles but extrapolted to all paricles) - #done: leo_PkPos_ch0 (calculated from c2c) --> position of the max amplitude calculated from the split and extrapolated c2c - #done: leo_FtEndPos_ch0 (from split) - + output_ds['leo_PkPos_ch4'] = output_ds['leo_PkPos_ch4'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") + return output_ds def leo_fit(my_binary,Globals=None): From ce4727bf63a9e54f5c256e12ffe1b3e46938502a Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 11 Sep 2025 11:04:49 +0300 Subject: [PATCH 19/43] bug --- pysp2/util/leo_fit.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index d15eb5e..8a66ace 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -62,12 +62,12 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): my_binary['FtPos_ch0'].values > Globals.ScatMinPeakPos)) scatter_low_gain_accepted = np.logical_and.reduce(( - my_binary['PkFWHM_ch0'].values > Globals.ScatMinWidth, - my_binary['PkFWHM_ch0'].values < Globals.ScatMaxWidth, - my_binary['PkHt_ch0'].values > Globals.ScatMinPeakHt2, - my_binary['PkHt_ch0'].values < Globals.ScatMaxPeakHt2, - my_binary['FtPos_ch0'].values < Globals.ScatMaxPeakPos, - my_binary['FtPos_ch0'].values > Globals.ScatMinPeakPos)) + my_binary['PkFWHM_ch4'].values > Globals.ScatMinWidth, + my_binary['PkFWHM_ch4'].values < Globals.ScatMaxWidth, + my_binary['PkHt_ch4'].values > Globals.ScatMinPeakHt2, + my_binary['PkHt_ch4'].values < Globals.ScatMaxPeakHt2, + my_binary['FtPos_ch4'].values < Globals.ScatMaxPeakPos, + my_binary['FtPos_ch4'].values > Globals.ScatMinPeakPos)) #no incandesence signal = True From c22cc3e1c0a43b0d327c98e5acaa0ac004ad99c2 Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 11 Sep 2025 11:14:06 +0300 Subject: [PATCH 20/43] bugfixes --- pysp2/util/leo_fit.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 8a66ace..6fb7d00 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -1,8 +1,8 @@ import numpy as np from scipy.optimize import curve_fit import xarray as xr -#from .peak_fit import _gaus, _do_fit_records -from pysp2.util.peak_fit import _do_fit_records, _gaus +from .peak_fit import _gaus, _do_fit_records +#from pysp2.util.peak_fit import _do_fit_records, _gaus def beam_shape(my_binary, beam_position_from='split point', Globals=None): """ @@ -83,7 +83,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #iloc = "event_index" and "index" with the scattering only particle events iloc_high_gain = np.argwhere(only_scattering_high_gain).flatten() - iloc_low_gain = np.argwhere(only_scattering_high_gain).flatten() + iloc_low_gain = np.argwhere(only_scattering_low_gain).flatten() print('High gain scattering particles for beam analysis :: ', np.sum(only_scattering_high_gain)) @@ -155,7 +155,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): my_high_gain_profiles[i, :] = profile_ch0 for i in my_low_gain_scatterers['event_index']: - data_ch4 = my_high_gain_scatterers['Data_ch4'].sel(event_index=i).values + data_ch4 = my_low_gain_scatterers['Data_ch4'].sel(event_index=i).values #base level base_ch4 = np.mean(data_ch4[0:num_base_pts_2_avg]) #peak height @@ -163,7 +163,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): #max peak position peak_pos_ch4 = data_ch4.argmax() #split position - split_pos_ch7 = my_high_gain_scatterers['PkSplitPos_ch7'].sel(event_index=i).values + split_pos_ch7 = my_low_gain_scatterers['PkSplitPos_ch7'].sel(event_index=i).values #normalize the profile to range [~0,1] profile_ch4 = (data_ch4 - base_ch4) / peak_height_ch4 #insert the profile as it was recorded (no shifting due to PEAK POSITION or PSD POSITION) From 01313aa9d554a5ce5f8d89b39abe7169c609850f Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 11 Sep 2025 12:29:21 +0300 Subject: [PATCH 21/43] added ch4 to leo_fit --- pysp2/util/leo_fit.py | 55 ++++++++++++++++++++++++++++++++----------- 1 file changed, 41 insertions(+), 14 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 6fb7d00..433c7ae 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -385,46 +385,73 @@ def leo_fit(my_binary,Globals=None): #split_to_peak_high_gain[~bl_only_scattering_particles_ch0] = np.nan #get the information about the gaussian fits - pos = my_binary['leo_PkPos_ch0'].values + pos_ch0 = my_binary['leo_PkPos_ch0'].values + pos_ch4 = my_binary['leo_PkPos_ch4'].values + #amplitude = my_binary['PkHt_ch0'].values - width = my_binary['leo_PkFWHM_ch0'].values / 2.35482 #2*sqrt(2*log(2)) + width_ch0 = my_binary['leo_PkFWHM_ch0'].values / 2.35482 #2*sqrt(2*log(2)) + width_ch4 = my_binary['leo_PkFWHM_ch4'].values / 2.35482 #2*sqrt(2*log(2)) data_ch0 = my_binary['Data_ch0'].values + data_ch4 = my_binary['Data_ch4'].values #mean of the first num_base_pts_2_avg points #leo_base_ch0 = np.mean(data_ch0[:, 0:num_base_pts_2_avg], axis=1) #leo_base_ch0 = my_binary['leo_Base_ch0'].values data_ch0_sorted = np.sort(data_ch0[:, 0:num_base_pts_2_avg], axis=1) + data_ch4_sorted = np.sort(data_ch4[:, 0:num_base_pts_2_avg], axis=1) leo_base_ch0 = np.min(data_ch0_sorted[:, 0:int(num_base_pts_2_avg)], axis=1) - - leo_fit_max_pos = my_binary['leo_EndPos_ch0'].astype(int).values + leo_base_ch4 = np.min(data_ch4_sorted[:, 0:int(num_base_pts_2_avg)], axis=1) + + leo_fit_max_pos_ch0 = my_binary['leo_EndPos_ch0'].astype(int).values + leo_fit_max_pos_ch4 = my_binary['leo_EndPos_ch4'].astype(int).values leo_AmpFactor_ch0 = my_binary['leo_AmpFactor_ch0'].values + leo_AmpFactor_ch4 = my_binary['leo_AmpFactor_ch4'].values leo_PkHt_ch0 = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan - leo_PkHt_ch0 = np.zeros_like(my_binary['PkHt_ch0'].values)*np.nan + leo_PkHt_ch4 = np.zeros_like(my_binary['PkHt_ch4'].values)*np.nan + #High gain for i in range(my_binary.sizes['event_index']): #max_value = data_ch0[i,:].max() - data_ch0[i,:].min() #bins_ = bins[num_base_pts_2_avg:leo_fit_max_pos[i]] - bins_ = bins[leo_fit_max_pos[i]-3:leo_fit_max_pos[i]] + bins_ = bins[leo_fit_max_pos_ch0[i]-3:leo_fit_max_pos_ch0[i]] if len(bins_) < 2: leo_PkHt_ch0[i] = np.nan continue #signals - data_ch0_ = data_ch0[i, num_base_pts_2_avg:leo_fit_max_pos[i]] - data_ch0_ = data_ch0[i, leo_fit_max_pos[i]-3:leo_fit_max_pos[i]] + data_ch0_ = data_ch0[i, num_base_pts_2_avg:leo_fit_max_pos_ch0[i]] + data_ch0_ = data_ch0[i, leo_fit_max_pos_ch0[i]-3:leo_fit_max_pos_ch0[i]] leo_coeff, var_matrix = curve_fit( - lambda x, a: _gaus(x, a, pos[i], width[i], leo_base_ch0[i]), + lambda x, a: _gaus(x, a, pos_ch0[i], width_ch0[i], leo_base_ch0[i]), bins_[:], data_ch0_[:], p0=[data_ch0[i,:].max()], maxfev=100, ftol=1e-5, method='lm' ) #, bounds=(0, 1e6)) #, method='lm' leo_PkHt_ch0[i] = leo_coeff[0] - leo_PkHt_ch0[i] = (data_ch0[i, leo_fit_max_pos[i]] - leo_base_ch0[i]) * leo_AmpFactor_ch0[i] - + leo_PkHt_ch0[i] = (data_ch0[i, leo_fit_max_pos_ch0[i]] - leo_base_ch0[i]) * leo_AmpFactor_ch0[i] + + #Low gain + for i in range(my_binary.sizes['event_index']): + #max_value = data_ch0[i,:].max() - data_ch0[i,:].min() + #bins_ = bins[num_base_pts_2_avg:leo_fit_max_pos[i]] + bins_ = bins[leo_fit_max_pos_ch4[i]-3:leo_fit_max_pos_ch4[i]] + if len(bins_) < 2: + leo_PkHt_ch4[i] = np.nan + continue + #signals + data_ch4_ = data_ch4[i, num_base_pts_2_avg:leo_fit_max_pos_ch4[i]] + data_ch4_ = data_ch4[i, leo_fit_max_pos_ch4[i]-3:leo_fit_max_pos_ch4[i]] + leo_coeff, var_matrix = curve_fit( + lambda x, a: _gaus(x, a, pos_ch4[i], width_ch4[i], leo_base_ch4[i]), + bins_[:], data_ch4_[:], p0=[data_ch4[i,:].max()], maxfev=100, + ftol=1e-5, method='lm' ) #, bounds=(0, 1e6)) #, method='lm' + leo_PkHt_ch4[i] = leo_coeff[0] + leo_PkHt_ch4[i] = (data_ch4[i, leo_fit_max_pos_ch4[i]] - leo_base_ch4[i]) * leo_AmpFactor_ch4[i] + + output_ds = my_binary.copy() output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0) - #output_ds['leo_FtAmp_ch0_'] = (('index'), leo_PkHt_ch0_) - #output_ds['leo_Base_ch0'] = (('index'), leo_base_ch0) + output_ds['leo_FtAmp_ch4'] = (('index'), leo_PkHt_ch4) output_ds['leo_Base_ch0'] = (('index'), leo_base_ch0) + output_ds['leo_Base_ch4'] = (('index'), leo_base_ch4) - return output_ds #my_binary = pysp2.util.beam_shape(my_binary, beam_position_from='peak maximum', Globals=global_settings) From 863e1df9f16571cfdc43a73739c566b7b68de51c Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 11 Sep 2025 20:55:58 +0300 Subject: [PATCH 22/43] bugfix for leo_AmpFactor_ch0 and leo_AmpFactor_ch4 --- pysp2/util/leo_fit.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 433c7ae..3924a83 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -28,7 +28,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): poistion. The maximum peak position is determied from the peak-height weighted average peak position. 'split point' = construct the beam profile around the split position. - The split position is taken from the split detector. Not working yet. + The split position is taken from the split detector. Use this if split + detector is installed. Returns ------- @@ -274,8 +275,10 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): moving_median_low_gain_c2c = np.nanmedian(moving_low_gain_c2c,axis=1) #JB LATER - #leo_FtMaxPosAmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - #leo_FtMaxPosAmpFactor_ch0[iloc[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor + leo_AmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan + leo_AmpFactor_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor + leo_AmpFactor_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan + leo_AmpFactor_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_avg_low_gain_max_leo_amplitude_factor leo_PkFWHM_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan @@ -312,14 +315,15 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): output_ds = my_binary.copy() - output_ds['leo_AmpFactor_ch0'] = (('event_index'), np.zeros_like(leo_PkFWHM_ch0) + - np.nanmean(moving_avg_high_gain_max_leo_amplitude_factor)) + output_ds['leo_AmpFactor_ch0'] = (('event_index'), leo_AmpFactor_ch0) + output_ds['leo_AmpFactor_ch0'] = output_ds['leo_AmpFactor_ch0'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) - output_ds['leo_AmpFactor_ch4'] = (('event_index'), np.zeros_like(leo_PkFWHM_ch4) + - np.nanmean(moving_avg_low_gain_max_leo_amplitude_factor)) + output_ds['leo_AmpFactor_ch4'] = (('event_index'), leo_AmpFactor_ch4) + output_ds['leo_AmpFactor_ch4'] = output_ds['leo_AmpFactor_ch4'].interpolate_na(dim="event_index", + method="nearest", fill_value="extrapolate") output_ds['leo_PkFWHM_ch4'] = (('event_index'), leo_PkFWHM_ch4) - #distance from cross-to-centre (split point to laser maximum intensity). #This comes from scattering only partilces output_ds['leo_PkPos_ch0'] = (('event_index'), leo_c2c_ch0) @@ -445,7 +449,6 @@ def leo_fit(my_binary,Globals=None): leo_PkHt_ch4[i] = leo_coeff[0] leo_PkHt_ch4[i] = (data_ch4[i, leo_fit_max_pos_ch4[i]] - leo_base_ch4[i]) * leo_AmpFactor_ch4[i] - output_ds = my_binary.copy() output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0) output_ds['leo_FtAmp_ch4'] = (('index'), leo_PkHt_ch4) From 355d03bd14091bd4ad0eeba7dd905b98d935f209 Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 11 Sep 2025 21:53:10 +0300 Subject: [PATCH 23/43] Changed description --- pysp2/util/leo_fit.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 3924a83..1929feb 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -33,12 +33,13 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): Returns ------- - coeff : numpy array with gaussian fit [amplitude, peakpos, width, base] - of the beam shape profile. - beam_profile : numpy array with the beam profile calculated from the mean - of all the profiles. The array has as many data points as - the input data. - + my_binary : xarray Dataset + Dataset with additional statistics and information about the + laser beam profile, splitpoint positions relative to the beam + profile etc. These are needed for the actual leo_fit() function. + All variables that are added to the xarray Dataset begin with + leo_ and are available for all particles. + """ num_base_pts_2_avg = 10 #take from globals @@ -457,4 +458,3 @@ def leo_fit(my_binary,Globals=None): return output_ds -#my_binary = pysp2.util.beam_shape(my_binary, beam_position_from='peak maximum', Globals=global_settings) From f0c13f906efb40d09e0044f479595eb95fca2877 Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 11 Sep 2025 21:54:48 +0300 Subject: [PATCH 24/43] Changed description --- pysp2/util/leo_fit.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 1929feb..fc79273 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -38,7 +38,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): laser beam profile, splitpoint positions relative to the beam profile etc. These are needed for the actual leo_fit() function. All variables that are added to the xarray Dataset begin with - leo_ and are available for all particles. + "leo_". These leo_ variables are available for all particles, hence + making the leo fit possible. """ From ac480dab172792fdda3520c13b79891a8ec81465 Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 11 Sep 2025 22:37:35 +0300 Subject: [PATCH 25/43] bugfix --- pysp2/util/particle_properties.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pysp2/util/particle_properties.py b/pysp2/util/particle_properties.py index b4a7b10..6a21a8f 100644 --- a/pysp2/util/particle_properties.py +++ b/pysp2/util/particle_properties.py @@ -66,7 +66,7 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None): 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 - rejectIncandWidthTotal = width < Globals.IncanMinWidth + rejectIncandWidthTotal = np.sum(width < Globals.IncanMinWidth) accepted_incand = np.logical_and(accepted_incand, input_ds['PkHt_ch2'].values >= Globals.IncanMinPeakHt1) accepted_incand = np.logical_and(accepted_incand, @@ -112,6 +112,8 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None): output_ds['ScatDiaBC50'] = xr.DataArray(1000*(0.013416 + 25.066*(Scatter**0.18057)), dims=['index']) + + sootMass_not_sat = factor * 1e-3 * ( Globals.c0Mass1 + Globals.c1Mass1*PkHt_ch1 + Globals.c2Mass1*PkHt_ch1**2) sootDiam_not_sat = (sootMass_not_sat/(0.5236e-9*Globals.densityBC))**(1./3.) From 69bd57eb149833b9dcbb5eb7da433d896aa10a7e Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 11 Sep 2025 22:37:58 +0300 Subject: [PATCH 26/43] cleaning the code --- pysp2/util/leo_fit.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index fc79273..fed6ad2 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -39,7 +39,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): profile etc. These are needed for the actual leo_fit() function. All variables that are added to the xarray Dataset begin with "leo_". These leo_ variables are available for all particles, hence - making the leo fit possible. + making the leo fit possible for incandesence particles as well. """ @@ -320,10 +320,10 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): output_ds['leo_AmpFactor_ch0'] = (('event_index'), leo_AmpFactor_ch0) output_ds['leo_AmpFactor_ch0'] = output_ds['leo_AmpFactor_ch0'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") - output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) output_ds['leo_AmpFactor_ch4'] = (('event_index'), leo_AmpFactor_ch4) output_ds['leo_AmpFactor_ch4'] = output_ds['leo_AmpFactor_ch4'].interpolate_na(dim="event_index", method="nearest", fill_value="extrapolate") + output_ds['leo_PkFWHM_ch0'] = (('event_index'), leo_PkFWHM_ch0) output_ds['leo_PkFWHM_ch4'] = (('event_index'), leo_PkFWHM_ch4) #distance from cross-to-centre (split point to laser maximum intensity). From 40b8cd6fce10b6bb201100207df3517225e30541 Mon Sep 17 00:00:00 2001 From: John Backman Date: Mon, 15 Sep 2025 09:35:34 +0300 Subject: [PATCH 27/43] cleaning up --- pysp2/util/leo_fit.py | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index fed6ad2..6a9c32f 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -376,20 +376,6 @@ def leo_fit(my_binary,Globals=None): num_base_pts_2_avg = 10 #max_amplitude_fraction = 0.03 - # bl_scattering_ok = my_binary['ScatRejectKey'].values == 0 - # bl_only_scattering_particles = np.logical_and(bl_scattering_ok, - # my_binary['PkHt_ch1'].values < Globals.IncanMinPeakHt1) - - # #Particles that only scatter light and ch0 not saturated - # bl_only_scattering_particles_ch0 = np.logical_and(my_binary['PkHt_ch0'].values < Globals.ScatMaxPeakHt1, - # bl_only_scattering_particles) - - #split to peak height difference (in bins) for scattering only particles - #split_to_peak_high_gain = my_binary['PkPos_ch0'].values - my_binary['PkSplitPos_ch3'].values - #For particles with inandesence signal, set to NaN since the peak needn't - #be where the laser intensity is the highest, so se to NaN - #split_to_peak_high_gain[~bl_only_scattering_particles_ch0] = np.nan - #get the information about the gaussian fits pos_ch0 = my_binary['leo_PkPos_ch0'].values pos_ch4 = my_binary['leo_PkPos_ch4'].values From f879c827d41ccb3c68531541682c4981a112cd80 Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 17 Sep 2025 11:32:01 +0300 Subject: [PATCH 28/43] do not use bad split positions for beam profile and no negative values --- pysp2/util/leo_fit.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 6a9c32f..c4c6dec 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -138,6 +138,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): peak_pos_ch0 = data_ch0.argmax() #split position split_pos_ch3 = my_high_gain_scatterers['PkSplitPos_ch3'].sel(event_index=i).values + if split_pos_ch3 > peak_pos_ch0: + continue #normalize the profile to range [~0,1] profile_ch0 = (data_ch0 - base_ch0) / peak_height_ch0 #insert the profile as it was recorded (no shifting due to PEAK POSITION or PSD POSITION) @@ -167,6 +169,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): peak_pos_ch4 = data_ch4.argmax() #split position split_pos_ch7 = my_low_gain_scatterers['PkSplitPos_ch7'].sel(event_index=i).values + if split_pos_ch7 > peak_pos_ch4: + continue #normalize the profile to range [~0,1] profile_ch4 = (data_ch4 - base_ch4) / peak_height_ch4 #insert the profile as it was recorded (no shifting due to PEAK POSITION or PSD POSITION) @@ -436,6 +440,10 @@ def leo_fit(my_binary,Globals=None): ftol=1e-5, method='lm' ) #, bounds=(0, 1e6)) #, method='lm' leo_PkHt_ch4[i] = leo_coeff[0] leo_PkHt_ch4[i] = (data_ch4[i, leo_fit_max_pos_ch4[i]] - leo_base_ch4[i]) * leo_AmpFactor_ch4[i] + + #Only positive values, negative values are set to nan. + leo_PkHt_ch0 = np.where(leo_PkHt_ch0>0, leo_PkHt_ch0, np.nan) + leo_PkHt_ch4 = np.where(leo_PkHt_ch4>0, leo_PkHt_ch4, np.nan) output_ds = my_binary.copy() output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0) From 6737df39edee6e930dd03eb030160cae0699219c Mon Sep 17 00:00:00 2001 From: John Backman Date: Wed, 17 Sep 2025 11:32:50 +0300 Subject: [PATCH 29/43] add leo diamteters if data is available --- pysp2/util/particle_properties.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/pysp2/util/particle_properties.py b/pysp2/util/particle_properties.py index 6a21a8f..7fcae2b 100644 --- a/pysp2/util/particle_properties.py +++ b/pysp2/util/particle_properties.py @@ -112,7 +112,20 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None): output_ds['ScatDiaBC50'] = xr.DataArray(1000*(0.013416 + 25.066*(Scatter**0.18057)), dims=['index']) - + #leo stuff + if 'leo_FtAmp_ch0' and 'leo_FtAmp_ch4' in list(input_ds.variables): + print('adding LEO information to the data.') + #accepted_leo_scat_not_sat = accepted_incand + leo_PkHt_ch0 = input_ds['leo_FtAmp_ch0'] + leo_PkHt_ch4 = input_ds['leo_FtAmp_ch4'] + leo_Scat_not_sat = 1e-18 * (Globals.c0Scat1 + Globals.c1Scat1*leo_PkHt_ch0 + Globals.c2Scat1*leo_PkHt_ch0**2) + leo_Scat_sat = 1e-18 * (Globals.c0Scat2 + Globals.c1Scat2*leo_PkHt_ch4 + Globals.c2Scat2*leo_PkHt_ch4**2) + leo_Scatter = np.where(PkHt_ch0 < Globals.ScatMaxPeakHt1, leo_Scat_not_sat, leo_Scat_sat) + #only for inandesence particles + #leo_Scatter = np.where(accepted_incand, leo_Scatter, np.nan) + output_ds['leo_ScatDiaSO4'] = xr.DataArray(1000 * (-0.015256 + 16.835 * leo_Scatter ** 0.15502), + dims=['index']) + sootMass_not_sat = factor * 1e-3 * ( Globals.c0Mass1 + Globals.c1Mass1*PkHt_ch1 + Globals.c2Mass1*PkHt_ch1**2) From 2ff285f0b7de1e91bcf3cfaca7ef894ff02d0587 Mon Sep 17 00:00:00 2001 From: John Backman Date: Thu, 18 Sep 2025 16:08:50 +0300 Subject: [PATCH 30/43] added leo fits to process_psds() which returns a 3d array called leo_IncandScatNumEnsemble with core and shell diameters. --- pysp2/util/particle_properties.py | 36 ++++++++++++++++++++++++++++--- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/pysp2/util/particle_properties.py b/pysp2/util/particle_properties.py index 7fcae2b..30b3b46 100644 --- a/pysp2/util/particle_properties.py +++ b/pysp2/util/particle_properties.py @@ -126,7 +126,6 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None): output_ds['leo_ScatDiaSO4'] = xr.DataArray(1000 * (-0.015256 + 16.835 * leo_Scatter ** 0.15502), dims=['index']) - sootMass_not_sat = factor * 1e-3 * ( Globals.c0Mass1 + Globals.c1Mass1*PkHt_ch1 + Globals.c2Mass1*PkHt_ch1**2) sootDiam_not_sat = (sootMass_not_sat/(0.5236e-9*Globals.densityBC))**(1./3.) @@ -148,7 +147,8 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None): return output_ds -def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, avg_interval=10,deadtime_correction=False): +def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, + avg_interval=10,deadtime_correction=False): """ Processes the Scattering and BC mass size distributions: @@ -245,7 +245,9 @@ def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, avg_ MassScat2total = np.zeros_like(time_bins[:-1]) MassIncand2 = np.zeros_like(time_bins[:-1]) MassIncand2Sat = np.zeros_like(time_bins[:-1]) - + #LEO + leo_IncandScatNumEnsemble = np.zeros((len(time_bins[:-1]), num_bins, num_bins)) + for t in range(len(time_bins) - 1): parts_time = np.logical_and( time_wave >= time_bins[t], time_wave <= time_bins[t + 1]) @@ -288,6 +290,27 @@ def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, avg_ np.add.at(IncanMassEnsemble[t,:], ind, np.multiply( OneOfEvery[the_particles_incan], sootMass[the_particles_incan])) + #LEO 2D ----> + leo_ScatDiaSO4 = particle_ds['leo_ScatDiaSO4'].values / 1000. + leo_scatter_accept = ~np.isnan(leo_ScatDiaSO4) + + the_particles = np.logical_and.reduce((parts_time, incand_accept, leo_scatter_accept)) + # Remove oversize particles + the_particles = np.logical_and.reduce( + (the_particles, SizeIncandOnly < SpecSizeBins[-1] + deltaSize / 2)) + # Remove oversize particles + the_particles = np.logical_and.reduce( + (the_particles, leo_ScatDiaSO4 < SpecSizeBins[-1] + deltaSize / 2)) + ind_leo = np.searchsorted(SpecSizeBins+deltaSize / 2, + leo_ScatDiaSO4[the_particles], side='right') + ind_incan = np.searchsorted(SpecSizeBins+deltaSize / 2, + SizeIncandOnly[the_particles], side='right') + + #axis 0 = time, axis 1 incandesence, axis 2 leo size + np.add.at(leo_IncandScatNumEnsemble[t,:,:], (ind_incan,ind_leo), OneOfEvery[the_particles]) + #LEO 2D <---- + + scat_parts = np.logical_and(scatter_accept, parts_time) incan_parts = np.logical_and(incand_accept, parts_time) #ConcIncanCycle = OneOfEvery * np.sum(incan_parts) @@ -361,6 +384,7 @@ def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, avg_ ScatMassEnsembleBC[t, :] = ScatMassEnsembleBC[t, :] / FlowCycle #Not in use, always zero ScatNumEnsemble[t, :] = ScatNumEnsemble[t, :] / FlowCycle ScatMassEnsemble[t, :] = ScatMassEnsemble[t, :] / FlowCycle + leo_IncandScatNumEnsemble[t,:,:] = leo_IncandScatNumEnsemble[t,:,:] / FlowCycle base_time = pd.Timestamp('1904-01-01') @@ -452,6 +476,11 @@ def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, avg_ ScatNumEnsembleBC.attrs["standard_name"] = "scattering_number_distribution (black carbon)" ScatNumEnsembleBC.attrs["units"] = "cm-3 per bin" + leo_IncandScatNumEnsemble = xr.DataArray(leo_IncandScatNumEnsemble, dims=('time', 'num_bins', 'num_bins')) + leo_IncandScatNumEnsemble.attrs["long_name"] = "leo test" + leo_IncandScatNumEnsemble.attrs["standard_name"] = "leo test" + leo_IncandScatNumEnsemble.attrs["units"] = "cm-3 per bin" + SpecSizeBins = xr.DataArray(SpecSizeBins, dims=('num_bins')) SpecSizeBins.attrs["long_name"] = "Spectra size bin centers" SpecSizeBins.attrs["standard_name"] = "particle_diameter" @@ -475,6 +504,7 @@ def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, avg_ 'IncanNumEnsemble': IncanNumEnsemble, 'IncanMassEnsemble': IncanMassEnsemble, 'ScatNumEnsembleBC': ScatNumEnsembleBC, + 'leo_IncandScatNumEnsemble': leo_IncandScatNumEnsemble, 'NumFracBC': NumFracBC}) return psd_ds From df77f416e745f1dc331a66995d4814878df1f306 Mon Sep 17 00:00:00 2001 From: John Backman Date: Fri, 19 Sep 2025 14:04:47 +0300 Subject: [PATCH 31/43] added argument leo_fit to calc_diam_masses() and process_psds(). Default is Flase. Set to True for leo fits. --- pysp2/util/particle_properties.py | 90 ++++++++++++++++++------------- 1 file changed, 54 insertions(+), 36 deletions(-) diff --git a/pysp2/util/particle_properties.py b/pysp2/util/particle_properties.py index 30b3b46..e470b7a 100644 --- a/pysp2/util/particle_properties.py +++ b/pysp2/util/particle_properties.py @@ -3,10 +3,12 @@ import datetime import pandas as pd import warnings +from .leo_fit import beam_shape, leo_fit + from .DMTGlobals import DMTGlobals -def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None): +def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None, leo_fits=False): """ Calculates the scattering and incandescence diameters/BC masses for each particle. @@ -34,6 +36,9 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None): rejectFtPosTotal = 0 if Globals is None: Globals = DMTGlobals() + if leo_fits: + input_ds = beam_shape(input_ds, beam_position_from='split point', Globals=Globals) + input_ds = leo_fit(input_ds, Globals=Globals) PkHt_ch0 = np.nanmax(np.stack([input_ds['PkHt_ch0'].values, input_ds['FtAmp_ch0'].values]), axis=0) PkHt_ch4 = np.nanmax(np.stack([input_ds['PkHt_ch4'].values, input_ds['FtAmp_ch4'].values]), axis=0) @@ -112,16 +117,12 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None): output_ds['ScatDiaBC50'] = xr.DataArray(1000*(0.013416 + 25.066*(Scatter**0.18057)), dims=['index']) - #leo stuff - if 'leo_FtAmp_ch0' and 'leo_FtAmp_ch4' in list(input_ds.variables): - print('adding LEO information to the data.') - #accepted_leo_scat_not_sat = accepted_incand + if leo_fits: leo_PkHt_ch0 = input_ds['leo_FtAmp_ch0'] leo_PkHt_ch4 = input_ds['leo_FtAmp_ch4'] leo_Scat_not_sat = 1e-18 * (Globals.c0Scat1 + Globals.c1Scat1*leo_PkHt_ch0 + Globals.c2Scat1*leo_PkHt_ch0**2) leo_Scat_sat = 1e-18 * (Globals.c0Scat2 + Globals.c1Scat2*leo_PkHt_ch4 + Globals.c2Scat2*leo_PkHt_ch4**2) leo_Scatter = np.where(PkHt_ch0 < Globals.ScatMaxPeakHt1, leo_Scat_not_sat, leo_Scat_sat) - #only for inandesence particles #leo_Scatter = np.where(accepted_incand, leo_Scatter, np.nan) output_ds['leo_ScatDiaSO4'] = xr.DataArray(1000 * (-0.015256 + 16.835 * leo_Scatter ** 0.15502), dims=['index']) @@ -148,7 +149,7 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None): def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, - avg_interval=10,deadtime_correction=False): + avg_interval=10,deadtime_correction=False, leo_fits=False): """ Processes the Scattering and BC mass size distributions: @@ -245,8 +246,8 @@ def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, MassScat2total = np.zeros_like(time_bins[:-1]) MassIncand2 = np.zeros_like(time_bins[:-1]) MassIncand2Sat = np.zeros_like(time_bins[:-1]) - #LEO - leo_IncandScatNumEnsemble = np.zeros((len(time_bins[:-1]), num_bins, num_bins)) + if leo_fits: + leo_IncandScatNumEnsemble = np.zeros((len(time_bins[:-1]), num_bins, num_bins)) for t in range(len(time_bins) - 1): parts_time = np.logical_and( @@ -290,27 +291,25 @@ def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, np.add.at(IncanMassEnsemble[t,:], ind, np.multiply( OneOfEvery[the_particles_incan], sootMass[the_particles_incan])) - #LEO 2D ----> - leo_ScatDiaSO4 = particle_ds['leo_ScatDiaSO4'].values / 1000. - leo_scatter_accept = ~np.isnan(leo_ScatDiaSO4) - - the_particles = np.logical_and.reduce((parts_time, incand_accept, leo_scatter_accept)) - # Remove oversize particles - the_particles = np.logical_and.reduce( - (the_particles, SizeIncandOnly < SpecSizeBins[-1] + deltaSize / 2)) - # Remove oversize particles - the_particles = np.logical_and.reduce( - (the_particles, leo_ScatDiaSO4 < SpecSizeBins[-1] + deltaSize / 2)) - ind_leo = np.searchsorted(SpecSizeBins+deltaSize / 2, - leo_ScatDiaSO4[the_particles], side='right') - ind_incan = np.searchsorted(SpecSizeBins+deltaSize / 2, - SizeIncandOnly[the_particles], side='right') - - #axis 0 = time, axis 1 incandesence, axis 2 leo size - np.add.at(leo_IncandScatNumEnsemble[t,:,:], (ind_incan,ind_leo), OneOfEvery[the_particles]) - #LEO 2D <---- + if leo_fits: + leo_ScatDiaSO4 = particle_ds['leo_ScatDiaSO4'].values / 1000. + leo_scatter_accept = ~np.isnan(leo_ScatDiaSO4) + + the_particles = np.logical_and.reduce((parts_time, incand_accept, leo_scatter_accept)) + # Remove oversize particles (incand) + the_particles = np.logical_and.reduce( + (the_particles, SizeIncandOnly < SpecSizeBins[-1] + deltaSize / 2)) + # Remove oversize particles (leo scatter size) + the_particles = np.logical_and.reduce( + (the_particles, leo_ScatDiaSO4 < SpecSizeBins[-1] + deltaSize / 2)) + ind_leo = np.searchsorted(SpecSizeBins+deltaSize / 2, + leo_ScatDiaSO4[the_particles], side='right') + ind_incan = np.searchsorted(SpecSizeBins+deltaSize / 2, + SizeIncandOnly[the_particles], side='right') + + #axis 0 = time, axis 1 incandesence, axis 2 leo size + np.add.at(leo_IncandScatNumEnsemble[t,:,:], (ind_incan,ind_leo), OneOfEvery[the_particles]) - scat_parts = np.logical_and(scatter_accept, parts_time) incan_parts = np.logical_and(incand_accept, parts_time) #ConcIncanCycle = OneOfEvery * np.sum(incan_parts) @@ -384,7 +383,8 @@ def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, ScatMassEnsembleBC[t, :] = ScatMassEnsembleBC[t, :] / FlowCycle #Not in use, always zero ScatNumEnsemble[t, :] = ScatNumEnsemble[t, :] / FlowCycle ScatMassEnsemble[t, :] = ScatMassEnsemble[t, :] / FlowCycle - leo_IncandScatNumEnsemble[t,:,:] = leo_IncandScatNumEnsemble[t,:,:] / FlowCycle + if leo_fits: + leo_IncandScatNumEnsemble[t,:,:] = leo_IncandScatNumEnsemble[t,:,:] / FlowCycle base_time = pd.Timestamp('1904-01-01') @@ -476,11 +476,6 @@ def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, ScatNumEnsembleBC.attrs["standard_name"] = "scattering_number_distribution (black carbon)" ScatNumEnsembleBC.attrs["units"] = "cm-3 per bin" - leo_IncandScatNumEnsemble = xr.DataArray(leo_IncandScatNumEnsemble, dims=('time', 'num_bins', 'num_bins')) - leo_IncandScatNumEnsemble.attrs["long_name"] = "leo test" - leo_IncandScatNumEnsemble.attrs["standard_name"] = "leo test" - leo_IncandScatNumEnsemble.attrs["units"] = "cm-3 per bin" - SpecSizeBins = xr.DataArray(SpecSizeBins, dims=('num_bins')) SpecSizeBins.attrs["long_name"] = "Spectra size bin centers" SpecSizeBins.attrs["standard_name"] = "particle_diameter" @@ -504,7 +499,30 @@ def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, 'IncanNumEnsemble': IncanNumEnsemble, 'IncanMassEnsemble': IncanMassEnsemble, 'ScatNumEnsembleBC': ScatNumEnsembleBC, - 'leo_IncandScatNumEnsemble': leo_IncandScatNumEnsemble, 'NumFracBC': NumFracBC}) + if leo_fit: + #Make the DataArray with differently named size dimensions + leo_IncandScatNumEnsemble = xr.DataArray(leo_IncandScatNumEnsemble, dims=('time', 'incand_bins', 'leo_bins')) + leo_IncandScatNumEnsemble.attrs["long_name"] = "2D incandesence size (black carbon core) and leo size (shell) \ + number size distributions. Dimensions are (time,incandesence size, leo size)" + leo_IncandScatNumEnsemble.attrs["standard_name"] = "2D core / shell size number size distributions" + leo_IncandScatNumEnsemble.attrs["units"] = "cm-3 per bin" + + #Add new dimensions (xarray want's different dimension names although the dimensions are the + #same) + IncandSizeBins = xr.DataArray(SpecSizeBins, dims=('incand_bins')) + IncandSizeBins.attrs["long_name"] = "Incandesence (core) size bin centers" + IncandSizeBins.attrs["standard_name"] = "particle_diameter" + IncandSizeBins.attrs["units"] = "um" + + LeoSizeBins = xr.DataArray(SpecSizeBins, dims=('leo_bins')) + LeoSizeBins.attrs["long_name"] = "Leo (shell) size bin centers" + LeoSizeBins.attrs["standard_name"] = "particle_diameter" + LeoSizeBins.attrs["units"] = "um" + + #add the dimensions to the Dataset + psd_ds = psd_ds.expand_dims(dim={'incand_bins': IncandSizeBins, 'leo_bins': LeoSizeBins}) + #add the 2D (core,shell) variable + psd_ds['leo_IncandScatNumEnsemble'] = leo_IncandScatNumEnsemble return psd_ds From 6e444b8b46fcd0b85726368d38e57418a2eeb29a Mon Sep 17 00:00:00 2001 From: John Backman Date: Fri, 19 Sep 2025 14:05:38 +0300 Subject: [PATCH 32/43] commented out calculations that aren't needed --- pysp2/util/leo_fit.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index c4c6dec..efcd36c 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -46,7 +46,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): num_base_pts_2_avg = 10 #take from globals moving_average_window = 5 #make an argument out of this max_amplitude_fraction = 0.033 #make and argument out of this - bins = my_binary['columns'] + #bins = my_binary['columns'] # median_peak_width = my_binary['PkFWHM_ch0'].median().values / 2.35482 # median_peak_pos = my_binary['FtPos_ch0'].median().values @@ -264,12 +264,12 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): moving_median_low_gain_beam_width = np.nanmedian(moving_low_gain_beam_width,axis=1) #Moving leo_Base - moving_high_gain_base = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['Base_ch0'].values, - moving_average_window, axis=0) - moving_low_gain_base = np.lib.stride_tricks.sliding_window_view(my_low_gain_scatterers['Base_ch4'].values, - moving_average_window, axis=0) - moving_median_high_gain_base = np.nanpercentile(moving_high_gain_base, 10,axis=1) - moving_median_low_gain_base = np.nanpercentile(moving_low_gain_base, 10,axis=1) + #moving_high_gain_base = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['Base_ch0'].values, + # moving_average_window, axis=0) + #moving_low_gain_base = np.lib.stride_tricks.sliding_window_view(my_low_gain_scatterers['Base_ch4'].values, + # moving_average_window, axis=0) + #moving_median_high_gain_base = np.nanpercentile(moving_high_gain_base, 10,axis=1) + #moving_median_low_gain_base = np.nanpercentile(moving_low_gain_base, 10,axis=1) #Moving cross to centre (c2c) moving_high_gain_c2c = np.lib.stride_tricks.sliding_window_view(high_gain_c2c, @@ -293,10 +293,10 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): leo_PkFWHM_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_median_low_gain_beam_width #TAKE FROM ACTUAL PARTICLE TRACE WITH INCANDESENCE (NOT NEEDED IN ACTUAL LEO FIT?) - leo_Base_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - leo_Base_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_base - leo_Base_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan - leo_Base_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_median_low_gain_base + #leo_Base_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan + #leo_Base_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_base + #leo_Base_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan + #leo_Base_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_median_low_gain_base #JB OK leo_c2c_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan From 81f6b8d842a4a3f8976a197579868e0d9182b4b9 Mon Sep 17 00:00:00 2001 From: John Backman Date: Tue, 14 Oct 2025 23:08:16 +0300 Subject: [PATCH 33/43] bug fixes --- pysp2/util/particle_properties.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysp2/util/particle_properties.py b/pysp2/util/particle_properties.py index e470b7a..b95136d 100644 --- a/pysp2/util/particle_properties.py +++ b/pysp2/util/particle_properties.py @@ -500,7 +500,7 @@ def process_psds(particle_ds, hk_ds, config, deltaSize=0.005, num_bins=199, 'IncanMassEnsemble': IncanMassEnsemble, 'ScatNumEnsembleBC': ScatNumEnsembleBC, 'NumFracBC': NumFracBC}) - if leo_fit: + if leo_fits: #Make the DataArray with differently named size dimensions leo_IncandScatNumEnsemble = xr.DataArray(leo_IncandScatNumEnsemble, dims=('time', 'incand_bins', 'leo_bins')) leo_IncandScatNumEnsemble.attrs["long_name"] = "2D incandesence size (black carbon core) and leo size (shell) \ From fefb49269814f1a9272bc99607199e596ec7f694 Mon Sep 17 00:00:00 2001 From: John Backman Date: Tue, 14 Oct 2025 23:09:17 +0300 Subject: [PATCH 34/43] removed index variables as they are not needed --- pysp2/util/leo_fit.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index efcd36c..1df473e 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -94,17 +94,17 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): np.sum(only_scattering_low_gain)) #make an xarray of the purely scattering particles - my_high_gain_scatterers = my_binary.sel(index = only_scattering_high_gain, - event_index = only_scattering_high_gain) - my_low_gain_scatterers = my_binary.sel(index = only_scattering_low_gain, - event_index = only_scattering_low_gain) + my_high_gain_scatterers = my_binary.sel(event_index = only_scattering_high_gain)#, + #index = only_scattering_high_gain() + my_low_gain_scatterers = my_binary.sel(event_index = only_scattering_low_gain)#, + #index = only_scattering_low_gain) #numpy array for the normalized beam profiels - my_high_gain_profiles = np.zeros((my_high_gain_scatterers.sizes['index'], + my_high_gain_profiles = np.zeros((my_high_gain_scatterers.sizes['event_index'], #was index my_high_gain_scatterers.sizes['columns'])) \ * np.nan - my_low_gain_profiles = np.zeros((my_low_gain_scatterers.sizes['index'], + my_low_gain_profiles = np.zeros((my_low_gain_scatterers.sizes['event_index'], #was index my_low_gain_scatterers.sizes['columns'])) \ * np.nan From 4e5f7f95c164589dd7c2672a65812647cf73772b Mon Sep 17 00:00:00 2001 From: John Backman Date: Tue, 14 Oct 2025 23:42:54 +0300 Subject: [PATCH 35/43] removed beam_shape tests for now --- tests/test_gfit.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_gfit.py b/tests/test_gfit.py index 61a8b94..fa4c2a9 100644 --- a/tests/test_gfit.py +++ b/tests/test_gfit.py @@ -34,7 +34,7 @@ def test_psds(): np.testing.assert_almost_equal(my_psds['ScatMassEnsemble'].sum(), 3.15026266) np.testing.assert_almost_equal(my_psds['IncanMassEnsemble'].sum(), 0.08280955) np.testing.assert_almost_equal(my_binary['DeadtimeRelativeBias'].mean(), -0.00023515, decimal=5) - coeff, beam_profile = pysp2.util.beam_shape( - my_binary, beam_position_from='peak maximum', Globals=pysp2.util.DMTGlobals()) - np.testing.assert_almost_equal(coeff, [9.83851858e-01, 4.64317390e+01, - 1.14337852e+01, 4.46761788e-03]) + #coeff, beam_profile = pysp2.util.beam_shape( + # my_binary, beam_position_from='peak maximum', Globals=pysp2.util.DMTGlobals()) + #np.testing.assert_almost_equal(coeff, [9.83851858e-01, 4.64317390e+01, + # 1.14337852e+01, 4.46761788e-03]) From 4ad88205de6deee1c0440a0232aad9c8d140e2f6 Mon Sep 17 00:00:00 2001 From: John Backman Date: Mon, 10 Nov 2025 14:04:55 +0200 Subject: [PATCH 36/43] updated documents --- pysp2/util/leo_fit.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 1df473e..e87c74b 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -8,10 +8,10 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): """ Calculates the beam shape needed to determine the laser intensity profile - used for leading-edge-only fits. The beam position is first determined - from the position of the peak. The beam profile is then calculated by - positioning all the peaks to that position. The beam shape is then the mean - of all the normalized profiles arround that position. + used for leading-edge-only fits. The beam position is by default determined + from the position of the peak relative to the split point. The beam profile + is then calculated by positioning all the peaks to that position. The beam + shape is then the mean of all the normalized profiles arround that position. Parameters ---------- @@ -24,12 +24,12 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): signal limits. beam_position_from : str - 'peak maximum' = construct the beam profile arround the maximum peak - poistion. The maximum peak position is determied from the peak-height - weighted average peak position. 'split point' = construct the beam profile around the split position. The split position is taken from the split detector. Use this if split detector is installed. + 'peak maximum' = construct the beam profile arround the maximum peak + poistion. The maximum peak position is determied from the peak-height + weighted average peak position. Returns ------- From be2363857bbf51932aebd4220369ff00ed94b3a0 Mon Sep 17 00:00:00 2001 From: John Backman Date: Mon, 10 Nov 2025 16:58:01 +0200 Subject: [PATCH 37/43] removed comments --- pysp2/util/leo_fit.py | 55 +++++-------------------------------------- 1 file changed, 6 insertions(+), 49 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index e87c74b..796d1c9 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -29,7 +29,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): detector is installed. 'peak maximum' = construct the beam profile arround the maximum peak poistion. The maximum peak position is determied from the peak-height - weighted average peak position. + position. Returns ------- @@ -47,15 +47,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): moving_average_window = 5 #make an argument out of this max_amplitude_fraction = 0.033 #make and argument out of this #bins = my_binary['columns'] + - # median_peak_width = my_binary['PkFWHM_ch0'].median().values / 2.35482 - # median_peak_pos = my_binary['FtPos_ch0'].median().values - # peak_pos_ok = np.logical_and(my_binary['FtPos_ch0'].values >= median_peak_pos - median_peak_width, - # my_binary['FtPos_ch0'].values <= median_peak_pos + median_peak_width) - - - #change this so that it uses scat_reject_key instead. - #scatterin signal ok = True scatter_high_gain_accepted = np.logical_and.reduce(( my_binary['PkFWHM_ch0'].values > Globals.ScatMinWidth, my_binary['PkFWHM_ch0'].values < Globals.ScatMaxWidth, @@ -72,7 +65,6 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): my_binary['FtPos_ch4'].values < Globals.ScatMaxPeakPos, my_binary['FtPos_ch4'].values > Globals.ScatMinPeakPos)) - #no incandesence signal = True no_incand_trigged = np.logical_and( my_binary['PkHt_ch1'].values < Globals.IncanMinPeakHt1, @@ -107,27 +99,19 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): my_low_gain_profiles = np.zeros((my_low_gain_scatterers.sizes['event_index'], #was index my_low_gain_scatterers.sizes['columns'])) \ * np.nan - - #my_high_gain_profiles_ = np.zeros_like(my_high_gain_profiles)*np.nan - + mean_high_gain_max_peak_pos = int(np.nanmean(my_high_gain_scatterers['PkPos_ch0'].values)) mean_high_gain_split_pos_float = np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values) mean_high_gain_split_pos = np.round(mean_high_gain_split_pos_float).astype(np.int32) - #mean_high_gain_split_pos = int(np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values)) mean_low_gain_max_peak_pos = int(np.nanmean(my_low_gain_scatterers['PkPos_ch4'].values)) mean_low_gain_split_pos_float = np.nanmean(my_low_gain_scatterers['PkSplitPos_ch7'].values) mean_low_gain_split_pos = np.round(mean_low_gain_split_pos_float).astype(np.int32) - #mean_low_gain_split_pos = int(np.nanmean(my_low_gain_scatterers['PkSplitPos_ch7'].values)) #cross to center high_gain_c2c = my_high_gain_scatterers['FtPos_ch0'].values - my_high_gain_scatterers['PkSplitPos_ch3'].values low_gain_c2c = my_low_gain_scatterers['FtPos_ch4'].values - my_low_gain_scatterers['PkSplitPos_ch7'].values - - #high_gain_split_position = my_high_gain_scatterers['PkSplitPos_ch3'].values - #low_gain_split_position = my_low_gain_scatterers['PkSplitPos_ch7'].values - - #loop through all particle events (THIS CAN BE MADE SMARTER WITH MUCH OF - #THE CALCULATIONS BEFORE THE LOOP) --> TBD later + + #loop through all particle events for i in my_high_gain_scatterers['event_index']: data_ch0 = my_high_gain_scatterers['Data_ch0'].sel(event_index=i).values #base level @@ -142,8 +126,6 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): continue #normalize the profile to range [~0,1] profile_ch0 = (data_ch0 - base_ch0) / peak_height_ch0 - #insert the profile as it was recorded (no shifting due to PEAK POSITION or PSD POSITION) - #my_high_gain_profiles_[i,:] = profile #distance to the mean beam peak position if beam_position_from == 'peak maximum': peak_difference_ch0 = mean_high_gain_max_peak_pos - peak_pos_ch0 @@ -262,15 +244,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): moving_average_window, axis=0) moving_median_high_gain_beam_width = np.nanmedian(moving_high_gain_beam_width,axis=1) moving_median_low_gain_beam_width = np.nanmedian(moving_low_gain_beam_width,axis=1) - - #Moving leo_Base - #moving_high_gain_base = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['Base_ch0'].values, - # moving_average_window, axis=0) - #moving_low_gain_base = np.lib.stride_tricks.sliding_window_view(my_low_gain_scatterers['Base_ch4'].values, - # moving_average_window, axis=0) - #moving_median_high_gain_base = np.nanpercentile(moving_high_gain_base, 10,axis=1) - #moving_median_low_gain_base = np.nanpercentile(moving_low_gain_base, 10,axis=1) - + #Moving cross to centre (c2c) moving_high_gain_c2c = np.lib.stride_tricks.sliding_window_view(high_gain_c2c, moving_average_window, axis=0) @@ -280,7 +254,6 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): moving_median_high_gain_c2c = np.nanmedian(moving_high_gain_c2c,axis=1) moving_median_low_gain_c2c = np.nanmedian(moving_low_gain_c2c,axis=1) - #JB LATER leo_AmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_AmpFactor_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor leo_AmpFactor_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan @@ -291,33 +264,17 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): leo_PkFWHM_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_beam_width leo_PkFWHM_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan leo_PkFWHM_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_median_low_gain_beam_width - - #TAKE FROM ACTUAL PARTICLE TRACE WITH INCANDESENCE (NOT NEEDED IN ACTUAL LEO FIT?) - #leo_Base_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - #leo_Base_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_base - #leo_Base_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan - #leo_Base_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_median_low_gain_base - #JB OK leo_c2c_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_c2c_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_c2c leo_c2c_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan leo_c2c_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_median_low_gain_c2c - #JB OK leo_split_to_leo_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan leo_split_to_leo_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_avg_high_gain_split_to_leo_pos leo_split_to_leo_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan leo_split_to_leo_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_avg_low_gain_split_to_leo_pos - """ - WHAT IS NEEDED: - * GAUSSIAN WIDTH FOR ALL PARTICLES - FROM PROFILES (DONE) - * PEAK POSITION FOR ALL PARTICLES - FROM C2C (DONE) - * LAST BIN TO USE FOR LEO FIT (from split detector) (DONE) - * do also for LG. - * TAKE BASELINE SCAT FROM ACTUAL TRACERS IN LEO_FIT - """ output_ds = my_binary.copy() From 66b5918da6172203945a5e83e69bd472b4ddd91e Mon Sep 17 00:00:00 2001 From: John Backman Date: Tue, 11 Nov 2025 10:15:51 +0200 Subject: [PATCH 38/43] added arguments and comments --- pysp2/util/leo_fit.py | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 796d1c9..9d69c70 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -2,9 +2,10 @@ from scipy.optimize import curve_fit import xarray as xr from .peak_fit import _gaus, _do_fit_records -#from pysp2.util.peak_fit import _do_fit_records, _gaus -def beam_shape(my_binary, beam_position_from='split point', Globals=None): + +def beam_shape(my_binary, beam_position_from='split point', Globals=None, + moving_average_window = 5, max_amplitude_fraction = 0.033): """ Calculates the beam shape needed to determine the laser intensity profile @@ -31,6 +32,15 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): poistion. The maximum peak position is determied from the peak-height position. + moving_average_window : int + Number of beam profiles to use to calculate the shape of the laser + beam. This is done using scattering only partilces. + + max_amplitude_fraction : float + How much of the signal amplitude to use for the leading edge. The + amplitude is normalized to be between 0 and 1 so a + max_amplitude_fraction of 0.03 means that 3% of the signal is used. + Returns ------- my_binary : xarray Dataset @@ -43,11 +53,8 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None): """ - num_base_pts_2_avg = 10 #take from globals - moving_average_window = 5 #make an argument out of this - max_amplitude_fraction = 0.033 #make and argument out of this - #bins = my_binary['columns'] - + num_base_pts_2_avg = 10 + scatter_high_gain_accepted = np.logical_and.reduce(( my_binary['PkFWHM_ch0'].values > Globals.ScatMinWidth, @@ -335,21 +342,17 @@ def leo_fit(my_binary,Globals=None): bins = my_binary['columns'].astype('float').values #number of points at the beginning to use for base line average num_base_pts_2_avg = 10 - #max_amplitude_fraction = 0.03 #get the information about the gaussian fits pos_ch0 = my_binary['leo_PkPos_ch0'].values pos_ch4 = my_binary['leo_PkPos_ch4'].values - #amplitude = my_binary['PkHt_ch0'].values width_ch0 = my_binary['leo_PkFWHM_ch0'].values / 2.35482 #2*sqrt(2*log(2)) width_ch4 = my_binary['leo_PkFWHM_ch4'].values / 2.35482 #2*sqrt(2*log(2)) data_ch0 = my_binary['Data_ch0'].values data_ch4 = my_binary['Data_ch4'].values #mean of the first num_base_pts_2_avg points - #leo_base_ch0 = np.mean(data_ch0[:, 0:num_base_pts_2_avg], axis=1) - #leo_base_ch0 = my_binary['leo_Base_ch0'].values data_ch0_sorted = np.sort(data_ch0[:, 0:num_base_pts_2_avg], axis=1) data_ch4_sorted = np.sort(data_ch4[:, 0:num_base_pts_2_avg], axis=1) leo_base_ch0 = np.min(data_ch0_sorted[:, 0:int(num_base_pts_2_avg)], axis=1) @@ -364,8 +367,6 @@ def leo_fit(my_binary,Globals=None): #High gain for i in range(my_binary.sizes['event_index']): - #max_value = data_ch0[i,:].max() - data_ch0[i,:].min() - #bins_ = bins[num_base_pts_2_avg:leo_fit_max_pos[i]] bins_ = bins[leo_fit_max_pos_ch0[i]-3:leo_fit_max_pos_ch0[i]] if len(bins_) < 2: leo_PkHt_ch0[i] = np.nan @@ -382,8 +383,6 @@ def leo_fit(my_binary,Globals=None): #Low gain for i in range(my_binary.sizes['event_index']): - #max_value = data_ch0[i,:].max() - data_ch0[i,:].min() - #bins_ = bins[num_base_pts_2_avg:leo_fit_max_pos[i]] bins_ = bins[leo_fit_max_pos_ch4[i]-3:leo_fit_max_pos_ch4[i]] if len(bins_) < 2: leo_PkHt_ch4[i] = np.nan From 7e2e09e1022c28d7c3b074054202272d92763b11 Mon Sep 17 00:00:00 2001 From: John Backman Date: Tue, 11 Nov 2025 10:42:20 +0200 Subject: [PATCH 39/43] added text on leo_fits to docstring --- pysp2/util/particle_properties.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pysp2/util/particle_properties.py b/pysp2/util/particle_properties.py index b95136d..c4a4a5c 100644 --- a/pysp2/util/particle_properties.py +++ b/pysp2/util/particle_properties.py @@ -24,6 +24,8 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None, leo_fits=F Globals: DMTGlobals structure or None DMTGlobals structure containing calibration coefficients. Set to None to use default values for MOSAiC. + leo_fits: boolean + Default is False. If True, LEO-fits will be calculated as well. Returns ------- @@ -98,8 +100,6 @@ def calc_diams_masses(input_ds, debug=True, factor=1.0, Globals=None, leo_fits=F print("Number of incandescence particles accepted = %d" % numIncandFlag) print("Number of incandescence particles rejected for min. peak height = %d" % rejectMinIncandTotal) print("Number of incandescence particles rejected for peak width = %d" % rejectIncandWidthTotal) - #print("Number of incandescence particles rejected for fat peak = %d" % rejectFatPeakTotal) - #print("Number of incandescence particles rejected for peak pos. = %d" % rejectFtPosTotal) Scat_not_sat = 1e-18 * (Globals.c0Scat1 + Globals.c1Scat1*PkHt_ch0 + Globals.c2Scat1*PkHt_ch0**2) From 4306b05bddab5f0a583b9563368f281e3893883c Mon Sep 17 00:00:00 2001 From: John Backman Date: Tue, 11 Nov 2025 11:27:04 +0200 Subject: [PATCH 40/43] removing inf values if they exist --- pysp2/util/leo_fit.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 9d69c70..93faf0d 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -55,7 +55,6 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None, num_base_pts_2_avg = 10 - scatter_high_gain_accepted = np.logical_and.reduce(( my_binary['PkFWHM_ch0'].values > Globals.ScatMinWidth, my_binary['PkFWHM_ch0'].values < Globals.ScatMaxWidth, @@ -397,9 +396,11 @@ def leo_fit(my_binary,Globals=None): leo_PkHt_ch4[i] = leo_coeff[0] leo_PkHt_ch4[i] = (data_ch4[i, leo_fit_max_pos_ch4[i]] - leo_base_ch4[i]) * leo_AmpFactor_ch4[i] - #Only positive values, negative values are set to nan. - leo_PkHt_ch0 = np.where(leo_PkHt_ch0>0, leo_PkHt_ch0, np.nan) - leo_PkHt_ch4 = np.where(leo_PkHt_ch4>0, leo_PkHt_ch4, np.nan) + #Only positive and finite values + leo_PkHt_ch0 = np.where(leo_PkHt_ch0 > 0, leo_PkHt_ch0, np.nan) + leo_PkHt_ch4 = np.where(leo_PkHt_ch4 > 0, leo_PkHt_ch4, np.nan) + leo_PkHt_ch0 = np.where(leo_PkHt_ch0 < np.inf, leo_PkHt_ch0, np.nan) + leo_PkHt_ch4 = np.where(leo_PkHt_ch4 < np.inf, leo_PkHt_ch4, np.nan) output_ds = my_binary.copy() output_ds['leo_FtAmp_ch0'] = (('index'), leo_PkHt_ch0) From 2a5728fb099c398e69bf4a4d97a3c46ca08d0906 Mon Sep 17 00:00:00 2001 From: John Backman Date: Tue, 11 Nov 2025 11:27:53 +0200 Subject: [PATCH 41/43] added test for leo fits as well --- tests/test_gfit.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/tests/test_gfit.py b/tests/test_gfit.py index fa4c2a9..567d2ef 100644 --- a/tests/test_gfit.py +++ b/tests/test_gfit.py @@ -34,7 +34,15 @@ def test_psds(): np.testing.assert_almost_equal(my_psds['ScatMassEnsemble'].sum(), 3.15026266) np.testing.assert_almost_equal(my_psds['IncanMassEnsemble'].sum(), 0.08280955) np.testing.assert_almost_equal(my_binary['DeadtimeRelativeBias'].mean(), -0.00023515, decimal=5) - #coeff, beam_profile = pysp2.util.beam_shape( - # my_binary, beam_position_from='peak maximum', Globals=pysp2.util.DMTGlobals()) - #np.testing.assert_almost_equal(coeff, [9.83851858e-01, 4.64317390e+01, - # 1.14337852e+01, 4.46761788e-03]) + + my_binary = pysp2.util.beam_shape(my_binary, + beam_position_from = 'peak maximum', + Globals=pysp2.util.DMTGlobals()) + np.testing.assert_almost_equal(my_binary['leo_PkPos_ch0'].mean(), 46.47738538) + + my_binary = pysp2.util.beam_shape(my_binary, + beam_position_from = 'split point', + Globals=pysp2.util.DMTGlobals()) + np.testing.assert_almost_equal(my_binary['leo_PkPos_ch4'].mean(), 46.75399774) + + \ No newline at end of file From 8f72046057a6b5b4b724366d18bedaf1c9797d28 Mon Sep 17 00:00:00 2001 From: John Backman Date: Tue, 11 Nov 2025 14:36:21 +0200 Subject: [PATCH 42/43] changed variable names 'high_gain' to 'ch0' and 'low_gain' to 'ch4' to match the rest of the code --- pysp2/util/leo_fit.py | 208 ++++++++++++++++++++---------------------- 1 file changed, 100 insertions(+), 108 deletions(-) diff --git a/pysp2/util/leo_fit.py b/pysp2/util/leo_fit.py index 93faf0d..c636733 100644 --- a/pysp2/util/leo_fit.py +++ b/pysp2/util/leo_fit.py @@ -55,7 +55,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None, num_base_pts_2_avg = 10 - scatter_high_gain_accepted = np.logical_and.reduce(( + scatter_ch0_accepted = np.logical_and.reduce(( my_binary['PkFWHM_ch0'].values > Globals.ScatMinWidth, my_binary['PkFWHM_ch0'].values < Globals.ScatMaxWidth, my_binary['PkHt_ch0'].values > Globals.ScatMinPeakHt1, @@ -63,7 +63,7 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None, my_binary['FtPos_ch0'].values < Globals.ScatMaxPeakPos, my_binary['FtPos_ch0'].values > Globals.ScatMinPeakPos)) - scatter_low_gain_accepted = np.logical_and.reduce(( + scatter_ch4_accepted = np.logical_and.reduce(( my_binary['PkFWHM_ch4'].values > Globals.ScatMinWidth, my_binary['PkFWHM_ch4'].values < Globals.ScatMaxWidth, my_binary['PkHt_ch4'].values > Globals.ScatMinPeakHt2, @@ -77,49 +77,45 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None, my_binary['PkHt_ch5'].values < Globals.IncanMinPeakHt2) #Particles that only scatter light - only_scattering_high_gain = np.logical_and.reduce((scatter_high_gain_accepted, + only_scattering_ch0 = np.logical_and.reduce((scatter_ch0_accepted, no_incand_trigged)) - only_scattering_low_gain = np.logical_and.reduce((scatter_low_gain_accepted, + only_scattering_ch4 = np.logical_and.reduce((scatter_ch4_accepted, no_incand_trigged)) - #iloc = "event_index" and "index" with the scattering only particle events - iloc_high_gain = np.argwhere(only_scattering_high_gain).flatten() - iloc_low_gain = np.argwhere(only_scattering_low_gain).flatten() + iloc_ch0 = np.argwhere(only_scattering_ch0).flatten() + iloc_ch4 = np.argwhere(only_scattering_ch4).flatten() print('High gain scattering particles for beam analysis :: ', - np.sum(only_scattering_high_gain)) + np.sum(only_scattering_ch0)) print('Low gain scattering particles for beam analysis :: ', - np.sum(only_scattering_low_gain)) + np.sum(only_scattering_ch4)) #make an xarray of the purely scattering particles - my_high_gain_scatterers = my_binary.sel(event_index = only_scattering_high_gain)#, - #index = only_scattering_high_gain() - my_low_gain_scatterers = my_binary.sel(event_index = only_scattering_low_gain)#, - #index = only_scattering_low_gain) - + my_ch0_scatterers = my_binary.sel(event_index = only_scattering_ch0) + my_ch4_scatterers = my_binary.sel(event_index = only_scattering_ch4) #numpy array for the normalized beam profiels - my_high_gain_profiles = np.zeros((my_high_gain_scatterers.sizes['event_index'], #was index - my_high_gain_scatterers.sizes['columns'])) \ + my_ch0_profiles = np.zeros((my_ch0_scatterers.sizes['event_index'], #was index + my_ch0_scatterers.sizes['columns'])) \ * np.nan - my_low_gain_profiles = np.zeros((my_low_gain_scatterers.sizes['event_index'], #was index - my_low_gain_scatterers.sizes['columns'])) \ + my_ch4_profiles = np.zeros((my_ch4_scatterers.sizes['event_index'], #was index + my_ch4_scatterers.sizes['columns'])) \ * np.nan - mean_high_gain_max_peak_pos = int(np.nanmean(my_high_gain_scatterers['PkPos_ch0'].values)) - mean_high_gain_split_pos_float = np.nanmean(my_high_gain_scatterers['PkSplitPos_ch3'].values) - mean_high_gain_split_pos = np.round(mean_high_gain_split_pos_float).astype(np.int32) - mean_low_gain_max_peak_pos = int(np.nanmean(my_low_gain_scatterers['PkPos_ch4'].values)) - mean_low_gain_split_pos_float = np.nanmean(my_low_gain_scatterers['PkSplitPos_ch7'].values) - mean_low_gain_split_pos = np.round(mean_low_gain_split_pos_float).astype(np.int32) + mean_ch0_max_peak_pos = int(np.nanmean(my_ch0_scatterers['PkPos_ch0'].values)) + mean_ch0_split_pos_float = np.nanmean(my_ch0_scatterers['PkSplitPos_ch3'].values) + mean_ch0_split_pos = np.round(mean_ch0_split_pos_float).astype(np.int32) + mean_ch4_max_peak_pos = int(np.nanmean(my_ch4_scatterers['PkPos_ch4'].values)) + mean_ch4_split_pos_float = np.nanmean(my_ch4_scatterers['PkSplitPos_ch7'].values) + mean_ch4_split_pos = np.round(mean_ch4_split_pos_float).astype(np.int32) #cross to center - high_gain_c2c = my_high_gain_scatterers['FtPos_ch0'].values - my_high_gain_scatterers['PkSplitPos_ch3'].values - low_gain_c2c = my_low_gain_scatterers['FtPos_ch4'].values - my_low_gain_scatterers['PkSplitPos_ch7'].values + ch0_c2c = my_ch0_scatterers['FtPos_ch0'].values - my_ch0_scatterers['PkSplitPos_ch3'].values + ch4_c2c = my_ch4_scatterers['FtPos_ch4'].values - my_ch4_scatterers['PkSplitPos_ch7'].values #loop through all particle events - for i in my_high_gain_scatterers['event_index']: - data_ch0 = my_high_gain_scatterers['Data_ch0'].sel(event_index=i).values + for i in my_ch0_scatterers['event_index']: + data_ch0 = my_ch0_scatterers['Data_ch0'].sel(event_index=i).values #base level base_ch0 = np.mean(data_ch0[0:num_base_pts_2_avg]) #peak height @@ -127,28 +123,28 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None, #max peak position peak_pos_ch0 = data_ch0.argmax() #split position - split_pos_ch3 = my_high_gain_scatterers['PkSplitPos_ch3'].sel(event_index=i).values + split_pos_ch3 = my_ch0_scatterers['PkSplitPos_ch3'].sel(event_index=i).values if split_pos_ch3 > peak_pos_ch0: continue #normalize the profile to range [~0,1] profile_ch0 = (data_ch0 - base_ch0) / peak_height_ch0 #distance to the mean beam peak position if beam_position_from == 'peak maximum': - peak_difference_ch0 = mean_high_gain_max_peak_pos - peak_pos_ch0 + peak_difference_ch0 = mean_ch0_max_peak_pos - peak_pos_ch0 elif beam_position_from == 'split point': - peak_difference_ch0 = mean_high_gain_split_pos - split_pos_ch3 + peak_difference_ch0 = mean_ch0_split_pos - split_pos_ch3 #insert so that the peak is at the right position (accounts for #particles travelling at different speeds) if peak_difference_ch0 > 0: - my_high_gain_profiles[i, peak_difference_ch0:] = profile_ch0[:len(data_ch0) - + my_ch0_profiles[i, peak_difference_ch0:] = profile_ch0[:len(data_ch0) - peak_difference_ch0] elif peak_difference_ch0 < 0: - my_high_gain_profiles[i, :len(data_ch0)+peak_difference_ch0] = profile_ch0[-peak_difference_ch0:] + my_ch0_profiles[i, :len(data_ch0)+peak_difference_ch0] = profile_ch0[-peak_difference_ch0:] else: - my_high_gain_profiles[i, :] = profile_ch0 + my_ch0_profiles[i, :] = profile_ch0 - for i in my_low_gain_scatterers['event_index']: - data_ch4 = my_low_gain_scatterers['Data_ch4'].sel(event_index=i).values + for i in my_ch4_scatterers['event_index']: + data_ch4 = my_ch4_scatterers['Data_ch4'].sel(event_index=i).values #base level base_ch4 = np.mean(data_ch4[0:num_base_pts_2_avg]) #peak height @@ -156,131 +152,127 @@ def beam_shape(my_binary, beam_position_from='split point', Globals=None, #max peak position peak_pos_ch4 = data_ch4.argmax() #split position - split_pos_ch7 = my_low_gain_scatterers['PkSplitPos_ch7'].sel(event_index=i).values + split_pos_ch7 = my_ch4_scatterers['PkSplitPos_ch7'].sel(event_index=i).values if split_pos_ch7 > peak_pos_ch4: continue #normalize the profile to range [~0,1] profile_ch4 = (data_ch4 - base_ch4) / peak_height_ch4 #insert the profile as it was recorded (no shifting due to PEAK POSITION or PSD POSITION) - #my_high_gain_profiles_[i,:] = profile #distance to the mean beam peak position if beam_position_from == 'peak maximum': - peak_difference_ch4 = mean_low_gain_max_peak_pos - peak_pos_ch4 + peak_difference_ch4 = mean_ch4_max_peak_pos - peak_pos_ch4 elif beam_position_from == 'split point': - peak_difference_ch4 = mean_low_gain_split_pos - split_pos_ch7 + peak_difference_ch4 = mean_ch4_split_pos - split_pos_ch7 #insert so that the peak is at the right position (accounts for #particles travelling at different speeds) if peak_difference_ch4 > 0: - my_low_gain_profiles[i, peak_difference_ch4:] = profile_ch4[:len(data_ch4) - + my_ch4_profiles[i, peak_difference_ch4:] = profile_ch4[:len(data_ch4) - peak_difference_ch4] elif peak_difference_ch4 < 0: - my_low_gain_profiles[i, :len(data_ch4)+peak_difference_ch4] = profile_ch4[-peak_difference_ch4:] + my_ch4_profiles[i, :len(data_ch4)+peak_difference_ch4] = profile_ch4[-peak_difference_ch4:] else: - my_low_gain_profiles[i, :] = profile_ch4 + my_ch4_profiles[i, :] = profile_ch4 #MOVING AVERAGE OF THE BEAM PROFILE TO FIND THE DISTANCE BETWEEN THE SPLIT POINT AND THE POINT IN THE #LASER BEAM WHERE PARTICLES CAN START TO EVAPORATE #moving average of the beam shape with a window of moving_average_window - moving_high_gain_profile_window = np.lib.stride_tricks.sliding_window_view(my_high_gain_profiles, + moving_ch0_profile_window = np.lib.stride_tricks.sliding_window_view(my_ch0_profiles, moving_average_window, axis=0) - moving_avg_high_gain_profiles_ = np.nanmedian(moving_high_gain_profile_window,axis=2) - moving_low_gain_profile_window = np.lib.stride_tricks.sliding_window_view(my_low_gain_profiles, + moving_avg_ch0_profiles_ = np.nanmedian(moving_ch0_profile_window,axis=2) + moving_ch4_profile_window = np.lib.stride_tricks.sliding_window_view(my_ch4_profiles, moving_average_window, axis=0) - moving_avg_low_gain_profiles_ = np.nanmedian(moving_low_gain_profile_window,axis=2) + moving_avg_ch4_profiles_ = np.nanmedian(moving_ch4_profile_window,axis=2) - moving_avg_high_gain_profiles = np.zeros_like(moving_avg_high_gain_profiles_) * np.nan - moving_avg_high_gain_split_to_leo_pos = np.zeros(moving_avg_high_gain_profiles_.shape[0]) * np.nan - moving_avg_high_gain_max_leo_pos = np.zeros(moving_avg_high_gain_profiles_.shape[0]) * np.nan + moving_avg_ch0_profiles = np.zeros_like(moving_avg_ch0_profiles_) * np.nan + moving_avg_ch0_split_to_leo_pos = np.zeros(moving_avg_ch0_profiles_.shape[0]) * np.nan + moving_avg_ch0_max_leo_pos = np.zeros(moving_avg_ch0_profiles_.shape[0]) * np.nan - moving_avg_low_gain_profiles = np.zeros_like(moving_avg_low_gain_profiles_) * np.nan - moving_avg_low_gain_split_to_leo_pos = np.zeros(moving_avg_low_gain_profiles_.shape[0]) * np.nan - moving_avg_low_gain_max_leo_pos = np.zeros(moving_avg_low_gain_profiles_.shape[0]) * np.nan + moving_avg_ch4_profiles = np.zeros_like(moving_avg_ch4_profiles_) * np.nan + moving_avg_ch4_split_to_leo_pos = np.zeros(moving_avg_ch4_profiles_.shape[0]) * np.nan + moving_avg_ch4_max_leo_pos = np.zeros(moving_avg_ch4_profiles_.shape[0]) * np.nan - moving_avg_high_gain_max_leo_amplitude_factor = np.zeros(moving_avg_high_gain_profiles_.shape[0]) * np.nan - moving_avg_low_gain_max_leo_amplitude_factor = np.zeros(moving_avg_low_gain_profiles_.shape[0]) * np.nan + moving_avg_ch0_max_leo_amplitude_factor = np.zeros(moving_avg_ch0_profiles_.shape[0]) * np.nan + moving_avg_ch4_max_leo_amplitude_factor = np.zeros(moving_avg_ch4_profiles_.shape[0]) * np.nan - for i in range(moving_avg_high_gain_profiles_.shape[0]): - i_profile = moving_avg_high_gain_profiles_[i,:] + for i in range(moving_avg_ch0_profiles_.shape[0]): + i_profile = moving_avg_ch0_profiles_[i,:] i_max = np.nanargmax(i_profile) i_range = i_profile[i_max] - np.nanmin(i_profile[:i_max]) - moving_avg_high_gain_profiles[i,:] = (i_profile - np.nanmin(i_profile[:i_max])) / i_range + moving_avg_ch0_profiles[i,:] = (i_profile - np.nanmin(i_profile[:i_max])) / i_range #interpolate here to get the exact position in fraction (not integer) :: which posiiton (float) is the 0.03 cross in #and skip if it is the where the baseline is calculated - above_max_leo_pos = np.ndarray.flatten(np.argwhere(moving_avg_high_gain_profiles[i,:] >= max_amplitude_fraction)) - moving_avg_high_gain_max_leo_pos[i] = above_max_leo_pos[above_max_leo_pos>num_base_pts_2_avg].min()-1 + above_max_leo_pos = np.ndarray.flatten(np.argwhere(moving_avg_ch0_profiles[i,:] >= max_amplitude_fraction)) + moving_avg_ch0_max_leo_pos[i] = above_max_leo_pos[above_max_leo_pos>num_base_pts_2_avg].min()-1 - moving_avg_high_gain_split_to_leo_pos[i] = moving_avg_high_gain_max_leo_pos[i] - mean_high_gain_split_pos - moving_avg_high_gain_max_leo_amplitude_factor[i] = 1./ moving_avg_high_gain_profiles[i, np.round(moving_avg_high_gain_max_leo_pos[i]).astype(int)] + moving_avg_ch0_split_to_leo_pos[i] = moving_avg_ch0_max_leo_pos[i] - mean_ch0_split_pos + moving_avg_ch0_max_leo_amplitude_factor[i] = 1./ moving_avg_ch0_profiles[i, np.round(moving_avg_ch0_max_leo_pos[i]).astype(int)] - for i in range(moving_avg_low_gain_profiles_.shape[0]): - i_profile = moving_avg_low_gain_profiles_[i,:] + for i in range(moving_avg_ch4_profiles_.shape[0]): + i_profile = moving_avg_ch4_profiles_[i,:] i_max = np.nanargmax(i_profile) i_range = i_profile[i_max] - np.nanmin(i_profile[:i_max]) - moving_avg_low_gain_profiles[i,:] = (i_profile - np.nanmin(i_profile[:i_max])) / i_range + moving_avg_ch4_profiles[i,:] = (i_profile - np.nanmin(i_profile[:i_max])) / i_range #interpolate here to get the exact position in fraction (not integer) :: which posiiton (float) is the 0.03 cross in #and skip if it is the where the baseline is calculated - above_max_leo_pos = np.ndarray.flatten(np.argwhere(moving_avg_low_gain_profiles[i,:] >= max_amplitude_fraction)) - moving_avg_low_gain_max_leo_pos[i] = above_max_leo_pos[above_max_leo_pos>num_base_pts_2_avg].min()-1 + above_max_leo_pos = np.ndarray.flatten(np.argwhere(moving_avg_ch4_profiles[i,:] >= max_amplitude_fraction)) + moving_avg_ch4_max_leo_pos[i] = above_max_leo_pos[above_max_leo_pos>num_base_pts_2_avg].min()-1 - moving_avg_low_gain_split_to_leo_pos[i] = moving_avg_low_gain_max_leo_pos[i] - mean_low_gain_split_pos - moving_avg_low_gain_max_leo_amplitude_factor[i] = 1./ moving_avg_low_gain_profiles[i, np.round(moving_avg_low_gain_max_leo_pos[i]).astype(int)] + moving_avg_ch4_split_to_leo_pos[i] = moving_avg_ch4_max_leo_pos[i] - mean_ch4_split_pos + moving_avg_ch4_max_leo_amplitude_factor[i] = 1./ moving_avg_ch4_profiles[i, np.round(moving_avg_ch4_max_leo_pos[i]).astype(int)] #cleaning up - moving_avg_high_gain_max_leo_pos = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, - np.nan, moving_avg_high_gain_max_leo_pos) - moving_avg_high_gain_split_to_leo_pos = np.where(moving_avg_high_gain_split_to_leo_pos < -30. , - np.nan, moving_avg_high_gain_split_to_leo_pos) - moving_avg_high_gain_max_leo_amplitude_factor = np.where(moving_avg_high_gain_max_leo_pos < num_base_pts_2_avg, - np.nan, moving_avg_high_gain_max_leo_amplitude_factor) + moving_avg_ch0_max_leo_pos = np.where(moving_avg_ch0_max_leo_pos < num_base_pts_2_avg, + np.nan, moving_avg_ch0_max_leo_pos) + moving_avg_ch0_split_to_leo_pos = np.where(moving_avg_ch0_split_to_leo_pos < -30. , + np.nan, moving_avg_ch0_split_to_leo_pos) + moving_avg_ch0_max_leo_amplitude_factor = np.where(moving_avg_ch0_max_leo_pos < num_base_pts_2_avg, + np.nan, moving_avg_ch0_max_leo_amplitude_factor) - moving_avg_low_gain_max_leo_pos = np.where(moving_avg_low_gain_max_leo_pos < num_base_pts_2_avg, - np.nan, moving_avg_low_gain_max_leo_pos) - moving_avg_low_gain_split_to_leo_pos = np.where(moving_avg_low_gain_split_to_leo_pos < -30. , - np.nan, moving_avg_low_gain_split_to_leo_pos) - moving_avg_low_gain_max_leo_amplitude_factor = np.where(moving_avg_low_gain_max_leo_pos < num_base_pts_2_avg, - np.nan, moving_avg_low_gain_max_leo_amplitude_factor) + moving_avg_ch4_max_leo_pos = np.where(moving_avg_ch4_max_leo_pos < num_base_pts_2_avg, + np.nan, moving_avg_ch4_max_leo_pos) + moving_avg_ch4_split_to_leo_pos = np.where(moving_avg_ch4_split_to_leo_pos < -30. , + np.nan, moving_avg_ch4_split_to_leo_pos) + moving_avg_ch4_max_leo_amplitude_factor = np.where(moving_avg_ch4_max_leo_pos < num_base_pts_2_avg, + np.nan, moving_avg_ch4_max_leo_amplitude_factor) - #moving average of beam width - moving_high_gain_beam_width = np.lib.stride_tricks.sliding_window_view(my_high_gain_scatterers['PkFWHM_ch0'].values, + moving_ch0_beam_width = np.lib.stride_tricks.sliding_window_view(my_ch0_scatterers['PkFWHM_ch0'].values, moving_average_window, axis=0) - moving_low_gain_beam_width = np.lib.stride_tricks.sliding_window_view(my_low_gain_scatterers['PkFWHM_ch4'].values, + moving_ch4_beam_width = np.lib.stride_tricks.sliding_window_view(my_ch4_scatterers['PkFWHM_ch4'].values, moving_average_window, axis=0) - moving_median_high_gain_beam_width = np.nanmedian(moving_high_gain_beam_width,axis=1) - moving_median_low_gain_beam_width = np.nanmedian(moving_low_gain_beam_width,axis=1) + moving_median_ch0_beam_width = np.nanmedian(moving_ch0_beam_width,axis=1) + moving_median_ch4_beam_width = np.nanmedian(moving_ch4_beam_width,axis=1) #Moving cross to centre (c2c) - moving_high_gain_c2c = np.lib.stride_tricks.sliding_window_view(high_gain_c2c, + moving_ch0_c2c = np.lib.stride_tricks.sliding_window_view(ch0_c2c, moving_average_window, axis=0) - moving_low_gain_c2c = np.lib.stride_tricks.sliding_window_view(low_gain_c2c, + moving_ch4_c2c = np.lib.stride_tricks.sliding_window_view(ch4_c2c, moving_average_window, axis=0) - moving_median_high_gain_c2c = np.nanmedian(moving_high_gain_c2c,axis=1) - moving_median_low_gain_c2c = np.nanmedian(moving_low_gain_c2c,axis=1) + moving_median_ch0_c2c = np.nanmedian(moving_ch0_c2c,axis=1) + moving_median_ch4_c2c = np.nanmedian(moving_ch4_c2c,axis=1) - leo_AmpFactor_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - leo_AmpFactor_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_avg_high_gain_max_leo_amplitude_factor - leo_AmpFactor_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan - leo_AmpFactor_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_avg_low_gain_max_leo_amplitude_factor + leo_AmpFactor_ch0 = np.zeros(scatter_ch0_accepted.shape)*np.nan + leo_AmpFactor_ch0[iloc_ch0[:-moving_average_window+1]] = moving_avg_ch0_max_leo_amplitude_factor + leo_AmpFactor_ch4 = np.zeros(scatter_ch4_accepted.shape)*np.nan + leo_AmpFactor_ch4[iloc_ch4[:-moving_average_window+1]] = moving_avg_ch4_max_leo_amplitude_factor + leo_PkFWHM_ch0 = np.zeros(scatter_ch0_accepted.shape)*np.nan + leo_PkFWHM_ch0[iloc_ch0[:-moving_average_window+1]] = moving_median_ch0_beam_width + leo_PkFWHM_ch4 = np.zeros(scatter_ch4_accepted.shape)*np.nan + leo_PkFWHM_ch4[iloc_ch4[:-moving_average_window+1]] = moving_median_ch4_beam_width - leo_PkFWHM_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - leo_PkFWHM_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_beam_width - leo_PkFWHM_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan - leo_PkFWHM_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_median_low_gain_beam_width - - leo_c2c_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - leo_c2c_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_median_high_gain_c2c - leo_c2c_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan - leo_c2c_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_median_low_gain_c2c - - leo_split_to_leo_ch0 = np.zeros(scatter_high_gain_accepted.shape)*np.nan - leo_split_to_leo_ch0[iloc_high_gain[:-moving_average_window+1]] = moving_avg_high_gain_split_to_leo_pos - leo_split_to_leo_ch4 = np.zeros(scatter_low_gain_accepted.shape)*np.nan - leo_split_to_leo_ch4[iloc_low_gain[:-moving_average_window+1]] = moving_avg_low_gain_split_to_leo_pos + leo_c2c_ch0 = np.zeros(scatter_ch0_accepted.shape)*np.nan + leo_c2c_ch0[iloc_ch0[:-moving_average_window+1]] = moving_median_ch0_c2c + leo_c2c_ch4 = np.zeros(scatter_ch4_accepted.shape)*np.nan + leo_c2c_ch4[iloc_ch4[:-moving_average_window+1]] = moving_median_ch4_c2c + leo_split_to_leo_ch0 = np.zeros(scatter_ch0_accepted.shape)*np.nan + leo_split_to_leo_ch0[iloc_ch0[:-moving_average_window+1]] = moving_avg_ch0_split_to_leo_pos + leo_split_to_leo_ch4 = np.zeros(scatter_ch4_accepted.shape)*np.nan + leo_split_to_leo_ch4[iloc_ch4[:-moving_average_window+1]] = moving_avg_ch4_split_to_leo_pos output_ds = my_binary.copy() From 7de31cda8806d015b4702f0faac23f15e466057f Mon Sep 17 00:00:00 2001 From: John Backman <61827655+johbackm@users.noreply.github.com> Date: Wed, 12 Nov 2025 14:54:17 +0200 Subject: [PATCH 43/43] Update processing_raw_file.rst Changed typo. --- docs/source/users_guide/processing_raw_file.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/users_guide/processing_raw_file.rst b/docs/source/users_guide/processing_raw_file.rst index dacdcac..9227170 100644 --- a/docs/source/users_guide/processing_raw_file.rst +++ b/docs/source/users_guide/processing_raw_file.rst @@ -13,7 +13,7 @@ Processing a raw .sp2b file is as easy as: .. code-block:: python - waveforms = pysp2.io.read_sp2b(file_name) + waveforms = pysp2.io.read_sp2(file_name) Where file_name is the name of the file you want to read. The output is an xarray Dataset containing the waveforms in each channel for each particle. The