Added uncertainities

ercid_uncertshow_twice.py

+#!/usr/bin/env python3
+from rawprasslib import load_raw
+from rawautoparams import load_params
+from matplotlib import pyplot as plt
+import errno
+import numpy as np
+import os
+import time
+import scipy.optimize as spopt
+
+# SCRIPT DOES NOT CHECK FOR SANITY OF THE OUTCOME !!!
+# THAT'S YOUR JOB !!!
+
+# AE == appearance energy from LCQ obtained by sigmoid fitting
+# E_tcid == threshold collision induced dissociation energy from Armentrout
+# red dots .. residual intensity (TIC -fragment - parent)
+# blue dots .. relative intensity of the fragment
+# orange vertial bars .. standard deviation (std)
+# black vertilcal bars .. std of the mean (usually too small to see)
+# black lines - sigmoid fit, tangent in the inflex point.
+
+# Ab-initio conversion ratios derivation
+ev = 1.602176634e-19
+mol = 6.02214076e23
+farad = ev*mol
+ev2kj = farad/1000
+kj2ev = 1/ev2kj
+
+# Energy-resolved CID reference dissociation energies are taken from:
+#
+# BH_HH, BH_HOMe, BH_MeMe:
+# ========================
+# https://pubs.acs.org/doi/pdf/10.1021/acs.analchem.9b02257
+# Benzhydrylpyridinium Ions: A New Class of Thermometer Ions for the
+# Characterization of Electrospray-Ionization Mass Spectrometers
+#
+# R. Rahrt, T. Auth, M. Demireva, P. B. Armentrout, K. Koszinowski, Anal. Chem.
+# 2019, 91, 11703–11711.
+# Took from the Abstract - the values are the same as the one in the Table2
+#
+# p-H, p-Me, p-OMe, p-NO2, p-NO2_alt:
+# ===================================
+# https://link.springer.com/content/pdf/10.1007/s13361-017-1693-0.pdf
+# How Hot are Your Ions Really? A Threshold Collision-Induced Dissociation
+# Study of Substituted Benzylpyridinium “Thermometer” Ions
+#
+# J. E. Carpenter, C. P. McNary, A. Furin, A. F. Sweeney, P. B. Armentrout, J.
+# Am. Soc. Mass Spectrom. 2017, 28, 1876–1888.
+# Took from the Table2, p-NO2_alt means the dissociation of the NO2
+#
+#
+# Due to the argumentation used by Armenrout in the abovementioned paper (2017)
+# we ommited the values presented in:
+# https://www.sciencedirect.com/science/article/pii/S1387380616303608
+# Experimental bond dissociation energies of benzylpyridinium thermometer ions
+# determined by threshold-CID and RRKM modeling
+#
+# D. Gatineau, A. Memboeuf, A.  Milet, R. B. Cole, H. Dossmann, Y. Gimbert, D.
+# Lesage, International Journal of Mass Spectrometry 2017, 417, 69–75.
+
+
+# energies, uncertainty (eV) - from Armentrout experimental values
+energs = {"p-H": (2.58, 0.15),
+          "p-Me": (2.26, 0.13),
+          "p-OMe": (1.93, 0.08),
+          "p-NO2": (3.04, 0.12),
+          # "p-NO2_alt": (2.68, 0.13), - Removed, commented below.
+          "BH_HH": (1.79, 0.11),
+          "BH_MeMe": (1.55, 0.13),
+          "BH_HOMe": (1.37, 0.14)}
+
+# This value is one from the Armenthrout paper - After discussion we decided to
+# remove it as a precaution as the dissociated bond differs from the others
+# (daughter is BnPy^+, not RBn^+). This can resolve into different kinetics and
+# thus maybe a slightly deviated BDE_{tcid}
+
+
+# masses to select - parent, fragment
+ions = {"p-H": (170, 91),
+        "p-Me": (184, 105),
+        "p-OMe": (200, 121),
+        "p-NO2": (215, 136),
+        # "p-NO2_alt": (215, 169),
+        "BH_HH": (246, 167),
+        "BH_MeMe": (274, 195),
+        "BH_HOMe": (276, 197)}
+
+# a bit generous peakwidth
+peakwidth = 1.5
+
+
+def boltzmann(x, bottom, top, xhalf, slope):
+    """boltzmann sigmoid curve function"""
+    y = bottom + (top - bottom) / (1 + np.exp((x-xhalf)/slope))
+    return y
+
+
+def fitline(x, a):
+    return x*a
+
+
+def get_tangentparams(parameters):
+    """generate a, b for y = ax + b from sigomid fit parameters"""
+    bottom, top, xhalf, slope = parameters[0]
+    a = -(top - bottom) / (4*slope)
+    yhalf = top + (bottom - top) / 2
+    b = yhalf - a * xhalf
+    return a, b
+
+
+def get_r_sqrd(function, xdata, ydata, popt):
+    """routine for getting fit probability"""
+    residuals = ydata - function(xdata, *popt)
+    ss_res = np.sum(residuals**2)
+    ss_tot = np.sum((ydata-np.mean(ydata))**2)
+    r_sqrd = 1 - (ss_res / ss_tot)
+    return r_sqrd
+
+
+figure = plt.figure(figsize=(8, 8), dpi=300, constrained_layout=True)
+grid = figure.add_gridspec(int(np.ceil(len(energs)/2))+1, 2)
+
+
+xAE, yAE, dates, yerrs = [], [], [], []
+for n, i in enumerate(energs.keys()):
+    print("processing %s..." % i)
+
+    # data readout and processing
+    filename = i+".RAW" if not i == "p-NO2_alt" else "p-NO2.RAW"
+    if not os.path.isfile(filename):
+        raise FileNotFoundError(
+            errno.ENOENT, os.strerror(errno.ENOENT), filename)
+    dates.append(os.stat(filename).st_mtime)
+    matrix = load_raw(filename)[0]
+    legends = load_params(filename)
+    normenergs = np.asarray(legends[1]).T[4]
+    times = matrix[0][0]
+    assert len(times) == len(normenergs)
+    masses = matrix[1]
+    parentargs = np.where((masses > (ions[i][0]-peakwidth/2)) &
+                          (masses < (ions[i][0]+peakwidth/2)))
+    daughterargs = np.where((masses > (ions[i][1]-peakwidth/2)) &
+                            (masses < (ions[i][1]+peakwidth/2)))
+    # tic == zero -> bad scan
+    tic = np.sum(matrix[2].T, axis=0)
+    assert len(times) == len(tic)
+    good = np.where(tic > 0)
+    parentints = np.divide(np.sum(matrix[2].T[parentargs], axis=0), tic)[good]
+    daughterints = np.divide(np.sum(
+        matrix[2].T[daughterargs], axis=0), tic)[good]
+    rest = 1 - parentints - daughterints
+
+    # fitting
+    x = normenergs[good]
+    y = daughterints
+    xaver = np.unique(x)
+    xarrs = [np.where(x == xval) for xval in xaver]
+    yaver = [np.average(y[arr]) for arr in xarrs]
+    stds = [np.std(y[arr], ddof=1) for arr in xarrs]
+    stdmeans = [stds[m]/np.sqrt(len(xarrs)) for m, arr in enumerate(xarrs)]
+    try:
+        guess = [min(y), max(y), x[np.argmax(y > (max(y)/2))], 0.1]
+        parameters = spopt.curve_fit(boltzmann, x, y, guess)
+        print(parameters[0])
+    except RuntimeError:
+        print("optimization failed, giving up")
+        raise RuntimeError
+
+    a, b = get_tangentparams(parameters)
+    zerocross = -b/a
+    yline = [0, max(y)]
+    xline = [zerocross, (yline[1]-b)/a]
+    rsqrd = get_r_sqrd(boltzmann, x, y, parameters[0])
+    xfit = np.linspace(x[0], x[-1], 500)
+    yfit = boltzmann(xfit, *parameters[0])
+
+    plot = figure.add_subplot(grid[int(n/2), n % 2],
+                              ylabel="$I_{frag}\ /\ \Sigma I$",
+                              xlim=(min(x), max(x)), ylim=(-0.1, 1.1))
+    plot.spines['top'].set_visible(False)
+    plot.spines['right'].set_visible(False)
+    plot.errorbar(xaver, yaver, yerr=stdmeans, zorder=10,
+                  color=[0, 0, 0, 0.5], linestyle='None')
+    plot.errorbar(xaver, yaver, yerr=stds, zorder=9,
+                  color=[1, 0.5, 0, 1], linestyle='None')
+    plot.plot(x, rest, '.', color=[1, 0.8, 0.8, 0.1])
+    plot.plot(x, daughterints, '.', color=[0.4, 0.7, 1, 0.1])
+    plot.plot(xfit, yfit, color=[0, 0, 0, 0.4], zorder=11)
+    plot.plot(xline, yline, color=[0, 0, 0, 0.7], zorder=12)
+    plot.plot(xaver, yaver, '.', color=[0, 0, 0, 1], zorder=13)
+    plot.text(0.01, 1,
+              "{}\n$BDE_{{AE}}={:.2f}$\n$BDE_{{tcid}}={}\ eV$\n$r^2={:.4f}$"
+              .format(i, zerocross, energs[i][0], rsqrd),
+              transform=plot.axes.transAxes, va='top')
+
+    xAE.append(zerocross)
+    yAE.append(energs[i][0])
+    yerrs.append(energs[i][1])
+    print(zerocross)
+    print("DONE! (%s)" % i)
+
+finalfit = spopt.curve_fit(fitline, xAE, yAE)
+k = finalfit[0][0]
+rsq = get_r_sqrd(fitline, np.asarray(xAE), np.asarray(yAE), finalfit[0])
+
+figure.suptitle("Calibration plots\n{} - {}".format(
+    time.ctime(min(dates)), time.ctime(max(dates))))
+
+print("Plotting the fitted values")
+plot = figure.add_subplot(grid[-1, :], xlabel="BDE_$AE$",
+                          ylabel="$BDE_{{tcid}}$")
+plot.errorbar(xAE, yAE, yerr=yerrs, color=[0, 0, 0, 0.4], linestyle='None')
+plot.plot(xAE, yAE, '.', color='black')
+plot.plot((min(xAE), max(xAE)), (min(xAE)*k, max(xAE)*k), color='red')
+for n, i in enumerate(energs.keys()):
+    plot.annotate(i, (xAE[n]+0.1, yAE[n]), rotation=0, va='top')
+plot.text(0.01, 1, "$k={:.4f}\ eV\ ({:.4f}\ kJ/mol),\ R^2 ={:.2f}$".format(
+    k, k*farad/1000, rsq), transform=plot.axes.transAxes, va='top')
+plt.savefig("uncertdoublecol_twice.png")