continued with covariance calculation

tools/test/cov.py

@@ -12,6 +12,17 @@ import numpy.linalg as la
 use_scipy = False
 
 
+# column indices in mc file
+ki_start_idx = 0
+kf_start_idx = 3
+wi_idx = 9
+wf_idx = 10
+
+
+# constants
+sig2hwhm = np.sqrt(2. * np.log(2.))
+sig2fwhm = 2.*sig2hwhm
+
 if use_scipy:
 	import scipy as sp
 	import scipy.constants as co
@@ -29,10 +40,10 @@ k2_to_E = 1./E_to_k2
 #
 def load_events(filename):
 	dat = np.loadtxt(filename)
-	ki = dat[:, 0:3]
-	kf = dat[:, 3:6]
-	wi = dat[:, 9]
-	wf = dat[:, 10]
+	ki = dat[:, ki_start_idx:ki_start_idx+3]
+	kf = dat[:, kf_start_idx:kf_start_idx+3]
+	wi = dat[:, wi_idx]
+	wf = dat[:, wf_idx]
 
 	w = wi * wf
 	Q = ki - kf
@@ -82,11 +93,62 @@ def calc_covar(Q, E, w):
 	print("Resolution matrix in (Qpara, Qperp, Qup, E) system:\n%s\n" % Qres_Q)
 
 
+	#[ evals, evecs ] = la.eig(Qcov_Q)
+	#print("Ellipsoid fwhm radii:\n%s\n" % (np.sqrt(evals) * sig2fwhm))
+
+	return Qres_Q
+
+
+
+#
+# calculate the characteristics of a given ellipse
+#
+def descr_ellipse(quadric):
+	[ evals, evecs ] = la.eig(quadric)
+	fwhms = 1./np.sqrt(evals) * sig2fwhm
+
+	angles = np.array([])
+	if len(quadric) == 2:
+		angles = np.array([ np.arctan2(evecs[1][0], evecs[0][0]) ])
+
+	return [fwhms, angles/np.pi*180.]
+
 
+
+#
+# describe the ellipsoid by a principal axis trafo
+#
+def calc_ellipses(Qres_Q):
+	[fwhms, angles] = descr_ellipse(Qres_Q)
+	print("4d resolution ellipsoid fwhm lengths:\n%s\n" % fwhms)
+
+
+	# 2d sliced ellipses
+	Qres_QxE = np.delete(np.delete(Qres_Q, 2, axis=0), 2, axis=1)
+	Qres_QxE = np.delete(np.delete(Qres_QxE, 1, axis=0), 1, axis=1)
+	[fwhms_QxE, angles_QxE] = descr_ellipse(Qres_QxE)
+	print("Qx/E ellipse fwhm lengths and slope angle:\n%s, %f\n" % (fwhms_QxE, angles_QxE[0]))
+
+	Qres_QyE = np.delete(np.delete(Qres_Q, 2, axis=0), 2, axis=1)
+	Qres_QyE = np.delete(np.delete(Qres_QyE, 0, axis=0), 0, axis=1)
+	[fwhms_QyE, angles_QyE] = descr_ellipse(Qres_QyE)
+	print("Qy/E ellipse fwhm lengths and slope angle:\n%s, %f\n" % (fwhms_QyE, angles_QyE[0]))
+
+	Qres_QzE = np.delete(np.delete(Qres_Q, 1, axis=0), 1, axis=1)
+	Qres_QzE = np.delete(np.delete(Qres_QzE, 0, axis=0), 0, axis=1)
+	[fwhms_QzE, angles_QzE] = descr_ellipse(Qres_QzE)
+	print("Qz/E ellipse fwhm lengths and slope angle:\n%s, %f\n" % (fwhms_QzE, angles_QzE[0]))
+
+
+
+#
+# main
+#
 if __name__ == "__main__":
 	if len(sys.argv) <= 1:
 		print("Please specify a filename.")
 		exit(-1)
 
 	[Q, E, w] = load_events(sys.argv[1])
-	calc_covar(Q, E, w)
+	Qres = calc_covar(Q, E, w)
+	calc_ellipses(Qres)