Repository: ocelot-collab/ocelot

Issue:
Matching paired quadrupoles
Hello!
Ocelot's `match` function is working well for me to find the quadrupole's strengths to match the Twiss parameters of the beam at a certain element location. 👍
In my script, `match` is varying the strength of 4 quadrupoles to find the optimum.
Let's call these quadrupoles `Q1`, `Q2`, `Q3` and `Q4`.
My question: is it possible to use `match` in a way that the strengths are paired, e.g. `Q1.k1 = Q2.k1` and `Q3.k1 = -Q4.k1`?
A trick that seems to work is to make, for example, `Q2 = Q1`, and remove Q2 from the list of varying elements.
However, this does not work like this for the case where `Q3.k1 = -Q4.k1`, because `Q4 = - Q3` is not allowed.
Any ideas?
Cheers.

Hints from issue comments:
@iagapov: Hi Alberto, 

I think the easiest way would be to implement such things as a constraint, this could be done quickly but is not very effective, since you would have an extra dimension in the optimization problem that is not needed.

What you describe is what is done in mad-x, where you would say something like Q1.k1 := k, Q2.k1 := -k; vary, k in matching etc.


But how to best do it in python without creating ugly syntax, that one should look into, I’ll see if I can come up with something…Or maybe someone else can do it faster.


Ilya




> On Jan 17, 2025, at 09:36, Alberto de la Ossa ***@***.***> wrote:
> 
> 
> Hello!
> Ocelot's match function is working well for me to find the quadrupoles strengths to match the Twiss parameters of the beam at a certain element location. 👍
> In my script, match is varying the strength of 4 quadrupoles to find the optimum.
> Let's call these quadrupoles Q1, Q2, Q3 and Q4.
> My question: is it possible to use match in a way that the strengths are paired, e.g. Q1.k1 = Q2.k1 and Q3.k1 = -Q4.k1?
> A trick that seems to work is to make, for example, Q2 = Q1, and remove Q2 from the list of varying elements.
> However, this does not work like this for the case where Q3.k1 = -Q4.k1, because Q4 = - Q3 is not allowed.
> Any ideas?
> Cheers.
> 
> —
> Reply to this email directly, view it on GitHub <https://github.com/ocelot-collab/ocelot/issues/270>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABVBP2PMXJTKNT4AJCUVXBT2LC6ILAVCNFSM6AAAAABVLMWL6GVHI2DSMVQWIX3LMV43ASLTON2WKOZSG44TIOBXHEZTAMA>.
> You are receiving this because you are subscribed to this thread.
> 

@delaossa: Hello Ilya,
Thanks for your answer and apologies for the delay in getting back to this.
I have just noticed that the `match` function actually allows for coupling a series of quads to the same strength value:
https://github.com/ocelot-collab/ocelot/blob/0d19cf5da2f353a0451f66dc97b1f7d24049c1fa/ocelot/cpbd/match.py#L122
One just needs to add the coupled quads as a tuple to the list of varying parameters for `match`.
However, the way is implemented does not allow for "anti coupling", i.e. quads of opposite strengths.

I think that a similar and more general approach would be to pass a dictionary with the quads as keys and the coupling parameters as values, e.g.:
```py
vars = [{Q1: 1, Q2: 1}, {Q3: 1, Q4: -1}]
```
I have just given it a try in this PR #275 and it seems to work well.




Solution patch:
diff --git a/demos/ebeam/triplet.py b/demos/ebeam/triplet.py
index 68d2cb31..61fc3101 100644
--- a/demos/ebeam/triplet.py
+++ b/demos/ebeam/triplet.py
@@ -98,4 +98,24 @@ plt.xlabel("S, m")
 plt.ylabel("X, m")
 plt.plot(s, x1,'r', s, x2,'b')
 plt.legend(["R", "R+T"])
+plt.show()
+
+# trajectory with transverse offset
+init_coordinates = [-0.01, -0.005, 0, 0.005, 0.01]
+p_list = [Particle(x=q) for q in init_coordinates]
+P = [copy.deepcopy(p_list)]
+navi2 = Navigator(lat2)
+dz = 0.01
+
+for i in range(int(lat1.totalLen/dz)):
+    tracking_step(lat2, p_list, dz=dz, navi=navi2)  # R + T
+    P.append(copy.deepcopy(p_list))
+
+fig, ax_xy = plot_API(lat2, fig_name="Trajectories with X offsets")
+for i in range(len(init_coordinates)):
+    s = [p[i].s for p in P]
+    x = [p[i].x for p in P]
+    ax_xy.plot(s, x)
+ax_xy.set_xlabel("S [m]")
+ax_xy.set_ylabel("X [m]")
 plt.show()
