updated fesn script

scripts/FeSn_film_HB1A.py

@@ -1,38 +1,145 @@
 import matplotlib.pyplot as plt
 
+from tavi.data.fit import Fit1D
 from tavi.data.scan import Scan
 from tavi.instrument.resolution.cooper_nathans import CN
 from tavi.plotter import Plot1D
 from tavi.sample.xtal import Xtal
 
-instrument_config_json_path = "./scripts/IPTS32816_HB1A_exp1034/hb1a.json"
+instrument_config_json_path = "./test_data/IPTS32816_HB1A_exp1034/hb1a.json"
 tas = CN(SPICE_CONVENTION=True)
 tas.load_instrument_params_from_json(instrument_config_json_path)
 
-sample_json_path = "./test_data/IPTS32816_HB1A_exp1034/fesn.json"
-sample = Xtal.from_json(sample_json_path)
-tas.mount_sample(sample)
 
 ei = 14.4503
 ef = 14.4503
-hkl = (0, 0, 0)
 R0 = False
 
+sample_json_path_1 = "./test_data/IPTS32816_HB1A_exp1034/fesn1.json"
+sample_1 = Xtal.from_json(sample_json_path_1)
+tas.mount_sample(sample_1)
+
+rez1_l = tas.cooper_nathans(hkl_list=(0, 0, 0.5), ei=ei, ef=ef, R0=R0)
+rez1_l.plot_ellipses()
+rez1_q = tas.cooper_nathans(hkl_list=(0, 0, 0.5), ei=ei, ef=ef, R0=R0, projection=None)
+
+sample_json_path_2 = "./test_data/IPTS32816_HB1A_exp1034/fesn2.json"
+sample_2 = Xtal.from_json(sample_json_path_2)
+tas.mount_sample(sample_2)
+
+rez2_l = tas.cooper_nathans(hkl_list=(0, 0, 6), ei=ei, ef=ef, R0=R0)
+rez2_l.plot_ellipses()
+rez2_q = tas.cooper_nathans(hkl_list=(0, 0, 6), ei=ei, ef=ef, R0=R0, projection=None)
 
 path_to_spice_folder = "./test_data/IPTS32816_HB1A_exp1034/exp1034/"
 scan35 = Scan.from_spice(path_to_spice_folder, scan_num=35)
+
 fesn000p5_lscan = scan35.get_data(norm_to=(120, "mcu"))
 
+# -----------------------------Fit 1 L-----------------------------
+
+f1_lscan = Fit1D(fesn000p5_lscan)
+f1_lscan.add_background(model="Linear")
+f1_lscan.add_signal(model="Gaussian")
+pars = f1_lscan.guess()
+result = f1_lscan.fit(pars)
+print(f1_lscan.result.fit_report())
+
+p1 = Plot1D()
+p1.add_scan(fesn000p5_lscan, fmt="o")
+p1.add_fit(
+    f1_lscan,
+    label=f"FWHM={f1_lscan.result.params["s1_fwhm"].value:.4f}+/-{f1_lscan.result.params["s1_fwhm"].stderr:.4f}",
+)
+x = f1_lscan.result.params["s1_center"].value
+components = result.eval_components(result.params, x=x)
+y = components["s1_"] / 2 + components["b1_"]
+p1.add_reso_bar(
+    pos=(x, y), fwhm=rez1_l.coh_fwhms(axis=2), c="C3", label=f"Resolution FWHM={rez1_l.coh_fwhms(axis=2):.04f}"
+)
+fig, ax = plt.subplots()
+p1.plot(ax)
+ax.set_title("scan0035")
+
+# -----------------------------Fit 1 Q-----------------------------
+fesn000p5_qscan = scan35.get_data(axes=("q", None), norm_to=(120, "mcu"))
+f1_qscan = Fit1D(fesn000p5_qscan)
+f1_qscan.add_background(model="Linear")
+f1_qscan.add_signal(model="Gaussian")
+pars = f1_qscan.guess()
+result = f1_qscan.fit(pars)
+print(f1_qscan.result.fit_report())
+
+p1_2 = Plot1D()
+p1_2.add_scan(fesn000p5_qscan, fmt="o")
+p1_2.add_fit(
+    f1_qscan,
+    label=f"FWHM={f1_qscan.result.params["s1_fwhm"].value:.4f}+/-{f1_qscan.result.params["s1_fwhm"].stderr:.4f}",
+)
+x = f1_qscan.result.params["s1_center"].value
+components = result.eval_components(result.params, x=x)
+y = components["s1_"] / 2 + components["b1_"]
+p1_2.add_reso_bar(
+    pos=(x, y), fwhm=rez1_q.coh_fwhms(axis=0), c="C3", label=f"Resolution FWHM={rez1_q.coh_fwhms(axis=0):.04f}"
+)
+
+fig, ax = plt.subplots()
+p1_2.plot(ax)
+ax.set_title("scan0035")
+# -----------------------------Fit 2 L-----------------------------
+
+scan50 = Scan.from_spice(path_to_spice_folder, scan_num=50)
+
+substrate006_lscan = scan50.get_data(norm_to=(1, "mcu"))
+
+f2_lscan = Fit1D(substrate006_lscan)
+# f2_lscan.add_background(model="Constant")
+f2_lscan.add_signal(model="Gaussian")
+pars = f2_lscan.guess()
+result = f2_lscan.fit(pars)
+print(f2_lscan.result.fit_report())
+
 p2 = Plot1D()
