Commit 7749b36f authored by Yan's avatar Yan
Browse files

Added uncertainities - single graph

parent 6debba87
#!/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
factor = 10.502
filename = "BH_MeMe"
outfilename = "777_uncertshow"
parent = (272.8, 274.6)
daughter = (194.4, 195.8)
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=(5, 1.5), dpi=300, constrained_layout=True)
print("processing %s..." % filename)
# data readout and processing
filename = filename+".RAW"
if not os.path.isfile(filename):
raise FileNotFoundError(
errno.ENOENT, os.strerror(errno.ENOENT), filename)
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 > parent[0]) & (masses < parent[1]))
daughterargs = np.where((masses > daughter[0]) & (masses < daughter[1]))
# 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(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={:.2f}\ eV$\n$BDE={:.2f}\ kJ/mol$\n$r^2={:.4f}$"
.format(filename, zerocross, (zerocross*factor/farad)*1000,
zerocross*factor, rsqrd),
transform=plot.axes.transAxes, va='top')
print(zerocross)
print("DONE! (%s)" % filename)
plt.savefig(outfilename+".jpg")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment