Merge pull request #173 from pycroscopy/use_dict3

multiple scattering

pyTEMlib/eels_tools.py

@@ -617,6 +617,100 @@ def energy_loss_function(E: np.ndarray, Ep: float, Ew: float, A: float) -> np.nd
     return fitted_dataset
 
 
+def energy_loss_function(energy_scale: np.ndarray, p: np.ndarray) -> np.ndarray:
+    eps = 1 - p[0]**2/(energy_scale**2+p[1]**2) + 1j * p[1] * p[0]**2/energy_scale/(energy_scale**2+p[1]**2)
+    elf = (-1/eps).imag
+    return elf*p[2]
+
+
+def multiple_scattering(energy_scale: np.ndarray, p: list, core_loss=False)-> np.ndarray:
+    p = np.abs(p)
+    tmfp = p[3]
+    if core_loss:
+        dif = 1
+    else:
+        dif = 16
+    LLene = np.linspace(1, 2048-1,2048)/dif
+
+    SSD = energy_loss_function(LLene, p)
+    ssd  = np.fft.fft(SSD)
+    ssd2 = ssd.copy()
+
+    ### sum contribution from each order of scattering:
+    PSD = np.zeros(len(LLene))
+    for order in range(15):
+        # This order convoluted spectum 
+        # convoluted SSD is SSD2
+        SSD2 = np.fft.ifft(ssd).real
+
+        # scale right (could be done better? GERD) 
+        # And add this order to final spectrum
+        PSD += SSD2*abs(sum(SSD)/sum(SSD2)) / scipy.special.factorial(order+1)*np.power(tmfp, (order+1))*np.exp(-tmfp) #using equation 4.1 of egerton ed2
+
+        # next order convolution
+        ssd = ssd * ssd2
+
+    PSD /=tmfp*np.exp(-tmfp)
+    BGDcoef = scipy.interpolate.splrep(LLene, PSD, s=0)    
+    return scipy.interpolate.splev(energy_scale, BGDcoef)
+
+def fit_multiple_scattering(dataset: Union[sidpy.Dataset, np.ndarray], startFitEnergy: float, endFitEnergy: float,pin=None, number_workers: int = 4, number_threads: int = 8) -> Union[sidpy.Dataset, np.ndarray]:
+    """
+    Fit multiple scattering of plasmon peak in a TEM dataset.
+
+    
+    Parameters:
+        dataset: sidpy.Dataset or numpy.ndarray
+            The dataset containing TEM spectral data.
+        startFitEnergy: float
+            The start energy of the fitting window.
+        endFitEnergy: float
+            The end energy of the fitting window.
+        number_workers: int, optional
+            The number of workers for parallel processing (default is 4).
+        number_threads: int, optional
+            The number of threads for parallel processing (default is 8).
+
+    Returns:
+        fitted_dataset: sidpy.Dataset or numpy.ndarray
+            The dataset with fitted plasmon peak parameters. The dimensions and 
+            format depend on the input dataset.
+
+    Raises:
+        ValueError: If the input dataset does not have the expected dimensions or format.
+
+    Notes:
+        - The function uses the Drude model to fit plasmon peaks.
+        - The fitting parameters are peak position (Ep), peak width (Ew), and amplitude (A).
+        - If `plot_result` is True, the function plots Ep, Ew, and A as separate subplots.
+    """
+
+
+    # define window for fitting
+    energy = dataset.get_spectral_dims(return_axis=True)[0].values
+    start_fit_pixel = np.searchsorted(energy, startFitEnergy)
+    end_fit_pixel = np.searchsorted(energy, endFitEnergy)
+
+    def errf_multi(p, y, x):
+        elf = multiple_scattering(x, p)
+        err = y - elf
+        #print (p,sum(np.abs(err)))
+        return np.abs(err) # /np.sqrt(y)
+
+    if pin is None:
+        pin = np.array([9,1,.7, 0.3])
+
+
+    fit_dset = np.array(dataset[start_fit_pixel:end_fit_pixel])
+    popt, lsq = leastsq(errf_multi, pin, args=(fit_dset, energy[start_fit_pixel:end_fit_pixel]), maxfev=2000)
+
+    multi = dataset.like_data(multiple_scattering(energy, popt))
+
+
+    multi.metadata['multiple_scattering'] = {'parameter': popt}
+    return multi
+
+
 
 def drude_simulation(dset, e, ep, ew, tnm, eb):
     """probabilities of dielectric function eps relative to zero-loss integral (i0 = 1)

pyTEMlib/info_widget.py

@@ -40,15 +40,15 @@ def get_image_sidebar() -> Any:
                                           layout=ipywidgets.Layout(width='auto', grid_area='header'),
                                           style=ipywidgets.ButtonStyle(button_color='lightblue'))
     row += 1
-    side_bar[row, :2] = ipywidgets.FloatText(value=7.5, description='Conv.Angle:', disabled=False, color='black',
+    side_bar[row, :2] = ipywidgets.FloatText(value=-1, description='Conv.Angle:', disabled=False, color='black',
                                              layout=ipywidgets.Layout(width='200px'))
     side_bar[row, 2] = ipywidgets.widgets.Label(value="mrad", layout=ipywidgets.Layout(width='100px'))
     row += 1
-    side_bar[row, :2] = ipywidgets.FloatText(value=0.1, description='Coll.Angle:', disabled=False, color='black',
+    side_bar[row, :2] = ipywidgets.FloatText(value=-0.1, description='Coll.Angle:', disabled=False, color='black',
                                              layout=ipywidgets.Layout(width='200px'))
     side_bar[row, 2] = ipywidgets.widgets.Label(value="mrad", layout=ipywidgets.Layout(width='100px'))
     row += 1
-    side_bar[row, :2] = ipywidgets.FloatText(value=0.1, description='Acc Voltage:', disabled=False, color='black',
+    side_bar[row, :2] = ipywidgets.FloatText(value=.1, description='Acc Voltage:', disabled=False, color='black',
                                              layout=ipywidgets.Layout(width='200px'))
     side_bar[row, 2] = ipywidgets.widgets.Label(value="keV", layout=ipywidgets.Layout(width='100px'))
 
@@ -123,11 +123,11 @@ def get_info_sidebar() -> Any:
                                           layout=ipywidgets.Layout(width='auto', grid_area='header'),
                                           style=ipywidgets.ButtonStyle(button_color='lightblue'))
     row += 1
-    side_bar[row, :2] = ipywidgets.FloatText(value=7.5, description='Conv.Angle:', disabled=False, color='black',
+    side_bar[row, :2] = ipywidgets.FloatText(value=-1, description='Conv.Angle:', disabled=False, color='black',
                                              layout=ipywidgets.Layout(width='200px'))
     side_bar[row, 2] = ipywidgets.widgets.Label(value="mrad", layout=ipywidgets.Layout(width='100px'))
     row += 1
-    side_bar[row, :2] = ipywidgets.FloatText(value=0.1, description='Coll.Angle:', disabled=False, color='black',
+    side_bar[row, :2] = ipywidgets.FloatText(value=-0.1, description='Coll.Angle:', disabled=False, color='black',
                                              layout=ipywidgets.Layout(width='200px'))
     side_bar[row, 2] = ipywidgets.widgets.Label(value="mrad", layout=ipywidgets.Layout(width='100px'))
     row += 1
@@ -251,6 +251,7 @@ def __init__(self, datasets, sidebar, tab_title=None):
         self.dir_name = file_tools.get_last_path()
 
         self.key = None
+        self.new_info = False
         self.image = 'Sum'
 
         self.save_path = True
@@ -632,6 +633,7 @@ def select_main(self, value=0):
 
         self.key = self.dataset_list[0].split(':')[0]
         self.dataset = self.datasets[self.key]
+        self.new_info = True
 
         self.selected_dataset = self.dataset
         if len(self.image_list) > 0:
@@ -822,9 +824,10 @@ def set_flux(self, value):
         self.info_tab[11, 0].value = np.round(self.parent.datasets[self.parent.key].metadata['experiment']['flux_ppm'], 2)
 
     def set_microscope_parameter(self, value):
-        self.parent.datasets[self.key].metadata['experiment']['convergence_angle'] = self.info_tab[5, 0].value
-        self.parent.datasets[self.key].metadata['experiment']['collection_angle'] = self.info_tab[6, 0].value
-        self.parent.datasets[self.key].metadata['experiment']['acceleration_voltage'] = self.info_tab[7, 0].value*1000
+        if not self.parent.new_info:
+            self.parent.datasets[self.key].metadata['experiment']['convergence_angle'] = self.info_tab[5, 0].value
+            self.parent.datasets[self.key].metadata['experiment']['collection_angle'] = self.info_tab[6, 0].value
+            self.parent.datasets[self.key].metadata['experiment']['acceleration_voltage'] = self.info_tab[7, 0].value*1000
 
     def cursor2energy_scale(self, value):
         self.energy_scale = self.parent.datasets[self.key].get_spectral_dims(return_axis=True)[0]
@@ -869,6 +872,7 @@ def update_sidebar(self):
         spectrum_list = ['None']
         reference_list = ['None']
         data_list = []
+
         self.key = self.info_key = self.parent.info_key
 
         spectrum_data = False
@@ -880,7 +884,6 @@ def update_sidebar(self):
                     if 'SPECTR' in self.parent.datasets[key].data_type.name:
                         spectrum_data = True
                         spectrum_list.append(f'{key}: {self.parent.datasets[key].title}')
-                        print( self.info_key, key)
                         if self.info_key == key:
                             info_index = len(spectrum_list)-1
                 reference_list.append(f'{key}: {self.parent.datasets[key].title}')
@@ -894,8 +897,9 @@ def update_sidebar(self):
         else:
             for i in range(15, 18):
                 self.info_tab[i, 0].layout.display = "flex"
-        
+
         if 'None' not in self.key:
+            self.parent.new_info = True
             energy_scale = self.parent.datasets[self.key].get_spectral_dims(return_axis=True)
             if len(energy_scale) == 0:
                 return
@@ -910,6 +914,7 @@ def update_sidebar(self):
             self.info_tab[5, 0].value = np.round(self.parent.datasets[self.key].metadata['experiment']['convergence_angle'], 1)
             self.info_tab[6, 0].value = np.round(self.parent.datasets[self.key].metadata['experiment']['collection_angle'], 1)
             self.info_tab[7, 0].value = np.round(self.parent.datasets[self.key].metadata['experiment']['acceleration_voltage']/1000, 1)
+            print(self.parent.datasets[self.key].metadata['experiment']['acceleration_voltage'])
             self.info_tab[10, 0].value = np.round(self.parent.datasets[self.key].metadata['experiment']['exposure_time'], 4)
             if 'flux_ppm' not in self.parent.datasets[self.key].metadata['experiment']:
                 self.parent.datasets[self.key].metadata['experiment']['flux_ppm'] = 0
@@ -927,6 +932,7 @@ def update_sidebar(self):
                     ll_key = f'{ll_key}: {self.parent.datasets[ll_key].title}'
                     self.lowloss_key = ll_key
             self.info_tab[9, 0].value = ll_key
+            self.parent.new_info = False
 
     def update_dataset(self, value=0):
         self.key = self.info_tab[0, 0].value.split(':')[0]