-p2.add_scan(fesn000p5_lscan, fmt="o")
+p2.add_scan(substrate006_lscan, fmt="o")
+p2.add_fit(
+    f2_lscan,
+    label=f"FWHM={f2_lscan.result.params["s1_fwhm"].value:.4f}+/-{f2_lscan.result.params["s1_fwhm"].stderr:.4f}",
+)
+x = f2_lscan.result.params["s1_center"].value
+components = result.eval_components(result.params, x=x)
+y = components["s1_"] / 2
+p2.add_reso_bar(
+    pos=(x, y), fwhm=rez2_l.coh_fwhms(axis=2), c="C3", label=f"Resolution FWHM={rez2_l.coh_fwhms(axis=2):.04f}"
+)
+
 fig, ax = plt.subplots()
 p2.plot(ax)
+ax.set_title("scan0060")
 
-scan60 = Scan.from_spice(path_to_spice_folder, scan_num=50)
-substrate006_lscan = scan60.get_data(norm_to=(1, "mcu"))
+# -----------------------------Fit 2 Q-----------------------------
+
+substrate006_qscan = scan50.get_data(axes=("q", None), norm_to=(1, "mcu"))
+f2_qscan = Fit1D(substrate006_qscan)
+# f2_lscan.add_background(model="Constant")
+f2_qscan.add_signal(model="Gaussian")
+pars = f2_qscan.guess()
+result = f2_qscan.fit(pars)
+print(f2_qscan.result.fit_report())
+
+p2_2 = Plot1D()
+p2_2.add_scan(substrate006_qscan, fmt="o")
+p2_2.add_fit(
+    f2_qscan,
+    label=f"FWHM={f2_qscan.result.params["s1_fwhm"].value:.4f}+/-{f2_qscan.result.params["s1_fwhm"].stderr:.4f}",
+)
+x = f2_qscan.result.params["s1_center"].value
+components = result.eval_components(result.params, x=x)
+y = components["s1_"] / 2
+p2_2.add_reso_bar(
+    pos=(x, y), fwhm=rez2_q.coh_fwhms(axis=0), c="C3", label=f"Resolution FWHM={rez2_q.coh_fwhms(axis=0):.04f}"
+)
 
-p1 = Plot1D()
-p1.add_scan(substrate006_lscan, fmt="o")
 fig, ax = plt.subplots()
-p1.plot(ax)
+p2_2.plot(ax)
+ax.set_title("scan0060")
 plt.show()

src/tavi/data/fit.py

