script to calculate new calibration constants from monochromator output taken with a Hg lamp

hg_calibration.py

+#!/usr/bin/python
+
+"""
+A small analysis program for data taken with the polychromator 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.
+
+The result is new calibration curve, including uncertainties
+"""
+
+import sys
+import matplotlib as mpl
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.cbook as cbook
+from scipy import optimize as opt
+from scipy.optimize import curve_fit
+
+file=str(sys.argv[1])
+print (file)      # input file must contain (only) lines of type: \twavelength\tPindiodeCurrent
+
+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   = 0      # 0: calibrate lambda_theor vs. lambda_meas   1: calibrate lambda_meas vs. lambda_theor
+
+data = np.loadtxt(file, delimiter='\t', converters = {0: lambda s: float(s.strip() or 0)})
+
+x= data[:,1]
+y= data[:,2]
+
+xxorg= [365.015, 404.656, 435.833, 546.074 ]   # from http://opticalengineering.spiedigitallibrary.org/article.aspx?articleid=1556406
+xx   = list(map (lambda a: a-wloffset,xxorg) ) # [365.015-wloffset, 404.656-wloffset, 435.833-wloffset, 546.074-wloffset ]
+xxerr=[0.001, 0.001, 0.001, 0.001 ]
+yy   =[]
+yyorg=[]
+yyerr=[]
+ww   =[]
+
+indices1 = (np.nonzero((x>xxorg[0]+shift-window)&(x<xxorg[0]+shift+window)))
+x1 = x[indices1]
+y1 = y[indices1]
+indices2 = (np.nonzero((x>xxorg[1]+shift-window)&(x<xxorg[1]+shift+window)))
+x2 = x[indices2]
+y2 = y[indices2]
+indices3 = (np.nonzero((x>xxorg[2]+shift-window)&(x<xxorg[2]+shift+window)))
+x3 = x[indices3]                                                  
+y3 = y[indices3]                                                  
+indices4 = (np.nonzero((x>xxorg[3]+shift-window)&(x<xxorg[3]+shift+window)))
+x4 = x[indices4]
+y4 = y[indices4]
+
+
+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')
+
+popt,pcov = curve_fit(gaus,x1,y1,p0=[2.5E-9,xxorg[0]+shift,2.])
+print (pcov)
+f2.plot(x1,gaus(x1,*popt),'ro:',label="1: %.1f nm (vs. %.1f nm)" % (popt[1],xxorg[0]))
+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.))
+ww.append(1./yyerr[-1]**2)
+
+popt,pcov = curve_fit(gaus,x2,y2,p0=[4.7E-9,xxorg[1]+shift,2.])
+f2.plot(x2,gaus(x2,*popt),'ro:',label="2: %.1f nm (vs. %.1f nm)" % (popt[1],xxorg[1]))
+yy.append(popt[1]-wloffset)
+yyorg.append(popt[1])
+perr = np.sqrt(np.diag(pcov))
+###print (perr)
+#yyerr.append(perr[1])
+yyerr.append(popt[2]/np.sqrt(4.))
+ww.append(1./yyerr[-1]**2)
+
+popt,pcov = curve_fit(gaus,x3,y3,p0=[5.2E-9,xxorg[2]+shift,2.])
+f2.plot(x3,gaus(x3,*popt),'ro:',label="3: %.1f nm (vs. %.1f nm)" % (popt[1],xxorg[2]))
+yy.append(popt[1]-wloffset)
+yyorg.append(popt[1])
+perr = np.sqrt(np.diag(pcov))
+#print (perr)
+#yyerr.append(perr[1])
+yyerr.append(popt[2]/np.sqrt(4.))
+ww.append(1./yyerr[-1]**2)
+
+popt,pcov = curve_fit(gaus,x4,y4,p0=[5.2E-9,xxorg[3]+shift,2.])
+f2.plot(x4,gaus(x4,*popt),'ro:',label="4: %.1f nm (vs. %.1f nm)" % (popt[1],xxorg[3]))
+yy.append(popt[1]-wloffset)
+yyorg.append(popt[1])
+perr = np.sqrt(np.diag(pcov))
+#print (perr)
+#yyerr.append(perr[1])
+yyerr.append(popt[2]/np.sqrt(4.))
+ww.append(1./yyerr[-1]**2)
+
+
+f2.legend()
+f2.set_title("Hg Lamp")    
+f2.set_xlabel('wavelength (nm)')
+f2.set_ylabel('Pin Diode current')
+
+f3 = plt.subplot(212)
+
+if (invert):
+    m,b = np.polyfit(yy , xx, 1) 
+else:
+    m,b = np.polyfit(xx , yy, 1) 
+
+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),
+    it is an exact copy of the function gsl_fit_wlinear.
+    """
+    # 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)
+    dx = xa-wm_x
+    dy = ya-wm_y
+    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
+    cov_00 = (1.0/W) * (1.0 + wm_x**2/wm_dx2)
+    cov_11 = 1.0 / (W*wm_dx2)
+    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
+
+
+aa,bb,cov_00,cov_11,cov_01 = wlinear_fit(np.array(xx),np.array(yy),np.array(ww))
+#print (xx, yy, ww, yyerr)
+
+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 
+
+if (invert):
+    f3.errorbar(np.array(yyorg),np.array(xxorg),np.array(xxerr),np.array(yyerr),'yo')
+else:
+    f3.errorbar(np.array(xxorg),np.array(yyorg),np.array(xxerr),np.array(yyerr),'yo')
+
+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.))  
+    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.))  
+    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.))  
+
+f3.legend(loc="lower right")
+#f3.set_title("Calibration")    
+if (invert):
+    f3.set_ylabel('wavelength expected (nm)')
+    f3.set_xlabel('wavelength measured (nm)')
+else:
+    f3.set_xlabel('wavelength expected (nm)')
+    f3.set_ylabel('wavelength measured (nm)')
+
+plt.show()
+
+
+