\ No newline at end of file
diff --git a/ocelot/cpbd/match.py b/ocelot/cpbd/match.py
index 5ca5c123..65ea4341 100644
--- a/ocelot/cpbd/match.py
+++ b/ocelot/cpbd/match.py
@@ -68,7 +68,12 @@ def match(lat, constr, vars, tw, verbose=True, max_iter=1000, method='simplex',
                         - difference between ELEM1 and ELEM2 of twiss parameter "muy" (can be any) == "val"
 
     :param vars: list of elements e.g. vars = [QF, QD] or it can be initial twiss parameters:
-                vars = [[tws0, 'beta_x'], [tws0, 'beta_y'], [tws0, 'alpha_x'], [tws0, 'alpha_y']]
+                vars = [[tws0, 'beta_x'], [tws0, 'beta_y'], [tws0, 'alpha_x'], [tws0, 'alpha_y']].
+                A tuple of quadrupoles can be passed as a variable to constrain their strengths
+                to the same value, e.g., vars = [(QF, QD)].
+                A dictionary with quadrupoles as keys and their relative strengths as values
+                can be used for more flexible constraints, e.g., vars = [{QF: 1.0, QD: -1.0}],
+                which constrains QD and QF to have equal strengths with opposite signs.
     :param tw: initial Twiss
     :param verbose: allow print output of minimization procedure
     :param max_iter:
@@ -116,12 +121,17 @@ def match(lat, constr, vars, tw, verbose=True, max_iter=1000, method='simplex',
                     vars[i].k1 = x[i]
             if isinstance(vars[i], list):
                 if isinstance(vars[i][0], Twiss) and isinstance(vars[i][1], str):
-
                     k = vars[i][1]
                     tw_loc.__dict__[k] = x[i]
-            if isinstance(vars[i], tuple):  # all quads strength in tuple varied simultaneously
+            if isinstance(vars[i], tuple):
+            # all quads strength in tuple varied simultaneously
                 for v in vars[i]:
                     v.k1 = x[i]
+            if isinstance(vars[i], dict):
+            # all quads strength in dict keys varied simultaneously
+            # with coupling parameters given as values.
+                for q in vars[i].keys():
+                    q.k1 = vars[i][q] * x[i]
 
         err = 0.0
         if "periodic" in constr.keys():
@@ -256,6 +266,9 @@ def match(lat, constr, vars, tw, verbose=True, max_iter=1000, method='simplex',
                 #     x[i] = 0.0
         if vars[i].__class__ == tuple:
             x[i] = vars[i][0].k1
+        if vars[i].__class__ == dict:
+            q = list(vars[i].keys())[0]
+            x[i] = q.k1 / vars[i][q]
         if vars[i].__class__ == Quadrupole:
             x[i] = vars[i].k1
         if vars[i].__class__ == Drift:
diff --git a/ocelot/gui/accelerator.py b/ocelot/gui/accelerator.py
index 43aaa261..49fd1c79 100644
--- a/ocelot/gui/accelerator.py
+++ b/ocelot/gui/accelerator.py
@@ -857,7 +857,7 @@ def show_e_beam(p_array, nparts_in_slice=5000, smooth_param=0.05, nbins_x=200, n
                 title=None, figsize=None, grid=True,
                 filename=None, headtail=True,
                 filter_base=2, filter_iter=2,
-                tau_units="mm", cmap="my_rainbow", **kwargs
+                tau_units="mm", p_units=None, cmap="my_rainbow", **kwargs
                 ):
     """
     Shows e-beam slice parameters (current, emittances, energy spread)
@@ -881,7 +881,7 @@ def show_e_beam(p_array, nparts_in_slice=5000, smooth_param=0.05, nbins_x=200, n
     :param headtail: True, shows where is the beam head is.
     :param filter_base: support of rectangle filter is 2*p+1
     :param filter_iter: the number of the filter iterations
-    :param tau_units: unints of longitudinal coordinate tau - 'm', 'mm' or 'um'
+    :param tau_units: unints of longitudinal coordinate tau - 'm', 'mm', 'um', "ps", "fs
     :param cmap: color map
     :return:
     """
@@ -891,10 +891,26 @@ def show_e_beam(p_array, nparts_in_slice=5000, smooth_param=0.05, nbins_x=200, n
     elif tau_units == "um":
         tau_factor = 1e6  # um
         tau_label = r"$s\,[\mathrm{\mu{}m}]$"
+    elif tau_units == "ps":
+        tau_factor = 1e12/speed_of_light  # pm
+        tau_label = r"$t\,[\mathrm{ps}]$"
+    elif tau_units == "fs":
+        tau_factor = 1e15/speed_of_light  # fm
+        tau_label = r"$t\,[\mathrm{fs}]$"
     else:
         tau_factor = 1e3  # mm
         tau_label = r"$s\,[\mathrm{mm}]$"
 
+    if p_units == "MeV":
+        pref = np.sqrt(p_array.E ** 2 / m_e_GeV ** 2 - 1) * m_e_GeV
+        #Enew = p_array.p()[0] * pref + p_array.E
+        p_factor = pref*1e3
+        p_label = r'$E\,[MeV]$'
+    else:
+        p_factor = 1e2
+        p_label = r'$\delta_E\,[\%]$'
+
+
     p_array_copy = deepcopy(p_array)
     if inverse_tau:
         p_array_copy.tau()[:] *= -1
@@ -963,8 +979,8 @@ def show_e_beam(p_array, nparts_in_slice=5000, smooth_param=0.05, nbins_x=200, n
 
     ax_ps = plt.subplot(322, sharex=ax_sp)
 
-    show_density(p_array_copy.tau() * tau_factor, p_array_copy.p() * 1e2, ax=ax_ps, nbins_x=nbins_x, nbins_y=nbins_y,
-                 interpolation=interpolation, ylabel=r'$\delta_E\,[\%]$',
+    show_density(p_array_copy.tau() * tau_factor, p_array_copy.p() * p_factor, ax=ax_ps, nbins_x=nbins_x, nbins_y=nbins_y,
+                 interpolation=interpolation, ylabel=p_label,
                  title="Longitudinal phase space", grid=grid, show_xtick_label=False, cmap=cmap)
 
     if inverse_tau:
diff --git a/ocelot/utils/image_analysis.py b/ocelot/utils/image_analysis.py
index 8320ee66..ae69d864 100644
--- a/ocelot/utils/image_analysis.py
+++ b/ocelot/utils/image_analysis.py
@@ -1,211 +1,177 @@
 import numpy as np
 from scipy import ndimage
 from scipy.signal import find_peaks
-from scipy.ndimage import gaussian_filter, uniform_filter
+from scipy.ndimage import gaussian_filter
 from scipy.optimize import curve_fit
 
 
 def simple_filter_and_mask(image, sigma=5, threshold=0.0):
     """
-    The function returns the same image but with the background cut (masked).
-    First, the function filters the original image with a Gaussian filter (sigma=5 by default).
-    Second, a mask is generated corresponding to areas where the intensity of the filtered image is less than the 'threshold'.
-    Lastly, the mask is applied to the original image by zeroing out the masked area.
-
-    :param image: 2D numpy array - particle density distribution
-    :param sigma: gaussian filter 5 pixes by default
-    :param threshold: below threshold the area will be set to zero
-    :return: masked image with the same sizes
-    """
+    Applies a Gaussian filter to the image and masks (sets to zero) pixels below a dynamic threshold.
+
+    Parameters:
+        image (np.ndarray): Input image.
+        sigma (float): Standard deviation for Gaussian kernel.
+        threshold (float): Initial relative threshold for masking.
 
+    Returns:
+        np.ndarray: Image with masked low-intensity background.
+    """
     Hg = ndimage.gaussian_filter(image, sigma=sigma, truncate=2)
     inds_mask_neg = ((Hg - np.max(Hg) * threshold) < 0).nonzero()
     for i in range(10):
-        img_neg = np.copy(image)
-        img_neg[:] = 1
+        img_neg = np.ones_like(image)
         img_neg[inds_mask_neg] = 0.0
         proj = np.sum(img_neg, axis=0)
-        if np.all(proj[0:20] == 0) and np.all(proj[-20:] == 0):
+        if np.all(proj[:20] == 0) and np.all(proj[-20:] == 0):
             break
         threshold += 0.01
-        print(f"threshold = {threshold}")
         inds_mask_neg = ((Hg - np.max(Hg) * threshold) < 0).nonzero()
-
     image[inds_mask_neg] = 0.0
-
     return image
 
 
 def simple_filter_and_mask_v2(image, sigma=5, threshold=0.0, max_iterations=10):
     """
-    The function returns the same image but with the background cut (masked).
-    1) The function filters the original image with a Gaussian filter (sigma=5 by default).
-    2) A mask is generated for areas where intensity < threshold * max(Hg).
-    3) threshold is iteratively adjusted so that the first and last 20 columns become fully masked.
-    4) The mask is then applied in-place (zeroing out masked areas) to the original image.
-
-    :param image: 2D numpy array - particle density distribution
-    :param sigma: Gaussian filter width (5 pixels by default)
-    :param threshold: Initial threshold fraction of the max intensity in the filtered image
-    :param max_iterations: Maximum number of increments to find a threshold that masks edges
-    :return: The modified (in-place) masked image
-    """
+    Enhanced version of the filtering and masking function with iterative threshold adjustment.
 
-    # 1) Filter the image
+    Parameters:
+        image (np.ndarray): Input image.
+        sigma (float): Standard deviation for Gaussian filter.
+        threshold (float): Initial relative threshold.
+        max_iterations (int): Number of iterations to adjust the threshold.
+
+    Returns:
+        np.ndarray: Image with background masked out.
+    """
     Hg = ndimage.gaussian_filter(image, sigma=sigma, truncate=2)
     max_val = np.max(Hg)
-
-    # 2) Iteratively adjust threshold until the edges are fully masked or we hit max_iterations
     for _ in range(max_iterations):
         mask = (Hg < threshold * max_val)
-        # Check if the first 20 and last 20 columns are fully masked
         if mask[:, :20].all() and mask[:, -20:].all():
             break
         threshold += 0.01
-        # Optionally: print the updated threshold if you want to track it
-        # print(f"threshold = {threshold}")
-
-    # 3) Create the final mask after adjustments
-    mask = (Hg < threshold * max_val)
-
-    # 4) Zero out the masked area in-place
-    image[mask] = 0.0
-
+    image[Hg < threshold * max_val] = 0.0
     return image
 
 
 def get_image_proj(image):
     """
-    function gets projection on horizontal asis
-    :param image: 2D numpy array
-    :return: 1D numpy array - projection on horizontal axis
+    Computes the horizontal projection (sum along vertical axis) of the image.
+
+    Parameters:
+        image (np.ndarray): Input image.
+
+    Returns:
+        np.ndarray: 1D array of summed pixel values along columns.
     """
-    proj = np.sum(image, axis=0)
-    return proj
+    return np.sum(image, axis=0)
 
 
 def find_peaks_regions(y):
     """
-    Finds peaks in the 1D array y and returns the indices of the regions around the peaks.
-    """
-    # Find peaks
-    peaks, properties = find_peaks(y, prominence=0.05 * np.max(y))
+    Finds prominent peaks and returns the base regions around them.
 
-    # You can adjust 'prominence' to suit your data
-    # Alternatively, use 'height', 'distance', 'width', etc.
+    Parameters:
+        y (np.ndarray): 1D input signal.
 
+    Returns:
+        tuple or list of tuples: Left and right indices of the most prominent peak or peaks.
+    """
+    peaks, properties = find_peaks(y, prominence=0.05 * np.max(y))
     regions = []
     for peak in peaks:
-        # Find the left and right bases of the peak
-        left_base = properties['left_bases'][np.where(peaks == peak)[0][0]]
-        right_base = properties['right_bases'][np.where(peaks == peak)[0][0]]
-        regions.append((left_base, right_base))
-    #print("find_peaks_regions ", regions, len(regions))
+        left = properties['left_bases'][np.where(peaks == peak)[0][0]]
+        right = properties['right_bases'][np.where(peaks == peak)[0][0]]
+        regions.append((left, right))
     if len(regions) == 1:
         return regions[0]
     elif len(regions) == 0:
         return regions
-    dists = [r[1] - r[0] for r in regions]
-    ind = np.argmax(dists)
-    return regions[ind]
+    return regions[np.argmax([r[1] - r[0] for r in regions])]
 
 
 def roi_1d(y, **kwargs):
-    threshold = kwargs.get("threshold", 0)
-    lmargin = kwargs.get("lmargin", 0)
-    rmargin = kwargs.get("rmargin", 0)
-    #regions = roi_1d_simple(y, threshold=threshold, lmargin=lmargin, rmargin=rmargin)
-    regions = find_peaks_regions(y)
-    return regions
+    """
+    Returns region of interest (ROI) from a 1D signal using peak detection.
+
+    Parameters:
+        y (np.ndarray): 1D signal.
+        **kwargs: Additional arguments passed to `find_peaks_regions`.
+
+    Returns:
+        tuple: Indices defining the ROI.
+    """
+    return find_peaks_regions(y)
 
 
 def crop_1d(x, y, threshold=0.01, lmargin=0, rmargin=0):
     """
-    function crops current distribution below thresholds + margins
-    :param x:
-    :param y:
-    :param threshold: find indices when yi == max(y)*threshold
-    :param lmargin: add number of pixes (indices) on the left side
-    :param rmargin: add number of pixes (indices) on the right side
-    :return: x[indx1:indx2], y[indx1:indx2]
+    Crops 1D data based on signal region of interest.
+
+    Parameters:
+        x (np.ndarray): X-axis values.
+        y (np.ndarray): Y-axis values (signal).
+        threshold (float): Threshold for ROI detection.
+        lmargin (int): Left margin to extend ROI.
+        rmargin (int): Right margin to extend ROI.
+
+    Returns:
+        tuple: Cropped (x, y) arrays.
     """
     indx1, indx2 = roi_1d(y, threshold=threshold, lmargin=lmargin, rmargin=rmargin)
-    x_new = x[indx1:indx2]
-    y_new = y[indx1:indx2]
-    return x_new, y_new
+    return x[indx1:indx2], y[indx1:indx2]
 
 
 def roi_2d(image, threshold=0.01, hor_margin=[5, 5], ver_margin=[5, 5]):
     """
-    function finds indices of image below threshold.
-    :param image:
-    :param threshold:
-    :param hor_margin: [left margin, right margin]
-    :param ver_margin: [bottoь margin, top margin]
-    :return:
-    """
-    h_proj = get_image_proj(image)
-    h_proj = gaussian_filter(h_proj, sigma=2)
-    #h_proj = uniform_filter(h_proj, size=2)
-    ix1, ix2 = roi_1d(y=h_proj, threshold=threshold, lmargin=hor_margin[0], rmargin=hor_margin[1])
+    Computes 2D region of interest using horizontal and vertical projections.
 
-    v_proj = get_image_proj(image.T)
-    v_proj = gaussian_filter(v_proj, sigma=2)
-    #v_proj = uniform_filter(v_proj, size=2)
-    iy1, iy2 = roi_1d(y=v_proj, threshold=threshold, lmargin=ver_margin[0], rmargin=ver_margin[1])
+    Parameters:
+        image (np.ndarray): 2D input image.
+        threshold (float): Threshold for ROI detection.
+        hor_margin (list): Left and right horizontal margins.
+        ver_margin (list): Top and bottom vertical margins.
 
+    Returns:
+        tuple: Indices (ix1, ix2, iy1, iy2) for cropping the image.
+    """
+    h_proj = gaussian_filter(get_image_proj(image), sigma=2)
+    ix1, ix2 = roi_1d(h_proj, threshold=threshold, lmargin=hor_margin[0], rmargin=hor_margin[1])
+    v_proj = gaussian_filter(get_image_proj(image.T), sigma=2)
+    iy1, iy2 = roi_1d(v_proj, threshold=threshold, lmargin=ver_margin[0], rmargin=ver_margin[1])
     return ix1, ix2, iy1, iy2
 
 
 def crop_2d(image, threshold=0.01, hor_margin=[5, 5], ver_margin=[5, 5]):
     """
-    function crop image
+    Crops 2D image to region of interest based on projections.
+
+    Parameters:
+        image (np.ndarray): Input image.
+        threshold (float): Threshold for ROI detection.
+        hor_margin (list): Horizontal margins for cropping.
+        ver_margin (list): Vertical margins for cropping.
 
-    :param image:
-    :param threshold:
-    :param hor_margin:
-    :param ver_margin:
-    :return:
+    Returns:
+        np.ndarray: Cropped image.
     """
     ix1, ix2, iy1, iy2 = roi_2d(image, threshold=threshold, hor_margin=hor_margin, ver_margin=ver_margin)
-    #print("crop = ", ix1, ix2, iy1, iy2)
     return image[iy1:iy2, ix1:ix2]
 
 
 def gauss_fit(x, y, sigma_estm=None, mu_estm=0.0):
-    # TODO: create unit tests
     """
-    Fit a Gaussian function to the given data.
-
-    Parameters
-    ----------
-    x : array_like
-        Independent variable data (e.g., positions).
-    y : array_like
-        Dependent variable data (e.g., intensities at positions x).
-    sigma_estm : float, optional
-        Initial estimate of the standard deviation (sigma) of the Gaussian.
-        If None, it is estimated as one-fourth of the x-range.
-    mu_estm : float, optional
-        Initial estimate of the mean (mu) of the Gaussian. Default is 0.0.
-
-    Returns
-    -------
-    A : float
-        Amplitude of the fitted Gaussian function.
-    mu : float
-        Mean (center position) of the fitted Gaussian function.
-    sigma : float
-        Standard deviation (spread or width) of the fitted Gaussian function.
-
-    Notes
-    -----
-    The Gaussian function is defined as:
-        f(x) = A * exp(- (x - mu)^2 / (2 * sigma^2))
-    where:
-        - A is the amplitude,
-        - mu is the mean,
-        - sigma is the standard deviation.
+    Fits a Gaussian function to 1D data.
+
+    Parameters:
+        x (np.ndarray): X values.
+        y (np.ndarray): Y values.
+        sigma_estm (float): Initial estimate for standard deviation.
+        mu_estm (float): Initial estimate for mean.
+
+    Returns:
+        list: Fitted parameters [A, mu, sigma].
     """
     if sigma_estm is None:
         sigma_estm = (np.max(x) - np.min(x)) / 4.0
@@ -213,95 +179,107 @@ def gauss_fit(x, y, sigma_estm=None, mu_estm=0.0):
     def gauss(x, A, mu, sigma):
         return A * np.exp(- (x - mu) ** 2 / (2.0 * sigma ** 2))
 
-    # Initial guess for the fitting coefficients: [A, mu, sigma]
     p0 = [np.max(y), mu_estm, sigma_estm]
-
     coeff, _ = curve_fit(gauss, x, y, p0=p0)
-    A, mu, sigma = coeff
-    return A, mu, sigma
+    return coeff
 
 
 def get_slice_parameters(image, gauss_fit_enabled=True):
-    # TODO: create unit tests
     """
-    Calculate slice parameters of an image spot, optionally fitting Gaussian profiles.
-
-    Parameters
-    ----------
-    image : 2D array_like
-        The input image containing the spot (intensity values).
-    gauss_fit_enabled : bool, optional
-        If True, perform Gaussian fitting on each column (slice) of the image.
-        Default is True.
-
-    Returns
-    -------
-    x : ndarray
-        Array of x-axis coordinates (column indices).
-    y : ndarray
-        Array of y-axis coordinates (row indices).
-    proj_x : ndarray
-        Projection of the image along the y-axis (sum over rows for each column).
-    proj_y : ndarray
-        Projection of the image along the x-axis (sum over columns for each row).
-    slice_energy : ndarray
-        Weighted mean position (energy) along the y-axis for each column.
-    slice_energy_spread : ndarray
-        Standard deviation (spread) of the energy along the y-axis for each column.
-    x_sigma : float
-        Standard deviation of the spot along the x-axis.
-
-    Notes
-    -----
-    - The function computes the centroid and spread of the intensity distribution
-      for each column (x-axis) of the image.
-    - If `gauss_fit_enabled` is True, it refines the centroid and spread by fitting
-      a Gaussian function to the intensity profile of each column.
+    Calculates image slice parameters like mean, sigma, and energy distribution.
+
+    Parameters:
+        image (np.ndarray): 2D input image.
+        gauss_fit_enabled (bool): If True, performs Gaussian fitting per column.
+
+    Returns:
+        tuple: x, y, proj_x, proj_y, slice_energy, slice_energy_spread, x_sigma
     """
     ny, nx = image.shape
     x = np.arange(nx)
     y = np.arange(ny)
-
-    # Projections along x and y axes
-    proj_x = np.sum(image, axis=0)  # Sum over rows (y-axis) for each column
-    proj_y = np.sum(image, axis=1)  # Sum over columns (x-axis) for each row
-
+    proj_x = np.sum(image, axis=0)
+    proj_y = np.sum(image, axis=1)
     total_intensity = np.sum(proj_y)
-
-    # Mean position along x-axis
     x_mean = np.sum(x * proj_x) / total_intensity
-
-    # Standard deviation along x-axis
+    print("X_MEAN = ", x_mean)
     x_variance = np.sum((x - x_mean) ** 2 * proj_x) / total_intensity
     x_sigma = np.sqrt(x_variance)
 
-    # Weighted mean position (slice energy) along y-axis for each column
     slice_energy = np.sum(image * y[:, np.newaxis], axis=0) / proj_x
-
-    # Compute the squared differences from the mean for each column
     y_diff_squared = (y[:, np.newaxis] - slice_energy) ** 2
-
-    # Weighted sum of squared differences (variance) for each column
     y_variance = np.sum(image * y_diff_squared, axis=0) / proj_x
-
-    # Standard deviation (spread) along y-axis for each column
     slice_energy_spread = np.sqrt(y_variance)
 
     if gauss_fit_enabled:
-        # Refine slice energy and spread by Gaussian fitting
         for i in range(nx):
             try:
-                A, mu, sigma = gauss_fit(
-                    y,
-                    image[:, i],
-                    sigma_estm=slice_energy_spread[i],
-                    mu_estm=slice_energy[i]
-                )
+                A, mu, sigma = gauss_fit(y, image[:, i], sigma_estm=slice_energy_spread[i], mu_estm=slice_energy[i])
                 slice_energy[i] = mu
                 slice_energy_spread[i] = sigma
-            except RuntimeError as e:
-                print(f"Gaussian fit failed for slice {i}: {e}")
             except Exception as e:
-                print(f"An error occurred for slice {i}: {e}")
+                print(f"Fit failed at column {i}: {e}")
 
     return x, y, proj_x, proj_y, slice_energy, slice_energy_spread, abs(x_sigma)
+
+
+def center_images(images):
+    """
+    Centers a list of images based on their center of mass.
+
+    Parameters:
+        images (list of np.ndarray): List of 2D images.
+
+    Returns:
+        list of np.ndarray: Centered images.
+    """
+    if len(images) <= 1:
+        return images
+
+    jj, ii = [], []
+    nyy, nxx = [], []
+    for img in images:
+        j, i = ndimage.center_of_mass(img)
+        ny, nx = img.shape
+        jj.append(int(j))
+        ii.append(int(i))
+        nyy.append(ny)
+        nxx.append(nx)
+
+    j0 = np.array(jj) - np.min(jj)
+    j1 = np.array(jj) + np.min(np.array(nyy) - np.array(jj))
+    i0 = np.array(ii) - np.min(ii)
+    i1 = np.array(ii) + np.min(np.array(nxx) - np.array(ii))
+
+    return [images[i][j0[i]:j1[i], i0[i]:i1[i]] for i in range(len(images))]
+
+
+def averaging_images(images):
+    """
+    Averages a list of images.
+
+    Parameters:
+        images (list of np.ndarray): List of 2D images.
+
+    Returns:
+        np.ndarray: Averaged image.
+    """
+    return np.mean(images, axis=0) if len(images) > 1 else images[0]
+
+
+def shearing(image, shear):
+    """
+    Applies horizontal shearing transformation to an image.
+
+    Parameters:
+        image (np.ndarray): Input image.
+        shear (float): Shear angle in radians.
+
+    Returns:
+        np.ndarray: Sheared image.
+    """
+    X, Y = np.meshgrid(np.arange(image.shape[1]), np.arange(image.shape[0]))
+    coords = np.vstack((X.flatten(), Y.flatten()))
+    transform = np.array([[1, np.tan(shear)], [0, 1]])
+    new_coords = np.dot(transform, coords)
+    return ndimage.map_coordinates(image, coordinates=[new_coords[1], new_coords[0]], order=1).reshape(image.shape)
\ No newline at end of file


Test patch:
diff --git a/unit_tests/ebeam_test/matching/match_test.py b/unit_tests/ebeam_test/matching/match_test.py
index 7a1582f0..bc60d39a 100644
--- a/unit_tests/ebeam_test/matching/match_test.py
+++ b/unit_tests/ebeam_test/matching/match_test.py
@@ -111,6 +111,37 @@ def test_quad_match(lattice, lattice_inj=None, update_ref_values=False):
     result = check_dict(tws, tws_ref, TOL, 'absotute', assert_info=' tws after matching - ')
     assert check_result(result)
 
+def test_quad_match_with_dict(lattice, lattice_inj=None, update_ref_values=False):
+    """After matching R maxtrix calculcation test"""
+    tws0 = Twiss()
+    tws0.beta_x = 8.4
+    tws0.beta_y = 8.4
+    tws0.alpha_x = -55.8
+    tws0.alpha_y = -55.8
+    tws0.E = 0.005071
+    q1_k1 = q1.k1
+    q2_k1 = q2.k1
+    print(q1_k1)
+    # match beta functions at the end
+    constr = {end: {"beta_x": 10, "beta_y": 10}}
+    vars = [{q1: 1, q2: -1}]
+
+    match(lattice, constr, vars, tws0, verbose=False)
+
+    tws = twiss(lattice, tws0, nPoints=20)
+
+    tws = obj2dict(tws)
+    q1.k1 = q1_k1
+    q2.k1 = q2_k1
+    lattice.update_transfer_maps()
+    if update_ref_values:
+        return tws
+
+    tws_ref = json_read(REF_RES_DIR + sys._getframe().f_code.co_name + '.json')
+
+    result = check_dict(tws, tws_ref, TOL, 'absotute', assert_info=' tws after matching - ')
+    assert check_result(result)
+
 def test_drift_match(lattice, lattice_inj=None, update_ref_values=False):
     """After matching R maxtrix calculcation test"""
     d = Drift(l=1)
@@ -326,6 +357,7 @@ def test_update_ref_values(lattice, lattice_inj, cmdopt):
     update_functions.append('test_twiss')
     update_functions.append('test_solenoid_match')
     update_functions.append('test_quad_match')
+    update_functions.append('test_quad_match_with_dict')
     update_functions.append('test_bend_k_match')
     update_functions.append('test_bend_angle_match')
     update_functions.append('test_inj_lattice')
diff --git a/unit_tests/ebeam_test/matching/ref_results/test_quad_match_with_dict.json b/unit_tests/ebeam_test/matching/ref_results/test_quad_match_with_dict.json
new file mode 100644
index 00000000..d58333a3
--- /dev/null
+++ b/unit_tests/ebeam_test/matching/ref_results/test_quad_match_with_dict.json
@@ -0,0 +1 @@
+[{"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 8.4, "beta_y": 8.4, "alpha_x": -55.8, "alpha_y": -55.8, "gamma_x": 370.79047619047617, "gamma_y": 370.79047619047617, "Dx": 0.0, "Dy": 0.0, "Dxp": 0.0, "Dyp": 0.0, "mux": 0.0, "muy": 0.0, "E": 0.005071, "s": 0.0, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 22.300501780767704, "beta_y": 22.300501780767704, "alpha_x": -90.92751879699247, "alpha_y": -90.92751879699247, "gamma_x": 370.79047619047617, "gamma_y": 370.79047619047617, "Dx": 0.0, "Dy": 0.0, "Dxp": 0.0, "Dyp": 0.0, "mux": 0.006921901446994267, "muy": 0.006921901446994267, "E": 0.005071, "s": 0.09473684210526315, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 37.52370631231519, "beta_y": 37.52370631231519, "alpha_x": -86.98785387916804, "alpha_y": -86.98785387916804, "gamma_x": 201.6828151119957, "gamma_y": 201.6828151119957, "Dx": 0.0, "Dy": 0.0, "Dxp": 0.0, "Dyp": 0.0, "mux": 0.010037647487170388, "muy": 0.010037647487170388, "E": 0.005071, "s": 0.18947368421052627, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 40.47777959300505, "beta_y": 40.47777959300505, "alpha_x": -47.61635583887656, "alpha_y": -47.61635583887656, "gamma_x": 56.038581319971925, "gamma_y": 56.038581319971925, "Dx": 0.0, "Dy": 0.0, "Dxp": 0.0, "Dyp": 0.0, "mux": 0.012005817134821491, "muy": 0.012005817134821491, "E": 0.005071, "s": 0.28421052631578947, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 45.968940693433616, "beta_y": 45.968940693433616, "alpha_x": -43.209882924477355, "alpha_y": -43.209882924477355, "gamma_x": 40.638177738428624, "gamma_y": 40.638177738428624, "Dx": 0.0, "Dy": 0.0, "Dxp": 0.0, "Dyp": 0.0, "mux": 0.014138526776139458, "muy": 0.014138526776139458, "E": 0.005071, "s": 0.3789473684210526, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 50.10446237975667, "beta_y": 59.1916497920312, "alpha_x": 13.786052567712044, "alpha_y": -114.90494267422494, "gamma_x": 3.813138317933743, "gamma_y": 223.07446907392256, "Dx": 0.0, "Dy": 0.0, "Dxp": 0.0, "Dyp": 0.0, "mux": 0.016071424501463507, "muy": 0.015992952737053978, "E": 0.005071, "s": 0.47368421052631565, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 40.90143207544426, "beta_y": 94.1399123765844, "alpha_x": 78.65490796389417, "alpha_y": -271.5104560914535, "gamma_x": 151.2806332892063, "gamma_y": 783.0783554598393, "Dx": 0.0, "Dy": 0.0, "Dxp": 0.0, "Dyp": 0.0, "mux": 0.01811434623470734, "muy": 0.017292709579838077, "E": 0.005071, "s": 0.5684210526315789, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 24.82697111319107, "beta_y": 162.78336311386937, "alpha_x": 76.87026521687012, "alpha_y": -422.6946586655822, "gamma_x": 238.04908168487097, "gamma_y": 1097.6046387457268, "Dx": 0.0, "Dy": 0.0, "Dxp": 0.0, "Dyp": 0.0, "mux": 0.021068850338526613, "muy": 0.018062763013096998, "E": 0.005071, "s": 0.6631578947368421, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 13.193021293691565, "beta_y": 240.4498483665515, "alpha_x": 42.53521910178664, "alpha_y": -300.69593471879267, "gamma_x": 137.21230518308872, "gamma_y": 376.0411818541468, "Dx": 0.0, "Dy": 0.0, "Dxp": 0.0, "Dyp": 0.0, "mux": 0.02638711628544802, "muy": 0.018534024479257897, "E": 0.005071, "s": 0.7578947368421051, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 7.620224270087071, "beta_y": 262.3786400087272, "alpha_x": 19.081995170493315, "alpha_y": 80.42951206555958, "gamma_x": 47.91493356960716, "gamma_y": 24.658662804597174, "Dx": 0.0, "Dy": 0.0, "Dxp": 0.0, "Dyp": 0.0, "mux": 0.03606387421681277, "muy": 0.01890222330963198, "E": 0.005071, "s": 0.8526315789473682, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 5.051782419678054, "beta_y": 220.941747106626, "alpha_x": 11.065216567909658, "alpha_y": 254.13946009964354, "gamma_x": 24.4347454898125, "gamma_y": 292.32983818385554, "Dx": 0.0, "Dy": 0.0, "Dxp": 0.0, "Dyp": 0.0, "mux": 0.05151705067888613, "muy": 0.01929099417030921, "E": 0.005071, "s": 0.9473684210526313, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 3.1684438564019235, "beta_y": 175.73513365169165, "alpha_x": 8.883556971489327, "alpha_y": 219.0573185382315, "gamma_x": 25.22297635295669, "gamma_y": 273.06496889956577, "Dx": 4.454148795067312e-05, "Dy": 0.0, "Dxp": 0.002115407315099699, "Dyp": 0.0, "mux": 0.07518401146541957, "muy": 0.01977207891779319, "E": 0.005071, "s": 1.0421052631578946, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 1.6925989695281154, "beta_y": 138.0208235466544, "alpha_x": 6.6480555316565475, "alpha_y": 180.2284432298867, "gamma_x": 26.702510852047677, "gamma_y": 235.35065879452853, "Dx": 0.00046980351033437814, "Dy": 0.0, "Dxp": 0.00685561403095432, "Dyp": 0.0, "mux": 0.11604306005449305, "muy": 0.020381288560347703, "E": 0.005071, "s": 1.1368421052631579, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 0.6650168981624808, "beta_y": 106.91882204763286, "alpha_x": 4.153259782621912, "alpha_y": 151.2656827627013, "gamma_x": 27.442260298001894, "gamma_y": 214.01570222567554, "Dx": 0.0013170499348952756, "Dy": 0.0, "Dxp": 0.009984121072529538, "Dyp": 0.0, "mux": 0.20535726343439098, "muy": 0.021162018149062184, "E": 0.005071, "s": 1.2315789473684209, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 0.12437965227380282, "beta_y": 80.17876158847056, "alpha_x": 1.5534667017585728, "alpha_y": 130.99051097290044, "gamma_x": 27.44226029800188, "gamma_y": 214.01570222567557, "Dx": 0.002262914036503338, "Dy": 0.0, "Dxp": 0.009984121072529538, "Dyp": 0.0, "mux": 0.5410287570627665, "muy": 0.022185222139589567, "E": 0.005071, "s": 1.3263157894736841, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 0.0763620386839955, "beta_y": 57.25424045688666, "alpha_x": -1.047844621269505, "alpha_y": 111.92020395704111, "gamma_x": 27.474100829148515, "gamma_y": 218.79832749189654, "Dx": 0.0031983190643804778, "Dy": 0.0, "Dxp": 0.008991208128262448, "Dyp": 0.0, "mux": 2.347868665002464, "muy": 0.02358321050024861, "E": 0.005071, "s": 1.4210526315789471, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 0.5240924800649228, "beta_y": 37.62049965663021, "alpha_x": -3.6922940251464986, "alpha_y": 94.70445366145444, "gamma_x": 27.920711944425996, "gamma_y": 238.432068292153, "Dx": 0.003840023655718392, "Dy": 0.0, "Dxp": 0.0045659941395693295, "Dyp": 0.0, "mux": 2.8420012747357077, "muy": 0.025621440853202268, "E": 0.005071, "s": 1.5157894736842104, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 1.4837540402200942, "beta_y": 21.583032211649716, "alpha_x": -6.453792108479431, "alpha_y": 73.871872471687, "gamma_x": 28.745621864082427, "gamma_y": 252.8863177773109, "Dx": 0.004067306013441289, "Dy": 0.0, "Dxp": 0.0006680359987834266, "Dyp": 0.0, "mux": 2.9497968626730655, "muy": 0.02894133046690596, "E": 0.005071, "s": 1.6105263157894736, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 2.9645717550666846, "beta_y": 9.855928589863296, "alpha_x": -9.177061548234605, "alpha_y": 49.914221313836535, "gamma_x": 28.74562186408243, "gamma_y": 252.8863177773109, "Dx": 0.004130593634378666, "Dy": 0.0, "Dxp": 0.0006680359987834266, "Dyp": 0.0, "mux": 2.9949829526552545, "muy": 0.03543689578930288, "E": 0.005071, "s": 1.7052631578947366, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}, {"emit_x": 0.0, "emit_y": 0.0, "emit_xn": 0.0, "emit_yn": 0.0, "eigemit_1": 0.0, "eigemit_2": 0.0, "beta_x": 4.961377363761635, "beta_y": 2.668169397985354, "alpha_x": -11.90033098798979, "alpha_y": 25.95657015598593, "gamma_x": 28.745621864082423, "gamma_y": 252.8863177773105, "Dx": 0.0041938812553160435, "Dy": 0.0, "Dxp": 0.0006680359987834266, "Dyp": 0.0, "mux": 3.0196877313718917, "muy": 0.053912052774346894, "E": 0.005071, "s": 1.8, "q": 0.0, "x": 0.0, "y": 0.0, "p": 0.0, "tau": 0.0, "xp": 0.0, "yp": 0.0, "xx": 0.0, "xpx": 0.0, "pxpx": 0.0, "yy": 0.0, "ypy": 0.0, "pypy": 0.0, "tautau": 0.0, "xy": 0.0, "pxpy": 0.0, "xpy": 0.0, "ypx": 0.0, "pp": 0}]
\ No newline at end of file