@@ -211,12 +211,18 @@ def x_to_plot(self, num_of_pts: Optional[int]):
             raise ValueError(f"num_of_points={num_of_pts} needs to be an integer.")
         return x_to_plot
 
-    def eval(self, pars: Parameters, num_of_pts: Optional[int] = 100) -> FitData1D:
+    def eval(self, pars: Optional[Parameters], num_of_pts: Optional[int] = 100, x=None) -> FitData1D:
+        if pars is None:
+            pars = self.result.params
 
-        x_to_plot = self.x_to_plot(num_of_pts)
-        y_to_plot = self.model.eval(pars, x=x_to_plot)
+        if x is not None:
+            return self.model.eval(pars, x=x)
 
-        return FitData1D(x_to_plot, y_to_plot)
+        else:
+            x_to_plot = self.x_to_plot(num_of_pts)
+            y_to_plot = self.model.eval(pars, x=x_to_plot)
+
+            return FitData1D(x_to_plot, y_to_plot)
 
     def fit(self, pars: Parameters) -> ModelResult:
 

src/tavi/instrument/resolution/ellipsoid.py

@@ -4,7 +4,6 @@
 import numpy as np
 from mpl_toolkits.axisartist import Axes
 
-from tavi.instrument.resolution.curve import ResoCurve
 from tavi.instrument.resolution.ellipse import ResoEllipse
 from tavi.plotter import Plot2D
 from tavi.sample.xtal import Xtal
@@ -121,43 +120,59 @@ def volume(self):
 
     def coh_fwhms(self, axis=None):
         """Coherent FWHM"""
-
-        curve = ResoCurve()
         idx = int(axis)
-        curve.fwhm = np.array(sig2fwhm / np.sqrt(self.mat[idx, idx]))
-        if idx == 3:
-            curve.cen = self.en
-        else:
-            curve.cen = self.q[idx]
-
-        curve.r0 = self.r0
-        curve.xlabel = self.axes_labels[idx]
-        curve.title = f"q={np.round(self.hkl,3)}, en={np.round(self.en,3)}"
-        curve.legend = f"coherent FWHM={np.round(curve.fwhm, 3)}"
-        return curve
+        return sig2fwhm / np.sqrt(self.mat[idx, idx])
 
     def incoh_fwhms(self, axis=None):
         """Incoherent FWHMs"""
-
-        curve = ResoCurve()
         idx = int(axis)
 
         reso = self.mat
         for i in (3, 2, 1, 0):
             if not i == idx:
                 reso = ResoEllipsoid.quadric_proj(reso, i)
 
-        curve.fwhm = (1.0 / np.sqrt(np.abs(reso[0, 0])) * sig2fwhm,)
-        if idx == 3:
-            curve.cen = self.en
-        else:
-            curve.cen = self.q[idx]
+        return 1.0 / np.sqrt(np.abs(reso[0, 0])) * sig2fwhm
+
+    # def coh_fwhms(self, axis=None):
+    #     """Coherent FWHM"""
+
+    #     curve = ResoCurve()
+    #     idx = int(axis)
+    #     curve.fwhm = np.array(sig2fwhm / np.sqrt(self.mat[idx, idx]))
+    #     if idx == 3:
+    #         curve.cen = self.en
+    #     else:
+    #         curve.cen = self.q[idx]
+
+    #     curve.r0 = self.r0
+    #     curve.xlabel = self.axes_labels[idx]
+    #     curve.title = f"q={np.round(self.hkl,3)}, en={np.round(self.en,3)}"
+    #     curve.legend = f"coherent FWHM={np.round(curve.fwhm, 3)}"
+    #     return curve
+
+    # def incoh_fwhms(self, axis=None):
+    #     """Incoherent FWHMs"""
+
+    #     curve = ResoCurve()
+    #     idx = int(axis)
+
+    #     reso = self.mat
+    #     for i in (3, 2, 1, 0):
+    #         if not i == idx:
+    #             reso = ResoEllipsoid.quadric_proj(reso, i)
+
+    #     curve.fwhm = (1.0 / np.sqrt(np.abs(reso[0, 0])) * sig2fwhm,)
+    #     if idx == 3:
+    #         curve.cen = self.en
+    #     else:
+    #         curve.cen = self.q[idx]
 
