diff --git a/ETA.py b/ETA.py index ed4f964..9d4ed13 100644 --- a/ETA.py +++ b/ETA.py @@ -384,39 +384,58 @@ else: # Calibration if calibration != None: - calib_peaks = np.loadtxt( calibration ) - mean_calib_up = np.mean( data[ + import wavelength_calibration as wave_calib + peaks_calib = np.loadtxt( calibration ) # sorted list + peaks_calib = np.sort( peaks_calib ) + mean_up = np.mean( data[ calibrations[ 'top' ][ 'y' ][ 'min' ] : calibrations[ 'top' ][ 'y' ][ 'max' ], calibrations[ 'top' ][ 'x' ][ 'min' ] : calibrations[ 'top' ][ 'x' ][ 'max' ] ] , axis = 0 ) - mean_calib_up -= np.min( mean_calib_up ) - mean_calib_up /= np.max( mean_calib_up ) + peaks_up = np.array( + utils.retrieve_peaks( + mean_up , + window_size = 1 , + max_window_size = 1, + ) + ).astype( int ) - peaks_up = find_peaks( mean_calib_up , height = 0.5 )[0] + calibrations[ 'top' ][ 'x' ][ 'min' ] - diff = np.array( [ np.inf ] ) - polyval = np.polyfit( peaks_up[ : len( calib_peaks ) ] , calib_peaks[ : len( peaks_up ) ] , 1 ) + peakless_up = wave_calib.remove_peaks( mean_up , peaks_up ) + calib = { # hard-coded for now + 'first': 3 , + 'last': 20, + } + first , last = wave_calib.get_extremities( peakless_up , peaks_up , mean_up ) + up = { + 'first': first, + 'last' : last , + } + peaks_up = peaks_up[ up[ 'first' ] : up[ 'last' ] + 1 ] + peaks_calib = peaks_calib[ calib[ 'first' ] : calib[ 'last' ] + 1 ] - while np.max( diff ) > 1000: # we do not have the exact same number of peak (and the same peaks) in calibration model - diff = abs( np.polyval( polyval , peaks_up[ : len( calib_peaks ) ] ) - calib_peaks[ : len( peaks_up ) ] ) - point_to_delete = np.argmax( diff ) # get the "worst" point - peaks_up_new = np.delete( peaks_up , [ point_to_delete ] ) # remove from data value ( other peak detected ) - calib_peaks_new = np.delete( calib_peaks , [ point_to_delete ] ) # remove from model ( peak not detected ) + peaks_up = wave_calib.only_keep_calib( peaks_up , peaks_calib ) - polyfull_up = np.polyfit( peaks_up_new[ : len( calib_peaks ) ] , calib_peaks[ : len( peaks_up_new ) ] , 1 , full = True ) - polyfull_calib = np.polyfit( peaks_up[ : len( calib_peaks_new ) ] , calib_peaks_new[ : len( peaks_up ) ] , 1 , full = True ) + polyval_wavelength = np.polyfit( # x = 0 begins at the start of + peaks_up , # calibartions[ 'top' ][ 'x' ][ 'min' ] + peaks_calib, + 1 , + ) - if polyfull_up[ 1 ][ 0 ] < polyfull_calib[ 1 ][ 0 ]: # which one is a better fit ? - peaks_up = peaks_up_new - polyval = polyfull_up[ 0 ] - else: - calib_peaks = calib_peaks_new - polyval = polyfull_calib[ 0 ] + wavelength = np.polyval( + polyval_wavelength, + np.arange( + 0 , + len( mean_up ), + 1 , + ) + ) - x = np.arange( calibrations[ 'top' ][ 'x' ][ 'min' ] , calibrations[ 'top' ][ 'x' ][ 'max' ] , 1 ) - x = np.polyval( polyval , x ) - - plt.margins( x = 0 ) - plt.plot( x , mean_calib_up ) + plt.plot( + peaks_up, + peaks_calib - np.polyval( + polyval_wavelength, + peaks_up , + ) , + ) plt.show() exit()