Add method to calibrate dataset/ refine spectral line poisition

[PauTorru]

A common issue when running espm is that the peak poisitions of the G matrix do not perfectly align with where the peaks are in the dataset, as well as their shape.

A good solution would be to perform a "precalibration" step, where the average spectrum is fit by the G matrix but alowing shifting the positions of the gaussians.

The idea would be to fit the average data with a set of free gaussians with initial positions and intensity ratios set by the kV_xrays.json tables and then fit with least squares, allowing small energy shifts and small changes in width. The result of this fit should then be used to set the energy scale, energy offset, gaussian width and intercept (as is called in espm) and small shifts in the line positions).

Ideally, this can set a proper G matrix that will improve the results of the subsequent NMF decompositon.

The addition of this functionality will be part of Perry's summer intership with us.

[adriente]

I think this is a crucial feature that needs to be implemented as soon as possible. I was planning to do this but I couldn't find the time.

I have a few suggestions :

• I don't recommend setting the ratios of the gaussian but to use fully free gaussians (except that they should be in a certain energy range). The current energy calibration method is not good but is not too far away either so it should be good starting point for the fit.
• Originally I was thinking of implementing the fit such that it produces a table similar (same format) to the kV_xrays.json table but with modified energies. It would then be stored internally in the .hspy file. The espm modelling would then use the generated tables for performing decomposition, etc ...
  It would then be easier for the user to come back to the result of the fit and even tune the table if needed. It would also prepare a next step : we could provide the user with functions to visualize or browse espm tables and even apply modifications (such as applying experimentally determined k factors).

[adriente]

By the way :

It could be interesting not to use gaussian to fit but 2nd order polynomials. It only works well over a limited number of points (less than 10) around the maximum of the peak. It is basically similar to performing a taylor development, that's why it is only valid very locally. The advantage are :

• It doesn't care about the actual correct theoretical model
• It is very fast to calculate
• It is very stable, if you have peaks around, a background, etc ... It is not an issue.

It is a well known tricks in Xray diffraction.

[perrymok]

in fact i have tried fitting a degree N-1 polynomial using N points, where N is the number after filtering out bad lines. the problem is that low energy lines can act quite funny.

[adriente]

And did you try 2nd order on N points ?

[perrymok]

yes actually i plotted it out vs different degrees -- all degrees of polynomial agree that the energy misalignment is linear (they start diverging after 15 keV since there are no sample lines after 15)

fitting peaks individually independent of each other (dashed lines are theoretical lines, and dotted lines are calibrated lines)

the problem with this approach is that, for example grey lines are Re, which has a very low concentration, thus its lines basically disappeared into other bell curves or the background, making its naive fitted lines unrealistically far away from where it should be.

calibrated peaks using polynomial of degree 3

using polynomial, the problem for low concentration or low cs lines is stablised, but the polynomial calibration is actually even worse then the theoretical table for large peaks.

i will use some other strategies. more coming next monday.

[adriente]

Could you, if possible, detail a bit more the last two graphs ?

What model did you choose to model the peak position (with equation) ? On how many points ? Did you put boundaries to the fit ? Could you show a list of results (peak name, peak energy (tabulated), peak energy (fitted), peak intensity, etc ...) ?

Shall we organise a remote meeting ? (if necessary)

[perrymok]

I am not sure what you mean by model, but i read the 200keV table and the dashed (- - - -) lines are the lines present in the table. The dotted (.....) lines are calibrated lines.

For the first plot, it was done by the following algorithm:

def gauss(A, x, x0, sigma, C, m):
    return A * np.exp(-((x - x0) ** 2) / (2 * sigma**2)) + C + m * (x - x0)

$A\cdot\mathcal N(x_0, \sigma) + (m (x-x_0)+C)$ where $m(x-x_0)+C$ is the background noise.

def fit_single_peak(E_theoretical, sigma_expected, window):
    mask = (energy_axis >= E_theoretical - window) & (
        energy_axis <= E_theoretical + window
    )
    xdata = energy_axis[mask]
    ydata = average_spectrum[mask]

    idx_max = np.argmax(ydata)
    x0_guess = xdata[idx_max]
    A_guess = np.max(ydata) - np.min(ydata)
    sigma_guess = sigma_expected
    C_guess = np.min(ydata)

    p0 = [A_guess, x0_guess, sigma_guess, C_guess, 0.0]
    bounds = (
        [0.0, E_theoretical - window, 0.9 * sigma_expected, 0.0, -1.0],
        [np.inf, E_theoretical + window, 1.1 * sigma_expected, np.inf, 1.0],
    )
    popt, _ = curve_fit(
        gauss,
        xdata,
        ydata,
        p0,
        bounds=bounds,
    )
    return popt

The above function is ran either for all lines in the table, or a selected subset of lines, with sigma_expected = model.width_slope * E_theoretical + model.width_intercept / 2.3548 and window = 3 * sigma_expected. It basically tries to fit gauss in the ball of radius window around E_theoretical, where E_theoretical is the line energy given by the table.

For the second plot, we collect the E_theoretical and the returned calibrated energies, poly fit a degree $d$ polynomial $p$, and set calibrated energy = p(E_theoretical).

Let's look at the polynomial $p$ fitted using different degrees:

I imagine we want the resulting polynomial to look linear, so any $p$ of degree $>8$ is not good. Each dot in this plot is the theoretical energy vs calibrated energy. Notice that there are a lot of pink and grey dots below the black line:

They are Pt and Re, which are of very low concentration in the sample, meaning they get fitted into another peak on the left, which in turns drag the fitted $p$ away from what it should be.

To combat this, I have added a weight for each line in the poly fit by the peak magnitude, the result is way better compared to plot 2 from the previous comment:

I cannot get over the fact that it can happen that the lines of an element with low concentration aligns closely to the lines of another abundant element, making it "very important" during the fitting. Also, I just don't see how this is better than the existing affine calibration using 2 lines.

In fact, if we know which element has the highest concentration, this problem would have been much easier.

[adriente]

Ok, I see one main issue here. I think that instead of fitting the peaks with gaussians we should fit each peak with a 2nd order polynomial on a severly reduced amount of points for each peak. I believe, that second order polynomial results can be sorted into correct ones :

• the maximum of the polynomial sits in the selected range
  and discard incorrect ones :
• The maximum of the polynomial sits outside (or on the edge) of the selected range.
  You don't need all peaks to participate to the calibration.
  You can also constraint the polynomial to negative 2nd degree coeff, so that it has a bell shape.

For the energy scale, so far we should keep the 1st order polynomial. For quite a few datasets there are peaks beyond 15keV.

[PauTorru]

You can also constraint the polynomial to negative 2nd degree coeff, so that it has a bell shape.

Is there a physical reason for this? I honestly also expect a 2nd degree polynomical to be largely enough, but why just negative coefficient?

[PauTorru]

In fact, if we know which element has the highest concentration, this problem would have been much easier.

Why can know by runing the quan spectrum function on the average spectrum.
Aside from that I'd like to know how (est.G@est.W@est.H).mean(1) compares to s.data.mean((0,1)) once you run decompositions on calibrated/uncalibrated data to see how well the calibration works.
@perrymok do you think you can produce this plots in on a couple of the datasets you have?

[perrymok]

refer to the calibration branch for the polyfit calibration and an example notebook in the notebook folder

as a side note, i think we can include the plot_table function in the new calibration.ipynb in the library as a utility function.