-        curve.r0 = self.r0
-        curve.xlabel = self.axes_labels[idx]
-        curve.title = f"q={np.round(self.hkl,3)}, en={np.round(self.en,3)}"
-        curve.legend = f"incoherent FWHM={np.round(curve.fwhm , 3)}"
-        return curve
+    #     curve.r0 = self.r0
+    #     curve.xlabel = self.axes_labels[idx]
+    #     curve.title = f"q={np.round(self.hkl,3)}, en={np.round(self.en,3)}"
+    #     curve.legend = f"incoherent FWHM={np.round(curve.fwhm , 3)}"
+    #     return curve
 
     # def principal_fwhms_calc(self):
     #     """FWHMs of principal axes in 4D"""
@@ -283,4 +298,5 @@ def plot_ellipses(self):
             )
             p.plot(ax)
 
+        fig.suptitle(f"Q={self.hkl}")
         fig.tight_layout(pad=2)

src/tavi/plotter.py

@@ -11,6 +11,14 @@
 from tavi.instrument.resolution.ellipse import ResoEllipse
 
 
+class ResoBar(object):
+    def __init__(self, pos: tuple[float, float], fwhm: float) -> None:
+
+        self.pos = pos
+        self.fwhm = fwhm
+        self.fmt: dict = {}
+
+
 def tr(x, y, angle):
     x, y = np.asarray(x), np.asarray(y)
     return x + y / np.tan(angle / 180 * np.pi), y
@@ -39,6 +47,7 @@ def __init__(self) -> None:
         # self.ax = None
         self.scan_data: list[ScanData1D] = []
         self.fit_data: list[FitData1D] = []
+        self.reso_data: list[ResoBar] = []
         self.title = ""
         self.xlabel = None
         self.ylabel = None
@@ -98,6 +107,12 @@ def add_fit_components(self, fit_data: Fit1D, num_of_pts: Optional[int] = 100, *
         else:
             raise ValueError(f"Invalid input fit_data={fit_data}")
 
+    def add_reso_bar(self, pos: tuple, fwhm: float, **kwargs):
+        reso_data = ResoBar(pos, fwhm)
+        for key, val in kwargs.items():
+            reso_data.fmt.update({key: val})
+        self.reso_data.append(reso_data)
+
     def plot(self, ax):
         for data in self.scan_data:
             if data.err is None:
@@ -108,6 +123,12 @@ def plot(self, ax):
         for fit in self.fit_data:
             ax.plot(fit.x, fit.y, **fit.fmt)
 
+        for reso in self.reso_data:
+            x, y = reso.pos
+            if "capsize" not in reso.fmt:
+                reso.fmt.update({"capsize": 3})
+            ax.errorbar(x, y, xerr=reso.fwhm / 2, **reso.fmt)
+
         if self.xlim is not None:
             ax.set_xlim(left=self.xlim[0], right=self.xlim[1])
         if self.ylim is not None:

test_data/IPTS32816_HB1A_exp1034/fesn1.json

@@ -0,0 +1,33 @@
+{
+    "type": "xtal",
+    "a": 4.765809,
+    "b": 4.765809,
+    "c": 4.439077,
+    "alpha": 90,
+    "beta": 90,
+    "gamma": 120,
+    "shape": "cylindrical",
+    "mosaic_h": 0,
+    "mosaic_v": 0,
+    "ub_matrix": [
+        -0.004704,
+        0.005858,
+        0.225057,
+        -0.208766,
+        -0.210872,
+        0.000705,
+        0.122875,
+        -0.119174,
+        0.009815
+    ],
+    "plane_normal": [
+        -0.043593,
+        0.008692,
+        0.999012
+    ],
+    "in_plane_ref": [
+        -0.000000,
+        -0.999962,
+        0.008700
+    ]
+}