I have this code that aims to fit the data in here DriveFilelink

The function read_file is utilized to extract information in a structured manner. However, I’m encountering challenges with the Gaussian fit, evident from the discrepancies observed in the Gaussian fit image. This issue appears to stem from the constraints imposed on certain parameters, such as fixing the mean of all the Gaussians and three out of the four amplitudes. Despite these constraints, they are necessary as they are based on available information.

Does anyone know how I can fix this problem?

The code is the following:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from scipy.optimize import curve_fit
from google.colab import drive
import os
# Mount Google Drive
drive.mount('/content/drive')


# Function to read numeric data from a file
def read_numeric_data(file_path):
    with open(file_path, 'r', encoding='latin1') as f:
        data = []
        live_time = None
        real_time = None
        for line in f:
            if 'LIVE_TIME' in line:
                live_time = line.strip()
            elif 'REAL_TIME' in line:
                real_time = line.strip()
            elif all(char.isdigit() or char in ".-+e" for char in line.strip()):
                row = [float(x) for x in line.split()]
                data.extend(row)
    return np.array(data), live_time, real_time



file_path = '/content/drive/MyDrive/ProjetoXRF_Manuel/April2024/20240402_In.mca'
data, _, _ = read_numeric_data(file_path)

a = -0.0188026396003431
b = 0.039549044037714
Data = data





# Función para convolucionar múltiples gaussianas
def multi_gaussian(x, *params):
    
    eps = 1e-10
    y = np.zeros_like(x)
    amplitude_relations = [1,0.5538, 0.1673 , 0.1673*0.5185 ]
    meanvalues= [24.210, 24.002 , 27.276 , 27.238]
    amplitude = params[0]

    for i in range(4):
        sigma = params[i * 3 + 2]  # Get the fitted standard deviation for this Gaussian
        
        y += amplitude*amplitude_relations[i]* np.exp(-(x - meanvalues[i])**2 / (2 * sigma**2 + eps))
    return y


sigma=[]
area=[]

# Función para graficar el espectro de energía convolucionado
def plot_convolved_spectrum(Data, a, b,i, ax=None):

    maxim = np.max(Data)
    workingarray = np.squeeze(Data)

    # Definir puntos máximos
    peaks, _ = find_peaks(workingarray / maxim, height=0.1)
    peak_values = workingarray[peaks] / maxim
    peak_indices = peaks

    # Calcular valores de energía correspondientes a los picos
    energy_spectrum = a + b * peak_indices

    # Preparar datos para convolución
    data = workingarray[:885] / maxim
    data_y = data / data.max()
    data_x = a + b * np.linspace(0, 885, num=len(data_y))

    # Ajustar múltiples gaussianas al espectro de energía
    peak_indices2, _ = find_peaks(data_y, height=0.1)
    peak_amplitudes = [1,0.5538, 0.1673 , 0.1673*0.5185 ]
    peak_means = [24.210, 24.002 , 27.276 , 27.238]
    peak_sigmas = [0.1] * 4

    params_init = list(zip(peak_amplitudes, peak_means, peak_sigmas))
    params_init = np.concatenate(params_init)

    # Ajustar curva
    params, params_cov = curve_fit(multi_gaussian, data_x, data_y, p0=params_init)

    # Obtener una interpolación de alta resolución del ajuste
    x_fine = np.linspace(data_x.min(), data_x.max(), num=20000)

    y_fit = multi_gaussian(x_fine, *params)

    #  Data Graphic
    ax.scatter(data_x, data_y, color='black', marker='o', s=20, label = 'Data' )
    y_data_error =np.sqrt(workingarray[:885])
    ax.plot(data_x, data_y + y_data_error/maxim, color='black',linestyle='--')
    ax.plot(data_x, data_y - y_data_error/maxim, color='black',linestyle='--')
    # Fit Graphic

    ax.plot(x_fine, y_fit, label="Fit", linewidth=1.5)



    # Extraer desviaciones estándar y amplitudes
    sigmas_array = params[2::3]
    # Calcular sigma promedio
    sigma.append(np.mean(sigmas_array))

    # Configuración del gráfico
    ax.set_xlabel('Energy (KeV)')
    ax.set_ylabel('Normalized Data')
    ax.legend()
    ax.set_title('Convolved Energy Spectrum')

    # Imprimir información
    print("Standard deviations:", sigmas_array)






fig, ax = plt.subplots()
plot_convolved_spectrum(Data, a, b,1,ax=ax)
ax.set_xlim(22, 28)
plt.show()