How to apply Kramers Kronig transformation on a NiO reflectivity or absorption spectrum?

Question

I have an experimental NiO absorption spectrum recorded in total electron yield across the M edge of Ni. I would like to extract the imaginary (and real) part of the refraction index (which also depends on the energy). Additionally, I have a reflectivity spectrum of NiO and I would also like to know how to extract the real and imaginary part of the refraction index.
Here is what I tried so far in Python:

1. I tried fitting the peaks in the reflectivity (or absorption) spectrum with a function I derived in the following manner: -I modelled the real and the imaginary part of the dielectric function with the Lorentz oscillator function -Using the real and imaginary part of the dielectric function, I derived the function of the real and imaginary part of the refractive index -Using the imaginary part of the refractive index, I derived the absorption coefficient form -Then, I used a superposition of these coefficients with different scaling factors to account for the peaks in the absorption spectra. The fitting parameters are the scaling factors and the parameters of the lorentzian model. Here is a bit of insight into my functions:

energy_scale=np.arange(0,100, 0.001)
def real_comp2(w, w_p, w_0, nu):
    vp = 8.854187817e-12
    return vp+ (vp*w_p**2*(-w**2+w_0**2))/((w**2-w_0**2)**2+w**2*nu**2)

def imag_comp2(w, w_p, w_0, nu):
    vp = 8.854187817e-12
    return -(vp*w_p**2*w*nu/((w**2-w_0**2)**2+w**2*nu**2))

def real_refr(w, w_p, w_0, nu):
    # vacuum_permittivity = 8.854187817e-12  # Farads per meter
    # srvp = np.sqrt(vacuum_permittivity)
    
    return 0.707* np.sqrt(real_comp2(w, w_p, w_0, nu)+np.sqrt((real_comp2(w, w_p, w_0, nu))**2+(imag_comp2(w, w_p, w_0, nu))**2))

def imag_refr(w, w_p, w_0, nu):
    # vacuum_permittivity = 8.854187817e-12  # Farads per meter
    # srvp = np.sqrt(vacuum_permittivity)
    return 0.707* np.sqrt(-real_comp2(w, w_p, w_0, nu)+np.sqrt((real_comp2(w, w_p, w_0, nu))**2+(imag_comp2(w, w_p, w_0, nu))**2))


def reflectivity2(w, w_p, w_0, nu):
    return ((real_refr(w, w_p, w_0, nu)-1)**2+imag_refr(w, w_p, w_0, nu)**2)/((real_refr(w, w_p, w_0, nu)+1)**2+imag_refr(w, w_p, w_0, nu)**2)

def reflectivity2_norm(w, w_p, w_0, nu, A,c):
    return c+-A*(reflectivity2(w, w_p, w_0, nu) - reflectivity2(w, w_p, w_0, nu).min())/(reflectivity2(w, w_p, w_0, nu).max() - reflectivity2(w, w_p, w_0, nu).min())


def absorption_coeff(w, w_p, w_0, nu):
    return 4*np.pi*0.08*1e7*w_0*imag_refr(w, w_p, w_0, nu)

def absorption_coeff_norm(w, w_p, w_0, nu, A):
    return A*(absorption_coeff(w, w_p, w_0, nu) - absorption_coeff(w, w_p, w_0, nu).min())/(absorption_coeff(w, w_p, w_0, nu).max() - absorption_coeff(w, w_p, w_0, nu).min())

def abs_superposition(w, w_p1, w_01, nu1, A1, w_p2, w_02, nu2, A2, w_p3, w_03, nu3, A3, c): #, w_p4, w_04, nu4, A4, c):
    return c+ absorption_coeff_norm(w, w_p1, w_01, nu1, A1)+ absorption_coeff_norm(w, w_p2, w_02, nu2, A2) + absorption_coeff_norm(w, w_p3, w_03, nu3, A3)# + absorption_coeff_norm(w, w_p4, w_04, nu4, A4)

def reflectivity_superposition(w, w_p1, w_01, nu1,  w_p2, w_02, nu2,  w_p3, w_03, nu3):
    return reflectivity2(w, w_p1, w_01, nu1) + reflectivity2(w, w_p2, w_02, nu2) + reflectivity2(w, w_p3, w_03, nu3)

Even if there are dimensional problems in my code, I think my approach is logically wrong. I don't think I can just eyeball which peaks (and how many) to fit in the absorption spectrum (or reflectivity spectrum) and make a random superposition for the model function. I should probably also subtract the background and take the geometry of the experiment into consideration as well (grazing incidence).

Is there a better way to do this? Maybe some source code or a program that takes everything into account? I have downloaded Larch, but I don't know how to do this with that program. Thank you very much for your help! If you need additional information, I am happy to provide it in the comments.