Commit 623670aa authored by Markus Gaug's avatar Markus Gaug
Browse files

added changes suggested in merge request discussion

parent a8342c0c
Pipeline #6011 passed with stages
in 4 minutes and 18 seconds
#!/usr/bin/python
"""
A small analysis program for data taken with the monochromator from a Hg spectral lamp
The program read the lvm-file from input, fits 4 Gaussian peaks around the expected lines
and derives a linear regression between fitted and expected intensity peak positions.
A small analysis program for data taken with the monochromator
from a Hg spectral lamp
The result is new calibration curve, including uncertainties
The program reads the lvm-file from input, fits 4 Gaussian peaks
around the expected lines and derives a linear regression between
fitted and expected intensity peak positions.
The result is a new calibration curve, including uncertainties
"""
import sys
......@@ -13,67 +16,81 @@ import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
filename=str(sys.argv[1])
print (filename) # input file must contain (only) lines of type: \twavelength\tPindiodeCurrent
filename = str(sys.argv[1])
# input file must contain (only) lines of type: \twavelength\tPindiodeCurrent
print(filename)
# 0: start calibration from lambda=0,
# otherwise choose wloffset as reference point
wloffset = 0.
# -3.5 for 1HgMarkus, -6.5 for 1HgCalcPere
shift = 0
# fit ranges: (central wavelength + shift) +- window
window = 5.
# False: calibrate lambda_theor vs. lambda_meas
# True: calibrate lambda_meas vs. lambda_theor
invert = False
wloffset = 0. # 0: start calibration from lambda=0, otherwise choose wloffset as reference point
shift = -6.5 # -3.5 for 1HgMarkus, -6.5 for 1HgCalcPere
window = 5. # fit ranges: (central wavelength + shift) +- window
invert = False # 0: calibrate lambda_theor vs. lambda_meas 1: calibrate lambda_meas vs. lambda_theor
data = np.loadtxt(filename, delimiter='\t', converters = {0: lambda s: float(s.strip() or 0)})
data = np.loadtxt(filename, delimiter='\t',
converters={0: lambda s: float(s.strip() or 0),
2: lambda s: float(s)*1.0e9}) # convert to nA
x= data[:,1]
y= data[:,2]
x = data[:, 1]
y = data[:, 2]
xxorg= np.array([365.015, 404.656, 435.833, 546.074 ]) # from http://opticalengineering.spiedigitallibrary.org/article.aspx?articleid=1556406
# http://opticalengineering.spiedigitallibrary.org/article.aspx?articleid=1556406
xxorg = np.array([365.015, 404.656, 435.833, 546.074])
xx = xxorg
xx -= wloffset # list(map (lambda a: a-wloffset,xxorg) ) # [365.015-wloffset, 404.656-wloffset, 435.833-wloffset, 546.074-wloffset ]
xxerr=np.array([0.001, 0.001, 0.001, 0.001 ])
yy =[]
yyorg=[]
yyerr=[]
ww =[]
# list(map (lambda a: a-wloffset,xxorg) )
xx -= wloffset
xxerr = np.array([0.001, 0.001, 0.001, 0.001])
yy = []
yyorg = []
yyerr = []
ww = []
def gaus(x,a,x0,sigma):
def gaus(x, a, x0, sigma):
return a*np.exp(-(x-x0)**2/(2*sigma**2))
f2 = plt.subplot(211)
f2.plot(x,y,'b+:',label='data')
f2.plot(x, y, 'b+:', label='data', color='red')
amp_guesses = [2.5E-9, 4.7E-9, 5.2E-9, 5.2E-9]
amp_guesses = [2.5, 4.7, 5.2, 5.2]
colors = ['violet', 'blue', 'green', 'orange']
for i, wavelength in enumerate(xxorg):
indices = np.nonzero((x > wavelength + shift - window) & (x < wavelength + shift + window))
indices = np.nonzero((x > wavelength+shift-window) & (x < wavelength+shift+window))
x_window = x[indices]
y_window = y[indices]
popt, pcov = curve_fit(gaus, x_window, y_window, p0=[amp_guesses[i], wavelength + shift, 2])
popt, pcov = curve_fit(gaus, x_window, y_window, p0=[amp_guesses[i], wavelength+shift, 2])
if i is 0:
print(pcov)
plot_label = "%d: %.1f nm (vs. %.1f nm)" % (i + 1, popt[1], wavelength)
f2.plot(x_window, gaus(x_window, *popt), 'ro:', label=plot_label)
f2.plot(x_window, gaus(x_window, *popt), 'ro:', label=plot_label, color=colors[i])
yy.append(popt[1] - wloffset)
yyorg.append(popt[1])
perr = np.sqrt(np.diag(pcov))
# yyerr.append(perr[1])
yyerr.append(popt[2]/np.sqrt(4)) # why 4? not len(x_window)?
yyerr.append(popt[2]/np.sqrt(len(x_window)))
ww.append(1./yyerr[-1]**2)
f2.legend()
f2.set_title("Hg Lamp")
f2.legend(fontsize=8)
f2.set_title("Hg Lamp")
f2.set_xlabel('wavelength (nm)')
f2.set_ylabel('Pin Diode current')
f2.set_ylabel('Pin Diode current (nA)')
f3 = plt.subplot(212)
if (invert):
m,b = np.polyfit(yy , xx, 1)
m, b = np.polyfit(yy, xx, 1)
else:
m,b = np.polyfit(xx , yy, 1)
m, b = np.polyfit(xx, yy, 1)
def wlinear_fit (xa,ya,w) :
def wlinear_fit(xa, ya, w):
"""
Fit (x,y,w) to a linear function, using exact formulae for weighted linear
regression. This code was translated from the GNU Scientific Library (GSL),
......@@ -82,12 +99,12 @@ def wlinear_fit (xa,ya,w) :
# compute the weighted means and weighted deviations from the means
# wm denotes a "weighted mean", wm(f) = (sum_i w_i f_i) / (sum_i w_i)
W = np.sum(w)
wm_x = np.average(xa,weights=w)
wm_y = np.average(ya,weights=w)
wm_x = np.average(xa, weights=w)
wm_y = np.average(ya, weights=w)
dx = xa-wm_x
dy = ya-wm_y
wm_dx2 = np.average(dx**2,weights=w)
wm_dxdy = np.average(dx*dy,weights=w)
wm_dx2 = np.average(dx**2, weights=w)
wm_dxdy = np.average(dx*dy, weights=w)
# In terms of y = a + b x
b = wm_dxdy / wm_dx2
a = wm_y - wm_x*b
......@@ -96,39 +113,62 @@ def wlinear_fit (xa,ya,w) :
cov_01 = -wm_x / (W*wm_dx2)
# Compute chi^2 = \sum w_i (y_i - (a + b * x_i))^2
# chi_2 = np.sum (w * ((ya-(a+b*xa))**2))
return a,b,cov_00,cov_11,cov_01
return a, b, cov_00, cov_11, cov_01
aa,bb,cov_00,cov_11,cov_01 = wlinear_fit(np.array(xx),np.array(yy),np.array(ww))
#print (xx, yy, ww, yyerr)
aa, bb, cov_00, cov_11, cov_01 = wlinear_fit(np.array(xx),
np.array(yy),
np.array(ww))
chi2 = 0
for i in range(len(xx)):
print (yy[i]-(aa+bb*xx[i]))
print (yyerr[i])
chi2 += ((yy[i]-(aa+bb*xx[i]))/yyerr[i])**2
# print (xx, yy, ww, yyerr)
print(yy-(aa+bb*xx))
print(yyerr)
chi2 = np.sum(((yy - (aa+bb*xx))/yyerr)**2)
if (invert):
f3.errorbar(np.array(yyorg),xxorg,xxerr,np.array(yyerr),'yo')
f3.errorbar(np.array(yyorg), xxorg, xxerr, np.array(yyerr), 'yo')
else:
f3.errorbar(xxorg,np.array(yyorg),xxerr,np.array(yyerr),'yo')
f3.errorbar(xxorg, np.array(yyorg), xxerr, np.array(yyerr), 'yo')
fit_param_dict = {'offset': wloffset,
'm': m,
'b': b,
'aa': aa,
'bb': bb,
'bb_inv': 1./bb,
'bb_cov': np.sqrt(cov_11),
'aa_inv': -aa/bb,
'aa_cov': np.sqrt(cov_00)/bb,
'chi2': chi2/2.}
if (wloffset != 0.):
f3.plot(x,m*x +b, label="(y-%.0f) = %.3f*(x-%.0f) + %.3f" % (wloffset,m,wloffset,b))
if (invert):
f3.plot(x,1/bb*x-aa/bb,label="(y-%.0f) = (%.4f$\pm$%.4f)*(x-%.0f) + (%.3f$\pm$%.3f), $\chi^{2}/NDF$=%.1f" % (wloffset,1./bb,np.sqrt(cov_11),wloffset,-aa/bb,np.sqrt(cov_00)/bb,chi2/2.))
# linear unweighted with offset
fit_label = "(y-{offset:.0f}) = {m:.3f}*(x-{offset:.0f})+{b:.3f}".format(**fit_param_dict)
f3.plot(x, m*x+b, label=fit_label)
if invert is True:
# linear weighted with offset, lambda_measured vs. lambda_theoretical
fit_label = r"(y-{offset:.0f}) = ({bb_inv:.4f}$\pm${bb_cov:.4f})*(x-{offset:.0f}) + ({aa_inv:.3f}$\pm${aa_cov:.3f}), $\chi^{{2}}/NDF$={chi2:.1f}".format(**fit_param_dict)
f3.plot(x, 1/bb*x - aa/bb, label=fit_label)
else:
f3.plot(x,bb*x+aa,label="(y-%.0f) = (%.4f$\pm$%.4f)*(x-%.0f) + (%.3f$\pm$%.3f), $\chi^{2}/NDF$=%.1f" % (wloffset,bb,np.sqrt(cov_11),wloffset,aa,np.sqrt(cov_00)/bb,chi2/2.))
else:
f3.plot(x,m*x+b,label="y = %.3f*x + %.3f" % (m,b))
if (invert):
f3.plot(x,1/bb*x-aa/bb,label="y = (%.4f$\pm$%.4f)*x + (%.3f$\pm$%.3f), $\chi^{2}/NDF$=%.1f" % (1./bb,np.sqrt(cov_11),-aa/bb,np.sqrt(cov_00)/bb,chi2/2.))
# linear weighted with offset, lambda_theoretical vs. lambda_measured
fit_label = r"(y-{offset:.0f}) = ({bb:.4f}$\pm${bb_cov:.4f})*(x-{offset:.0f}) + ({aa:.3f}$\pm${aa_cov:.3f}), $\chi^{{2}}/NDF$={chi2:.1f}".format(**fit_param_dict)
f3.plot(x, bb*x + aa, label=fit_label)
else:
# linear unweighted without offset
fit_label = "y = {m:.3f}*x+{b:.3f}".format(**fit_param_dict)
f3.plot(x, m*x+b, label=fit_label)
if invert is True:
# linear weighted without offset, lambda_measured vs. lambda_theoretical
fit_label = r"y = ({bb_inv:.4f}$\pm${bb_cov:.4f})*x + ({aa_inv:.3f}$\pm${aa_cov:.3f}), $\chi^{{2}}/NDF$={chi2:.1f}".format(**fit_param_dict)
f3.plot(x, 1/bb*x-aa/bb, label=fit_label)
else:
f3.plot(x,bb*x+aa,label="y = (%.4f$\pm$%.4f)*x + (%.3f$\pm$%.3f), $\chi^{2}/NDF$=%.1f" % (bb,np.sqrt(cov_11),aa,np.sqrt(cov_00)/bb,chi2/2.))
# linear weighted without offset, lambda_theoretical vs. lambda_measured
fit_label = r"y = ({bb:.4f}$\pm${bb_cov:.4f})*x + ({aa:.3f}$\pm${aa_cov:.3f}), $\chi^{{2}}/NDF$={chi2:.1f}".format(**fit_param_dict)
f3.plot(x, bb*x+aa, label=fit_label)
f3.legend(loc="lower right")
#f3.set_title("Calibration")
if (invert):
f3.legend(loc="lower right", fontsize=9)
# f3.set_title("Calibration")
if invert is True:
f3.set_ylabel('wavelength expected (nm)')
f3.set_xlabel('wavelength measured (nm)')
else:
......@@ -137,5 +177,7 @@ else:
plt.show()
#if invert is True:
# return 1./bb, -aa/bb, np.sqrt(cov_11), np.sqrt(cov_00)/bb
#else:
# return bb, aa, np.sqrt(cov_11), np.sqrt(cov_00)/bb
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